carmichael:=
proc(k)
// find all Carmichael numbers up to k 
/* use the following theorem: n is a Carmichael number if and only if
  n is squarefree, not prime, and p-1 divides n-1 for any prime divisor p of n 
*/
local a,i,b;

begin
for i from 2 to k do
	a:= stdlib::ifactor(i);
	if nops(a)>3 then
		if { op(a,2*l+1) $ l=1..nops(a)-1 div 2 } = {1} then 
			b:={ op(a,2*l)-1 $l=1..nops(a) div 2 };
			if map(b, n->  i-1 mod n)  = { 0} then
				 print(i)
         end_if;
		end_if;
	end_if;
end_for;
end_proc:

/* Running time: a bit slow. Reason: even if p-1 does not divide n-1
  for a certain p, the other prime factors of n are tested nevertheless. */

// time(carmichael(10000)) = 8930