```  1 /*==============================================================================
2 |
3 |  File Name:
4 |
5 |     ppArith.c
6 |
7 |  Description:
8 |
9 |     Routines for integer arithmetic modulo p.
10 |
11 |  Functions:
12 |
13 |      mod
14 |      power
15 |      power_mod
16 |      is_primitive_root
17 |
18 |  LEGAL
19 |
20 |     Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.
22 |
23 |     This program is free software: you can redistribute it and/or modify
25 |     the Free Software Foundation, either version 3 of the License, or
26 |     (at your option) any later version.
27 |
28 |     This program is distributed in the hope that it will be useful,
29 |     but WITHOUT ANY WARRANTY; without even the implied warranty of
30 |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
31 |     GNU General Public License for more details.
32 |
33 |     You should have received a copy of the GNU General Public License
34 |     along with this program.  If not, see <http://www.gnu.org/licenses/>.
35 |
36 |     The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
37 |     with !DOT! replaced by . and the !AT! replaced by @
38 |
39 ==============================================================================*/
40
41 /*------------------------------------------------------------------------------
42 |                                Include Files                                 |
43 ------------------------------------------------------------------------------*/
44
45 #include "Primpoly.h"
46
47 /*==============================================================================
48 |                                  mod                                         |
49 ================================================================================
50
51 DESCRIPTION
52
53      A mod function for both positive and negative arguments.
54
55 INPUT
56
57     n (int, -infinity < n < infinity)
58     p (int, p > 0)
59
60 RETURNS
61
62     n mod p (0 <= n mod p < p)
63
64 EXAMPLE
65
66     The C language gives -5 % 3 = - ( -(-5) mod 3 ) = -( 5 mod 3 ) = -2.
67     The result is nonzero, so we add p=3 to give 1.
68
69     C computes -3 % 3 = - ( -(-3) mod 3  ) = 0, which we leave unchanged.
70
71 METHOD
72
73      For n >= 0, the C language defines r = n mod p by the equation
74
75          n = kp + r    0 <= r < p
76
77      but when n < 0, C returns the quantity
78
79          r = - ( (-n) mod p )
80
81      To put the result into the correct range 0 to p-1, add p to r if
82      r is non-zero.
83
84      By the way, FORTRAN's MOD function does the same thing.
85
86 BUGS
87
88      None.
89
90 --------------------------------------------------------------------------------
91 |                                Function Call                                 |
92 ------------------------------------------------------------------------------*/
93
94 int
95     mod( int n, int p )
96 {
97
98 /*------------------------------------------------------------------------------
99 |                               Local Variables                                |
100 ------------------------------------------------------------------------------*/
101
102 register int
103     raw_mod ;      /*  The value n % p computed by C's mod function.  */
104
105 /*------------------------------------------------------------------------------
106 |                                Function Body                                 |
107 ------------------------------------------------------------------------------*/
108
109 raw_mod = n % p ;
110
111 if (raw_mod == 0)
112
113     return( 0 ) ;
114
115 else if (n >= 0)             /*  mod is not 0.  n >= 0. */
116
117     return( raw_mod ) ;
118
119 else
120
121     return( raw_mod + p ) ;  /* mod is not 0.  n < 0. */
122
123 } /* =========================== end of function mod ======================== */
124
125
126
127 /*==============================================================================
128 |                                       power                                  |
129 ================================================================================
130
131 DESCRIPTION
132
133      Raise a positive integer to a positive power exactly.
134
135 INPUT
136
137      x (int, x > 0)   The base.
138      y (int, y > 0)   The exponent.
139
140 RETURNS
141
142      y
143     x
144
145 EXAMPLE CALLING SEQUENCE
146
147      power( 2, 3 ) => 8
148
149 METHOD
150
151      Multiply x by itself y times.
152
153 BUGS
154
155      None.
156
157 --------------------------------------------------------------------------------
158 |                                Function Call                                 |
159 ------------------------------------------------------------------------------*/
160
161 bigint power( int x, int y )
162 {
163
164 /*------------------------------------------------------------------------------
165 |                               Local Variables                                |
166 ------------------------------------------------------------------------------*/
167
168 bigint
169     pow ;   /*  x to the y.  */
170
171 /*------------------------------------------------------------------------------
172 |                                Function Body                                 |
173 ------------------------------------------------------------------------------*/
174
175 for (pow = 1 ;  y > 0 ;  --y)
176
177     pow *= x ;
178
179 return( pow ) ;
180
181 } /* ========================= end of function power ======================== */
182
183
184
185 /*==============================================================================
186 |                                power_mod                                     |
187 ================================================================================
188
189 DESCRIPTION
190
191      Raise a positive integer to a positive power modulo an integer.
192
193 INPUT
194
195      a (int, a >= 0)   The base.
196      n (int, n >= 0)   The exponent.
197      p (int, p >= 2)   The modulus.
198
199 RETURNS
200
201       n
202      a  (modulo p)
203     -1 if a < 0, n < 0, or p <= 1, or 0^0 case.
204
205 EXAMPLE CALLING SEQUENCE
206
207      power_mod( 2, 3, 7 ) => 1
208
209 METHOD
210
211      Multiplication by repeated squaring.
212
213 BUGS
214
215      None.
216
217 --------------------------------------------------------------------------------
218 |                                Function Call                                 |
219 ------------------------------------------------------------------------------*/
220
221 int power_mod( int a, int n, int p )
222 {
223
224 /*------------------------------------------------------------------------------
225 |                               Local Variables                                |
226 ------------------------------------------------------------------------------*/
227
228 /*
229     mask = 1000 ... 000.  That is, all bits of mask are zero except for the
230     most significant bit of the computer word holding its value.  Signed or
231     unsigned type for bigint shouldn't matter here, since we just mask and
232     compare bits.
233 */
234
235 int
236     mask = (1 << (8 * sizeof( int ) - 1)),
237     bit_count = 0 ;  /* Number of bits in exponent to the right of the leading bit. */
238
239 /*  Use extra precision since we need to square an integer within the
240     calculations below. */
241 long int product = (int) a ;
242
243
244 /*------------------------------------------------------------------------------
245 |                                Function Body                                 |
246 ------------------------------------------------------------------------------*/
247
248 /*
249     Out of range conditions.
250 */
251
252 if (a < 0 || n < 0 || p <= 1 || (a == 0 && n == 0))
253
254     return -1 ;
255
256 /*  Quick return for 0 ^ n, a^0 and a^1.  */
257 if (a == 0)
258     return 0 ;
259
260 if (n == 0)
261     return 1 ;
262
263 if (n == 1)
264     return a % p ;
265
266 /*
267     Advance the leading bit of the exponent up to the word's left
268     hand boundary.  Count how many bits were to the right of the
270 */
271
272 while (! (n & mask))
273 {
274     n <<= 1 ;
275     ++bit_count ;
276 }
277
278 bit_count = (8 * sizeof( int )) - bit_count ;
279
280
281 /*
283     Thereafter, square for every 0 bit;  square and multiply by x for
284     every 1 bit.
285 */
286
287 while ( --bit_count > 0 )
288 {
289     /*  Expose the next bit. */
290     n <<= 1 ;
291
292     /* Square modulo p. */
293     product = (product * product) % p ;
294
295     /*  Leading bit is 1: multiply by a modulo p. */
297         product = (a * product) % p ;
298 }
299
300 /*  We don't lose precision because product is always reduced modulo p. */
301 return (int) product ;
302
303 } /* ========================= end of function power_mod ==================== */
304
305
306
307 /*==============================================================================
308 |                             is_primitive_root                                |
309 ================================================================================
310
311 DESCRIPTION
312
313      Test whether the number a is a primitive root of the prime p.
314
315 INPUT
316
317      a      (int, a >= 1)     The number to be tested.
318      p      (int, p >= 2)     A prime number.
319
320 RETURNS
321
322      YES    If a is a primitive root of p.
323      NO     If a isn't a primitive root of p or if inputs are out of range.
324
325 EXAMPLE
326                                                      1      2      3
327      3 is a primitive root of the prime p = 7 since 3 = 3, 3 = 2, 3 = 6,
328       4      5                                       p-1    6
329      3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3    = 3 = 1 (mod 7).
330
331      So 3 is a primitive root of 7 because it has maximal period.  a = 2
332                                                       3
333      isn't a primitive root of p = 7 because already 2 = 1 (mod 7).
334
335 METHOD
336
337     From number theory, a is a primitive root of the prime p iff
338      (p-1)/q
339     a        != 1 (mod p) for all prime divisors q of (p-1).
340                             (p-1)
341     Don't need to check if a     = 1 (mod p):  Fermat's little
342     theorem guarantees it.
343
344 BUGS
345
346     We don't check if p is not a prime.
347
348 --------------------------------------------------------------------------------
349 |                                Function Call                                 |
350 ------------------------------------------------------------------------------*/
351
352 int is_primitive_root( int a, int p )
353 {
354
355 /*------------------------------------------------------------------------------
356 |                               Local Variables                                |
357 ------------------------------------------------------------------------------*/
358
359 bigint
360     primes[ MAXNUMPRIMEFACTORS ] ; /* The distinct prime factors of p-1. */
361
362 int
363     count[ MAXNUMPRIMEFACTORS ],   /* ... and their multiplicities.      */
364     t,                             /*  Primes are indexed from 0 to t.   */
365     i ;                            /*  Loop index.                       */
366
367 /*------------------------------------------------------------------------------
368 |                                Function Body                                 |
369 ------------------------------------------------------------------------------*/
370
371 /*  Error for out of range inputs, including p
372     being an even number greater than 2.
373 */
374 if (p < 2 || a < 1 || (p > 2 && (p % 2 == 0)))
375     return NO ;
376
377 /*  Special cases:
378     1 is the only primitive root of 2;  i.e. 1 generates the unit elements
379     of GF( 2 );  2 is the only primitive root of 3, and 2 and 3 are the only
380     primitive roots of 5, etc.  See http://mathworld.wolfram.com/PrimitiveRoot.html
381 */
382 if ( (p == 2  && a == 1) ||
383      (p == 3  && a == 2) ||
384      (p == 5  && (a == 2 || a == 3)) ||
385      (p == 7  && (a == 3 || a == 5)) ||
386      (p == 11 && (a == 2 || a == 6 || a == 7 || a ==  8)) ||
387      (p == 13 && (a == 2 || a == 6 || a == 7 || a == 11)) )
388     return YES ;
389
390
391 /*  Reduce a down modulo p. */
392 a = a % p ;
393
394
395 /*  a = 0 (mod p);  Zero can't be a primitive root of p. */
396 if (a == 0)
397     return NO ;
398
399 /*  Factor p-1 into primes. */
400 t = factor( (bigint)(p - 1), primes, count ) ;
401
402
403 /*  Check if a^(p-1)/q <> 1 (mod p) for all primes q | (p-1).
404     If so, we have a primitive root of p, otherwise, we exit early.
405  */
406 for (i = 0 ;  i <= t ;  ++i)
407 {
408     /*  The prime factors can be safely cast from bigint to int
409         because p-1 needs only integer precision.  */
410     if (power_mod( a, (p-1) / (int)primes[ i ], p ) == 1)
411         return NO ;
412 }
413
414 return YES ;
415
416 } /* =============== end of function is_primitive_root ===================== */
417
418
419 /* Divide modulo p.
420
421 Do extended Euclid's algorithm on u and v to find u1, u2, and u3 such that
422
423     u u1 + v u2 = u3 = gcd( u, v ).
424
425 Now let v = p = prime number, so gcd = u3 = 1.  Then we get
426
427   u u1 + p ? = 1
428
429 or u u1 = 0 (mod p) which makes u (mod p) the unique multiplicative
430 inverse of u.
431
432
433 Assume 0 <= u < p, p is prime >= 2.
434
435 */
436 int inverse_mod_p( int u, int p )
437 {
438     /*  Do Euclid's algorithm to find the quantities */
439     int t1 = 0 ;
440     int t3 = 0 ;
441     int q  = 0 ;
442
443     int u1 = 1 ;
444     int u3 = u ;
445     int v1 = 0 ;
446     int v3 = p ;
447
448     int inv_v = 0 ;
449
450     while( v3 != 0)
451     {
452         q = (int)(u3 / v3) ;
453
454         t1 = u1 - v1 * q ;
455         t3 = u3 - v3 * q ;
456
457         u1 = v1 ;
458         u3 = v3 ;
459
460         v1 = t1 ;
461         v3 = t3 ;
462     }
463
464     inv_v = mod( u1, p ) ;
465
466     /* Self check. */
467     if ( mod( u * inv_v, p ) != 1)
468     {
469         return 0 ;
470     }
471
472     return inv_v ;
473 }
```