Date: Wed, 27 Oct 2004 08:47:32 -0400
Reply-To: "Richard A. DeVenezia" <radevenz@IX.NETCOM.COM>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: "Richard A. DeVenezia" <radevenz@IX.NETCOM.COM>
Subject: Re: SAS and Amicable Numbers
Martin O'Connell wrote:
> Hi all,
>
> This probably falls under the heading of "if you have any code and/or
> thoughts laying around to help, please pass on". Otherwise delete
> this one.
>
> Anyway, 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).
>
> Not required, but he also wants to see how many he can find
> numerically. I thought it was a harmless enough problem to attack
> with SAS and we came up with the following to generate the pairs <=
> 100000:
>
> data numbers;
> do i = 1 to 100000;
> factors = 0;
> do j = 1 to int(i/2) + 1;
> if i/j eq int(i/j) then factors + j;
> end;
> number = i;
> output;
> end;
> keep number factors;
> run;
>
> Then I merge/purge to get the necessary pairs.
>
> This is fine as written but this factoization will not extend much
> higher before things get miserably slower.
>
> My questions are:
> 1) Are there any better factorization methods that SAS has, and
> 2) (Math Question) When I am looking for a amicable pair for n, is it
> necessary to consider 1, 2, ..., n-1.
>
> Thanks for any help.
>
> Martin
Here is a program that factors numbers, computes primes, and combines prime
factors into individual factors for summation. On a P4 3ghz it might take
45 minutes to check the first million integers.
I also found these pages interesting and helpful
http://xraysgi.ims.uconn.edu:8080/amicable.html
http://en.wikipedia.org/wiki/Amicable_numbers
To see the puts as the step runs, be sure to activate the LOG window first,
and
issue command AUTOSCROLL 1. Then submit the program
==================
dm 'clear log' log;
/* Richard A. DeVenezia
* Oct. 26, 2004
*
* Amicable numbers
*
* Per SAS-L Post by Martin OConnell
Anyway, 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).
...
My factoization code will not extend much higher
before things get miserably slower.
My questions are:
1) Are there any better factorization methods that SAS has, and
2) (Math Question) When I am looking for a amicable pair for n, is it
necessary to consider 1, 2, ..., n-1.
*/
%let N_NATURALS = 1e6;
data Amicable (keep=N);
link PREP_PRIME;
* Note: at N=669688 the factorization of other in the aliquot will
cause a warning;
do N = 1 to &N_NATURALS;
number_to_factor = N;
link FACTOR;
link FACTORSUM;
number_to_factor = factor_sum;
link FACTOR;
link FACTORSUM;
if N = factor_sum then do;
no_time = datetime() - no_start;
if no_time > 0 then
no_rate = no_count / no_time;
* report amicable numbers, the number of factors
* and the average number of numbers examined per second since
* the last amicable number reported;
put N '-> ' number_to_factor '-> ' factor_sum _nf= no_rate;
output;
no_count = 0;
no_start = datetime();
end;
else do;
no_count + 1;
end;
end;
put N=;
stop;
PREP_PRIME:
array prime [ 100000 ] _temporary_;
prime [ 1 ] = 2;
prime [ 2 ] = 3;
__ix = 2;
return;
NEXT_PRIME:
__N = prime [ __ix ] ;
NEXT_CANDIDATE:
* test each candidate for divisibility
* by a prior known prime (upto sqrt of candidate)
* when a divisibility testing fails, the candidate
* is stored in the next entry of the prime array
*;
__N = __N + 2;
__sqrt = sqrt ( __N );
__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:
* SAS numeric can store every integer from 1..2**53.
* After 2**53, some integers can not be stored in a SAS variable
* The product of the first 14 primes exceeds 2**53.
* Thus the most number of distinct prime factors that a number
* (which is also a SAS precise integer) can have is 13.
*;
array PF[13,2] _temporary_; * ,1 will store prime factor : ,2
power;
_n = number_to_factor;
_nf = 0;
_ix = 0;
_power = 0;
* clear factorization array;
do _i = 1 to hbound(pf,1);
do _j = 1 to hbound(pf,2);
PF [_i,_j] = .;
end;
end;
do while ( _n > 1 );
* select next prime for divisibility testing;
if _ix < hbound(prime) 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;
* determine power of this prime factor;
do while (mod (_n, _factor) = 0);
_power + 1;
_n = _n / _factor;
end;
* record this prime factor and its power;
if _power then do;
_nf + 1;
PF [ _nf , 1 ] = _factor;
PF [ _nf , 2 ] = _power;
_power = 0;
end;
end;
return; * comment return to show factorization;
do _i = 1 to _nf;
do _j = 1 to hbound(pf,2);
put PF [_i,_j] = @;
end;
put;
end;
return;
FACTORSUM:
* This gets a little hairy, partly because there are are no
built-in
* programmatic statements for doing recursion in datastep.
* (
* recursion is possible, but special considerations are needed.
* See Chang Chung SAS-L post for more info
*
http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0401C&L=sas-l&P=R2603
* )
*
* 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
* unfortunately 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 ii;
%do ii = 1 %to &n;
%let i = %eval(&n-&ii+1);
end;
end;
%end;
%mend;
options mprint;
factor_sum = 1 - number_to_factor;
%comboDos (n=13)
%comboEnds (n=13)
return;
run;
==================
--
Richard A. DeVenezia
http://www.devenezia.com/
