48 lines
1.2 KiB
Text
48 lines
1.2 KiB
Text
program duffinian_numbers;
|
|
init sigma := divisor_sum_table(20000);
|
|
|
|
print("First 50 Duffinian numbers:");
|
|
loop for n in first(50, routine is_duffinian) do
|
|
nprint(lpad(str n, 6));
|
|
if (col +:= 1) mod 10 = 0 then print; end if;
|
|
end loop;
|
|
|
|
print;
|
|
print("First 15 Duffinian triplets:");
|
|
loop for n in first(15, routine is_duffinian_triplet) do
|
|
print(+/[lpad(str i, 6) : i in [n,n+1,n+2]]);
|
|
end loop;
|
|
|
|
proc first(num, pred);
|
|
ls := [];
|
|
loop while #ls < num do
|
|
loop until call(pred, n) do
|
|
n +:= 1;
|
|
end loop;
|
|
ls with:= n;
|
|
end loop;
|
|
return ls;
|
|
end proc;
|
|
|
|
proc is_duffinian_triplet(n);
|
|
return and/[is_duffinian(i) : i in [n,n+1,n+2]];
|
|
end proc;
|
|
|
|
proc is_duffinian(n);
|
|
return sigma(n) > n+1 and gcd(n, sigma(n)) = 1;
|
|
end proc;
|
|
|
|
proc gcd(a,b);
|
|
return if b=0 then a else gcd(b, a mod b) end;
|
|
end proc;
|
|
|
|
proc divisor_sum_table(sz);
|
|
ds := [0] * sz;
|
|
loop for i in [1..sz] do
|
|
loop for j in [i,i+i..sz] do
|
|
ds(j) +:= i;
|
|
end loop;
|
|
end loop;
|
|
return ds;
|
|
end proc;
|
|
end program;
|