# After Rémy Sigrist's PARI program. rad := n -> NumberTheory:-Radical(n): bitset := (m, v) -> Bits:-GetBits(m, v) = 1: best := []: explore := proc(n, seq, m, r, w) global best; local v, rv; if w = 1 and nops(best) < nops(seq) then best := seq fi; for v from 1 to n do if not bitset(m, v) and igcd(r, v) = w then rv := rad(v); explore(n, [op(seq), v], m + 2^v, rv, rv/w) fi; od; end: row := proc(n) global best; local r, v; r := 1; best := []; for v from 1 to n do r := rad(v); explore(n, [v], 2^v, r, r) od; return best; end: PrintTriangle := proc(upto) local n, r; for n from 1 to upto do if not isprime(n) then r := row(n) fi; lprint([n], r); od; end: PrintTriangle(13);