From 4acdb92b6c4b0a52f84be62ee70663a6e751ade4 Mon Sep 17 00:00:00 2001 From: Steve Linton Date: Tue, 13 Nov 2018 17:24:38 +0000 Subject: [PATCH] Fix some really stupid performance issues which become much more visible with bigger fields --- src/finfield.c | 186 +++++++++++++++++++++++++------------------------ 1 file changed, 96 insertions(+), 90 deletions(-) diff --git a/src/finfield.c b/src/finfield.c index 037781890a..2de6d3e67f 100644 --- a/src/finfield.c +++ b/src/finfield.c @@ -93,54 +93,69 @@ static Obj TYPE_FFE0; // used for successor bags static Obj TYPE_KERNEL_OBJECT; + +// used to call out to GAP static Obj PrimitiveRootMod; +/**************************************************************************** + ** + *F lookupPrimePower ( q ) . . . . . .search for a prime power in tables + ** + ** searches the tables of prime powers from ffdata.c for q, returns the index + ** of q in SizeFF if it is present and 0 if not (the 0 position of SizeFF is + ** unused). + */ + +static FF lookupPrimePower(UInt q) +{ + UInt l,n; + FF ff; + UInt e; + /* search through the finite field table */ + l = 1; + n = NUM_SHORT_FINITE_FIELDS; + ff = 0; + while (l <= n && SizeFF[l] <= q && q <= SizeFF[n]) { + /* interpolation search */ + /* cuts iterations roughly in half compared to binary search at + * the expense of additional divisions. */ + e = (q - SizeFF[l] + 1) * (n - l) / (SizeFF[n] - SizeFF[l] + 1); + ff = l + e; + if (SizeFF[ff] == q) + break; + if (SizeFF[ff] < q) + l = ff + 1; + else + n = ff - 1; + } + if (ff < 1 || ff > NUM_SHORT_FINITE_FIELDS) + return 0; + if (SizeFF[ff] != q) + return 0; + return ff; +} + /**************************************************************************** ** *F FiniteField(

,) . . . make the small finite field with elements ** ** 'FiniteField' returns the small finite field with

