59 lines
1.4 KiB
Plaintext
59 lines
1.4 KiB
Plaintext
duffinianNumbers: procedure options(main);
|
|
%replace MAXSIGMA by 10000;
|
|
declare sigma (1:MAXSIGMA) fixed;
|
|
|
|
calculateSigmaTable: procedure;
|
|
declare (i, j) fixed;
|
|
do i=1 to MAXSIGMA;
|
|
sigma(i) = 0;
|
|
end;
|
|
do i=1 to MAXSIGMA;
|
|
do j=i to MAXSIGMA by i;
|
|
sigma(j) = sigma(j) + i;
|
|
end;
|
|
end;
|
|
end calculateSigmaTable;
|
|
|
|
gcd: procedure(aa,bb) returns(fixed);
|
|
declare (a, aa, b, bb, c) fixed;
|
|
a = aa;
|
|
b = bb;
|
|
do while(b > 0);
|
|
c = mod(a,b);
|
|
a = b;
|
|
b = c;
|
|
end;
|
|
return(a);
|
|
end gcd;
|
|
|
|
duffinian: procedure(n) returns(bit);
|
|
declare n fixed;
|
|
return(sigma(n) > n+1 & gcd(n, sigma(n)) = 1);
|
|
end duffinian;
|
|
|
|
triplet: procedure(n) returns(bit);
|
|
declare n fixed;
|
|
return(duffinian(n) & duffinian(n+1) & duffinian(n+2));
|
|
end triplet;
|
|
|
|
declare (i, n) fixed;
|
|
|
|
call calculateSigmaTable;
|
|
put skip list('First 50 Duffinian numbers:');
|
|
put skip;
|
|
n=0;
|
|
do i=1 to 50;
|
|
do n=n+1 repeat(n+1) while(^duffinian(n)); end;
|
|
put edit(n) (F(5));
|
|
if mod(i,10) = 0 then put skip;
|
|
end;
|
|
|
|
put skip;
|
|
put skip list('First 15 Duffinian triplets:');
|
|
n=0;
|
|
do i=1 to 15;
|
|
do n=n+1 repeat(n+1) while(^triplet(n)); end;
|
|
put skip edit(n, n+1, n+2) (F(7),F(7),F(7));
|
|
end;
|
|
end duffinianNumbers;
|