1 /*
2  * $Id: gcd.i,v 1.1 2005-09-18 22:06:00 dhmunro Exp $
3  * GCD, LCM, and prime factorization routines.
4  */
5 /* Copyright (c) 2005, The Regents of the University of California.
6  * All rights reserved.
7  * This file is part of yorick (http://yorick.sourceforge.net).
8  * Read the accompanying LICENSE file for details.
9  */
10 
gcd(a,b)11 func gcd(a, b)
12 /* DOCUMENT gcd(a,b)
13      returns the GCD (greatest common divisor) of A and B, which must
14      be one of the integer data types.  A and B may be conformable
15      arrays; the semantics of the gcd call are the same as any other
16      binary operation.  Uses Euclid's celebrated algorithm.
17      The absolute values of A and B are taken before the operation
18      commences; if either A or B is 0, the return value will be 0.
19    SEE ALSO: lcm, is_prime, factorize
20  */
21 {
22   a= abs(a);
23   b= abs(b);
24   c= min(a, b);
25   a= max(a, b);
26   b= c;          /* simplifies c=0 case */
27 
28   if (dimsof(a)(1)) {
29     /* array case */
30     for (list=where(c) ; numberof(list) ; list=where(c)) {
31       b(list)= bl= c(list);
32       c(list)= a(list) % bl;
33       a(list)= bl;
34     }
35 
36   } else {
37     /* scalar case can be less baroque */
38     while (c) {
39       b= c;
40       c= a % b;
41       a= b;
42     }
43   }
44 
45   return b;
46 }
47 
lcm(a,b)48 func lcm(a, b)
49 /* DOCUMENT lcm(a,b)
50      returns the LCM (least common multiple) of A and B, which must
51      be one of the integer data types.  A and B may be conformable
52      arrays; the semantics of the lcm call are the same as any other
53      binary operation.
54      The absolute values of A and B are taken before the operation
55      commences; if either A or B is 0, the return value will be 0.
56    SEE ALSO: gcd, is_prime, factorize
57  */
58 {
59   d= gcd(a, b);
60   /* two potential problems: zero divide and overflow - handle the
61      first but not the second */
62   return abs(a*b)/(d+!d);
63 }
64 
is_prime(x)65 func is_prime(x)
66 /* DOCUMENT is_prime(x)
67      return non-zero if and only if X (which must be a scalar integer)
68      is prime.  May return a false positive if X is greater than about
69      3e9, since at most 20000 candidate factors are checked.
70      The absolute value of X is taken first; zero is not prime, but 1 is.
71    SEE ALSO: gcd, lcm, factorize
72  */
73 {
74   x= long(abs(x));
75   if (x<2) return x==1;
76   /* make a list of factors which includes 2, 3, and all larger
77      odd numbers not divisible by 3 less or equal to sqrt(x) */
78   top= min(long((sqrt(x)+8.)/3.+0.5), 20000);
79   factors= ((3*indgen(2:top))/2)*2 - 7;
80   factors(1:2)= [2,3];
81   return allof(x%factors);
82 }
83 
factorize(x)84 func factorize(x)
85 /* DOCUMENT factorize(x)
86      return list of prime factors of X and their powers as an n-by-2
87      array.  May include a large non-prime factor if X exceeds 3e9.
88      In any event, product(result(,1)^result(,2)) will equal abs(X).
89      X must be a scalar integer type.
90    SEE ALSO: gcd, lcm, is_prime
91  */
92 {
93   x= long(abs(x));
94   if (x<2) return [[x],[1]];
95 
96   /* first get rid of any prime factors less than 102 */
97   primes= [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,
98            71,73,79,83,89,97,101];
99   primes= primes(where(!(x%primes)));
100   powers= _fact_extract(primes, x);   /* returns "deflated" x */
101 
102   if (x>1) {
103     /* large prime factors require a less direct approach */
104     top= min(long((sqrt(x)+2.)/3.+0.5), 20000);
105     if (top>=35) {
106       /* trial divisors are all odd numbers 103 or greater which are
107          not divisible by three and do not exceed sqrt(x) */
108       trial= ((3*indgen(35:top))/2)*2 - 1;
109       /* discard all trial divisors which do not divide x */
110       trial= trial(where(!(x%trial)));
111       /* the smallest remaining divisor must be prime - remove it and
112          all its subsequent multiples to find the next prime divisor
113          and so on until the list contains only primes */
114       list= trial;
115       for (n=0 ; numberof(trial) ; ++n) {
116         list(n+1)= trial(1);
117         trial= trial(where(trial%trial(1)));
118       }
119       if (n) {
120         trial= list(1:n);
121         grow, primes, trial;
122         grow, powers, _fact_extract(trial,x);
123       }
124     }
125     if (x>1) {
126       grow, primes, [x];
127       grow, powers, [1];
128     }
129   }
130 
131   return [primes, powers];
132 }
133 
134 func _fact_extract(primes, &x)
135 {
136   if (is_void(primes)) return [];
137   /* first get largest power of each prime less than or equal x */
138   powers= long(log(x)/log(primes)+1.e-6);
139   factors= gcd(primes^powers, x);
140   x/= long(exp(sum(log(factors)))+0.5);
141   return long(log(factors)/log(primes)+0.5);
142 }
143