1 /*============================================================================== 2 | 3 | File Name: 4 | 5 | ppPolyArith.c 6 | 7 | Description: 8 | 9 | Polynomial arithmetic and exponentiation routines. 10 | 11 | Functions: 12 | 13 | eval_poly 14 | linear_factor 15 | is_integer 16 | construct_power_table 17 | auto_convolve 18 | convolve 19 | coeff_of_square 20 | coeff_of_product 21 | square 22 | product 23 | times_x 24 | x_to_power 25 | 26 | LEGAL 27 | 28 | Primpoly Version 16.1 - A Program for Computing Primitive Polynomials. 29 | Copyright (C) 1999-2021 by Sean Erik O'Connor. All Rights Reserved. 30 | 31 | This program is free software: you can redistribute it and/or modify 32 | it under the terms of the GNU General Public License as published by 33 | the Free Software Foundation, either version 3 of the License, or 34 | (at your option) any later version. 35 | 36 | This program is distributed in the hope that it will be useful, 37 | but WITHOUT ANY WARRANTY; without even the implied warranty of 38 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39 | GNU General Public License for more details. 40 | 41 | You should have received a copy of the GNU General Public License 42 | along with this program. If not, see <http://www.gnu.org/licenses/>. 43 | 44 | The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com 45 | with !DOT! replaced by . and the !AT! replaced by @ 46 | 47 ==============================================================================*/ 48 49 /*------------------------------------------------------------------------------ 50 | Include Files | 51 ------------------------------------------------------------------------------*/ 52 53 #include <stdio.h> 54 #include "Primpoly.h" 55 56 57 /*============================================================================== 58 | eval_poly | 59 ================================================================================ 60 61 DESCRIPTION 62 63 Evaluate the polynomial f( x ) with modulo p arithmetic. 64 65 INPUT 66 67 f (int *) nth degree mod p polynomial 68 69 n n-1 70 f( x ) = x + a x + ... + a 0 <= a < p 71 n-1 0 i 72 x (int) 0 <= x < p 73 74 n (int) n >= 1 75 76 p (int) 77 78 RETURNS 79 80 f( x ) (int) 81 82 EXAMPLE 83 4 84 Let n = 4, p = 5 and f(x) = x + 3x + 3. Then f(x) = 85 (((x + 0)x + 0)x + 3)x + 3. f(2) = (((2 + 0)2 + 0)2 + 3) = 86 (8 + 3)2 + 3 = 1 + 2 + 3 (mod 5) = 0. 87 88 METHOD 89 90 By Horner's rule, f(x) = ( ... ( (x + a )x + ... )x + a . 91 n-1 0 92 We evaluate recursively, f := f * x + a (mod p), starting 93 i 94 with f = 1 and i = n-1. 95 96 BUGS 97 98 None. 99 100 -------------------------------------------------------------------------------- 101 | Function Call | 102 ------------------------------------------------------------------------------*/ 103 104 int 105 eval_poly( int * f, int x, int n, int p ) 106 { 107 108 /*------------------------------------------------------------------------------ 109 | Local Variables | 110 ------------------------------------------------------------------------------*/ 111 112 int 113 val = 1 ; /* The value of f(x) which is returned. */ 114 115 /*------------------------------------------------------------------------------ 116 | Function Body | 117 ------------------------------------------------------------------------------*/ 118 119 while ( --n >= 0 ) 120 121 val = mod( val * x + f[ n ], p ) ; 122 123 return( val ) ; 124 125 126 } /* ========================= end of function eval_poly ==================== */ 127 128 129 /*============================================================================== 130 | linear_factor | 131 ================================================================================ 132 133 DESCRIPTION 134 135 Check if the polynomial f(x) has any linear factors. 136 137 INPUT 138 139 f (int *) nth degree monic mod p polynomial f(x). 140 141 n (int) Its degree. 142 143 p (int) Modulus for coefficient arithmetic. 144 145 RETURNS 146 147 YES if f( a ) = 0 (mod p) for a = 1, 2, ... p-1. 148 NO otherwise 149 150 i.e. check if f(x) contains a linear factor (x - a). We don't need to test 151 for the root a = 0 because f(x) was chosen in main to have a non-zero 152 constant term, hence no zero root. 153 154 EXAMPLE 155 4 156 Let n = 4, p = 5 and f(x) = x + 3x + 3. Then f(1) = 2 (mod 5), but 157 f(2) = 0 (mod 5), so we exit immediately with a yes. 158 159 METHOD 160 161 Evaluate f(x) at x = 1, ..., p-1 by Horner's rule. Return instantly the 162 moment f(x) evaluates to 0. 163 164 BUGS 165 166 None. 167 168 -------------------------------------------------------------------------------- 169 | Function Call | 170 ------------------------------------------------------------------------------*/ 171 172 int 173 linear_factor( int * f, int n, int p ) 174 { 175 176 /*------------------------------------------------------------------------------ 177 | Local Variables | 178 ------------------------------------------------------------------------------*/ 179 180 int i ; /* Loop counter. */ 181 182 /*------------------------------------------------------------------------------ 183 | Function Body | 184 ------------------------------------------------------------------------------*/ 185 186 for (i = 1 ; i <= p-1 ; ++i) 187 188 if (eval_poly( f, i, n, p ) == 0) 189 190 return( YES ) ; 191 192 return( NO ) ; 193 194 } /* ====================== end of function linear_factor =================== */ 195 196 197 /*============================================================================== 198 | is_integer | 199 ================================================================================ 200 201 DESCRIPTION 202 203 Test if a polynomial is a constant. 204 205 INPUT 206 207 t (int *) mod p polynomial t(x). 208 n (int) Its degree. 209 210 RETURNS 211 212 YES if t(x) is a constant polynomial. 213 NO otherwise 214 215 EXAMPLE 216 2 217 Let n = 3. Then t(x) = 2 x is not constant but t(x) = 2 is. 218 219 METHOD 220 221 A constant polynomial is zero in its first through n th degree 222 terms. Return immediately with NO if any such term is non-zero. 223 224 BUGS 225 226 None. 227 228 -------------------------------------------------------------------------------- 229 | Function Call | 230 ------------------------------------------------------------------------------*/ 231 232 int 233 is_integer( int * t, int n ) 234 { 235 236 /*------------------------------------------------------------------------------ 237 | Function Body | 238 ------------------------------------------------------------------------------*/ 239 240 while ( n >= 1 ) 241 242 if (t[ n-- ] != 0) 243 244 return( NO ) ; 245 246 return( YES ) ; 247 248 } /* ====================== end of function is_integer ====================== */ 249 250 251 /*============================================================================== 252 | construct_power_table | 253 ================================================================================ 254 255 DESCRIPTION 256 257 Construct a table of powers of x: 258 259 n 2n-2 260 x (mod f(x), p) ... x (mod f(x), p) 261 262 INPUT 263 264 f (int *) Coefficients of f(x), a monic polynomial of degree n. 265 n (int, -infinity < n < infinity) 266 p (int, p > 0) 267 268 RETURNS 269 270 power_table (int *) power_table[i][j] is the coefficient of 271 j n+i 272 x in x (mod f(x), p) where 0 <= i <= n-2 and 0 <= j <= n-1. 273 274 EXAMPLE 275 4 2 4 276 Let n = 4, p = 5 and f(x) = x + x + 2x + 3. Then x = 277 278 2 2 279 -( x + 2x + 3) = 4 x + 3 x + 2 (mod f(x), 5), and we get 280 281 4 2 282 x (mod f(x), 5) = 4 x + 3 x + 2 = power_table[0]. 283 284 5 2 3 2 285 x (mod f(x), 5) = x( 4 x + 3 x + 2) = 4 x + 3 x + 2x 286 = power_table[1]. 287 288 6 3 2 4 3 2 289 x (mod f(x), 5) = x( 4 x + 3 x + 2 x) = 4x + 3 x + 2 x 290 291 2 3 2 292 = 4 ( 4x + 3 x + 2) + 3 x + 2 x = 293 294 3 2 295 = 3 x + 3 x + 2 x + 3 = power_table[2]. 296 297 j 298 power_table[i][j]: | 0 1 2 3 299 ---+------------- 300 0 | 2 3 4 0 301 i 1 | 0 2 3 4 302 2 | 3 2 3 3 303 304 METHOD 305 n-1 306 Beginning with t( x ) = x, compute the next power of x from the last 307 n 308 one by multiplying by x. If necessary, substitute x = 309 n-1 310 -( a x + ... + a ) to reduce the degree. This formula comes from 311 n-1 0 312 n n-1 313 the observation f( x ) = x + a x + ... + a = 0 (mod f(x), p). 314 n-1 0 315 BUGS 316 317 None. 318 319 -------------------------------------------------------------------------------- 320 | Function Call | 321 ------------------------------------------------------------------------------*/ 322 323 void 324 construct_power_table( int power_table[][ MAXDEGPOLY ], int * f, 325 int n, int p ) 326 { 327 328 /*------------------------------------------------------------------------------ 329 | Local Variables | 330 ------------------------------------------------------------------------------*/ 331 332 int 333 i, j, /* Loop counters. */ 334 coeff, /* Coefficient of x ^ n in t(x) */ 335 t[ MAXDEGPOLY + 1 ] ; /* t(x) is temporary storage for x ^ k (mod f(x),p) 336 n <= k <= 2n-2. Its degree can go as high as 337 n before it is reduced again. */ 338 339 /*------------------------------------------------------------------------------ 340 | Function Body | 341 ------------------------------------------------------------------------------*/ 342 343 /* 344 n-1 345 Initialize t( x ) = x 346 */ 347 348 for (i = 0 ; i <= n-2 ; ++i) 349 350 t[ i ] = 0 ; 351 352 t[ n-1 ] = 1 ; 353 354 355 356 /* 357 i+n 358 Fill the ith row of the table with x (mod f(x), p). 359 */ 360 361 for (i = 0 ; i <= n-2 ; ++i) 362 { 363 364 /* Compute t(x) = x t(x) */ 365 for (j = n ; j >= 1 ; --j) 366 367 t[ j ] = t[ j-1 ] ; 368 369 t[ 0 ] = 0 ; 370 371 372 /* Coefficient of the nth degree term of t(x). */ 373 if ( (coeff = t[ n ]) != 0) 374 { 375 t[ n ] = 0 ; /* Zero out the x ^ n th term. */ 376 377 /* 378 n n n-1 379 Replace x with x (mod f(x), p) = -(a x + ... + a ) 380 n-1 0 381 */ 382 383 for (j = 0 ; j <= n-1 ; ++j) 384 385 t[ j ] = mod( t[ j ] + 386 mod( -coeff * f[ j ], p ), 387 p ) ; 388 } /* end if */ 389 390 391 /* Copy t(x) into row i of power_table. */ 392 for (j = 0 ; j <= n-1 ; ++j) 393 394 power_table[i][j] = t[ j ] ; 395 396 } /* end for */ 397 398 return ; 399 400 } /* ================== end of function construct_power_table =============== */ 401 402 403 404 /*============================================================================== 405 | autoconvolve | 406 ================================================================================ 407 408 DESCRIPTION 409 410 Compute a convolution-like sum for use in function coeff_of_square. 411 412 INPUT 413 n-1 414 t (int *) Coefficients of t(x) = t x + ... + t x + t 415 n-1 1 0 416 k (int), 0 <= k <= 2n - 2) 417 418 lower (int, 0 <= lower <= n-1) 419 420 uppper (int, 0 <= upper <= n-1) 421 422 p (int, p > 0) 423 424 RETURNS 425 426 upper 427 --- 428 \ t t But define the sum to be 0 when lower > upper to catch 429 / i k-i the special cases that come up in function coeff_of_square. 430 --- 431 i=lower 432 433 EXAMPLE 434 3 2 435 Suppose t(x) = 4 x + x + 3 x + 3, lower = 1, upper = 3, n = 3, 436 437 and p = 5. For k = 3, auto_convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] + 438 439 t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3. For lower = 0, 440 441 upper = -1, or for lower = 3 and upper = 2, auto_convolve = 0, no matter what 442 443 k is. 444 445 METHOD 446 447 A "for" loop in the C language is not executed when its lower limit 448 449 exceeds its upper limit, leaving sum = 0. 450 451 BUGS 452 453 None. 454 455 -------------------------------------------------------------------------------- 456 | Function Call | 457 ------------------------------------------------------------------------------*/ 458 459 int 460 auto_convolve( int * t, int k, int lower, int upper, int p ) 461 { 462 463 /*------------------------------------------------------------------------------ 464 | Local Variables | 465 ------------------------------------------------------------------------------*/ 466 467 int 468 sum = 0, /* Convolution sum. */ 469 i ; /* Loop counter. */ 470 471 /*------------------------------------------------------------------------------ 472 | Function Body | 473 ------------------------------------------------------------------------------*/ 474 475 for (i = lower ; i <= upper ; ++i) 476 477 sum = mod( sum + mod( t[ i ] * t[ k - i ], p ), p ) ; 478 479 return( sum ) ; 480 481 } /* ==================== end of function auto_convolve ===================== */ 482 483 484 /*============================================================================== 485 | convolve | 486 ================================================================================ 487 488 DESCRIPTION 489 490 Compute a convolution-like sum. 491 492 INPUT 493 n-1 494 s (int *) Coefficients of s(x) = s x + ... + s x + s 495 n-1 1 0 496 n-1 497 t (int *) Coefficients of t(x) = t x + ... + t x + t 498 n-1 1 0 499 500 k (int), 0 <= k <= 2n - 2) 501 502 lower (int, 0 <= lower <= n-1) 503 504 uppper (int, 0 <= upper <= n-1) 505 506 p (int, p > 0) 507 508 RETURNS 509 510 upper 511 --- 512 \ s t But define the sum to be 0 when lower > upper to catch 513 / i k-i the special cases 514 --- 515 i=lower 516 517 EXAMPLE 518 3 2 519 Suppose s(x) = 4 x + x + 3 x + 3, 520 521 3 2 522 Suppose t(x) = 4 x + x + 3 x + 3, 523 524 525 lower = 1, upper = 3, n = 3, 526 527 and p = 5. For k = 3, convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] + 528 529 t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3. For lower = 0, 530 531 upper = -1, or for lower = 3 and upper = 2, convolve = 0, no matter what 532 533 k is. 534 535 METHOD 536 537 A "for" loop in the C language is not executed when its lower limit 538 539 exceeds its upper limit, leaving sum = 0. 540 541 BUGS 542 543 None. 544 545 -------------------------------------------------------------------------------- 546 | Function Call | 547 ------------------------------------------------------------------------------*/ 548 549 int 550 convolve( int * s, int * t, int k, int lower, int upper, int p ) 551 { 552 553 /*------------------------------------------------------------------------------ 554 | Local Variables | 555 ------------------------------------------------------------------------------*/ 556 557 int 558 sum = 0, /* Convolution sum. */ 559 i ; /* Loop counter. */ 560 561 /*------------------------------------------------------------------------------ 562 | Function Body | 563 ------------------------------------------------------------------------------*/ 564 565 for (i = lower ; i <= upper ; ++i) 566 567 sum = mod( sum + mod( s[ i ] * t[ k - i ], p ), p ) ; 568 569 return( sum ) ; 570 571 } /* ======================= end of function convolve ======================= */ 572 573 574 /*============================================================================== 575 | coeff_of_square | 576 ================================================================================ 577 578 DESCRIPTION 579 th 2 580 Return the coefficient of the k power of x in t( x ) modulo p. 581 582 INPUT 583 584 t (int *) Coefficients of t(x), of degree <= n-1. 585 586 k (int, 0 <= k <= 2n-2) 587 588 n (int, -infinity < n < infinity) 589 590 p (int, p > 0) 591 592 RETURNS 593 th 2 594 coefficient of the k power of x in t(x) mod p. 595 596 EXAMPLE 597 3 2 2 598 Let n = 4, p = 5, and t(x) = 4 x + x + 4. t(x) = 599 600 0 1 2 601 (4 * 4) x + (2 * 1 * 4) x + (2 * 1 * 4 + 1 * 1) x + 602 603 3 4 604 (2 * 1 * 1 + 2 * 4 * 4) x + (2 * 4 * 1 + 1 * 1) x + 605 606 5 6 607 (2 * 4 * 1) x + (4 * 4) x = 608 609 6 5 4 3 2 610 x + 3 x + 4 x + 4 x + 4 x + 3 x + 1, all modulo 5. 611 612 k | 0 1 2 3 4 5 6 613 ----------------+--------------------- 614 coeff_of_square | 1 3 4 4 4 3 1 615 616 METHOD 617 2 618 The formulas were gotten by writing out the product t (x) explicitly. 619 620 The sum is 0 in two cases: 621 622 (1) when k = 0 and the limits of summation are 0 to -1 623 624 (2) k = 2n - 2, when the limits of summation are n to n-1. 625 626 To derive the formulas, let 627 628 n-1 629 Let t (x) = t x + ... + t x + t 630 n-1 1 0 631 632 Look at the formulas in coeff_of_product for each power of x, 633 replacing s with t, and observe that half of the terms are 634 duplicates, so we can save computation time. 635 636 Inspection yields the formulas, 637 638 639 for 0 <= k <= n-1, even k, 640 641 k/2-1 642 --- 2 643 2 \ t t + t 644 / i k-i k/2 645 --- 646 i=0 647 648 for 0 <= k <= n-1, odd k, 649 650 (k-1)/2 651 --- 652 2 \ t t 653 / i k-i 654 --- 655 i=0 656 657 658 and for n <= k <= 2n-2, even k, 659 660 n-1 661 --- 2 662 2 \ t t + t 663 / i k-i k/2 664 --- 665 i=k/2+1 666 667 and for n <= k <= 2n-2, odd k, 668 669 n-1 670 --- 671 2 \ t t 672 / i k-i 673 --- 674 i=(k+1)/2 675 676 677 678 BUGS 679 680 None. 681 682 -------------------------------------------------------------------------------- 683 | Function Call | 684 ------------------------------------------------------------------------------*/ 685 686 int 687 coeff_of_square( int * t, int k, int n, int p ) 688 { 689 690 /*------------------------------------------------------------------------------ 691 | Local Variables | 692 ------------------------------------------------------------------------------*/ 693 694 int 695 sum = 0 ; /* kth coefficient of t(x) ^ 2. */ 696 697 /*------------------------------------------------------------------------------ 698 | Function Body | 699 ------------------------------------------------------------------------------*/ 700 701 if (0 <= k && k <= n-1) 702 { 703 704 if (k % 2 == 0) /* Even k */ 705 706 sum = mod( mod( 2 * auto_convolve( t, k, 0, k/2 - 1, p ), p) + 707 t[ k/2 ] * t[ k/2 ], p ) ; 708 709 else /* Odd k */ 710 711 sum = mod( 2 * auto_convolve( t, k, 0, (k-1)/2, p ), p) ; 712 } 713 else if (n <= k && k <= 2 * n - 2) 714 { 715 716 if (k % 2 == 0) /* Even k */ 717 718 sum = mod( mod( 2 * auto_convolve( t, k, k/2 + 1, n-1, p ), p) + 719 t[ k/2 ] * t[ k/2 ], p ) ; 720 721 else /* Odd k */ 722 723 sum = mod( 2 * auto_convolve( t, k, (k+1)/2, n-1, p ), p) ; 724 } 725 726 return( sum ) ; 727 728 } /* ================== end of function coeff_of_square ===================== */ 729 730 731 /*============================================================================== 732 | coeff_of_product | 733 ================================================================================ 734 735 DESCRIPTION 736 th 737 Return the coefficient of the k power of x in s( x ) t( x ) modulo p. 738 739 INPUT 740 741 t (int *) Coefficients of t(x), of degree <= n-1. 742 743 k (int, 0 <= k <= 2n-2) 744 745 n (int, -infinity < n < infinity) 746 747 p (int, p > 0) 748 749 RETURNS 750 th 751 coefficient of the k power of x in s(x) t(x) mod p. 752 753 EXAMPLE 754 3 2 2 755 Let n = 4, p = 5, t(x) = 4 x + x + 4, s( x ) = 3 x + x + 2 756 757 We'll do the case k=3, 758 759 t3 s0 + t2 s1 + t1 s2 + t0 s3 = 4 * 2 + 1 * 1 + 0 * 3 + 4 * 0 = 9 = 4 (mod 5). 760 761 762 k | 0 1 2 3 4 5 6 763 -----------------+--------------------- 764 coeff_of_product | 3 4 4 4 2 2 0 765 766 METHOD 767 768 The formulas were gotten by writing out the product s(x) t (x) explicitly. 769 770 The sum is 0 in two cases: 771 772 (1) when k = 0 and the limits of summation are 0 to -1 773 774 (2) k = 2n - 2, when the limits of summation are n to n-1. 775 776 777 To derive the formulas, let 778 779 n-1 780 Let s (x) = s x + ... + s x + s and 781 n-1 1 0 782 783 n-1 784 t (x) = t x + ... + t x + t 785 n-1 1 0 786 787 and multiply out the terms, collecting like powers of x: 788 789 790 Power of x Coefficient 791 ========================== 792 2n-2 793 x s t 794 n-1 n-1 795 796 2n-3 797 x s t + s t 798 n-2 n-1 n-1 n-2 799 800 2n-4 801 x s t + s t + s t 802 n-3 n-1 n-2 n-2 n-3 n-1 803 804 2n-5 805 x s t + s t + s t + s t 806 n-4 n-1 n-3 n-2 n-2 n-3 n-1 n-4 807 808 . . . 809 810 n 811 x s t + s t + ... + s t 812 1 n-1 2 n-2 n-1 1 813 814 n-1 815 x s t + s t + ... + s t 816 0 n-1 1 n-2 n-1 0 817 818 . . . 819 820 3 821 x s t + s t + s t + s t 822 0 3 1 2 2 1 3 0 823 824 2 825 x s t + s t + s t 826 0 2 1 1 2 0 827 828 829 x s t + s t 830 0 1 1 0 831 832 1 s t 833 0 0 834 835 836 Inspection yields the formulas, 837 838 839 for 0 <= k <= n-1, 840 841 k 842 --- 843 \ s t 844 / i k-i 845 --- 846 i=0 847 848 849 and for n <= k <= 2n-2, 850 851 n-1 852 --- 853 \ s t 854 / i k-i 855 --- 856 i=k-n+1 857 858 BUGS 859 860 None. 861 862 -------------------------------------------------------------------------------- 863 | Function Call | 864 ------------------------------------------------------------------------------*/ 865 866 int 867 coeff_of_product( int * s, int * t, int k, int n, int p ) 868 { 869 870 /*------------------------------------------------------------------------------ 871 | Local Variables | 872 ------------------------------------------------------------------------------*/ 873 874 int 875 sum = 0 ; /* kth coefficient of t(x) ^ 2. */ 876 877 /*------------------------------------------------------------------------------ 878 | Function Body | 879 ------------------------------------------------------------------------------*/ 880 881 if (0 <= k && k <= n-1) 882 { 883 sum = convolve( s, t, k, 0, k, p ) ; 884 } 885 else if (n <= k && k <= 2 * n - 2) 886 { 887 sum = convolve( s, t, k, k - n + 1, n - 1, p ) ; 888 } 889 890 return( sum ) ; 891 892 } /* ================== end of function coeff_of_product ==================== */ 893 894 895 /*============================================================================== 896 | square | 897 ================================================================================ 898 899 DESCRIPTION 900 2 901 Compute t ( x ) (mod f(x), p) 902 903 Uses a precomputed table of powers of x. 904 905 INPUT 906 907 t (int *) Coefficients of t (x), of degree <= n-1. 908 909 power_table (int **) x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic. 910 911 n (int, -infinity < n < infinity) Degree of f(x). 912 913 p (int, p > 0) Mod p coefficient arithmetic. 914 915 OUTPUT 916 2 917 t (int *) Overwritten with t (x) (mod f(x), p) 918 919 EXAMPLE 920 3 2 921 Let n = 4, p = 5, and t(x) = 4 x + x + 4. Let f(x) = 922 923 4 2 2 6 5 4 3 2 924 x + x + 2 x + 3. Then t (x) = x + 3 x + 4 x + 4 x + 4 x + 3 x 925 926 2 3 2 927 + 1. t (x) (mod f(x), p) = 1 * (3 x + 3 x + 2 x + 3) + 928 929 3 2 2 3 2 930 3 * (4 x + 3 x + 2 x) + 4 * (4 x + 3 x + 2) + 4 x + 4 x + 3 x + 1 931 932 3 2 933 = 4 x + 2 x + 3 x + 2. 934 935 936 METHOD 937 2 2n-2 n n-1 938 Let t (x) = t x + ... + t x + t x + ... + t . 939 2n-2 n n-1 0 940 941 Compute the coefficients t using the function coeff_of_square. 942 k 943 944 2 945 The next step is to reduce t (x) modulo f(x). To do so, replace 946 947 k k 948 each non-zero term t x, n <= k <= 2n-2, by the term t * [ x mod f(x), p)] 949 k k 950 951 which we get from the array power_table. 952 953 BUGS 954 955 None. 956 957 -------------------------------------------------------------------------------- 958 | Function Call | 959 ------------------------------------------------------------------------------*/ 960 961 void 962 square( int * t, int power_table[][ MAXDEGPOLY ], int n, int p ) 963 { 964 965 /*------------------------------------------------------------------------------ 966 | Local Variables | 967 ------------------------------------------------------------------------------*/ 968 969 int 970 i, j, /* Loop counters. */ 971 coeff, /* Coefficient of x ^ k term of t(x) ^2 */ 972 temp[ MAXDEGPOLY + 1 ] ; /* Temporary storage for the new t(x). */ 973 974 /*------------------------------------------------------------------------------ 975 | Function Body | 976 ------------------------------------------------------------------------------*/ 977 978 /* 979 0 n-1 980 Compute the coefficients of x , ..., x. These terms do not require 981 982 reduction mod f(x) because their degree is less than n. 983 */ 984 985 for (i = 0 ; i <= n ; ++i) 986 987 temp[ i ] = coeff_of_square( t, i, n, p ) ; 988 989 /* 990 n 2n-2 k 991 Compute the coefficients of x , ..., x. Replace t x with 992 k k 993 t * [ x (mod f(x), p) ] from array power_table when t is 994 k k 995 non-zero. 996 997 */ 998 999 for (i = n ; i <= 2 * n - 2 ; ++i) 1000 1001 if ( (coeff = coeff_of_square( t, i, n, p )) != 0 ) 1002 1003 for (j = 0 ; j <= n - 1 ; ++j) 1004 1005 temp[ j ] = mod( temp[ j ] + 1006 mod( coeff * power_table[ i - n ] [ j ], p ), 1007 p ) ; 1008 1009 for (i = 0 ; i <= n - 1 ; ++i) 1010 1011 t[ i ] = temp[ i ] ; 1012 1013 } /* ======================== end of function square ======================== */ 1014 1015 1016 1017 /*============================================================================== 1018 | product | 1019 ================================================================================ 1020 1021 DESCRIPTION 1022 1023 Compute s( x ) t( x ) (mod f(x), p) 1024 1025 Uses a precomputed table of powers of x. 1026 1027 INPUT 1028 1029 s (int *) Coefficients of s(x), of degree <= n-1. 1030 1031 t (int *) Coefficients of t(x), of degree <= n-1. 1032 1033 power_table (int **) x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic. 1034 1035 n (int, -infinity < n < infinity) Degree of f(x). 1036 1037 p (int, p > 0) Mod p coefficient arithmetic. 1038 1039 OUTPUT 1040 1041 s (int *) Overwritten with s( x ) t( x ) (mod f(x), p) 1042 1043 EXAMPLE 1044 3 2 2 1045 Let n = 4 and p = 5, t( x ) = 4 x + x + 4, s( x ) = 3 x + x + 2 1046 1047 5 4 3 2 1048 Then s( x ) t( x ) = 2 x + 2 x + 4 x + 4 x + 4 x + 3, modulo 5, 1049 1050 4 2 1051 and after reduction modulo f( x ) = x + x + 2 x + 3, using the power 1052 1053 4 2 5 3 2 1054 table entries for x = 4 x + 3 x + 2 and x = 4 x + 3 x + 2 x, we get 1055 1056 3 2 1057 s( x ) t( x ) (mod f( x ), p) = 2 x + 3 x + 4 x + 2 1058 1059 1060 METHOD 1061 1062 Compute the coefficients using the function coeff_of_product. 1063 1064 The next step is to reduce s(x) t(x) modulo f(x) and p. To do so, replace 1065 1066 k k 1067 each non-zero term t x, n <= k <= 2n-2, by the term t * [ x mod f(x), p)] 1068 k k 1069 1070 which we get from the array power_table. 1071 1072 BUGS 1073 1074 None. 1075 1076 -------------------------------------------------------------------------------- 1077 | Function Call | 1078 ------------------------------------------------------------------------------*/ 1079 1080 void 1081 product( int * s, int * t, int power_table[][ MAXDEGPOLY ], int n, int p ) 1082 { 1083 1084 /*------------------------------------------------------------------------------ 1085 | Local Variables | 1086 ------------------------------------------------------------------------------*/ 1087 1088 int 1089 i, j, /* Loop counters. */ 1090 coeff, /* Coefficient of x ^ k term of t(x) ^2 */ 1091 temp[ MAXDEGPOLY + 1 ] ; /* Temporary storage for the new t(x). */ 1092 1093 /*------------------------------------------------------------------------------ 1094 | Function Body | 1095 ------------------------------------------------------------------------------*/ 1096 1097 /* 1098 0 n-1 1099 Compute the coefficients of x , ..., x. These terms do not require 1100 1101 reduction mod f(x) because their degree is less than n. 1102 */ 1103 1104 for (i = 0 ; i <= n ; ++i) 1105 1106 temp[ i ] = coeff_of_product( s, t, i, n, p ) ; 1107 1108 /* 1109 n 2n-2 k 1110 Compute the coefficients of x , ..., x. Replace t x with 1111 k k 1112 t * [ x (mod f(x), p) ] from array power_table when t is 1113 k k 1114 non-zero. 1115 1116 */ 1117 1118 for (i = n ; i <= 2 * n - 2 ; ++i) 1119 1120 if ( (coeff = coeff_of_product( s, t, i, n, p )) != 0 ) 1121 1122 for (j = 0 ; j <= n - 1 ; ++j) 1123 1124 temp[ j ] = mod( temp[ j ] + 1125 mod( coeff * power_table[ i - n ] [ j ], p ), 1126 p ) ; 1127 1128 for (i = 0 ; i <= n - 1 ; ++i) 1129 1130 s[ i ] = temp[ i ] ; 1131 1132 } /* ======================== end of function product ======================= */ 1133 1134 1135 /*============================================================================== 1136 | times_x | 1137 ================================================================================ 1138 1139 DESCRIPTION 1140 1141 Compute x t(x) (mod f(x), p) 1142 1143 Uses a precomputed table of powers of x. 1144 1145 INPUT 1146 2 1147 t (int *) Coefficients of t (x), of degree <= n-1. 1148 1149 power_table (int **) x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic. 1150 1151 n (int, -infinity < n < infinity) Degree of f(x). 1152 1153 p (int, p > 0) Mod p coefficient arithmetic. 1154 1155 OUTPUT 1156 1157 t (int *) Overwritten with x t(x) (mod f(x), p) 1158 1159 EXAMPLE 1160 3 2 1161 Let n = 4, p = 5, and t(x) = 2 x + 4 x + 3 x. Let f(x) = 1162 4 2 4 3 2 1163 x + x + 2 x + 3. Then x t (x) = 2 x + 4 x + 3 x and 1164 2 3 2 1165 x t(x) (mod f(x), p) = 2 * (4 x + 3 x + 2) + 4 x + 3 x = 1166 3 2 1167 4 x + x + x + 4. 1168 3 2 1169 = 4 x + 2 x + 3 x + 2. 1170 1171 METHOD 1172 n-1 1173 Multiply t(x) = t x + ... + t by shifting the coefficients 1174 n-1 0 1175 n 1176 to the left. If an x term appears, eliminate it by 1177 1178 substitution using power_table. 1179 1180 BUGS 1181 1182 None. 1183 1184 -------------------------------------------------------------------------------- 1185 | Function Call | 1186 ------------------------------------------------------------------------------*/ 1187 1188 void 1189 times_x( int * t, int power_table[][ MAXDEGPOLY ], int n, int p ) 1190 { 1191 1192 /*------------------------------------------------------------------------------ 1193 | Local Variables | 1194 ------------------------------------------------------------------------------*/ 1195 1196 int 1197 i, /* Loop counter. */ 1198 coeff ; /* Coefficient of x ^ n term of x t(x). */ 1199 1200 /*------------------------------------------------------------------------------ 1201 | Function Body | 1202 ------------------------------------------------------------------------------*/ 1203 1204 /* 1205 Multiply t(x) by x. Do it by shifting the coefficients left in the array. 1206 */ 1207 1208 coeff = t[ n - 1 ] ; 1209 1210 for (i = n-1 ; i >= 1 ; --i) 1211 1212 t[ i ] = t[ i-1 ] ; 1213 1214 t[ 0 ] = 0 ; 1215 1216 /* 1217 n n 1218 Zero out t x . Replace it with t * [ x (mod f(x), p) ] from array 1219 n 1220 power_table. 1221 */ 1222 1223 if (coeff != 0) 1224 { 1225 for (i = 0 ; i <= n - 1 ; ++i) 1226 1227 t[ i ] = mod( t[ i ] + 1228 mod( coeff * power_table[ 0 ] [ i ], p ), 1229 p ) ; 1230 } 1231 1232 } /* ======================== end of function times_x ======================= */ 1233 1234 1235 1236 /*============================================================================== 1237 | x_to_power | 1238 ================================================================================ 1239 1240 DESCRIPTION 1241 m 1242 Compute g(x) = x (mod f(x), p). 1243 1244 INPUT 1245 1246 m (int) 1 <= m <= r. The exponent. 1247 1248 power_table (int **) x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic. 1249 1250 n (int, n >= 1) Degree of monic polynomial f(x). 1251 1252 p (int, p >= 2) Modulo p coefficient arithmetic. 1253 1254 OUTPUT 1255 1256 g (int *) Polynomial of degree <= n-1. 1257 1258 EXAMPLE 1259 4 2 1260 Let n = 4, p = 5, f(x) = x + x + 2 x + 3, and m = 156. 1261 1262 156 = 0 . . . 0 1 0 0 1 1 1 0 0 (binary representation) 1263 |<- ignore ->| S S SX SX SX S S (exponentiation rule, 1264 S = square, X = multiply by x) 1265 m 1266 x (mod f(x), p) = 1267 1268 2 2 1269 S x = x 1270 1271 4 2 1272 S x = 4 x + 3 x + 2 1273 1274 8 3 2 1275 S x = 4 x + 4 x + 1 1276 1277 9 3 2 1278 X x = 4 x + x + 3 x + 3 1279 1280 18 2 1281 S x = 2 x + x + 2 1282 1283 19 3 2 1284 X x = 2 x + x + 2 x 1285 1286 38 3 2 1287 S x = 2 x + 4 x + 3 x 1288 1289 39 3 2 1290 X x = 4 x + 2 x + x + 4 1291 1292 78 3 2 1293 S x = 4 x + 2 x + 3 x + 2 1294 1295 156 1296 S x = 3 1297 1298 METHOD 1299 1300 Exponentiation by repeated squaring, using precomputed table of 1301 powers. See ART OF COMPUTER PROGRAMMING, vol. 2, 2nd id, 1302 D. E. Knuth, pgs 441-443. 1303 1304 n 2n-2 1305 x, ... , x (mod f(x), p) 1306 1307 m 1308 to find g(x) = x (mod f(x), p), expand m into binary, 1309 1310 k k-1 1311 m = a 2 + a 2 + . . . + a 2 + a 1312 k k-1 1 0 1313 1314 m 1315 where a = 1, and split x apart into 1316 k 1317 1318 k k 1319 2 2 a 2 a a 1320 m k-1 1 0 1321 x = x x . . . x x 1322 1323 1324 Then to raise x to the mth power, do 1325 1326 1327 t( x ) = x 1328 1329 return if m = 1 1330 1331 1332 for i = k-1 downto 0 do 1333 1334 2 1335 t(x) = t(x) (mod f(x), p) // Square each time. 1336 1337 if a = 1 then 1338 i 1339 1340 t(x) = x t(x) (mod f(x), p) // Times x only if current bit is 1 1341 endif 1342 1343 endfor 1344 k 1345 2 1346 The initial t(x) = x gets squared k times to give x . If a = 1 for 1347 i 1348 0 <= i <= k-1, we multiply by x which then gets squared i times more 1349 1350 i 1351 2 1352 to give x . On a binary computer, we use bit shifting and masking to 1353 1354 identify the k bits { a . . . a } to the right of the leading 1 1355 k-1 0 1356 1357 bit. There are log ( m ) - 1 squarings and (number of 1 bits) - 1 1358 2 1359 multiplies. 1360 1361 BUGS 1362 1363 None. 1364 1365 -------------------------------------------------------------------------------- 1366 | Function Call | 1367 ------------------------------------------------------------------------------*/ 1368 1369 void 1370 x_to_power( bigint m, int * g, int power_table[][ MAXDEGPOLY ], int n, int p ) 1371 { 1372 1373 /*------------------------------------------------------------------------------ 1374 | Local Variables | 1375 ------------------------------------------------------------------------------*/ 1376 1377 /* 1378 mask = 1000 ... 000. That is, all bits of mask are zero except for the 1379 most significant bit of the computer word holding its value. Signed or 1380 unsigned type for bigint shouldn't matter here, since we just mask and 1381 compare bits. 1382 */ 1383 1384 bigint 1385 mask = ((bigint)1 << (NUMBITS - 1)) ; 1386 1387 int 1388 bit_count = 0, /* Number of bits in m to the right of the leading bit. */ 1389 i ; /* Loop counter. */ 1390 1391 1392 /*------------------------------------------------------------------------------ 1393 | Function Body | 1394 ------------------------------------------------------------------------------*/ 1395 1396 #ifdef DEBUG_PP_PRIMPOLY 1397 printf( "x to power = x ^ %lld\n", m ) ; 1398 #endif 1399 1400 /* 1401 Initialize g(x) to x. Exit right away if m = 1. 1402 */ 1403 1404 for (i = 0 ; i <= n-1 ; ++i) 1405 1406 g[ i ] = 0 ; 1407 1408 g[ 1 ] = 1 ; 1409 1410 if (m == 1) 1411 1412 return ; 1413 1414 /* 1415 Advance the leading bit of m up to the word's left hand boundary. 1416 Count how many bits were to the right of the leading bit. 1417 */ 1418 1419 while (! (m & mask)) 1420 { 1421 m <<= 1 ; 1422 ++bit_count ; 1423 } 1424 1425 bit_count = NUMBITS - bit_count ; 1426 1427 1428 1429 /* 1430 Exponentiation by repeated squaring. Discard the leading 1 bit. 1431 Thereafter, square for every 0 bit; square and multiply by x for 1432 every 1 bit. 1433 */ 1434 1435 while ( --bit_count > 0 ) 1436 { 1437 1438 m <<= 1 ; /* Expose the next bit. */ 1439 1440 square( g, power_table, n, p ) ; 1441 1442 #ifdef DEBUG_PP_PRIMPOLY 1443 printf( " after squaring, poly = \n" ) ; 1444 write_poly( g, n-1 ) ; 1445 printf( "\n\n" ); 1446 #endif 1447 1448 if (m & mask) /* Leading bit is 1. */ 1449 { 1450 times_x( g, power_table, n, p ) ; 1451 1452 #ifdef DEBUG_PP_PRIMPOLY 1453 printf( " after additional times x, poly = \n" ) ; 1454 write_poly( g, n-1 ) ; 1455 printf( "\n\n" ); 1456 #endif 1457 } 1458 } 1459 1460 } /* ===================== end of function x_to_power ======================= */