upto(n) = {
	pr = primes(primepi(n)); pr = concat(pr, 2*n);
	prp = vector(n); 
	t = 0; 
	prind = 1;
	for(i = 1, n, 
		if(i >= pr[prind],
			prind++; 
		);
		prp[i] = prind-1;
	);
	
	my(res = List()); 
	inc = 10000; 
	printat = inc; 
	for(i = 2, #pr, 
		pr[i]+=pr[i-1]
	);
	for(i = 1, n, 
		if(i > printat, 
			print1(i\inc", "); 
			printat += inc; 
		);
		if(is(i), 
			listput(res, i); 
			print(); print(); print(res); print(); print(); 
		)
	)
}

addhelp(is, "See if n is a term. Call this from upto above as it uses precomputed stuff from there.")
is(n) = {
	my(d = divisors(n), s = 0);
	for(i = 1, #d-1, 
		s += d[i] * pr[prp[n/d[i]]];
		f = factor(n / d[i]);
		for(j = 1, #f~, 
			s-=(d[i] * f[j, 1])
		)
	);
	s % n == 0
}