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.4 - A Program for Computing Primitive Polynomials.
29| Copyright (C) 1999-2025 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 seanerikoconnor!AT!gmail!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
61DESCRIPTION
62
63 Evaluate the polynomial f( x ) with modulo p arithmetic.
64
65INPUT
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
78RETURNS
79
80 f( x ) (int)
81
82EXAMPLE
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
88METHOD
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
96BUGS
97
98 None.
99
100--------------------------------------------------------------------------------
101| Function Call |
102------------------------------------------------------------------------------*/
103
104int
105 eval_poly( int * f, int x, int n, int p )
106{
107
108/*------------------------------------------------------------------------------
109| Local Variables |
110------------------------------------------------------------------------------*/
111
112int
113 val = 1 ; /* The value of f(x) which is returned. */
114
115/*------------------------------------------------------------------------------
116| Function Body |
117------------------------------------------------------------------------------*/
118
119while ( --n >= 0 )
120
121 val = mod( val * x + f[ n ], p ) ;
122
123return( val ) ;
124
125
126} /* ========================= end of function eval_poly ==================== */
127
128
129/*==============================================================================
130| linear_factor |
131================================================================================
132
133DESCRIPTION
134
135 Check if the polynomial f(x) has any linear factors.
136
137INPUT
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
145RETURNS
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
154EXAMPLE
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
159METHOD
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
164BUGS
165
166 None.
167
168--------------------------------------------------------------------------------
169| Function Call |
170------------------------------------------------------------------------------*/
171
172int
173 linear_factor( int * f, int n, int p )
174{
175
176/*------------------------------------------------------------------------------
177| Local Variables |
178------------------------------------------------------------------------------*/
179
180int i ; /* Loop counter. */
181
182/*------------------------------------------------------------------------------
183| Function Body |
184------------------------------------------------------------------------------*/
185
186for (i = 1 ; i <= p-1 ; ++i)
187
188 if (eval_poly( f, i, n, p ) == 0)
189
190 return( YES ) ;
191
192return( NO ) ;
193
194} /* ====================== end of function linear_factor =================== */
195
196
197/*==============================================================================
198| is_integer |
199================================================================================
200
201DESCRIPTION
202
203 Test if a polynomial is a constant.
204
205INPUT
206
207 t (int *) mod p polynomial t(x).
208 n (int) Its degree.
209
210RETURNS
211
212 YES if t(x) is a constant polynomial.
213 NO otherwise
214
215EXAMPLE
216 2
217 Let n = 3. Then t(x) = 2 x is not constant but t(x) = 2 is.
218
219METHOD
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
224BUGS
225
226 None.
227
228--------------------------------------------------------------------------------
229| Function Call |
230------------------------------------------------------------------------------*/
231
232int
233 is_integer( int * t, int n )
234{
235
236/*------------------------------------------------------------------------------
237| Function Body |
238------------------------------------------------------------------------------*/
239
240while ( n >= 1 )
241
242 if (t[ n-- ] != 0)
243
244 return( NO ) ;
245
246return( YES ) ;
247
248} /* ====================== end of function is_integer ====================== */
249
250
251/*==============================================================================
252| construct_power_table |
253================================================================================
254
255DESCRIPTION
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
262INPUT
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
268RETURNS
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
274EXAMPLE
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
304METHOD
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
315BUGS
316
317 None.
318
319--------------------------------------------------------------------------------
320| Function Call |
321------------------------------------------------------------------------------*/
322
323void
324 construct_power_table( int power_table[][ MAXDEGPOLY ], int * f,
325 int n, int p )
326{
327
328/*------------------------------------------------------------------------------
329| Local Variables |
330------------------------------------------------------------------------------*/
331
332int
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
348for (i = 0 ; i <= n-2 ; ++i)
349
350 t[ i ] = 0 ;
351
352t[ 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
361for (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
398return ;
399
400} /* ================== end of function construct_power_table =============== */
401
402
403
404/*==============================================================================
405| autoconvolve |
406================================================================================
407
408DESCRIPTION
409
410 Compute a convolution-like sum for use in function coeff_of_square.
411
412INPUT
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
424RETURNS
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
433EXAMPLE
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
445METHOD
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
451BUGS
452
453 None.
454
455--------------------------------------------------------------------------------
456| Function Call |
457------------------------------------------------------------------------------*/
458
459int
460 auto_convolve( int * t, int k, int lower, int upper, int p )
461{
462
463/*------------------------------------------------------------------------------
464| Local Variables |
465------------------------------------------------------------------------------*/
466
467int
468 sum = 0, /* Convolution sum. */
469 i ; /* Loop counter. */
470
471/*------------------------------------------------------------------------------
472| Function Body |
473------------------------------------------------------------------------------*/
474
475for (i = lower ; i <= upper ; ++i)
476
477 sum = mod( sum + mod( t[ i ] * t[ k - i ], p ), p ) ;
478
479return( sum ) ;
480
481} /* ==================== end of function auto_convolve ===================== */
482
483
484/*==============================================================================
485| convolve |
486================================================================================
487
488DESCRIPTION
489
490 Compute a convolution-like sum.
491
492INPUT
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
508RETURNS
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
517EXAMPLE
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
535METHOD
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
541BUGS
542
543 None.
544
545--------------------------------------------------------------------------------
546| Function Call |
547------------------------------------------------------------------------------*/
548
549int
550 convolve( int * s, int * t, int k, int lower, int upper, int p )
551{
552
553/*------------------------------------------------------------------------------
554| Local Variables |
555------------------------------------------------------------------------------*/
556
557int
558 sum = 0, /* Convolution sum. */
559 i ; /* Loop counter. */
560
561/*------------------------------------------------------------------------------
562| Function Body |
563------------------------------------------------------------------------------*/
564
565for (i = lower ; i <= upper ; ++i)
566
567 sum = mod( sum + mod( s[ i ] * t[ k - i ], p ), p ) ;
568
569return( sum ) ;
570
571} /* ======================= end of function convolve ======================= */
572
573
574/*==============================================================================
575| coeff_of_square |
576================================================================================
577
578DESCRIPTION
579 th 2
580 Return the coefficient of the k power of x in t( x ) modulo p.
581
582INPUT
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
592RETURNS
593 th 2
594 coefficient of the k power of x in t(x) mod p.
595
596EXAMPLE
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
616METHOD
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
678BUGS
679
680 None.
681
682--------------------------------------------------------------------------------
683| Function Call |
684------------------------------------------------------------------------------*/
685
686int
687 coeff_of_square( int * t, int k, int n, int p )
688{
689
690/*------------------------------------------------------------------------------
691| Local Variables |
692------------------------------------------------------------------------------*/
693
694int
695 sum = 0 ; /* kth coefficient of t(x) ^ 2. */
696
697/*------------------------------------------------------------------------------
698| Function Body |
699------------------------------------------------------------------------------*/
700
701if (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}
713else 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
726return( sum ) ;
727
728} /* ================== end of function coeff_of_square ===================== */
729
730
731/*==============================================================================
732| coeff_of_product |
733================================================================================
734
735DESCRIPTION
736 th
737 Return the coefficient of the k power of x in s( x ) t( x ) modulo p.
738
739INPUT
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
749RETURNS
750 th
751 coefficient of the k power of x in s(x) t(x) mod p.
752
753EXAMPLE
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
766METHOD
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
858BUGS
859
860 None.
861
862--------------------------------------------------------------------------------
863| Function Call |
864------------------------------------------------------------------------------*/
865
866int
867 coeff_of_product( int * s, int * t, int k, int n, int p )
868{
869
870/*------------------------------------------------------------------------------
871| Local Variables |
872------------------------------------------------------------------------------*/
873
874int
875 sum = 0 ; /* kth coefficient of t(x) ^ 2. */
876
877/*------------------------------------------------------------------------------
878| Function Body |
879------------------------------------------------------------------------------*/
880
881if (0 <= k && k <= n-1)
882{
883 sum = convolve( s, t, k, 0, k, p ) ;
884}
885else if (n <= k && k <= 2 * n - 2)
886{
887 sum = convolve( s, t, k, k - n + 1, n - 1, p ) ;
888}
889
890return( sum ) ;
891
892} /* ================== end of function coeff_of_product ==================== */
893
894
895/*==============================================================================
896| square |
897================================================================================
898
899DESCRIPTION
900 2
901 Compute t ( x ) (mod f(x), p)
902
903 Uses a precomputed table of powers of x.
904
905INPUT
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
915OUTPUT
916 2
917 t (int *) Overwritten with t (x) (mod f(x), p)
918
919EXAMPLE
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
936METHOD
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
953BUGS
954
955 None.
956
957--------------------------------------------------------------------------------
958| Function Call |
959------------------------------------------------------------------------------*/
960
961void
962 square( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
963{
964
965/*------------------------------------------------------------------------------
966| Local Variables |
967------------------------------------------------------------------------------*/
968
969int
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
985for (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
999for (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
1009for (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
1021DESCRIPTION
1022
1023 Compute s( x ) t( x ) (mod f(x), p)
1024
1025 Uses a precomputed table of powers of x.
1026
1027INPUT
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
1039OUTPUT
1040
1041 s (int *) Overwritten with s( x ) t( x ) (mod f(x), p)
1042
1043EXAMPLE
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
1060METHOD
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
1072BUGS
1073
1074 None.
1075
1076--------------------------------------------------------------------------------
1077| Function Call |
1078------------------------------------------------------------------------------*/
1079
1080void
1081product( int * s, int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1082{
1083
1084/*------------------------------------------------------------------------------
1085| Local Variables |
1086------------------------------------------------------------------------------*/
1087
1088int
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
1104for (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
1118for (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
1128for (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
1139DESCRIPTION
1140
1141 Compute x t(x) (mod f(x), p)
1142
1143 Uses a precomputed table of powers of x.
1144
1145INPUT
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
1155OUTPUT
1156
1157 t (int *) Overwritten with x t(x) (mod f(x), p)
1158
1159EXAMPLE
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
1171METHOD
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
1180BUGS
1181
1182 None.
1183
1184--------------------------------------------------------------------------------
1185| Function Call |
1186------------------------------------------------------------------------------*/
1187
1188void
1189 times_x( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1190{
1191
1192/*------------------------------------------------------------------------------
1193| Local Variables |
1194------------------------------------------------------------------------------*/
1195
1196int
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
1208coeff = t[ n - 1 ] ;
1209
1210for (i = n-1 ; i >= 1 ; --i)
1211
1212 t[ i ] = t[ i-1 ] ;
1213
1214t[ 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
1223if (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
1240DESCRIPTION
1241 m
1242 Compute g(x) = x (mod f(x), p).
1243
1244INPUT
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
1254OUTPUT
1255
1256 g (int *) Polynomial of degree <= n-1.
1257
1258EXAMPLE
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
1298METHOD
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
1361BUGS
1362
1363 None.
1364
1365--------------------------------------------------------------------------------
1366| Function Call |
1367------------------------------------------------------------------------------*/
1368
1369void
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
1384bigint
1385 mask = ((bigint)1 << (NUMBITS - 1)) ;
1386
1387int
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
1397printf( "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
1404for (i = 0 ; i <= n-1 ; ++i)
1405
1406 g[ i ] = 0 ;
1407
1408g[ 1 ] = 1 ;
1409
1410if (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
1419while (! (m & mask))
1420{
1421 m <<= 1 ;
1422 ++bit_count ;
1423}
1424
1425bit_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
1435while ( --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 ======================= */