60 lines
1.6 KiB
Text
60 lines
1.6 KiB
Text
smith: procedure options(main);
|
|
/* find the digit sum of N */
|
|
digitSum: procedure(nn) returns(fixed);
|
|
declare (n, nn, s) fixed;
|
|
s = 0;
|
|
do n=nn repeat(n/10) while(n>0);
|
|
s = s + mod(n,10);
|
|
end;
|
|
return(s);
|
|
end digitSum;
|
|
|
|
/* find and count factors of N */
|
|
factors: procedure(nn, facs) returns(fixed);
|
|
declare (n, nn, cnt, fac, facs(16)) fixed;
|
|
cnt = 0;
|
|
if nn<=1 then return(0);
|
|
|
|
/* factors of two */
|
|
do n=nn repeat(n/2) while(mod(n,2)=0);
|
|
cnt = cnt + 1;
|
|
facs(cnt) = 2;
|
|
end;
|
|
|
|
/* take out odd factors */
|
|
do fac=3 repeat(fac+2) while(fac <= n);
|
|
do n=n repeat(n/fac) while(mod(n,fac) = 0);
|
|
cnt = cnt + 1;
|
|
facs(cnt) = fac;
|
|
end;
|
|
end;
|
|
|
|
return(cnt);
|
|
end factors;
|
|
|
|
/* see if a number is a Smith number */
|
|
smith: procedure(n) returns(bit);
|
|
declare (n, nfacs, facsum, i, facs(16)) fixed;
|
|
nfacs = factors(n, facs);
|
|
if nfacs <= 1 then
|
|
return('0'b); /* primes are not Smith numbers */
|
|
|
|
facsum = 0;
|
|
do i=1 to nfacs;
|
|
facsum = facsum + digitSum(facs(i));
|
|
end;
|
|
return(facsum = digitSum(n));
|
|
end smith;
|
|
|
|
/* print all Smith numbers up to 10000 */
|
|
declare (i, cnt) fixed;
|
|
cnt = 0;
|
|
do i=2 to 9999;
|
|
if smith(i) then do;
|
|
put edit(i) (F(5));
|
|
cnt = cnt + 1;
|
|
if mod(cnt,16) = 0 then put skip;
|
|
end;
|
|
end;
|
|
put skip list('Found', cnt, 'Smith numbers.');
|
|
end smith;
|