^ elements. */ -FF FiniteField ( - UInt p, - UInt d ) + +FF FiniteFieldBySize ( UInt q) { FF ff; /* finite field, result */ Obj tmp; /* temporary bag */ Obj succBag; /* successor table bag */ FFV * succ; /* successor table */ FFV * indx; /* index table */ - UInt q; /* size of finite field */ UInt poly; /* Conway polynomial of extension */ UInt i, l, f, n, e; /* loop variables */ Obj root; /* will be a primitive root mod p */ - /* calculate size of field */ - q = 1; - for ( i = 1; i <= d; i++ ) q *= p; + + ff = lookupPrimePower(q); - /* search through the finite field table */ - l = 1; n = NUM_SHORT_FINITE_FIELDS; - ff = 0; - while (l <= n && SizeFF[l] <= q && q <= SizeFF[n]) { - /* interpolation search */ - /* cuts iterations roughly in half compared to binary search at - * the expense of additional divisions. */ - e = (q - SizeFF[l]+1) * (n-l) / (SizeFF[n]-SizeFF[l]+1); - ff = l + e; - if (SizeFF[ff] == q) - break; - if (SizeFF[ff] < q) - l = ff+1; - else - n = ff-1; - } - if (ff < 1 || ff > NUM_SHORT_FINITE_FIELDS) - return 0; - if (CharFF[ff] != p) - return 0; - if (SizeFF[ff] != q) - return 0; #ifdef HPCGAP /* Important correctness concern here: * @@ -177,21 +192,22 @@ FF FiniteField ( /* if q is a prime find the smallest primitive root $e$, use $x - e$ */ - if ( d == 1 ) { - if ( p < 65537 ) { + if ( DEGR_FF(ff) == 1 ) { + if ( q < 65537 ) { /* for smaller primes we do this in the kernel for performance and bootstrapping - reasons */ - for ( e = 1, i = 1; i != p-1; ++e ) { + reasons + TODO -- review the threshold */ + for ( e = 1, i = 1; i != q-1; ++e ) { for ( f = e, i = 1; f != 1; ++i ) - f = (f * e) % p; + f = (f * e) % q; } } else { /* Otherwise we ask the library */ - root=CALL_1ARGS(PrimitiveRootMod, INTOBJ_INT(p)); + root=CALL_1ARGS(PrimitiveRootMod, INTOBJ_INT(q)); e = INT_INTOBJ(root)+1; } - poly = p-(e-1); + poly = q-(e-1); } /* otherwise look up the polynomial used to construct this field */ @@ -202,6 +218,7 @@ FF FiniteField ( /* construct 'indx' such that 'e = x^(indx[e]-1) % poly' for every e */ indx[ 0 ] = 0; + UInt p = CHAR_FF(ff); for ( e = 1, n = 0; n < q-1; ++n ) { indx[ e ] = n + 1; /* e =p*e mod poly =x*e mod poly =x*x^n mod poly =x^{n+1} mod poly */ @@ -254,6 +271,20 @@ FF FiniteField ( return ff; } +FF FiniteField ( + UInt p, + UInt d ) +{ + UInt q,i; + /* calculate size of field */ + q = 1; + for ( i = 1; i <= d; i++ ) + q *= p; + return FiniteFieldBySize(q); +} + + + /**************************************************************************** ** @@ -1487,9 +1518,6 @@ Obj FuncZ ( Obj q ) { FF ff; /* the finite field */ - UInt p; /* characteristic */ - UInt d; /* degree */ - UInt r; /* temporary */ /* check the argument */ if ( (IS_INTOBJ(q) && (INT_INTOBJ(q) > SIZE_LARGEST_INTERNAL_FF)) || @@ -1504,66 +1532,44 @@ Obj FuncZ ( return FuncZ( self, q ); } - /* compute the prime and check that is a prime power */ - if ( INT_INTOBJ(q) % 2 == 0 ) { - p = 2; - } - else { - p = 3; - while ( INT_INTOBJ(q) % p != 0 ) { - p += 2; - } - } - d = 1; - r = p; - while ( r < INT_INTOBJ(q) ) { - d = d + 1; - r = r * p; - } - if ( r != INT_INTOBJ(q) ) { - q = ErrorReturnObj( - "Z: must be a positive prime power (not a %s)", - (Int)TNAM_OBJ(q), 0L, - "you can replace via 'return ;'" ); - return FuncZ( self, q ); + ff = FiniteFieldBySize(INT_INTOBJ(q)); + + if (!ff) { + ErrorMayQuit("Z: must be a positive prime power (not a %s)", + (Int)TNAM_OBJ(q), 0L); } - /* get the finite field */ - ff = FiniteField( p, d ); - /* make the root */ - return NEW_FFE( ff, (p == 2 && d == 1 ? 1 : 2) ); + return NEW_FFE( ff, (q == INTOBJ_INT(2)) ? 1 : 2 ); } -Obj FuncZ2 ( Obj self, Obj p, Obj d) +Obj FuncZ2(Obj self, Obj p, Obj d) { - FF ff; - Int ip,id,id1; - UInt q; - if (ARE_INTOBJS(p,d)) - { - ip = INT_INTOBJ(p); - id = INT_INTOBJ(d); - if (ip > 1 && id > 0 && id <= DEGREE_LARGEST_INTERNAL_FF && ip <= SIZE_LARGEST_INTERNAL_FF) - { - id1 = id; - q = ip; - while (--id1 > 0 && q <= SIZE_LARGEST_INTERNAL_FF) - q *= ip; - if (q <= SIZE_LARGEST_INTERNAL_FF) - { - /* get the finite field */ - ff = FiniteField( ip, id ); - - if ( ff == 0 || CHAR_FF(ff) != ip ) - ErrorMayQuit("Z:

must be a prime", 0, 0); - - /* make the root */ - return NEW_FFE( ff, (ip == 2 && id == 1 ? 1 : 2) ); + FF ff; + Int ip, id, id1; + UInt q; + if (ARE_INTOBJS(p, d)) { + ip = INT_INTOBJ(p); + id = INT_INTOBJ(d); + if (ip > 1 && id > 0 && id <= DEGREE_LARGEST_INTERNAL_FF && + ip <= SIZE_LARGEST_INTERNAL_FF) { + id1 = id; + q = ip; + while (--id1 > 0 && q <= SIZE_LARGEST_INTERNAL_FF) + q *= ip; + if (q <= SIZE_LARGEST_INTERNAL_FF) { + /* get the finite field */ + ff = FiniteFieldBySize(q); + + if (ff == 0 || CHAR_FF(ff) != ip) + ErrorMayQuit("Z:

must be a prime", 0, 0); + + /* make the root */ + return NEW_FFE( ff, (q == 2) ? 1 : 2 ); } } } - return CALL_2ARGS(ZOp, p, d); + return CALL_2ARGS(ZOp, p, d); }