It's fun to translate terse APL to somewhat terse numpy. The result still can be very compact and you can parse it easily if you're used to looking at numpy:
s = arange(2, 50); p = outer(s, s).ravel(); sorted(set(s) - set(p)) /:~ s -. p [ p =: s*/s [ s=: 2+i.48
test =: (* = *. * +.) & ?
*./ test~ 1000 # 100 def d(x): N = arange(1, x + 1); return N[x % N == 0]
def m(x, n): return x * arange(1, n + 1)
def gcd(x, y): return max(set(d(x)) & set(d(y)))
def lcm(x, y): return min(set(m(x, y)) & set(m(y, x)))
def test(x, y): return gcd(x, y) * lcm(x, y) == x * y
all([test(x, y) for (x, y) in randint(1, 100, (1000, 2))]) # True
⎕io←0 ⍝ to match J
test ← (× = ∨ × ∧) ∘ ?
∧⌿ test⍨ 1000⍴100