Martin O'Connell wrote: > 1) Are there any better factorization methods that SAS has, and

My earlier post ran quite slowly due to a dumb boundary check.

While factorizing there is a while (_n_ > 1) which means stop factoring when last prime divisor is checked. It also means that when the original number is prime, it is checked for divisibility by all known prior primes (of which all will fail) causing a huge amount of useless checking.

A much smarter bound while (_factor > _sqrt_n) causes the loop to stop when no larger prime could possibly divide the subfactor candidate.

The dumb way took 25min to scan 1e6 integers, the smarter way only 25sec ! A 60x improvment.

==================== dm 'clear log' log;

/* Richard A. DeVenezia * Oct. 26, 2004 * * Amicable numbers * * Per SAS-L Post by Martin O'Connell:

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; %let N_PRIMES = 1e5;

%let N_PRIMES = %sysevalf (&N_PRIMES + 0); %put &N_PRIMES;

data Amicable (keep=N); * / debug;

link PREP_PRIME;

* Note: at N=669688 the other in the aliquot will blow the prime stack;

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 12.-L; output; no_count = 0; no_start = datetime(); end; else do; no_count + 1; end; end;

put N=;

stop;

PREP_PRIME: array prime [ &N_PRIMES ] _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;

_sqrt_n = sqrt(_n);

do until ( _factor > _sqrt_n );

* 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;

* 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;

_sqrt_n = sqrt(_n); end; end;

if _n > 1 and _nf then do; _nf + 1; PF [ _nf , 1 ] = _n; PF [ _nf , 2 ] = 1; end;

return; *comment return to show factorization; do _i = 1 to _nf; do _j = 1 to dim(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;

if _nf = 0 then do; factor_sum = 1; return; end;

options mprint;

factor_sum = 1 - number_to_factor;

%comboDos (n=13) %comboEnds (n=13) return;

run; ====================

-- Richard A. DeVenezia http://www.devenezia.com/