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.4 - A Program for Computing Primitive Polynomials.
21| Copyright (C) 1999-2025 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 seanerikoconnor!AT!gmail!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
51DESCRIPTION
52
53 A mod function for both positive and negative arguments.
54
55INPUT
56
57 n (int, -infinity < n < infinity)
58 p (int, p > 0)
59
60RETURNS
61
62 n mod p (0 <= n mod p < p)
63
64EXAMPLE
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
71METHOD
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
86BUGS
87
88 None.
89
90--------------------------------------------------------------------------------
91| Function Call |
92------------------------------------------------------------------------------*/
93
94int
95 mod( int n, int p )
96{
97
98/*------------------------------------------------------------------------------
99| Local Variables |
100------------------------------------------------------------------------------*/
101
102register int
103 raw_mod ; /* The value n % p computed by C's mod function. */
104
105/*------------------------------------------------------------------------------
106| Function Body |
107------------------------------------------------------------------------------*/
108
109raw_mod = n % p ;
110
111if (raw_mod == 0)
112
113 return( 0 ) ;
114
115else if (n >= 0) /* mod is not 0. n >= 0. */
116
117 return( raw_mod ) ;
118
119else
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
131DESCRIPTION
132
133 Raise a positive integer to a positive power exactly.
134
135INPUT
136
137 x (int, x > 0) The base.
138 y (int, y > 0) The exponent.
139
140RETURNS
141
142 y
143 x
144
145EXAMPLE CALLING SEQUENCE
146
147 power( 2, 3 ) => 8
148
149METHOD
150
151 Multiply x by itself y times.
152
153BUGS
154
155 None.
156
157--------------------------------------------------------------------------------
158| Function Call |
159------------------------------------------------------------------------------*/
160
161bigint power( int x, int y )
162{
163
164/*------------------------------------------------------------------------------
165| Local Variables |
166------------------------------------------------------------------------------*/
167
168bigint
169 pow ; /* x to the y. */
170
171/*------------------------------------------------------------------------------
172| Function Body |
173------------------------------------------------------------------------------*/
174
175for (pow = 1 ; y > 0 ; --y)
176
177 pow *= x ;
178
179return( pow ) ;
180
181} /* ========================= end of function power ======================== */
182
183
184
185/*==============================================================================
186| power_mod |
187================================================================================
188
189DESCRIPTION
190
191 Raise a positive integer to a positive power modulo an integer.
192
193INPUT
194
195 a (int, a >= 0) The base.
196 n (int, n >= 0) The exponent.
197 p (int, p >= 2) The modulus.
198
199RETURNS
200
201 n
202 a (modulo p)
203 -1 if a < 0, n < 0, or p <= 1, or 0^0 case.
204
205EXAMPLE CALLING SEQUENCE
206
207 power_mod( 2, 3, 7 ) => 1
208
209METHOD
210
211 Multiplication by repeated squaring.
212
213BUGS
214
215 None.
216
217--------------------------------------------------------------------------------
218| Function Call |
219------------------------------------------------------------------------------*/
220
221int 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
235int
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. */
241long int product = (int) a ;
242
243
244/*------------------------------------------------------------------------------
245| Function Body |
246------------------------------------------------------------------------------*/
247
248/*
249 Out of range conditions.
250*/
251
252if (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. */
257if (a == 0)
258 return 0 ;
259
260if (n == 0)
261 return 1 ;
262
263if (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
272while (! (n & mask))
273{
274 n <<= 1 ;
275 ++bit_count ;
276}
277
278bit_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
287while ( --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. */
301return (int) product ;
302
303} /* ========================= end of function power_mod ==================== */
304
305
306
307/*==============================================================================
308| is_primitive_root |
309================================================================================
310
311DESCRIPTION
312
313 Test whether the number a is a primitive root of the prime p.
314
315INPUT
316
317 a (int, a >= 1) The number to be tested.
318 p (int, p >= 2) A prime number.
319
320RETURNS
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
325EXAMPLE
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
335METHOD
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
344BUGS
345
346 We don't check if p is not a prime.
347
348--------------------------------------------------------------------------------
349| Function Call |
350------------------------------------------------------------------------------*/
351
352int is_primitive_root( int a, int p )
353{
354
355/*------------------------------------------------------------------------------
356| Local Variables |
357------------------------------------------------------------------------------*/
358
359bigint
360 primes[ MAXNUMPRIMEFACTORS ] ; /* The distinct prime factors of p-1. */
361
362int
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*/
374if (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*/
382if ( (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. */
392a = a % p ;
393
394
395/* a = 0 (mod p); Zero can't be a primitive root of p. */
396if (a == 0)
397 return NO ;
398
399/* Factor p-1 into primes. */
400t = 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 */
406for (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
414return YES ;
415
416} /* =============== end of function is_primitive_root ===================== */
417
418
419/* Divide modulo p.
420
421Do 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
425Now let v = p = prime number, so gcd = u3 = 1. Then we get
426
427 u u1 + p ? = 1
428
429or u u1 = 0 (mod p) which makes u (mod p) the unique multiplicative
430inverse of u.
431
432
433Assume 0 <= u < p, p is prime >= 2.
434
435*/
436int 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}
