addhelp(upto, "all positive integers k where k/A007954(k) <= f and is integer.");
upto(f) = {
	ints = vector(floor(f));
	fracs = vector(floor(f)); 
	maxn = f; 
	for(i = 1, 9, 
		process([[[i], i]]);
	); 
	res = vector(vecsum(ints));
	t = 0; 
	ind = 0; 
	for(i = 1, #ints,
		for(j = 1, ints[i], 
			t++; 
			ind++; 
			res[ind] = t;
		); 
		t+=fracs[i];
	);
	res
}

process(v) = {
	my(c);
	for(i = 1, #v, 
		processperms(v[i][1]);
		c = childrenbounded(v[i][1], v[i][2]);
		process(c);
	)	
}

processperms(v) = {
vp = vecprod(v);	
forperm(v, p, 		
	c = fromdigits(Vec(p))/vp; 
	if(c > maxn, 
		return
	);
	if(denominator(c) > 1, 
		fracs[floor(c)]++
	,
		ints[floor(c)]++	
	)
)	
}

childrenbounded(v, p) = {
	my(res = List(), v = concat(0, v)); 
	for(i = 1, v[2],
		v[1] = i; 
		fr = fromdigits(v); 
		if(fr\(p*i) <= maxn, 
			listput(res, [v, p*i]);
		)
	); res
}