options nosource nomprint;
dm 'log; autoscroll 1' log;
/* Richard A. DeVenezia
* Oct. 26, 2004
*
* Amicable numbers
*
* Per SAS-L Post by Martin OConnell
My son is working on a math report on Euler and one of the things
he is highlighting is Euler's work on "Amicable Numbers". These are
non "Perfect") numbers whose factors add up the the other number. The
smallest pair being (220,284), since
1+2+4+5+10+11+20+22+44+55+110 = 284 (factors of 220)
1+2+4+71+142 = 220 (factors of 284).
* Let SF(N) be the sum of the factors of N (excluding N).
* Two numbers are amicable if
* N1 = SF(N2) and
* N2 = SF(N1)
*
* Notes:
* A SAS numeric can represent every integer from 1..2**53 (~9e15 or nine quadrillion).
* After 2**53, some integers can not be represented. (see IEEE 754)
* The product of the first 14 primes [2x3x5x7x11x13x17x19x23x29x31x37x41x43]
* exceeds 2**53.
* Thus, in this program, 13 is the most number of distinct prime factors
* that a number (which is also assuredly represented in SAS) can have.
*
*/
%let N_SEARCH_FROM = 0; * at which number should searches for amicable numbers begin;
%let DURATION = 3; *60; *3600; * how long should the program run (seconds) ;
%let N_PRIMES = 1e5; * how many prime numbers should be tracked? ;
%let N_PRIMES = %sysevalf (&N_PRIMES + 0);
data Amicable (keep=N number_to_factor elapsed rate rename=(number_to_factor=N2))
; */ debug;
* Prep prime array;
array prime [ &N_PRIMES ] _temporary_;
prime [ 1 ] = 2;
prime [ 2 ] = 3;
__ix = 2; * the index of the last computed prime;
__sqrt_of = 3;
__sqrt = 1;
t0 = datetime();
N = &N_SEARCH_FROM - 1;
do until (datetime() - t0 > &DURATION);
N + 1;
number_to_factor = N;
link FACTOR;
link FACTORSUM;
* Only list candidates N1 < N2;
if N >= factor_sum then continue;
number_to_factor = factor_sum;
link FACTOR;
link FACTORSUM;
* Amicable check;
if N ne factor_sum then continue;
* Compute global throughput at this point;
elapsed = datetime() - t0;
if elapsed then
rate = N / elapsed;
* Log amicable numbers;
put N '-> ' number_to_factor '-> ' factor_sum ;
output;
end;
* After duration has elapsed report last
* number tested and last prime computed;
put N= prime[__ix]=;
stop;
NEXT_PRIME:
* Find the next prime after the last computed one;
__N = prime [ __ix ] ;
NEXT_CANDIDATE:
* Test each candidate for divisibility
* by each prior known prime (upto sqrt of candidate).
* When divisibility testing fails, the candidate
* is known to be prime and is stored in the next
* entry of the prime array
*;
__N = __N + 2;
if __N >= __sqrt_of then do;
__sqrt + 1;
__sqrt_of = __sqrt ** 2;
end;
__p = prime [ 1 ];
do __j = 1 to __ix while ( __p <= __sqrt );
if mod (__N,__p) = 0 then goto NEXT_CANDIDATE;
__p = prime [ __j ];
end;
__ix + 1;
prime [ __ix ] = __N;
* put prime [ __ix ] =;
return;
FACTOR:
* A SAS numeric can represent every integer from 1..2**53 (~9e15 or nine quadrillion).
* After 2**53, some integers can not be represented. (see IEEE 754)
* The product of the first 14 primes [2x3x5x7x11x13x17x19x23x29x31x37x41x47]
* exceeds 2**53.
* Thus, in this program, 13 is the most number of distinct prime factors
* that a number (which is also assuredly represented in SAS) can have.
*;
array PF[13,2] _temporary_; * ,1 will store prime factor : ,2 power;
_n = number_to_factor;
_nf = 0;
_ix = 1;
_factor = prime [ _ix ];
_sqrt_n = sqrt(_n);
do while ( _factor <= _sqrt_n );
* determine power of this prime factor;
_power = 0;
do while (mod (_n, _factor) = 0);
_power + 1;
_n = _n / _factor;
end;
* record prime factor and its power;
if _power then do;
_nf + 1;
PF [ _nf , 1 ] = _factor;
PF [ _nf , 2 ] = _power;
if _n = 1 then leave;
_sqrt_n = sqrt(_n);
end;
* select next prime for divisibility testing;
if _ix < &N_PRIMES then do;
_ix + 1;
if _ix > __ix then link NEXT_PRIME;
_factor = prime [ _ix ];
end;
else do;
* this will happen when pi(N) > dim(primes);
put 'WARNING: Out of prime real estate. Skipping factorization of ' number_to_factor;
return;
end;
end;
if _n > 1 and _nf then do;
_nf + 1;
PF [ _nf , 1 ] = _n;
PF [ _nf , 2 ] = 1;
end;
return; * comment this line to show factorization;
put / number_to_factor ':';
do _i = 1 to _nf;
do _j = 1 to dim(pf,2);
put PF [_i,_j] = @;
end;
put;
end;
return;
FACTORSUM:
* Determine each combination of the prime factors, furthermore
* iterate from 1 to the power of the prime factor to determine
* individual factors for use in summation
*
* The combinatoric aspect is generated using nested loops,
* which must be done using macro.
*;
%macro comboDos (n=, debug=no);
%local i;
ix0 = 0;
array ix ix1-ix&n;
array ixp ixp1-ixp&n;
%do i = 1 %to &n;
do ix&i = ix%eval(&i-1)+1 to &n while (ix&i<=_nf);
* iterate from 1 to the power of the prime factor;
do ixp&i = 1 to PF[ix&i,2];
* compute individual factor from a combinatoric perspective;
factor = 1;
do j = 1 to &i;
factor = factor * ( PF[ix[j],1] ** ixp[j] );
end;
* put ix1-ix&i @8 ': ' ixp1-ixp&i @18 factor=;
factor_sum + factor;
%end;
%mend;
%macro comboEnds (n=);
%local i ;
%do i = 1 %to &n;
end;
end;
%end;
%mend;
if _nf = 0 then do;
factor_sum = 1;
return;
end;
options nomprint;
factor_sum = 1 - number_to_factor;
%comboDos (n=13)
%comboEnds (n=13)
return;
run ;
/*
goptions reset=all;
symbol1 i=rQ v=dot c=black;
data amy;
set amicable;
logn = log(n);
lograte = log(rate);
run;
proc gplot data=amy;
plot lograte * logn ;
where logn > 10;
run;
quit;
data _null_;
N = 1e9;
logn = log(N);
lograte = 11.89636 + 0.043*logn - 0.008567*logn**2; * 2minute;
rate = exp (lograte);
time = N / rate;
put lograte= rate= time= time8.;
lograte = 12.02125 + 0.027565*logn - 0.008026*logn**2; *5minute;
rate = exp (lograte);
time = N / rate;
put lograte= rate= time= time8.;
lograte = 12.00645 + 0.028099*logn - 0.008029*logn**2; *10minute;
rate = exp (lograte);
time = N / rate;
put lograte= rate= time= time8.;
lograte = 12.95519 - 0.120549*logn - 0.002202*logn**2; *10s;
rate = exp (lograte);
time = N / rate;
put lograte= rate= time= time8.;
run;
*/
options source;