1 /****************************************************************************
2 **
3 ** This file is part of GAP, a system for computational discrete algebra.
4 **
5 ** Copyright of GAP belongs to its developers, whose names are too numerous
6 ** to list here. Please refer to the COPYRIGHT file for details.
7 **
8 ** SPDX-License-Identifier: GPL-2.0-or-later
9 **
10 ** This file contains the functions to compute with elements from small
11 ** finite fields.
12 **
13 ** The concepts of this kernel module are documented in finfield.h
14 */
15
16 #include "finfield.h"
17
18 #include "ariths.h"
19 #include "bool.h"
20 #include "calls.h"
21 #include "error.h"
22 #include "finfield_conway.h"
23 #include "gvars.h"
24 #include "io.h"
25 #include "lists.h"
26 #include "modules.h"
27 #include "opers.h"
28 #include "plist.h"
29
30 #ifdef HPCGAP
31 #include "hpc/aobjects.h"
32 #include "hpc/thread.h"
33 #endif
34
35 /****************************************************************************
36 **
37 *V SuccFF . . . . . Tables for finite fields which are computed on demand
38 *V TypeFF
39 *V TypeFF0
40 **
41 ** SuccFF holds a plain list of successor lists.
42 ** TypeFF holds the types of typical elements of the finite fields.
43 ** TypeFF0 holds the types of the zero elements of the finite fields.
44 */
45
46 Obj SuccFF;
47 static Obj TypeFF;
48 static Obj TypeFF0;
49
50
51 /****************************************************************************
52 **
53 *V TYPE_FFE . . . . . kernel copy of GAP function TYPE_FFE
54 *V TYPE_FFE0 . . . . . kernel copy of GAP function TYPE_FFE0
55 **
56 ** These GAP functions are called to compute types of finite field elemnents
57 */
58 static Obj TYPE_FFE;
59 static Obj TYPE_FFE0;
60
61 /****************************************************************************
62 **
63 *V PrimitiveRootMod
64 **
65 ** Local copy of GAP function PrimitiveRootMod, used when initializing new
66 ** fields.
67 */
68 static Obj PrimitiveRootMod;
69
70
71 /****************************************************************************
72 **
73 *F LookupPrimePower(<q>) . . . . . . . . search for a prime power in tables
74 **
75 ** Searches the tables of prime powers from ffdata.c for q, returns the
76 ** index of q in SizeFF if it is present and 0 if not (the 0 position of
77 ** SizeFF is unused).
78 */
LookupPrimePower(UInt q)79 static FF LookupPrimePower(UInt q)
80 {
81 UInt l, n;
82 FF ff;
83 UInt e;
84
85 /* search through the finite field table */
86 l = 1;
87 n = NUM_SHORT_FINITE_FIELDS;
88 ff = 0;
89 while (l <= n && SizeFF[l] <= q && q <= SizeFF[n]) {
90 /* interpolation search */
91 /* cuts iterations roughly in half compared to binary search at
92 * the expense of additional divisions. */
93 e = (q - SizeFF[l] + 1) * (n - l) / (SizeFF[n] - SizeFF[l] + 1);
94 ff = l + e;
95 if (SizeFF[ff] == q)
96 break;
97 if (SizeFF[ff] < q)
98 l = ff + 1;
99 else
100 n = ff - 1;
101 }
102 if (ff < 1 || ff > NUM_SHORT_FINITE_FIELDS)
103 return 0;
104 if (SizeFF[ff] != q)
105 return 0;
106 return ff;
107 }
108
109
110 /****************************************************************************
111 **
112 *F FiniteField(<p>,<d>) . make the small finite field with <p>^<d> elements
113 *F FiniteFieldBySize(<q>) . . make the small finite field with <q> elements
114 **
115 ** The work is done in the Lookup function above, and in FiniteFieldBySize
116 ** where the successor tables are computed.
117 */
118
FiniteFieldBySize(UInt q)119 FF FiniteFieldBySize(UInt q)
120 {
121 FF ff; /* finite field, result */
122 Obj tmp; /* temporary bag */
123 Obj succBag; /* successor table bag */
124 FFV * succ; /* successor table */
125 FFV * indx; /* index table */
126 UInt p; /* characteristic of the field */
127 UInt poly; /* Conway polynomial of extension */
128 UInt i, l, f, n, e; /* loop variables */
129 Obj root; /* will be a primitive root mod p */
130
131 ff = LookupPrimePower(q);
132 if (!ff)
133 return 0;
134
135 #ifdef HPCGAP
136 /* Important correctness concern here:
137 *
138 * The values of SuccFF, TypeFF, and TypeFF0 are set in that
139 * order, separated by write barriers. This can happen concurrently
140 * in a different thread.
141 *
142 * Thus, after observing that TypeFF0 has been set, we can be sure
143 * that the thread also sees the values of SuccFF and TypeFF.
144 * This is ensured by the read barrier in ATOMIC_ELM_PLIST().
145 *
146 * In the worst case, we may do the following calculations once per
147 * thread and throw them away for all but one thread. Correctness
148 * is still ensured through the use of ATOMIC_SET_ELM_PLIST_ONCE(),
149 * which results in all threads sharing the same types and successor
150 * tables.
151 */
152 if (ATOMIC_ELM_PLIST(TypeFF0, ff))
153 return ff;
154 #else
155 if (ELM_PLIST(TypeFF0, ff))
156 return ff;
157 #endif
158
159 // determine the characteristic of the field
160 p = CHAR_FF(ff);
161
162 /* allocate a bag for the successor table and one for a temporary */
163 tmp = NewKernelBuffer(sizeof(Obj) + q * sizeof(FFV));
164 succBag = NewKernelBuffer(sizeof(Obj) + q * sizeof(FFV));
165
166 indx = (FFV *)(1 + ADDR_OBJ(tmp));
167 succ = (FFV *)(1 + ADDR_OBJ(succBag));
168
169 /* if q is a prime find the smallest primitive root $e$, use $x - e$ */
170
171 if (DEGR_FF(ff) == 1) {
172 if (p < 65537) {
173 /* for smaller primes we do this in the kernel for performance and
174 bootstrapping reasons
175 TODO -- review the threshold */
176 for (e = 1, i = 1; i != p - 1; ++e) {
177 for (f = e, i = 1; f != 1; ++i)
178 f = (f * e) % p;
179 }
180 }
181 else {
182 /* Otherwise we ask the library */
183 root = CALL_1ARGS(PrimitiveRootMod, INTOBJ_INT(p));
184 e = INT_INTOBJ(root) + 1;
185 }
186 poly = p - (e - 1);
187 }
188
189 /* otherwise look up the polynomial used to construct this field */
190 else {
191 for (i = 0; PolsFF[i] != q; i += 2)
192 ;
193 poly = PolsFF[i + 1];
194 }
195
196 /* construct 'indx' such that 'e = x^(indx[e]-1) % poly' for every e */
197 indx[0] = 0;
198 for (e = 1, n = 0; n < q - 1; ++n) {
199 indx[e] = n + 1;
200 /* e =p*e mod poly =x*e mod poly =x*x^n mod poly =x^{n+1} mod poly */
201 if (p != 2) {
202 f = p * (e % (q / p));
203 l = ((p - 1) * (e / (q / p))) % p;
204 e = 0;
205 for (i = 1; i < q; i *= p)
206 e = e + i * ((f / i + l * (poly / i)) % p);
207 }
208 else {
209 if (2 * e & q)
210 e = 2 * e ^ poly ^ q;
211 else
212 e = 2 * e;
213 }
214 }
215
216 /* construct 'succ' such that 'x^(n-1)+1 = x^(succ[n]-1)' for every n */
217 succ[0] = q - 1;
218 for (e = 1, f = p - 1; e < q; e++) {
219 if (e < f) {
220 succ[indx[e]] = indx[e + 1];
221 }
222 else {
223 succ[indx[e]] = indx[e + 1 - p];
224 f += p;
225 }
226 }
227
228 /* enter the finite field in the tables */
229 #ifdef HPCGAP
230 MakeBagReadOnly(succBag);
231 ATOMIC_SET_ELM_PLIST_ONCE(SuccFF, ff, succBag);
232 CHANGED_BAG(SuccFF);
233 tmp = CALL_1ARGS(TYPE_FFE, INTOBJ_INT(p));
234 ATOMIC_SET_ELM_PLIST_ONCE(TypeFF, ff, tmp);
235 CHANGED_BAG(TypeFF);
236 tmp = CALL_1ARGS(TYPE_FFE0, INTOBJ_INT(p));
237 ATOMIC_SET_ELM_PLIST_ONCE(TypeFF0, ff, tmp);
238 CHANGED_BAG(TypeFF0);
239 #else
240 ASS_LIST(SuccFF, ff, succBag);
241 CHANGED_BAG(SuccFF);
242 tmp = CALL_1ARGS(TYPE_FFE, INTOBJ_INT(p));
243 ASS_LIST(TypeFF, ff, tmp);
244 CHANGED_BAG(TypeFF);
245 tmp = CALL_1ARGS(TYPE_FFE0, INTOBJ_INT(p));
246 ASS_LIST(TypeFF0, ff, tmp);
247 CHANGED_BAG(TypeFF0);
248 #endif
249
250 /* return the finite field */
251 return ff;
252 }
253
FiniteField(UInt p,UInt d)254 FF FiniteField(UInt p, UInt d)
255 {
256 UInt q, i;
257 FF ff;
258
259 q = 1;
260 for (i = 1; i <= d; i++)
261 q *= p;
262
263 ff = FiniteFieldBySize(q);
264 if (ff != 0 && CHAR_FF(ff) != p)
265 return 0;
266 return ff;
267 }
268
269
270 /****************************************************************************
271 **
272 *F CommonFF(<f1>,<d1>,<f2>,<d2>) . . . . . . . . . . . . find a common field
273 **
274 ** 'CommonFF' returns a small finite field that can represent elements of
275 ** degree <d1> from the small finite field <f1> and elements of degree <d2>
276 ** from the small finite field <f2>. Note that this is not guaranteed to be
277 ** the smallest such field. If <f1> and <f2> have different characteristic
278 ** or the smallest common field, is too large, 'CommonFF' returns 0.
279 */
CommonFF(FF f1,UInt d1,FF f2,UInt d2)280 FF CommonFF (
281 FF f1,
282 UInt d1,
283 FF f2,
284 UInt d2 )
285 {
286 UInt p; /* characteristic */
287 UInt d; /* degree */
288
289 /* trivial case first */
290 if ( f1 == f2 ) {
291 return f1;
292 }
293
294 /* get and check the characteristics */
295 p = CHAR_FF( f1 );
296 if ( p != CHAR_FF( f2 ) ) {
297 return 0;
298 }
299
300 /* check whether one of the fields will do */
301 if ( DEGR_FF(f1) % d2 == 0 ) {
302 return f1;
303 }
304 if ( DEGR_FF(f2) % d1 == 0 ) {
305 return f2;
306 }
307
308 /* compute the neccessary degree */
309 d = d1;
310 while ( d % d2 != 0 ) {
311 d += d1;
312 }
313
314 /* try to build the field */
315 return FiniteField( p, d );
316 }
317
318
319 /****************************************************************************
320 **
321 *F CharFFE(<ffe>) . . . . . . . . . characteristic of a small finite field
322 **
323 ** 'CharFFE' returns the characteristic of the small finite field in which
324 ** the element <ffe> lies.
325 */
CharFFE(Obj ffe)326 UInt CharFFE (
327 Obj ffe )
328 {
329 return CHAR_FF( FLD_FFE(ffe) );
330 }
331
FuncCHAR_FFE_DEFAULT(Obj self,Obj ffe)332 static Obj FuncCHAR_FFE_DEFAULT(Obj self, Obj ffe)
333 {
334 return INTOBJ_INT( CHAR_FF( FLD_FFE(ffe) ) );
335 }
336
337
338 /****************************************************************************
339 **
340 *F DegreeFFE(<ffe>) . . . . . . . . . . . . degree of a small finite field
341 **
342 ** 'DegreeFFE' returns the degree of the smallest finite field in which the
343 ** element <ffe> lies.
344 */
DegreeFFE(Obj ffe)345 UInt DegreeFFE (
346 Obj ffe )
347 {
348 UInt d; /* degree, result */
349 FFV val; /* value of element */
350 FF fld; /* field of element */
351 UInt q; /* size of field */
352 UInt p; /* char. of field */
353 UInt m; /* size of minimal field */
354
355 /* get the value, the field, the size, and the characteristic */
356 val = VAL_FFE( ffe );
357 fld = FLD_FFE( ffe );
358 q = SIZE_FF( fld );
359 p = CHAR_FF( fld );
360
361 /* the zero element has a degree of one */
362 if ( val == 0 ) {
363 return 1L;
364 }
365
366 /* compute the degree */
367 m = p;
368 d = 1;
369 while ( (q-1) % (m-1) != 0 || (val-1) % ((q-1)/(m-1)) != 0 ) {
370 m *= p;
371 d += 1;
372 }
373
374 /* return the degree */
375 return d;
376 }
377
FuncDEGREE_FFE_DEFAULT(Obj self,Obj ffe)378 static Obj FuncDEGREE_FFE_DEFAULT(Obj self, Obj ffe)
379 {
380 return INTOBJ_INT( DegreeFFE( ffe ) );
381 }
382
383
384 /****************************************************************************
385 **
386 *F TypeFFE(<ffe>) . . . . . . . . . . type of element of small finite field
387 **
388 ** 'TypeFFE' returns the type of the element <ffe> of a small finite field.
389 **
390 ** 'TypeFFE' is the function in 'TypeObjFuncs' for elements in small finite
391 ** fields.
392 */
TypeFFE(Obj ffe)393 Obj TypeFFE(Obj ffe)
394 {
395 Obj types = (VAL_FFE(ffe) == 0) ? TypeFF0 : TypeFF;
396 #ifdef HPCGAP
397 return ATOMIC_ELM_PLIST(types, FLD_FFE(ffe));
398 #else
399 return ELM_PLIST(types, FLD_FFE(ffe));
400 #endif
401 }
402
403
404 /****************************************************************************
405 **
406 *F EqFFE(<opL>,<opR>) . . . . . . . test if finite field elements are equal
407 **
408 ** 'EqFFE' returns 'True' if the two finite field elements <opL> and <opR>
409 ** are equal and 'False' othwise.
410 **
411 ** This is complicated because it must account for the following situation.
412 ** Suppose 'a' is 'Z(3)', 'b' is 'Z(3^2)^4' and finally 'c' is 'Z(3^3)^13'.
413 ** Mathematically 'a' is equal to 'b', so we want 'a = b' to be 'true' and
414 ** since 'a' is represented over a subfield of 'b' this is no big problem.
415 ** Again 'a' is equal to 'c', and again we want 'a = c' to be 'true' and
416 ** again this is no problem since 'a' is represented over a subfield of 'c'.
417 ** Since '=' ought to be transitive we also want 'b = c' to be 'true' and
418 ** this is a problem, because they are represented over incompatible fields.
419 */
EqFFE(Obj opL,Obj opR)420 static Int EqFFE(Obj opL, Obj opR)
421 {
422 FFV vL, vR; /* value of left and right */
423 FF fL, fR; /* field of left and right */
424 UInt pL, pR; /* char. of left and right */
425 UInt qL, qR; /* size of left and right */
426 UInt mL, mR; /* size of minimal field */
427
428 /* get the values and the fields over which they are represented */
429 vL = VAL_FFE( opL );
430 vR = VAL_FFE( opR );
431 fL = FLD_FFE( opL );
432 fR = FLD_FFE( opR );
433
434 /* if the elements are represented over the same field, it is easy */
435 if ( fL == fR ) {
436 return (vL == vR);
437 }
438
439 /* elements in fields of different characteristic are different too */
440 pL = CHAR_FF( fL );
441 pR = CHAR_FF( fR );
442 if ( pL != pR ) {
443 return 0L;
444 }
445
446 /* the zero element is not equal to any other element */
447 if ( vL == 0 || vR == 0 ) {
448 return (vL == 0 && vR == 0);
449 }
450
451 /* compute the sizes of the minimal fields in which the elements lie */
452 qL = SIZE_FF( fL );
453 mL = pL;
454 while ( (qL-1) % (mL-1) != 0 || (vL-1) % ((qL-1)/(mL-1)) != 0 ) mL *= pL;
455 qR = SIZE_FF( fR );
456 mR = pR;
457 while ( (qR-1) % (mR-1) != 0 || (vR-1) % ((qR-1)/(mR-1)) != 0 ) mR *= pR;
458
459 /* elements in different fields are different too */
460 if ( mL != mR ) {
461 return 0L;
462 }
463
464 /* otherwise compare the elements in the common minimal field */
465 return ((vL-1)/((qL-1)/(mL-1)) == (vR-1)/((qR-1)/(mR-1)));
466 }
467
468
469 /****************************************************************************
470 **
471 *F LtFFE(<opL>,<opR>) . . . . . . test if finite field elements is smaller
472 **
473 ** 'LtFFEFFE' returns 'True' if the finite field element <opL> is strictly
474 ** less than the finite field element <opR> and 'False' otherwise.
475 */
LtFFE(Obj opL,Obj opR)476 static Int LtFFE(Obj opL, Obj opR)
477 {
478 FFV vL, vR; /* value of left and right */
479 FF fL, fR; /* field of left and right */
480 UInt pL, pR; /* char. of left and right */
481 UInt qL, qR; /* size of left and right */
482 UInt mL, mR; /* size of minimal field */
483
484 /* get the values and the fields over which they are represented */
485 vL = VAL_FFE( opL );
486 vR = VAL_FFE( opR );
487 fL = FLD_FFE( opL );
488 fR = FLD_FFE( opR );
489
490 /* elements in fields of different characteristic are not comparable */
491 pL = CHAR_FF( fL );
492 pR = CHAR_FF( fR );
493 if ( pL != pR ) {
494 return (DoOperation2Args( LtOper, opL, opR ) == True);
495 }
496
497 /* the zero element is smaller than any other element */
498 if ( vL == 0 || vR == 0 ) {
499 return (vL == 0 && vR != 0);
500 }
501
502 /* get the sizes of the fields over which the elements are written */
503 qL = SIZE_FF( fL );
504 qR = SIZE_FF( fR );
505
506 /* Deal quickly with the case where both elements are written over the ground field */
507 if (qL ==pL && qR == pR)
508 return vL < vR;
509
510 /* compute the sizes of the minimal fields in which the elements lie */
511 mL = pL;
512 while ( (qL-1) % (mL-1) != 0 || (vL-1) % ((qL-1)/(mL-1)) != 0 ) mL *= pL;
513 mR = pR;
514 while ( (qR-1) % (mR-1) != 0 || (vR-1) % ((qR-1)/(mR-1)) != 0 ) mR *= pR;
515
516 /* elements in smaller fields are smaller too */
517 if ( mL != mR ) {
518 return (mL < mR);
519 }
520
521 /* otherwise compare the elements in the common minimal field */
522 return ((vL-1)/((qL-1)/(mL-1)) < (vR-1)/((qR-1)/(mR-1)));
523 }
524
525
526 /****************************************************************************
527 **
528 *F PrFFV(<fld>,<val>) . . . . . . . . . . . . . print a finite field value
529 **
530 ** 'PrFFV' prints the value <val> from the finite field <fld>.
531 **
532 */
PrFFV(FF fld,FFV val)533 static void PrFFV(FF fld, FFV val)
534 {
535 UInt q; /* size of finite field */
536 UInt p; /* char. of finite field */
537 UInt m; /* size of minimal field */
538 UInt d; /* degree of minimal field */
539
540 /* get the characteristic, order of the minimal field and the degree */
541 q = SIZE_FF( fld );
542 p = CHAR_FF( fld );
543
544 /* print the zero */
545 if ( val == 0 ) {
546 Pr( "%>0*Z(%>%d%2<)", (Int)p, 0L );
547 }
548
549 /* print a nonzero element as power of the primitive root */
550 else {
551
552 /* find the degree of the minimal field in that the element lies */
553 d = 1; m = p;
554 while ( (q-1) % (m-1) != 0 || (val-1) % ((q-1)/(m-1)) != 0 ) {
555 d++; m *= p;
556 }
557 val = (val-1) / ((q-1)/(m-1)) + 1;
558
559 /* print the element */
560 Pr( "%>Z(%>%d%<", (Int)p, 0L );
561 if ( d == 1 ) {
562 Pr( "%<)", 0L, 0L );
563 }
564 else {
565 Pr( "^%>%d%2<)", (Int)d, 0L );
566 }
567 if ( val != 2 ) {
568 Pr( "^%>%d%<", (Int)(val-1), 0L );
569 }
570 }
571
572 }
573
574
575 /****************************************************************************
576 **
577 *F PrFFE(<ffe>) . . . . . . . . . . . . . . . print a finite field element
578 **
579 ** 'PrFFE' prints the finite field element <ffe>.
580 */
PrFFE(Obj ffe)581 static void PrFFE(Obj ffe)
582 {
583 PrFFV( FLD_FFE(ffe), VAL_FFE(ffe) );
584 }
585
586
587 /****************************************************************************
588 **
589 *F SumFFEFFE(<opL>,<opR>) . . . . . . . . . . sum of finite field elements
590 **
591 ** 'SumFFEFFE' returns the sum of the two finite field elements <opL> and
592 ** <opR>. The sum is represented over the field over which the operands are
593 ** represented, even if it lies in a much smaller field.
594 **
595 ** If one of the elements is represented over a subfield of the field over
596 ** which the other element is represented, it is lifted into the larger
597 ** field before the addition.
598 **
599 ** 'SumFFEFFE' just does the conversions mentioned above and then calls the
600 ** macro 'SUM_FFV' to do the actual addition.
601 */
602 static Obj SUM_FFE_LARGE;
603
SumFFEFFE(Obj opL,Obj opR)604 static Obj SumFFEFFE(Obj opL, Obj opR)
605 {
606 FFV vL, vR, vX; /* value of left, right, result */
607 FF fL, fR, fX; /* field of left, right, result */
608 UInt qL, qR, qX; /* size of left, right, result */
609
610 /* get the values, handle trivial cases */
611 vL = VAL_FFE( opL );
612 vR = VAL_FFE( opR );
613
614 /* bring the two operands into a common field <fX> */
615 fL = FLD_FFE( opL );
616 qL = SIZE_FF( fL );
617 fR = FLD_FFE( opR );
618 qR = SIZE_FF( fR );
619
620 if ( qL == qR ) {
621 fX = fL;
622 }
623 else if ( qL % qR == 0 && (qL-1) % (qR-1) == 0 ) {
624 fX = fL;
625 if ( vR != 0 ) vR = (qL-1) / (qR-1) * (vR-1) + 1;
626 }
627 else if ( qR % qL == 0 && (qR-1) % (qL-1) == 0 ) {
628 fX = fR;
629 if ( vL != 0 ) vL = (qR-1) / (qL-1) * (vL-1) + 1;
630 }
631 else {
632 fX = CommonFF( fL, DegreeFFE(opL), fR, DegreeFFE(opR) );
633 if ( fX == 0 ) return CALL_2ARGS( SUM_FFE_LARGE, opL, opR );
634 qX = SIZE_FF( fX );
635 /* if ( vL != 0 ) vL = (qX-1) / (qL-1) * (vL-1) + 1; */
636 if ( vL != 0 ) vL = ((qX-1) * (vL-1)) / (qL-1) + 1;
637 /* if ( vR != 0 ) vR = (qX-1) / (qR-1) * (vR-1) + 1; */
638 if ( vR != 0 ) vR = ((qX-1) * (vR-1)) / (qR-1) + 1;
639 }
640
641 /* compute and return the result */
642 vX = SUM_FFV( vL, vR, SUCC_FF(fX) );
643 return NEW_FFE( fX, vX );
644 }
645
SumFFEInt(Obj opL,Obj opR)646 static Obj SumFFEInt(Obj opL, Obj opR)
647 {
648 FFV vL, vR, vX; /* value of left, right, result */
649 FF fX; /* field of result */
650 Int pX; /* char. of result */
651 const FFV* sX; /* successor table of result field */
652
653 /* get the field for the result */
654 fX = FLD_FFE( opL );
655 pX = CHAR_FF( fX );
656 sX = SUCC_FF( fX );
657
658 /* get the right operand */
659 vX = ((INT_INTOBJ( opR ) % pX) + pX) % pX;
660 if ( vX == 0 ) {
661 vR = 0;
662 }
663 else {
664 vR = 1;
665 for ( ; 1 < vX; vX-- ) vR = sX[vR];
666 }
667
668 /* get the left operand */
669 vL = VAL_FFE( opL );
670
671 /* compute and return the result */
672 vX = SUM_FFV( vL, vR, sX );
673 return NEW_FFE( fX, vX );
674 }
675
SumIntFFE(Obj opL,Obj opR)676 static Obj SumIntFFE(Obj opL, Obj opR)
677 {
678 FFV vL, vR, vX; /* value of left, right, result */
679 FF fX; /* field of result */
680 Int pX; /* char. of result */
681 const FFV* sX; /* successor table of result field */
682
683 /* get the field for the result */
684 fX = FLD_FFE( opR );
685 pX = CHAR_FF( fX );
686 sX = SUCC_FF( fX );
687
688 /* get the left operand */
689 vX = ((INT_INTOBJ( opL ) % pX) + pX) % pX;
690 if ( vX == 0 ) {
691 vL = 0;
692 }
693 else {
694 vL = 1;
695 for ( ; 1 < vX; vX-- ) vL = sX[vL];
696 }
697
698 /* get the right operand */
699 vR = VAL_FFE( opR );
700
701 /* compute and return the result */
702 vX = SUM_FFV( vL, vR, sX );
703 return NEW_FFE( fX, vX );
704 }
705
706
707 /****************************************************************************
708 **
709 *F ZeroFFE(<op>) . . . . . . . . . . . . . . zero of a finite field element
710 */
ZeroFFE(Obj op)711 static Obj ZeroFFE(Obj op)
712 {
713 FF fX; /* field of result */
714
715 /* get the field for the result */
716 fX = FLD_FFE( op );
717
718 /* return the result */
719 return NEW_FFE( fX, 0 );
720 }
721
722
723 /****************************************************************************
724 **
725 *F AInvFFE(<op>) . . . . . . . . . . additive inverse of finite field element
726 */
AInvFFE(Obj op)727 static Obj AInvFFE(Obj op)
728 {
729 FFV v, vX; /* value of operand, result */
730 FF fX; /* field of result */
731 const FFV* sX; /* successor table of result field */
732
733 /* get the field for the result */
734 fX = FLD_FFE( op );
735 sX = SUCC_FF( fX );
736
737 /* get the operand */
738 v = VAL_FFE( op );
739
740 /* compute and return the result */
741 vX = NEG_FFV( v, sX );
742 return NEW_FFE( fX, vX );
743 }
744
745
746 /****************************************************************************
747 **
748 *F DiffFFEFFE(<opL>,<opR>) . . . . . . . difference of finite field elements
749 **
750 ** 'DiffFFEFFE' returns the difference of the two finite field elements
751 ** <opL> and <opR>. The difference is represented over the field over which
752 ** the operands are represented, even if it lies in a much smaller field.
753 **
754 ** If one of the elements is represented over a subfield of the field over
755 ** which the other element is represented, it is lifted into the larger
756 ** field before the subtraction.
757 **
758 ** 'DiffFFEFFE' just does the conversions mentioned above and then calls the
759 ** macros 'NEG_FFV' and 'SUM_FFV' to do the actual subtraction.
760 */
761 static Obj DIFF_FFE_LARGE;
762
DiffFFEFFE(Obj opL,Obj opR)763 static Obj DiffFFEFFE(Obj opL, Obj opR)
764 {
765 FFV vL, vR, vX; /* value of left, right, result */
766 FF fL, fR, fX; /* field of left, right, result */
767 UInt qL, qR, qX; /* size of left, right, result */
768
769 /* get the values, handle trivial cases */
770 vL = VAL_FFE( opL );
771 vR = VAL_FFE( opR );
772
773 /* bring the two operands into a common field <fX> */
774 fL = FLD_FFE( opL );
775 qL = SIZE_FF( fL );
776 fR = FLD_FFE( opR );
777 qR = SIZE_FF( fR );
778
779 if ( qL == qR ) {
780 fX = fL;
781 }
782 else if ( qL % qR == 0 && (qL-1) % (qR-1) == 0 ) {
783 fX = fL;
784 if ( vR != 0 ) vR = (qL-1) / (qR-1) * (vR-1) + 1;
785 }
786 else if ( qR % qL == 0 && (qR-1) % (qL-1) == 0 ) {
787 fX = fR;
788 if ( vL != 0 ) vL = (qR-1) / (qL-1) * (vL-1) + 1;
789 }
790 else {
791 fX = CommonFF( fL, DegreeFFE(opL), fR, DegreeFFE(opR) );
792 if ( fX == 0 ) return CALL_2ARGS( DIFF_FFE_LARGE, opL, opR );
793 qX = SIZE_FF( fX );
794 /* if ( vL != 0 ) vL = (qX-1) / (qL-1) * (vL-1) + 1; */
795 if ( vL != 0 ) vL = ((qX-1) * (vL-1)) / (qL-1) + 1;
796 /* if ( vR != 0 ) vR = (qX-1) / (qR-1) * (vR-1) + 1; */
797 if ( vR != 0 ) vR = ((qX-1) * (vR-1)) / (qR-1) + 1;
798 }
799
800 /* compute and return the result */
801 vR = NEG_FFV( vR, SUCC_FF(fX) );
802 vX = SUM_FFV( vL, vR, SUCC_FF(fX) );
803 return NEW_FFE( fX, vX );
804 }
805
DiffFFEInt(Obj opL,Obj opR)806 static Obj DiffFFEInt(Obj opL, Obj opR)
807 {
808 FFV vL, vR, vX; /* value of left, right, result */
809 FF fX; /* field of result */
810 Int pX; /* char. of result */
811 const FFV* sX; /* successor table of result field */
812
813 /* get the field for the result */
814 fX = FLD_FFE( opL );
815 pX = CHAR_FF( fX );
816 sX = SUCC_FF( fX );
817
818 /* get the right operand */
819 vX = ((INT_INTOBJ( opR ) % pX) + pX) % pX;
820 if ( vX == 0 ) {
821 vR = 0;
822 }
823 else {
824 vR = 1;
825 for ( ; 1 < vX; vX-- ) vR = sX[vR];
826 }
827
828 /* get the left operand */
829 vL = VAL_FFE( opL );
830
831 /* compute and return the result */
832 vR = NEG_FFV( vR, sX );
833 vX = SUM_FFV( vL, vR, sX );
834 return NEW_FFE( fX, vX );
835 }
836
DiffIntFFE(Obj opL,Obj opR)837 static Obj DiffIntFFE(Obj opL, Obj opR)
838 {
839 FFV vL, vR, vX; /* value of left, right, result */
840 FF fX; /* field of result */
841 Int pX; /* char. of result */
842 const FFV* sX; /* successor table of result field */
843
844 /* get the field for the result */
845 fX = FLD_FFE( opR );
846 pX = CHAR_FF( fX );
847 sX = SUCC_FF( fX );
848
849 /* get the left operand */
850 vX = ((INT_INTOBJ( opL ) % pX) + pX) % pX;
851 if ( vX == 0 ) {
852 vL = 0;
853 }
854 else {
855 vL = 1;
856 for ( ; 1 < vX; vX-- ) vL = sX[vL];
857 }
858
859 /* get the right operand */
860 vR = VAL_FFE( opR );
861
862 /* compute and return the result */
863 vR = NEG_FFV( vR, sX );
864 vX = SUM_FFV( vL, vR, sX );
865 return NEW_FFE( fX, vX );
866 }
867
868
869 /****************************************************************************
870 **
871 *F ProdFFEFFE(<opL>,<opR>) . . . . . . . . product of finite field elements
872 **
873 ** 'ProdFFEFFE' returns the product of the two finite field elements <opL>
874 ** and <opR>. The product is represented over the field over which the
875 ** operands are represented, even if it lies in a much smaller field.
876 **
877 ** If one of the elements is represented over a subfield of the field over
878 ** which the other element is represented, it is lifted into the larger
879 ** field before the multiplication.
880 **
881 ** 'ProdFFEFFE' just does the conversions mentioned above and then calls the
882 ** macro 'PROD_FFV' to do the actual multiplication.
883 */
884 static Obj PROD_FFE_LARGE;
885
ProdFFEFFE(Obj opL,Obj opR)886 static Obj ProdFFEFFE(Obj opL, Obj opR)
887 {
888 FFV vL, vR, vX; /* value of left, right, result */
889 FF fL, fR, fX; /* field of left, right, result */
890 UInt qL, qR, qX; /* size of left, right, result */
891
892 /* get the values, handle trivial cases */
893 vL = VAL_FFE( opL );
894 vR = VAL_FFE( opR );
895
896 /* bring the two operands into a common field <fX> */
897 fL = FLD_FFE( opL );
898 qL = SIZE_FF( fL );
899 fR = FLD_FFE( opR );
900 qR = SIZE_FF( fR );
901
902 if ( qL == qR ) {
903 fX = fL;
904 }
905 else if ( qL % qR == 0 && (qL-1) % (qR-1) == 0 ) {
906 fX = fL;
907 if ( vR != 0 ) vR = (qL-1) / (qR-1) * (vR-1) + 1;
908 }
909 else if ( qR % qL == 0 && (qR-1) % (qL-1) == 0 ) {
910 fX = fR;
911 if ( vL != 0 ) vL = (qR-1) / (qL-1) * (vL-1) + 1;
912 }
913 else {
914 fX = CommonFF( fL, DegreeFFE(opL), fR, DegreeFFE(opR) );
915 if ( fX == 0 ) return CALL_2ARGS( PROD_FFE_LARGE, opL, opR );
916 qX = SIZE_FF( fX );
917 /* if ( vL != 0 ) vL = (qX-1) / (qL-1) * (vL-1) + 1; */
918 if ( vL != 0 ) vL = ((qX-1) * (vL-1)) / (qL-1) + 1;
919 /* if ( vR != 0 ) vR = (qX-1) / (qR-1) * (vR-1) + 1; */
920 if ( vR != 0 ) vR = ((qX-1) * (vR-1)) / (qR-1) + 1;
921 }
922
923 /* compute and return the result */
924 vX = PROD_FFV( vL, vR, SUCC_FF(fX) );
925 return NEW_FFE( fX, vX );
926 }
927
ProdFFEInt(Obj opL,Obj opR)928 static Obj ProdFFEInt(Obj opL, Obj opR)
929 {
930 FFV vL, vR, vX; /* value of left, right, result */
931 FF fX; /* field of result */
932 Int pX; /* char. of result */
933 const FFV* sX; /* successor table of result field */
934
935 /* get the field for the result */
936 fX = FLD_FFE( opL );
937 pX = CHAR_FF( fX );
938 sX = SUCC_FF( fX );
939
940 /* get the right operand */
941 vX = ((INT_INTOBJ( opR ) % pX) + pX) % pX;
942 if ( vX == 0 ) {
943 vR = 0;
944 }
945 else {
946 vR = 1;
947 for ( ; 1 < vX; vX-- ) vR = sX[vR];
948 }
949
950 /* get the left operand */
951 vL = VAL_FFE( opL );
952
953 /* compute and return the result */
954 vX = PROD_FFV( vL, vR, sX );
955 return NEW_FFE( fX, vX );
956 }
957
ProdIntFFE(Obj opL,Obj opR)958 static Obj ProdIntFFE(Obj opL, Obj opR)
959 {
960 FFV vL, vR, vX; /* value of left, right, result */
961 FF fX; /* field of result */
962 Int pX; /* char. of result */
963 const FFV* sX; /* successor table of result field */
964
965 /* get the field for the result */
966 fX = FLD_FFE( opR );
967 pX = CHAR_FF( fX );
968 sX = SUCC_FF( fX );
969
970 /* get the left operand */
971 vX = ((INT_INTOBJ( opL ) % pX) + pX) % pX;
972 if ( vX == 0 ) {
973 vL = 0;
974 }
975 else {
976 vL = 1;
977 for ( ; 1 < vX; vX-- ) vL = sX[vL];
978 }
979
980 /* get the right operand */
981 vR = VAL_FFE( opR );
982
983 /* compute and return the result */
984 vX = PROD_FFV( vL, vR, sX );
985 return NEW_FFE( fX, vX );
986 }
987
988
989 /****************************************************************************
990 **
991 *F OneFFE(<op>) . . . . . . . . . . . . . . . one of a finite field element
992 */
OneFFE(Obj op)993 static Obj OneFFE(Obj op)
994 {
995 FF fX; /* field of result */
996
997 /* get the field for the result */
998 fX = FLD_FFE( op );
999
1000 /* return the result */
1001 return NEW_FFE( fX, 1 );
1002 }
1003
1004
1005 /****************************************************************************
1006 **
1007 *F InvFFE(<op>) . . . . . . . . . . . . . . inverse of finite field element
1008 */
InvFFE(Obj op)1009 static Obj InvFFE(Obj op)
1010 {
1011 FFV v, vX; /* value of operand, result */
1012 FF fX; /* field of result */
1013 const FFV* sX; /* successor table of result field */
1014
1015 /* get the field for the result */
1016 fX = FLD_FFE( op );
1017 sX = SUCC_FF( fX );
1018
1019 /* get the operand */
1020 v = VAL_FFE( op );
1021 if ( v == 0 ) return Fail;
1022
1023 /* compute and return the result */
1024 vX = QUO_FFV( 1, v, sX );
1025 return NEW_FFE( fX, vX );
1026 }
1027
1028
1029 /****************************************************************************
1030 **
1031 *F QuoFFEFFE(<opL>,<opR>) . . . . . . . . quotient of finite field elements
1032 **
1033 ** 'QuoFFEFFE' returns the quotient of the two finite field elements <opL>
1034 ** and <opR>. The quotient is represented over the field over which the
1035 ** operands are represented, even if it lies in a much smaller field.
1036 **
1037 ** If one of the elements is represented over a subfield of the field over
1038 ** which the other element is represented, it is lifted into the larger
1039 ** field before the division.
1040 **
1041 ** 'QuoFFEFFE' just does the conversions mentioned above and then calls the
1042 ** macro 'QUO_FFV' to do the actual division.
1043 */
1044 static Obj QUO_FFE_LARGE;
1045
QuoFFEFFE(Obj opL,Obj opR)1046 static Obj QuoFFEFFE(Obj opL, Obj opR)
1047 {
1048 FFV vL, vR, vX; /* value of left, right, result */
1049 FF fL, fR, fX; /* field of left, right, result */
1050 UInt qL, qR, qX; /* size of left, right, result */
1051
1052 /* get the values, handle trivial cases */
1053 vL = VAL_FFE( opL );
1054 vR = VAL_FFE( opR );
1055
1056 /* bring the two operands into a common field <fX> */
1057 fL = FLD_FFE( opL );
1058 qL = SIZE_FF( fL );
1059 fR = FLD_FFE( opR );
1060 qR = SIZE_FF( fR );
1061
1062 if ( qL == qR ) {
1063 fX = fL;
1064 }
1065 else if ( qL % qR == 0 && (qL-1) % (qR-1) == 0 ) {
1066 fX = fL;
1067 if ( vR != 0 ) vR = (qL-1) / (qR-1) * (vR-1) + 1;
1068 }
1069 else if ( qR % qL == 0 && (qR-1) % (qL-1) == 0 ) {
1070 fX = fR;
1071 if ( vL != 0 ) vL = (qR-1) / (qL-1) * (vL-1) + 1;
1072 }
1073 else {
1074 fX = CommonFF( fL, DegreeFFE(opL), fR, DegreeFFE(opR) );
1075 if ( fX == 0 ) return CALL_2ARGS( QUO_FFE_LARGE, opL, opR );
1076 qX = SIZE_FF( fX );
1077 /* if ( vL != 0 ) vL = (qX-1) / (qL-1) * (vL-1) + 1; */
1078 if ( vL != 0 ) vL = ((qX-1) * (vL-1)) / (qL-1) + 1;
1079 /* if ( vR != 0 ) vR = (qX-1) / (qR-1) * (vR-1) + 1; */
1080 if ( vR != 0 ) vR = ((qX-1) * (vR-1)) / (qR-1) + 1;
1081 }
1082
1083 /* compute and return the result */
1084 if ( vR == 0 ) {
1085 ErrorMayQuit("FFE operations: <divisor> must not be zero", 0, 0);
1086 }
1087 vX = QUO_FFV( vL, vR, SUCC_FF(fX) );
1088 return NEW_FFE( fX, vX );
1089 }
1090
QuoFFEInt(Obj opL,Obj opR)1091 static Obj QuoFFEInt(Obj opL, Obj opR)
1092 {
1093 FFV vL, vR, vX; /* value of left, right, result */
1094 FF fX; /* field of result */
1095 Int pX; /* char. of result */
1096 const FFV* sX; /* successor table of result field */
1097
1098 /* get the field for the result */
1099 fX = FLD_FFE( opL );
1100 pX = CHAR_FF( fX );
1101 sX = SUCC_FF( fX );
1102
1103 /* get the right operand */
1104 vX = ((INT_INTOBJ( opR ) % pX) + pX) % pX;
1105 if ( vX == 0 ) {
1106 vR = 0;
1107 }
1108 else {
1109 vR = 1;
1110 for ( ; 1 < vX; vX-- ) vR = sX[vR];
1111 }
1112
1113 /* get the left operand */
1114 vL = VAL_FFE( opL );
1115
1116 /* compute and return the result */
1117 if ( vR == 0 ) {
1118 ErrorMayQuit("FFE operations: <divisor> must not be zero", 0, 0);
1119 }
1120 vX = QUO_FFV( vL, vR, sX );
1121 return NEW_FFE( fX, vX );
1122 }
1123
QuoIntFFE(Obj opL,Obj opR)1124 static Obj QuoIntFFE(Obj opL, Obj opR)
1125 {
1126 FFV vL, vR, vX; /* value of left, right, result */
1127 FF fX; /* field of result */
1128 Int pX; /* char. of result */
1129 const FFV* sX; /* successor table of result field */
1130
1131 /* get the field for the result */
1132 fX = FLD_FFE( opR );
1133 pX = CHAR_FF( fX );
1134 sX = SUCC_FF( fX );
1135
1136 /* get the left operand */
1137 vX = ((INT_INTOBJ( opL ) % pX) + pX) % pX;
1138 if ( vX == 0 ) {
1139 vL = 0;
1140 }
1141 else {
1142 vL = 1;
1143 for ( ; 1 < vX; vX-- ) vL = sX[vL];
1144 }
1145
1146 /* get the right operand */
1147 vR = VAL_FFE( opR );
1148
1149 /* compute and return the result */
1150 if ( vR == 0 ) {
1151 ErrorMayQuit("FFE operations: <divisor> must not be zero", 0, 0);
1152 }
1153 vX = QUO_FFV( vL, vR, sX );
1154 return NEW_FFE( fX, vX );
1155 }
1156
1157
1158 /****************************************************************************
1159 **
1160 *F PowFFEInt(<opL>,<opR>) . . . . . . . . . power of a finite field element
1161 **
1162 ** 'PowFFEInt' returns the power of the finite field element <opL> and the
1163 ** integer <opR>. The power is represented over the field over which the
1164 ** left operand is represented, even if it lies in a much smaller field.
1165 **
1166 ** 'PowFFEInt' just does the conversions mentioned above and then calls the
1167 ** macro 'POW_FFV' to do the actual exponentiation.
1168 */
PowFFEInt(Obj opL,Obj opR)1169 static Obj PowFFEInt(Obj opL, Obj opR)
1170 {
1171 FFV vL, vX; /* value of left, result */
1172 Int vR; /* value of right */
1173 FF fX; /* field of result */
1174 const FFV* sX; /* successor table of result field */
1175
1176 /* get the field for the result */
1177 fX = FLD_FFE( opL );
1178 sX = SUCC_FF( fX );
1179
1180 /* get the right operand */
1181 vR = INT_INTOBJ( opR );
1182
1183 /* get the left operand */
1184 vL = VAL_FFE( opL );
1185
1186 /* if the exponent is negative, invert the left operand */
1187 if ( vR < 0 ) {
1188 if ( vL == 0 ) {
1189 ErrorMayQuit("FFE operations: <divisor> must not be zero", 0, 0);
1190 }
1191 vL = QUO_FFV( 1, vL, sX );
1192 vR = -vR;
1193 }
1194
1195 /* catch the case when vL is zero. */
1196 if( vL == 0 ) return NEW_FFE( fX, (vR == 0 ? 1 : 0 ) );
1197
1198 /* reduce vR modulo the order of the multiplicative group first. */
1199 vR %= *sX;
1200
1201 /* compute and return the result */
1202 vX = POW_FFV( vL, vR, sX );
1203 return NEW_FFE( fX, vX );
1204 }
1205
1206
1207 /****************************************************************************
1208 **
1209 *F PowFFEFFE( <opL>, <opR> ) . . . . . . conjugate of a finite field element
1210 */
PowFFEFFE(Obj opL,Obj opR)1211 static Obj PowFFEFFE(Obj opL, Obj opR)
1212 {
1213 /* get the field for the result */
1214 if ( CHAR_FF( FLD_FFE(opL) ) != CHAR_FF( FLD_FFE(opR) ) ) {
1215 ErrorMayQuit("<x> and <y> have different characteristic", 0, 0);
1216 }
1217
1218 /* compute and return the result */
1219 return opL;
1220 }
1221
1222
1223 /****************************************************************************
1224 **
1225 *F FiltIS_FFE( <self>, <obj> ) . . . . . . . test for finite field elements
1226 **
1227 ** 'FuncIsFFE' implements the internal function 'IsFFE( <obj> )'.
1228 **
1229 ** 'IsFFE' returns 'true' if its argument <obj> is a finite field element
1230 ** and 'false' otherwise. 'IsFFE' will cause an error if called with an
1231 ** unbound variable.
1232 */
1233 static Obj IsFFEFilt;
1234
FiltIS_FFE(Obj self,Obj obj)1235 static Obj FiltIS_FFE(Obj self, Obj obj)
1236 {
1237 /* return 'true' if <obj> is a finite field element */
1238 if ( IS_FFE(obj) ) {
1239 return True;
1240 }
1241 else if ( TNUM_OBJ(obj) < FIRST_EXTERNAL_TNUM ) {
1242 return False;
1243 }
1244 else {
1245 return DoFilter( self, obj );
1246 }
1247 }
1248
1249
1250 /****************************************************************************
1251 **
1252 *F FuncLOG_FFE_DEFAULT( <self>, <opZ>, <opR> ) . logarithm of a ff constant
1253 **
1254 ** 'FuncLOG_FFE_DEFAULT' implements the function 'LogFFE( <z>, <r> )'.
1255 **
1256 ** 'LogFFE' returns the logarithm of the nonzero finite field element <z>
1257 ** with respect to the root <r> which must lie in the same field like <z>.
1258 */
1259 static Obj LOG_FFE_LARGE;
1260
FuncLOG_FFE_DEFAULT(Obj self,Obj opZ,Obj opR)1261 static Obj FuncLOG_FFE_DEFAULT(Obj self, Obj opZ, Obj opR)
1262 {
1263 FFV vZ, vR; /* value of left, right */
1264 FF fZ, fR, fX; /* field of left, right, common */
1265 UInt qZ, qR, qX; /* size of left, right, common */
1266 Int a, b, c, d, t; /* temporaries */
1267
1268 /* check the arguments */
1269 if (!IS_FFE(opZ) || VAL_FFE(opZ) == 0) {
1270 ErrorMayQuit("LogFFE: <z> must be a nonzero finite field element", 0,
1271 0);
1272 }
1273 if (!IS_FFE(opR) || VAL_FFE(opR) == 0) {
1274 ErrorMayQuit("LogFFE: <r> must be a nonzero finite field element", 0,
1275 0);
1276 }
1277
1278 /* get the values, handle trivial cases */
1279 vZ = VAL_FFE( opZ );
1280 vR = VAL_FFE( opR );
1281
1282 /* bring the two operands into a common field <fX> */
1283 fZ = FLD_FFE( opZ );
1284 qZ = SIZE_FF( fZ );
1285 fR = FLD_FFE( opR );
1286 qR = SIZE_FF( fR );
1287
1288 if ( qZ == qR ) {
1289 fX = fZ;
1290 qX = qZ;
1291 }
1292 else if ( qZ % qR == 0 && (qZ-1) % (qR-1) == 0 ) {
1293 fX = fZ;
1294 qX = qZ;
1295 if ( vR != 0 ) vR = (qZ-1) / (qR-1) * (vR-1) + 1;
1296 }
1297 else if ( qR % qZ == 0 && (qR-1) % (qZ-1) == 0 ) {
1298 fX = fR;
1299 qX = qR;
1300 if ( vZ != 0 ) vZ = (qR-1) / (qZ-1) * (vZ-1) + 1;
1301 }
1302 else {
1303 fX = CommonFF( fZ, DegreeFFE(opZ), fR, DegreeFFE(opR) );
1304 if ( fX == 0 ) return CALL_2ARGS( LOG_FFE_LARGE, opZ, opR );
1305 qX = SIZE_FF( fX );
1306 /* if ( vZ != 0 ) vZ = (qX-1) / (qZ-1) * (vZ-1) + 1; */
1307 if ( vZ != 0 ) vZ = ((qX-1) * (vZ-1)) / (qZ-1) + 1;
1308 /* if ( vR != 0 ) vR = (qX-1) / (qR-1) * (vR-1) + 1; */
1309 if ( vR != 0 ) vR = ((qX-1) * (vR-1)) / (qR-1) + 1;
1310 }
1311
1312 /* now solve <l> * (<vR>-1) = (<vZ>-1) % (<qX>-1) */
1313 a = 1; b = 0;
1314 c = (Int) (vR-1); d = (Int) (qX-1);
1315 while ( d != 0 ) {
1316 t = b; b = a - (c/d) * b; a = t;
1317 t = d; d = c - (c/d) * d; c = t;
1318 }
1319 if ( ((Int) (vZ-1)) % c != 0 ) {
1320 return Fail;
1321 }
1322
1323
1324 while (a < 0)
1325 a+= (qX -1)/c;
1326
1327 // return the logarithm
1328 return INTOBJ_INT( (((UInt) (vZ-1) / c) * a) % ((UInt) (qX-1)) );
1329 }
1330
1331
1332 /****************************************************************************
1333 **
1334 *F FuncINT_FFE_DEFAULT( <self>, <z> ) . . . . convert a ffe to an integer
1335 **
1336 ** 'FuncINT_FFE_DEFAULT' implements the internal function 'IntFFE( <z> )'.
1337 **
1338 ** 'IntFFE' returns the integer that corresponds to the finite field
1339 ** element <z>, which must of course be an element of a prime field, i.e.,
1340 ** the smallest integer <i> such that '<i> * <z>^0 = <z>'.
1341 */
1342 static Obj IntFF;
1343 #ifdef HPCGAP
1344 static Int NumFF;
1345 #endif
1346
INT_FF(FF ff)1347 static Obj INT_FF(FF ff)
1348 {
1349 Obj conv; /* conversion table, result */
1350 Int q; /* size of finite field */
1351 Int p; /* char of finite field */
1352 const FFV * succ; /* successor table of finite field */
1353 FFV z; /* one element of finite field */
1354 UInt i; /* loop variable */
1355
1356 /* if the conversion table is not already known, construct it */
1357 #ifdef HPCGAP
1358 if ( NumFF < ff || (MEMBAR_READ(), ATOMIC_ELM_PLIST(IntFF, ff) == 0)) {
1359 HashLock(&IntFF);
1360 #else
1361 if ( LEN_PLIST(IntFF) < ff || ELM_PLIST(IntFF,ff) == 0 ) {
1362 #endif
1363 q = SIZE_FF( ff );
1364 p = CHAR_FF( ff );
1365 conv = NEW_PLIST_IMM( T_PLIST, p-1 );
1366 succ = SUCC_FF( ff );
1367 SET_LEN_PLIST( conv, p-1 );
1368 z = 1;
1369 for ( i = 1; i < p; i++ ) {
1370 SET_ELM_PLIST( conv, (z-1)/((q-1)/(p-1))+1, INTOBJ_INT(i) );
1371 z = succ[ z ];
1372 }
1373 #ifdef HPCGAP
1374 GrowPlist(IntFF, ff);
1375 ATOMIC_SET_ELM_PLIST( IntFF, ff, conv );
1376 MEMBAR_WRITE();
1377 NumFF = LEN_PLIST(IntFF);
1378 HashUnlock(&IntFF);
1379 #else
1380 AssPlist( IntFF, ff, conv );
1381 #endif
1382 }
1383
1384 /* return the conversion table */
1385 #ifdef HPCGAP
1386 return ATOMIC_ELM_PLIST( IntFF, ff);
1387 #else
1388 return ELM_PLIST( IntFF, ff );
1389 #endif
1390 }
1391
1392
1393 static Obj FuncINT_FFE_DEFAULT(Obj self, Obj z)
1394 {
1395 FFV v; /* value of finite field element */
1396 FF ff; /* finite field */
1397 Int q; /* size of finite field */
1398 Int p; /* char of finite field */
1399 Obj conv; /* conversion table */
1400
1401 /* get the value */
1402 v = VAL_FFE( z );
1403
1404 /* special case for 0 */
1405 if ( v == 0 ) {
1406 return INTOBJ_INT( 0 );
1407 }
1408
1409 /* get the field, size, characteristic, and conversion table */
1410 ff = FLD_FFE( z );
1411 q = SIZE_FF( ff );
1412 p = CHAR_FF( ff );
1413 conv = INT_FF( ff );
1414
1415 /* check the argument */
1416 if ( (v-1) % ((q-1)/(p-1)) != 0 ) {
1417 ErrorMayQuit("IntFFE: <z> must lie in prime field", 0, 0);
1418 }
1419
1420 /* convert the value into the prime field */
1421 v = (v-1) / ((q-1)/(p-1)) + 1;
1422
1423 /* return the integer value */
1424 return ELM_PLIST( conv, v );
1425 }
1426
1427
1428 /****************************************************************************
1429 **
1430 *F FuncZ( <self>, <q> ) . . . return the generator of a small finite field
1431 **
1432 ** 'FuncZ' implements the internal function 'Z( <q> )'.
1433 **
1434 ** 'Z' returns the generators of the small finite field with <q> elements.
1435 ** <q> must be a positive prime power.
1436 */
1437 static Obj ZOp;
1438
1439 static Obj FuncZ(Obj self, Obj q)
1440 {
1441 FF ff; /* the finite field */
1442
1443 /* check the argument */
1444 if ( (IS_INTOBJ(q) && (INT_INTOBJ(q) > 65536)) ||
1445 (TNUM_OBJ(q) == T_INTPOS))
1446 return CALL_1ARGS(ZOp, q);
1447
1448 if ( !IS_INTOBJ(q) || INT_INTOBJ(q)<=1 ) {
1449 RequireArgument("Z", q, "must be a positive prime power");
1450 }
1451
1452 ff = FiniteFieldBySize(INT_INTOBJ(q));
1453
1454 if (!ff) {
1455 RequireArgument("Z", q, "must be a positive prime power");
1456 }
1457
1458 /* make the root */
1459 return NEW_FFE(ff, (q == INTOBJ_INT(2)) ? 1 : 2);
1460 }
1461
1462 static Obj FuncZ2(Obj self, Obj p, Obj d)
1463 {
1464 FF ff;
1465 Int ip, id, id1;
1466 UInt q;
1467 if (ARE_INTOBJS(p, d)) {
1468 ip = INT_INTOBJ(p);
1469 id = INT_INTOBJ(d);
1470 if (ip > 1 && id > 0 && id <= 16 && ip < 65536) {
1471 id1 = id;
1472 q = ip;
1473 while (--id1 > 0 && q <= 65536)
1474 q *= ip;
1475 if (q <= 65536) {
1476 /* get the finite field */
1477 ff = FiniteField(ip, id);
1478
1479 if (ff == 0 || CHAR_FF(ff) != ip)
1480 RequireArgument("Z", p, "must be a prime");
1481
1482 /* make the root */
1483 return NEW_FFE(ff, (ip == 2 && id == 1 ? 1 : 2));
1484 }
1485 }
1486 }
1487 return CALL_2ARGS(ZOp, p, d);
1488 }
1489
1490
1491 /****************************************************************************
1492 **
1493 *F * * * * * * * * * * * * * initialize module * * * * * * * * * * * * * * *
1494 */
1495
1496
1497 /****************************************************************************
1498 **
1499 *V BagNames . . . . . . . . . . . . . . . . . . . . . . . list of bag names
1500 */
1501 static StructBagNames BagNames[] = {
1502 { T_FFE, "ffe" },
1503 { -1, "" }
1504 };
1505
1506
1507 /****************************************************************************
1508 **
1509 *V GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export
1510 */
1511 static StructGVarFilt GVarFilts [] = {
1512
1513 GVAR_FILT(IS_FFE, "obj", &IsFFEFilt),
1514 { 0, 0, 0, 0, 0 }
1515
1516 };
1517
1518
1519 /****************************************************************************
1520 **
1521 *V GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
1522 */
1523 static StructGVarFunc GVarFuncs [] = {
1524
1525 GVAR_FUNC(CHAR_FFE_DEFAULT, 1, "z"),
1526 GVAR_FUNC(DEGREE_FFE_DEFAULT, 1, "z"),
1527 GVAR_FUNC(LOG_FFE_DEFAULT, 2, "z, root"),
1528 GVAR_FUNC(INT_FFE_DEFAULT, 1, "z"),
1529 GVAR_FUNC(Z, 1, "q"),
1530 { 0, 0, 0, 0, 0 }
1531
1532 };
1533
1534
1535 /****************************************************************************
1536 **
1537 *F InitKernel( <module> ) . . . . . . . . initialise kernel data structures
1538 */
1539 static Int InitKernel (
1540 StructInitInfo * module )
1541 {
1542 // set the bag type names (for error messages and debugging)
1543 InitBagNamesFromTable( BagNames );
1544
1545 /* install the type functions */
1546 ImportFuncFromLibrary( "TYPE_FFE", &TYPE_FFE );
1547 ImportFuncFromLibrary( "TYPE_FFE0", &TYPE_FFE0 );
1548 ImportFuncFromLibrary( "ZOp", &ZOp );
1549 InitFopyGVar( "PrimitiveRootMod", &PrimitiveRootMod );
1550 TypeObjFuncs[ T_FFE ] = TypeFFE;
1551
1552 /* create the fields and integer conversion bags */
1553 InitGlobalBag( &SuccFF, "src/finfield.c:SuccFF" );
1554 InitGlobalBag( &TypeFF, "src/finfield.c:TypeFF" );
1555 InitGlobalBag( &TypeFF0, "src/finfield.c:TypeFF0" );
1556 InitGlobalBag( &IntFF, "src/finfield.c:IntFF" );
1557
1558 /* install the functions that handle overflow */
1559 ImportFuncFromLibrary( "SUM_FFE_LARGE", &SUM_FFE_LARGE );
1560 ImportFuncFromLibrary( "DIFF_FFE_LARGE", &DIFF_FFE_LARGE );
1561 ImportFuncFromLibrary( "PROD_FFE_LARGE", &PROD_FFE_LARGE );
1562 ImportFuncFromLibrary( "QUO_FFE_LARGE", &QUO_FFE_LARGE );
1563 ImportFuncFromLibrary( "LOG_FFE_LARGE", &LOG_FFE_LARGE );
1564
1565 /* init filters and functions */
1566 InitHdlrFiltsFromTable( GVarFilts );
1567 InitHdlrFuncsFromTable( GVarFuncs );
1568 InitHandlerFunc( FuncZ2, "src/finfield.c: Z (2 args)");
1569
1570
1571 /* install the printing method */
1572 PrintObjFuncs[ T_FFE ] = PrFFE;
1573
1574 /* install the comparison methods */
1575 EqFuncs[ T_FFE ][ T_FFE ] = EqFFE;
1576 LtFuncs[ T_FFE ][ T_FFE ] = LtFFE;
1577
1578 /* install the arithmetic methods */
1579 ZeroFuncs[ T_FFE ] = ZeroFFE;
1580 ZeroMutFuncs[ T_FFE ] = ZeroFFE;
1581 AInvFuncs[ T_FFE ] = AInvFFE;
1582 AInvMutFuncs[ T_FFE ] = AInvFFE;
1583 OneFuncs [ T_FFE ] = OneFFE;
1584 OneMutFuncs [ T_FFE ] = OneFFE;
1585 InvFuncs [ T_FFE ] = InvFFE;
1586 InvMutFuncs [ T_FFE ] = InvFFE;
1587 SumFuncs[ T_FFE ][ T_FFE ] = SumFFEFFE;
1588 SumFuncs[ T_FFE ][ T_INT ] = SumFFEInt;
1589 SumFuncs[ T_INT ][ T_FFE ] = SumIntFFE;
1590 DiffFuncs[ T_FFE ][ T_FFE ] = DiffFFEFFE;
1591 DiffFuncs[ T_FFE ][ T_INT ] = DiffFFEInt;
1592 DiffFuncs[ T_INT ][ T_FFE ] = DiffIntFFE;
1593 ProdFuncs[ T_FFE ][ T_FFE ] = ProdFFEFFE;
1594 ProdFuncs[ T_FFE ][ T_INT ] = ProdFFEInt;
1595 ProdFuncs[ T_INT ][ T_FFE ] = ProdIntFFE;
1596 QuoFuncs[ T_FFE ][ T_FFE ] = QuoFFEFFE;
1597 QuoFuncs[ T_FFE ][ T_INT ] = QuoFFEInt;
1598 QuoFuncs[ T_INT ][ T_FFE ] = QuoIntFFE;
1599 PowFuncs[ T_FFE ][ T_INT ] = PowFFEInt;
1600 PowFuncs[ T_FFE ][ T_FFE ] = PowFFEFFE;
1601
1602 /* return success */
1603 return 0;
1604 }
1605
1606
1607 /****************************************************************************
1608 **
1609 *F InitLibrary( <module> ) . . . . . . . initialise library data structures
1610 */
1611 static Int InitLibrary (
1612 StructInitInfo * module )
1613 {
1614 /* create the fields and integer conversion bags */
1615 SuccFF = NEW_PLIST( T_PLIST, NUM_SHORT_FINITE_FIELDS );
1616 SET_LEN_PLIST( SuccFF, NUM_SHORT_FINITE_FIELDS );
1617
1618 TypeFF = NEW_PLIST( T_PLIST, NUM_SHORT_FINITE_FIELDS );
1619 SET_LEN_PLIST( TypeFF, NUM_SHORT_FINITE_FIELDS );
1620
1621 TypeFF0 = NEW_PLIST( T_PLIST, NUM_SHORT_FINITE_FIELDS );
1622 SET_LEN_PLIST( TypeFF0, NUM_SHORT_FINITE_FIELDS );
1623
1624 IntFF = NEW_PLIST( T_PLIST, NUM_SHORT_FINITE_FIELDS );
1625 SET_LEN_PLIST( IntFF, NUM_SHORT_FINITE_FIELDS );
1626
1627 #ifdef HPCGAP
1628 MakeBagPublic(SuccFF);
1629 MakeBagPublic(TypeFF);
1630 MakeBagPublic(TypeFF0);
1631 MakeBagPublic(IntFF);
1632 #endif
1633
1634 /* init filters and functions */
1635 InitGVarFiltsFromTable( GVarFilts );
1636 InitGVarFuncsFromTable( GVarFuncs );
1637 SET_HDLR_FUNC(ValGVar(GVarName("Z")), 2, FuncZ2);
1638
1639 /* return success */
1640 return 0;
1641 }
1642
1643
1644 /****************************************************************************
1645 **
1646 *F InitInfoFinfield() . . . . . . . . . . . . . . . table of init functions
1647 */
1648 static StructInitInfo module = {
1649 // init struct using C99 designated initializers; for a full list of
1650 // fields, please refer to the definition of StructInitInfo
1651 .type = MODULE_BUILTIN,
1652 .name = "finfield",
1653 .initKernel = InitKernel,
1654 .initLibrary = InitLibrary,
1655 };
1656
1657 StructInitInfo * InitInfoFinfield ( void )
1658 {
1659 return &module;
1660 }
1661