Constrained Cartesian Product

def product_with(base, *rest, &constraint) 
  return base.map{|a|a} if rest.empty?
  rest = product_with(*rest, &constraint)
  base.inject([]){|rs, a| 
  rest.inject(rs){|rs, b|
    if constraint.call a, b.kind_of?(Array)? b[0]: b
      rs<<[a, *b]
    else
      rs
    end
  }}		
end

p product_with(1..3, 1..3, 1..3, &proc{|a, b| true})  # no constraint
p product_with(1..3, 1..3, 1..3, &proc{|a, b| a < b}) # with constraint
p product_with(1..3, 1..4, 1..5, 1..5, &proc{|a, b| a < b}) # more args

Output:

[[1, 1, 1], [1, 1, 2], [1, 1, 3],
 [1, 2, 1], [1, 2, 2], [1, 2, 3],
 [1, 3, 1], [1, 3, 2], [1, 3, 3],
 [2, 1, 1], [2, 1, 2], [2, 1, 3],
 [2, 2, 1], [2, 2, 2], [2, 2, 3],
 [2, 3, 1], [2, 3, 2], [2, 3, 3],
 [3, 1, 1], [3, 1, 2], [3, 1, 3],
 [3, 2, 1], [3, 2, 2], [3, 2, 3],
 [3, 3, 1], [3, 3, 2], [3, 3, 3]]
[[1, 2, 3]]
[[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]

An Application to the Pythagorean Triples

# Pythagorean Triples
r = 1..100
ts = product_with(r, r, r, &proc{|a, b| a < b}
).select{|x, y, z|
    x**2 + y**2 == z**2
}

# reduction and print
ts.inject(ts){|ps, t|
  ps.select{|p|
    t[0]==p[0] && t[1]==p[1] && t[2]==p[2] || 
    p[0]%t[0]!=0 || p[1]%t[1]!=0 || p[2]%t[2]!=0 
  }
}.each {|p| 
  puts "#{p[0]}^2 + #{p[1]}^2 == #{p[2]}^2"
}

3^2 + 4^2 == 5^2
5^2 + 12^2 == 13^2
7^2 + 24^2 == 25^2
8^2 + 15^2 == 17^2
9^2 + 40^2 == 41^2
11^2 + 60^2 == 61^2
12^2 + 35^2 == 37^2
13^2 + 84^2 == 85^2
16^2 + 63^2 == 65^2
20^2 + 21^2 == 29^2
28^2 + 45^2 == 53^2
36^2 + 77^2 == 85^2
39^2 + 80^2 == 89^2
48^2 + 55^2 == 73^2
65^2 + 72^2 == 97^2