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.
 21 |     Copyright (C) 1999-2021 by Sean Erik O'Connor.  All Rights Reserved.
 22 |
 23 |     This program is free software: you can redistribute it and/or modify
 24 |     it under the terms of the GNU General Public License as published by
 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
269     leading bit.
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 /*
282     Exponentiation by repeated squaring.  Discard the leading 1 bit.
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. */
296     if (n & mask)
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 }