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