upto(n) = { fibs = List(); ulim = n; for(i = 2, oo, c = fibonacci(i); if(c <= n, listput(fibs, c) , break ); ); res = vector(#fibs + 1, i, oo); process(1, 0, 1); res = Set(res); my(recs = List(), r = 0); for(i = 1, #res, if(res[i] < oo, c = qfibdivsforupto(res[i]); if(c > r, r = c; listput(recs, res[i]); ); ); ); recs } qfibdivsforupto(n) = 1 + sum(i = 2, #fibs, n%fibs[i] == 0) process(n, ind, qd) = { c = qfibdivsforupto(n); res[c] = min(res[c], n); my(cn, i, cq); for(i = ind + 1, #fibs, cn = lcm(n, fibs[i]); if(cn <= ulim, cq = qfibdivsforupto(cn); process(cn, i, cq); ); ); }