/* GCD.I GCD, LCM, and prime factorization routines. $Id$ */ /* Copyright (c) 1995. The Regents of the University of California. All rights reserved. */ func gcd (a, b) /* DOCUMENT gcd(a,b) returns the GCD (greatest common divisor) of A and B, which must be one of the integer data types. A and B may be conformable arrays; the semantics of the gcd call are the same as any other binary operation. Uses Euclid's celebrated algorithm. The absolute values of A and B are taken before the operation commences; if either A or B is 0, the return value will be 0. SEE ALSO: lcm, is_prime, factorize */ { a= abs(a); b= abs(b); c= min(a, b); a= max(a, b); b= c; /* simplifies c=0 case */ if (dimsof(a)(1)) { /* array case */ for (list=where(c) ; numberof(list) ; list=where(c)) { b(list)= bl= c(list); c(list)= a(list) % bl; a(list)= bl; } } else { /* scalar case can be less baroque */ while (c) { b= c; c= a % b; a= b; } } return b; } func lcm (a, b) /* DOCUMENT lcm(a,b) returns the LCM (least common multiple) of A and B, which must be one of the integer data types. A and B may be conformable arrays; the semantics of the lcm call are the same as any other binary operation. The absolute values of A and B are taken before the operation commences; if either A or B is 0, the return value will be 0. SEE ALSO: gcd, is_prime, factorize */ { d= gcd(a, b); /* two potential problems: zero divide and overflow - handle the first but not the second */ return abs(a*b)/(d+!d); } func is_prime (x) /* DOCUMENT is_prime(x) return non-zero if and only if X (which must be a scalar integer) is prime. May return a false positive if X is greater than about 3e9, since at most 20000 candidate factors are checked. The absolute value of X is taken first; zero is not prime, but 1 is. SEE ALSO: gcd, lcm, factorize */ { x= long(abs(x)); if (x<2) return x==1; /* make a list of factors which includes 2, 3, and all larger odd numbers not divisible by 3 less or equal to sqrt(x) */ top= min(long((sqrt(x)+8.)/3.+0.5), 20000); factors= ((3*indgen(2:top))/2)*2 - 7; factors(1:2)= [2,3]; return allof(x%factors); } func factorize (x) /* DOCUMENT factorize(x) return list of prime factors of X and their powers as an n-by-2 array. May include a large non-prime factor if X exceeds 3e9. In any event, product(result(,1)^result(,2)) will equal abs(X). X must be a scalar integer type. SEE ALSO: gcd, lcm, is_prime */ { x= long(abs(x)); if (x<2) return [[x],[1]]; /* first get rid of any prime factors less than 102 */ primes= [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67, 71,73,79,83,89,97,101]; primes= primes(where(!(x%primes))); powers= _fact_extract(primes, x); /* returns "deflated" x */ if (x>1) { /* large prime factors require a less direct approach */ top= min(long((sqrt(x)+2.)/3.+0.5), 20000); if (top>=35) { /* trial divisors are all odd numbers 103 or greater which are not divisible by three and do not exceed sqrt(x) */ trial= ((3*indgen(35:top))/2)*2 - 1; /* discard all trial divisors which do not divide x */ trial= trial(where(!(x%trial))); /* the smallest remaining divisor must be prime - remove it and all its subsequent multiples to find the next prime divisor and so on until the list contains only primes */ list= trial; for (n=0 ; numberof(trial) ; ++n) { list(n+1)= trial(1); trial= trial(where(trial%trial(1))); } if (n) { trial= list(1:n); grow, primes, trial; grow, powers, _fact_extract(trial,x); } } if (x>1) { grow, primes, [x]; grow, powers, [1]; } } return [primes, powers]; } func _fact_extract(primes, &x) { if (is_void(primes)) return []; /* first get largest power of each prime less than or equal x */ powers= long(log(x)/log(primes)+1.e-6); factors= gcd(primes^powers, x); x/= long(exp(sum(log(factors)))+0.5); return long(log(factors)/log(primes)+0.5); }