1 /* File factor.c. */
2
3 #include "group.h"
4 #include "errmesg.h"
5
6 CHECK( factor)
7
8 extern Unsigned primeList[];
9
10
11 /*-------------------------- gcd ------------------------------------------*/
12
13 /* The function gcd( a, b) returns the greatest common divisor of log
14 integers a and b. */
15
gcd(unsigned long a,unsigned long b)16 unsigned long gcd(
17 unsigned long a,
18 unsigned long b)
19 {
20 unsigned long temp_a;
21
22 if (a < b)
23 {temp_a = a; a = b; b = temp_a;}
24 while (b != 0) {
25 temp_a = a;
26 a = b;
27 b = temp_a % b;
28 }
29 return a;
30 }
31
32
33 /*-------------------------- factMultiply ---------------------------------*/
34
35 /* The function factMultiply multiplies two factored integers. Specifically,
36 fmultiply( a, b) replaces a by the product a*b. (b remains unchanged.) */
37
factMultiply(FactoredInt * const a,FactoredInt * const b)38 void factMultiply(
39 FactoredInt *const a,
40 FactoredInt *const b) /* a = a * b */
41 {
42 Unsigned aIndex = 0, bIndex = 0, i;
43 a->prime[a->noOfFactors] = b->prime[b->noOfFactors] = MAX_INT;
44 while ( a->prime[aIndex] < MAX_INT || b->prime[bIndex] < MAX_INT )
45 if ( a->prime[aIndex] < b->prime[bIndex] )
46 ++aIndex;
47 else if ( a->prime[aIndex] == b->prime[bIndex] )
48 a->exponent[aIndex++] += b->exponent[bIndex++];
49 else
50 if ( a->noOfFactors < MAX_PRIME_FACTORS ) {
51 for ( i = ++a->noOfFactors ; i > aIndex ; --i ) {
52 a->prime[i] = a->prime[i-1];
53 a->exponent[i] = a->exponent[i-1];
54 }
55 a->prime[aIndex] = b->prime[bIndex];
56 a->exponent[aIndex++] = b->exponent[bIndex++];
57 }
58 else
59 ERROR1i( "fMultiply", "Number of prime factored exceeded bound "
60 "of ", MAX_PRIME_FACTORS, ".")
61 }
62
63
64
65 /*-------------------------- factDivide -----------------------------------*/
66
67 /* The function factDivide divides two factored integers. Specifically,
68 fmultiply( a, b) replaces a by the quotient a/b. (b remains unchanged.)
69 If b does not divide a, an error occurs. */
70
factDivide(FactoredInt * const a,FactoredInt * const b)71 void factDivide(
72 FactoredInt *const a,
73 FactoredInt *const b) /* a = a / b */
74 {
75 Unsigned aIndex = 0, bIndex = 0, i;
76 a->prime[a->noOfFactors] = b->prime[b->noOfFactors] = MAX_INT;
77 while ( a->prime[aIndex] < MAX_INT && b->prime[bIndex] < MAX_INT )
78 if ( a->prime[aIndex] < b->prime[bIndex] )
79 ++aIndex;
80 else if ( a->prime[aIndex] == b->prime[bIndex] )
81 if ( a->exponent[aIndex] > b->exponent[bIndex] )
82 a->exponent[aIndex++] -= b->exponent[bIndex++];
83 else if ( a->exponent[aIndex] == b->exponent[bIndex] ) {
84 --a->noOfFactors;
85 for ( i = aIndex ; i <= a->noOfFactors ; ++i ) {
86 a->prime[i] = a->prime[i+1];
87 a->exponent[i] = a->exponent[i+1];
88 }
89 ++bIndex;
90 }
91 else
92 ERROR( "factDivide", "Divisor does not divide dividend in factored "
93 "integer division.")
94 else
95 ERROR( "factDivide", "Divisor does not divide dividend in factored "
96 "integer division.")
97 }
98
99
100 /*-------------------------- factorize ------------------------------------*/
101
factorize(Unsigned n)102 FactoredInt factorize(
103 Unsigned n)
104 {
105 FactoredInt nFactored;
106 Unsigned i = 0, lastPrime = 0;
107
108 nFactored.noOfFactors = 0;
109 while ( primeList[i] * primeList[i] <= n && primeList[i] )
110 if ( n % primeList[i] == 0 ) {
111 if ( primeList[i] == lastPrime )
112 ++nFactored.exponent[nFactored.noOfFactors-1];
113 else if ( nFactored.noOfFactors < MAX_PRIME_FACTORS ) {
114 nFactored.prime[nFactored.noOfFactors] = primeList[i];
115 nFactored.exponent[nFactored.noOfFactors++] = 1;
116 lastPrime = primeList[i];
117 }
118 else
119 ERROR1i( "factorize", "Number of prime factors exceeded "
120 "maximum of ", MAX_PRIME_FACTORS, ".")
121 n /= primeList[i];
122 }
123 else
124 ++i;
125
126 if ( primeList[i] == 0 )
127 ERROR( "factorize", "Prime number list overflow.")
128 else if ( n > 1)
129 if ( n == lastPrime )
130 ++nFactored.exponent[nFactored.noOfFactors-1];
131 else if ( nFactored.noOfFactors < MAX_PRIME_FACTORS ) {
132 nFactored.prime[nFactored.noOfFactors] = n;
133 nFactored.exponent[nFactored.noOfFactors++] = 1;
134 }
135 else
136 ERROR1i( "factorize", "Number of prime factors exceeded "
137 "maximum of ", MAX_PRIME_FACTORS, ".")
138
139 return nFactored;
140 }
141
142
143 /*-------------------------- factEqual ------------------------------------*/
144
factEqual(FactoredInt * a,FactoredInt * b)145 BOOLEAN factEqual(
146 FactoredInt *a,
147 FactoredInt *b)
148 {
149 Unsigned i;
150
151 if ( a->noOfFactors != b->noOfFactors )
152 return FALSE;
153 for ( i = 0 ; i < a->noOfFactors ; ++i )
154 if ( a->prime[i] != b->prime[i] || a->exponent[i] != b->exponent[i] )
155 return FALSE;
156
157 return TRUE;
158 }
159