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 }
