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);
		);
	);
	
}