share/contrib/gf/gf.mac

   1 /*
   2 ------------------------------------------------------------------------
   3 Efficient Galois Fields in Maxima
   4
   5 by Alasdair McAndrew
   6 and later extended by Fabrizio Caruso and Jacopo Daurizio
   7
   8 it is distributed together with gf_roots by Jacopo Daurizio
   9
  10 Authors:
  11
  12 Fabrizio Caruso   (optimizations, testing)
  13 Jacopo D'Aurizio   (optimizations, modular roots)
  14 Alasdair McAndrew (original version of the package, pohlig-helman log, etc... )
  15 ------------------------------------------------------------------------*/
  16
  17 /* Released under terms of the GNU General Public License, version 2,
  18  * by permission of the authors to Robert Dodier circa 2007-12-02.
  19  */
  20
  21 /* Revision by Volker van Nek, 2012
  22  *
  23  * A lot of functions moved to src/numth.lisp.
  24  * lucas moved to src/combin.lisp next to fib.
  25  *
  26  * Remaining contents:
  27  *  1. Functions for Z_p
  28  *  2. Galois Fields (Itoh-Tsujii algorithm)
  29  *  3. Collected comments and algorithms for number theory
  30  *  4. Roots in finite fields (documented in gf_manual.pdf)
  31  */
  32
  33 /* *** Functions for Z_p, p prime ******************************************* */
  34
  35 /* Creates an array of couples [a!, 1/(a!)] mod p
  36 */
  37 zp_fact_array(p) := block([v2],
  38    v2 : make_array('fixnum, p, 2),
  39    for i:0 thru p-2 do (v2[i,0] : 1, v2[i,1] : 1),
  40    v2[p-1, 0] : -1,
  41    v2[p-1, 1] : -1,
  42    if p > 3 then (
  43       for i:2 thru p-1 do v2[i,0] : mod(i * v2[i-1,0], p),
  44       for i:2 thru p-2 do v2[i,1] : mod((-1)^(p-i) * v2[p-i-1,0], p) ),
  45    v2 )$
  46
  47
  48 /* Calculates the product a(a+1)(a+2)...(b-1)b in Z_p
  49 */
  50 zp_prod(a, b, p) :=
  51    if a = 1 and b = 0 then 1
  52    else if a = b then mod(a, p)
  53    else if (a = 1 or a = 2) and b = p-1 then b
  54    else block([s:1],
  55       for i:a thru b do s : mod(s*i, p),
  56       s )$
  57
  58 /* Calculates (n!)^(-1) in Z_p using Wilson Theorem
  59    This algorithm DOES NOT require explicit inversion through gcd()
  60 */
  61 zp_inv_fact(n, p) :=
  62    if n = 0 or n = 1 then 1
  63    else if n = p-1 then n
  64    else mod(- zp_prod(n+1, p-1, p), p)$
  65
  66
  67 /* Calculates the maximum exponent s for which p^s divides n!
  68 */
  69 zp_max_pow(n, p) :=
  70    if n < p then 0
  71    else block([k, s:0],
  72       k : length_base(p, n) - 1,
  73       for i:1 thru k do s : s + floor(n/p^i),
  74       s )$
  75
  76 length_base(b, n) := block([len:0],
  77    while n # 0 do (
  78       n : quotient(n, b),
  79       len : len + 1 ),
  80    len )$
  81
  82
  83 /* Calculates binomial(m,n) in Z_p using Wilson Theorem and fast inversions of k!
  84 */
  85 zp_binomial(m, n, p) :=
  86    if m < n then 0
  87    else (
  88       if m < 2*n then n : m-n,
  89       if n = 0 then 1
  90       else if n = 1 then mod(m, p)
  91       else if zp_max_pow(m, p) - zp_max_pow(m-n, p) > zp_max_pow(n,p) then 0
  92       else block([a,b,c, na,nb,nc],
  93          a : mod(n, p),   na : mod(floor(n/p), 2),
  94          b : mod(m, p),   nb : mod(floor(m/p), 2),
  95          c : mod(m-n, p), nc : mod(floor((m-n)/p), 2),
  96          mod(zp_prod(1, b, p) * zp_inv_fact(a, p) * zp_inv_fact(c, p) * (-1)^(na+nb+nc), p) ))$
  97
  98
  99
 100 /* *** Galois Fields ******************************************************** */
 101
 102 /* Inversion through Itoh-Tsujii algorithm -- Highly Experimental
 103 */
 104 gf_inv_ITA(poly) := block([p, r, s, t, u, gf_rat:false],
 105    p : gf_characteristic(),
 106    r : gf_order()/(p - 1),
 107    s : gf_exp(poly, r - 1),
 108    t : gf_mult(poly, s),
 109    u : inv_mod(t, p),
 110    gf_mult(s, u) )$
 111
 112 /* Remark (VvN):
 113 A possible Lisp-implementation is listed below.
 114 However it is much slower than $gf_inv in src (e.g. 30 times in F3^102).
 115
 116 (defmfun $gf_inv_ita (a)
 117   (gf-red? "gf_inv_ita")
 118   (let ((*ef-arith?*))
 119     (setq a (gf-inv-ita (gf-p2x a) *gf-red*))
 120     (when a (gf-x2p a)) ))
 121
 122 (defun gf-inv-ita (y red)
 123   #+ (or ccl ecl gcl) (declare (optimize (speed 3) (safety 0)))
 124   (when (null y)
 125     (gf-merror (intl:gettext "gf arithmetic: Quotient by zero")) )
 126   (let (($gf_rat) r x c)
 127     (setq r (truncate *gf-ord* (1- *gf-char*))
 128           x (gf-pow y (1- r) red)
 129           y (gf-times x y red) ) ;; y is (0 n) or nil
 130     (unless y (return-from gf-inv-ita))
 131     (setq c (inv-mod (cadr y) *gf-char*))
 132     (when c (gf-xctimes x c)) ))
 133 */
 134
 135
 136 /*
 137 gf_dotprod(lst1, lst2) - cut out
 138 use
 139 matrix_element_add : gf_add
 140 matrix_element_mult : gf_mult
 141 and the dot-operator itself
 142
 143 gf_mateval(M) - cut out
 144 use matrixmap('gf_eval, M)
 145 */
 146
 147
 148 gf_eval_poly(poly, var) :=
 149   sum(gf_eval(ratcoeff(rat(poly, var), var, i)) * var^i, i, 0, hipow(poly, var) )$
 150
 151
 152
 153 /* *** Collected comments and algorithms ************************************ */
 154
 155 /* Complessita' : detta C la cardinalita' del campo, si hanno in media Osoft(C^(1/t)) operazioni,  */
 156 /* dove t=log(C)/log(q), dove q e' il piu' grande primo che divide C                               */
 157 /* diff1,diff2,ll,  --- allocare o non allocare? non cambia niente, sembra                         */
 158
 159
 160 /* Berlekamp-Massey Algorithm
 161 */
 162 berlekamp_massey(seq) := block([x, m, r,r0,r1, v,v0,v1, q, h, d, p],
 163    m : length(seq),
 164    v0 : rat(0, x),
 165    v1 : rat(1, x),
 166    r0 : rat(x^m, x),
 167    r1 : rat(sum(seq[i] * x^(i-1), i, 1, m), x),
 168    h : hipow(r1, x),
 169    while 2*h >= m do (
 170       [q, r] : divide(r0, r1, x),
 171       [v0, v1] : [v1, rat(v0 - q*v1, x)],
 172       [r0, r1] : [r1, r],
 173       h : hipow(r1, x) ),
 174    d : max(1 + h, hipow(v1, x)),
 175    p : rat((x^d) * ratsubst(1/x, x, v1), x),
 176    h : hipow(p, x),
 177    rat(p/coeff(p, x, h), x) )$
 178
 179
 180
 181 /* Binary search in a sorted [by rising values] binary list, still using
 182    ?integer\-length(length(list)) for accessing list elements and therefore
 183    less efficient than it could be but still utilizabile for speeding up
 184    the original gf_ind  */
 185 binsearch(elem, list) := block([l:length(list), l2, pos, k:1, check],
 186    l2 : ?integer\-length(l) - 1,
 187    pos : ?ash(2, l2 - 1),
 188    while k < l2 + 2 do (
 189       if pos > l then (
 190          pos : pos - ?ash(2, l2 - k) + ?ash(2, l2 - k - 1),
 191          k : k + 1 )
 192    else (
 193       check : list[pos],
 194       if check = elem then return(pos)
 195       else if check > elem then
 196           pos : pos - ?ash(2, l2 - k) + ?ash(2, l2 - k - 1)
 197       else if check < elem then
 198          pos : pos + ?ash(2, l2 - k - 1),
 199       k : k+1 )),
 200   if k = l2 + 2 then 0
 201   else pos )$
 202
 203
 204 /* Bijection maps between N and N^2 */
 205 n2n2(n) := block([a,b],
 206    a : ceiling((ceiling(sqrt(8*n + 9)) - 3) / 2),
 207    b : n - ?ash(a*(a + 1), -1),
 208    [b, a - b] )$
 209
 210 n22n(n2) := block([a, x:first(n2), y:last(n2)],
 211    a : x + y,
 212    ?ash(a*(a + 1), -1) + x )$
 213
 214 giretto(n, p) := block([start, list],
 215    start : [0,0],
 216    list : [start],
 217    for i:1 thru n do (
 218       start : map(lambda([x], mod(x, p)), n2n2(n22n(start) + 1)),
 219       list : cons(start, list),
 220       print(i, " : ", start) ),
 221    reverse(list) )$
 222
 223
 224 /* Fast square detection               */
 225 /* first(args(factor(x^2-n)))=x^2      */
 226 issquare(n) := block([prime:2, prod:2, k],
 227    if mod(n, 4) = 3 or mod(n, 8) = 5 then (
 228       print("Invalid residue class modulo 4 or 8"),
 229       return() ),
 230    while prod < ?ash(n, 4) do (
 231       if mod(n, prime) = 0 then (
 232          k : 0,
 233          while mod(n, prime) = 0 do (
 234             n : n/prime,
 235             k : k+1 ),
 236          if oddp(k) then (
 237             print("The maximum power of ", prime, " dividing ", n, " is odd."),
 238             return() )),
 239       if jacobi(n, prime) # 1 then (
 240          print(n, "fails Legendre test for p=", prime),
 241          return() ),
 242       prime : next_prime(prime),
 243       prod : prod*prime ),
 244    is( prod > ?ash(n,4) ))$
 245
 246
 247
 248 /* dopo quanto possono occorrere k numeri interi consecutivi tutti composti? */
 249 /* 2*3*5+1 = 31  32,33,34,35,36                                              */
 250 /* 2*3*5*7+1 = 211  212,213,214,215,216,217,218,219,220                      */
 251 /* da M=1+prod_{j=1..n} p_j  a   M+p_{n+1}-1                                 */
 252 /* p_k interi consecutivi compaiono almeno dopo prod_{j=1..k-1} j log(j)     */
 253 /* k log(k) dopo k^k e^{-k}                                                  */
 254 /* m = k log(k)  log(m) = log(k) + log(log(k))                               */
 255 /* log(k) = log( m/ log(k) ) ~ log(m) - log(log(m))                          */
 256 /* M interi consecutivi dopo almeno (M/(e log M))^M                          */
 257
 258 /* base next_prime(2^30), generatore 2, facile log(16384)                    */
 259
 260
 261
 262 /* ************************************************************************** */
 263 /* Roots in finite fields                                                     */
 264 /*                                                                            */
 265 /* Coded by Jacopo D'Aurizio                                                  */
 266
 267
 268 /* ------------------------------------------- */
 269 /*  Some number theoretic subroutines for F_p  */
 270 /* ------------------------------------------- */
 271
 272
 273 /* Low-complexity Lucas-Muller method for square roots
 274 */
 275 lmp(a, p, m) := block([d1:a, d2:a^2 - 2, j],
 276    j : ?integer\-length(m),
 277    for i:2 thru j-1 do (
 278       if ?logbitp(j - i, m) then (
 279          d1 : mod(d1*d2 - a, p),
 280          d2 : mod(d2^2 - 2, p) )
 281       else (
 282          d2 : mod(d1*d2 - a, p),
 283          d1 : mod(d1^2 - 2, p) )),
 284    if ?logbitp(0, m) then mod(d1*d2 - a, p)
 285    else mod(d1^2 - 2, p) )$
 286
 287
 288 gf_sqrt([ap]) :=
 289    if length(ap) = 1 then gsqrt(ap[1], gf_characteristic())
 290    else gsqrt(ap[1], ap[2]);
 291
 292
 293 gsqrt(a, p) := block([t:0, s:-4, term],
 294    a : mod(a, p),
 295    if a = 0 then return([0,0]),
 296    /* if not primep(p) then error("ERROR: Second argument must be a prime number."), */
 297    if jacobi(a, p) # 1 then
 298       error("ERROR: First argument must be a quadratic residue."),
 299    if p = 2 then return(mod(a, 2)),
 300    if mod(p, 4) = 3 then
 301       term : power_mod(a, ?ash(p + 1, -2), p)
 302    else (
 303       do (
 304          s : mod(s + a*(?ash(t, 1) + 1), p),
 305          t : t + 1,
 306          if jacobi(s, p) = -1 then return() ),
 307       term : mod(lmp(s + 2, p, ?ash(p - 1, -2)) * inv_mod(t, p), p) ),
 308    [term, p - term] )$
 309
 310
 311 gf_cbrt([ap]) :=
 312    if length(ap) = 1 then gcbrt(ap[1], gf_characteristic())
 313    else gcbrt(ap[1], ap[2]);
 314
 315
 316 /* Generalized Shank's algorithm for cube roots
 317 */
 318 gcbrt(a, p) := block(
 319    [c,c_1,c_2, s:0, q:p-1, b:2, bq,bq0, r, t,t0, j:0, i, omega],
 320    if p = 2 then return(mod(a, 2)),
 321    if mod(p, 3) = 2 then return(power_mod(a, (?ash(p, 1) - 1)/3, p)),
 322    if mod(p, 9) = 4 then (
 323       omega : mod(?ash(p + 1, -1) * (gsqrt(-3, p)[1] -1), p),
 324       c : power_mod(a, (?ash(p, 1) + 1)/9, p))
 325    else if mod(p,9) = 7 then (
 326       omega : mod(?ash(p + 1, -1) * (gsqrt(-3, p)[1] - 1), p),
 327       c : power_mod(a, (p + 2)/9, p) )
 328    else (
 329       while mod(q, 3) = 0 do (
 330          s : s + 1,
 331          q : q/3 ),
 332       omega : power_mod(b, (p - 1)/3, p),
 333       while omega = 1 do (
 334          b : next_prime(b),
 335          omega : power_mod(b, (p - 1)/3, p) ),
 336       /* b:next_prime(b) e' sensato, tanto e' inutile provare */
 337       /* i composti se quello che li ha preceduti ha fallito. */
 338       /* Magari implementiamo un crivello in piccolo ad hoc.  */
 339       bq : power_mod(b, q, p),
 340       if mod(q, 3) = 1 then (
 341          r : power_mod(a, (q - 1)/3, p),
 342          t : mod(a*r*r*r, p),
 343          t0 : t,
 344          while t # 1  do (
 345             t : mod(t^3, p),
 346             j : j + 1 ),
 347          for i:1 thru s - j do (
 348             bq0 : bq,
 349             bq : mod(bq0^3, p) ),
 350          if mod(t0 - bq, p) # 0 then bq0 : mod(bq0*bq0, p),
 351          c : mod(bq0 * inv_mod(r, p), p) )
 352       else (
 353          r : power_mod(a, (q + 1)/3, p),
 354          t : mod(inv_mod(a, p) * r*r*r, p),
 355          t0 : t,
 356          while t # 1 do (
 357             t : mod(t^3, p),
 358             j : j + 1 ),
 359          for i:1 thru s - j do (
 360             bq0 : bq,
 361             bq : mod(bq0^3, p) ),
 362          if mod(t0 - bq, p) # 0 then bq0 : mod(bq0*bq0, p),
 363          c : mod(inv_mod(bq0, p) * r, p) )),
 364    c : mod(c, p),
 365    c_1 : mod(c*omega, p),
 366    c_2 : mod(-c - c_1, p),
 367    [c, c_1, c_2] )$
 368
 369
 370 /* Fast exponentiation of the Frobenius (companion) matrix       */
 371 /*     (0  0  -a)   corresponding to the linear map              */
 372 /* F = (1  0  -b)   p(x) -> x*p(x) (mod x^3+bx+a)                */
 373 /*     (0  1   0)   It uses only 3 parameters (F_22,F_23,F_32)   */
 374 /* that are sufficient to completely describe F^j for any j;     */
 375 /* it returns (F^m)_22, corresponding to the evaluation of       */
 376 /*                  x^m (mod x^3+bx+a)                           */
 377 /* in x=0.                                                       */
 378 f3power(a, b, p, m) := block([v_22:1, v_23:0, v_32:0,
 379                               w_22:1, w_23:0, w_32:0, c, j, k, l, n, index],
 380    if m = 0 then return(1),
 381    b : mod(b, p),
 382    a : mod(a, p),
 383    c : inv_mod(-a, p),
 384    j : ?integer\-length(m),
 385    for i:1 thru j do (
 386       k : mod(v_23 + b*v_32, p),
 387       l : mod(?ash(v_22, 1), p),
 388       if ?logbitp(j-i, m)  then (
 389          n : mod(a*v_32, p),
 390          w_22 : mod(l*v_23 - n*v_32, p),
 391          w_32 : mod(v_32*(k + v_23) + v_22*v_22, p),
 392          w_23 : mod(k*k - l*n - b*w_32, p) )
 393       else (
 394          w_22 : mod(v_32*(k + v_23) + v_22*v_22, p),
 395          w_23 : mod(l*v_23 - a*v_32*v_32, p),
 396          w_32 : mod(l*v_32 + c*k*k, p) ),
 397       v_22 : w_22,
 398       v_23 : w_23,
 399       v_32 : w_32 ),
 400     v_22 )$
 401
 402
 403 /* This subroutine returns (1 0) A^(m) (1 0) mod p, where     */
 404 /*          ( 0  -a )                                         */
 405 /*     A =  ( 1  -b )                                         */
 406 /* i.e the m-th term of the recurrence                        */
 407 /* A_0 = 1, A_1 = 0, A_{n+2} = -b A_{n+1} - a A_{n}           */
 408 /* with characteristic polynomial x^2 + b x + a               */
 409 f2power(a, b, p, m) := block([v_11:0, v_21:1, w_11, w_21, u, v, index],
 410    /* if m<2 then return(1-m), */
 411    /* b:mod(b,p), a:mod(a,p),  */
 412    index : ?integer\-length(m),
 413    for i:2 thru index do (
 414       if  ?logbitp(index - i, m) then (
 415          u : - a*v_21,
 416          v : v_11 - b*v_21,
 417          w_11 : mod(u*(v + v_11), p),
 418          w_21 : mod(v*v + v_21*u, p) )
 419       else (
 420          u : v_11*v_11,
 421          v : v_21*v_21,
 422          w_11 : mod(u - a*v, p),
 423          w_21 : mod((v_11 + v_21)^2 - (u + (b + 1)*v), p) ),
 424       v_11 : w_11,
 425       v_21 : w_21 ),
 426    mod(v_11, p) )$
 427
 428
 429 /* Stickelberger-Redei irreducibility test for cubic polynomial f(x)=x^3+bx+a    */
 430 /* f(x) is irreducible over F_p iff its discriminant D=-4b^3-27a^2               */
 431 /* is a quadratic residue AND 2*(sqrt(D)+sqrt(-27a^2))^2 is NOT a cubic residue. */
 432 redei(a, b, p, ome, p3) := block([d, g, y, r],
 433    d : mod(- 4*b^3 - 27*a^2, p),
 434    if jacobi(d, p) # 1 then return(),
 435    g : mod(msqrt(d, p)[1] + 3*a*ome, p),
 436    y : mod(?ash(g*g, 1), p),
 437    r : power_mod(y, p3, p),
 438    is( mod(r, p) # 1 ))$
 439
 440
 441 /* For VERY large p, it might be useful to factor (p-1) into prime factors */
 442 /* (p-1) = 2^(a_1) 3^(a_2) ... q_l^(a_l)                                   */
 443 /* then use a factor-addition-chain instead of the binary one:             */
 444 /* r:y, for i:l downto 1 do (                                              */
 445 /*                   r : power_mod(r,q_l^(a_l),p),                         */
 446 /*                   if r=1 then return(false)                             */
 447 /*                        ), if r=1 then return(true) else return(false)   */
 448 /* to have great chances to shorten the cubic-residue test.                */
 449
 450
 451 /* If a is a quadratic residue in F_p, this routine returns    */
 452 /* two integers in [0, p-1] such that their squares            */
 453 /* are congruent to a (mod p).                                 */
 454 /*                                                             */
 455 /* If b^2-4a is a quadratic non-residue mod p, the polynomial  */
 456 /* x^2 - bx + a                                                */
 457 /* has two roots in F_(p^2) = (Z/pZ[x])/(x^2-bx+a)             */
 458 /* the former being the Frobenius conjugate of the latter      */
 459 /* y = y^p                                                     */
 460 /* so that sqrt(a)=y^((p+1)/2) for Viete's theorem.            */
 461 /* The point is to compute                                     */
 462 /* mu = (b/2 + sqrt(b^2-4a)/2)^((p+1)/2) = J + K sqrt(b^2-4a)  */
 463 /* in a fast way. We notes that if a is a quadratic residue    */
 464 /* mu can't be a proper element of the field extension         */
 465 /* F_p (sqrt(b^2-4a))                                          */
 466 /* so we must have K=0, or, in the same way,                   */
 467 /* mu = (1/2) { (b/2 + sqrt(delta)/2)^((p+1)/2)                */
 468 /*             +(b/2 - sqrt(delta)/2)^((p+1)/2 }               */
 469 /* that is the (p+1)/2-th term of the sequence                 */
 470 /* A_0 = 1,  A_1 = b/2,  A_{n+2}=b*A_{n+1}-a*A_{n}             */
 471 /* All we need is the (p+1)/2-power of the matrix              */
 472 /*       ( b  -a )                                             */
 473 /*   A = ( 1   0 )    over Z/pZ                                */
 474 /* The subroutine "mpower" returns it with about 8 log(p)      */
 475 /* multiplications and 4 log(p) mod-reductions in Z/pZ,        */
 476 /* using the standard "repeat-squaring" algorithm.             */
 477
 478 /* Field extension algorithm for square-root extraction through companion matrix powers */
 479 /* returns two integers in [0,p-1] whose square are congruent to a (mod p)              */
 480 msqrt(a ,p) := block([b:0, term],
 481    if mod(a, p) = 0 then return([0,0]),
 482    /* if not(primep(p)) then error("ERROR: Second argument must be a prime number."),   */
 483    if jacobi(a, p) # 1 then
 484       error("ERROR: First argument must be a quadratic residue."),
 485    if p = 2 then return(mod(a, 2)),
 486    do (
 487       b : b+1,
 488       if jacobi(b*b - 4*a, p) = -1 then return() ),
 489    /* print(b,"^2 - 4*",a," is a quadratic non-residue"), */
 490    term : mod(f2power(a, b, p, ?ash(p + 1, -1)), p),
 491    [term, p - term] )$
 492
 493
 494 /* Shank's algorithm for square-root extraction                                         */
 495 /* returns two integers in [0,p-1] whose square are congruent to a (mod p)              */
 496 ssqrt(a, p) := block([nrflag:false, b:0, s:0, q:p-1, t,alpha,c,r, d],
 497    if mod(a, p) = 0 then return([0,0]),
 498    /* if not(primep(p)) then error("ERROR: Second argument must be a prime number."),   */
 499    if jacobi(a, p) # 1 then
 500       error("ERROR: First argument must be a quadratic residue."),
 501    if p = 2 then return(mod(a, 2)),
 502    do (
 503       b : b+1,
 504       if jacobi(b, p) = -1 then return() ),
 505    unless ?logbitp(0, q) do (
 506       s : s+1,
 507       q : ?ash(q, -1) ),
 508    t : q,
 509    alpha : inv_mod(a, p),
 510    c : power_mod(b, t, p),
 511    r : power_mod(a, ?ash((t + 1), -1), p),
 512    for i:1 thru s-1 do (
 513       d : mod(r*r*alpha, p),
 514       for j:1 thru s-i-1 do d : mod(d*d, p),
 515       if mod(d + 1, p) = 0 then r : mod(r*c, p),
 516       c : mod(c*c, p) ),
 517    r : mod(r, p),
 518    [r, p - r] )$
 519
 520
 521 /* Field extension algorithm for cubic-root extraction        */
 522 /* For a cubic residue a and a prime p=1(mod 3) mcbrt finds   */
 523 /* an irreducible polynomial of the form x^3+bx+a over F_p    */
 524 /* so that -a=x^(p^2+p+1) and returns x^((p^2+p+1)/3)         */
 525 /* computed over F_(p^3) = (F_p)/(x^3+bx+a)                   */
 526 /* This subroutine returns three integers in [0,p-1]          */
 527 /* whose third-power are congruent to a (mod p)               */
 528 scbrt(a, p) := block( [x, b:1, h1,h2,h3, c,c_1,c_2, om,omega_1, p3],
 529    if mod(a, p) = 0 then return(0),
 530    /* if not primep(p) then error("ERROR: Second argument must be a prime number."), */
 531    if p = 2 then return(mod(a, 2)),
 532    if mod(p, 3) = 2 then return(power_mod(a, (?ash(p, 1) - 1)/3, p)),
 533    om : msqrt(-3, p),
 534    p3 : (p - 1)/3,
 535    if power_mod(a, p3, p) # 1 then
 536       error("ERROR: First argument must be a cubic residue."),
 537    unless redei(a, b, p, om[1], p3) do b : b+1,
 538
 539    gf_minimal_set(p, x^3 + b*x + a),
 540    h1 : gf_exp(x, p), /* get the benefit from gf_exp in numth.lisp (VvN) */
 541    h2 : gf_exp(h1, p3),
 542    h3 : gf_exp(x, (?ash(p, 1) + 1)/3),
 543    c : mod(- gf_mult(h2, h3), p),
 544    gf_unset(),
 545
 546    omega_1 : mod((-1 + om[1])*(?ash(p + 1, -1)), p),
 547    c_1 : mod(omega_1*c, p),
 548    c_2 : mod(- c - c_1, p),
 549    [c, c_1, c_2] )$
 550
 551
 552 /* Field extension algorithm for cubic-root extraction                     */
 553 /* through companion matrix powers. This subroutine returns three integers */
 554 /* in [0,p-1] whose third-power are congruent to a (mod p)                 */
 555 mcbrt(a, p) := block([b:1, c,c_1,c_2, om,omega_1, p3],
 556    if mod(a, p) = 0 then return(0),
 557    /* if not(primep(p)) then error("ERROR: Second argument must be a prime number."), */
 558    if p = 2 then return(mod(a, 2)),
 559    if mod(p, 3) = 2 then return(power_mod(a, (?ash(p, 1) - 1)/3, p)),
 560    om : msqrt(-3, p),
 561    p3 : (p - 1)/3,
 562    if power_mod(a, p3, p) # 1 then
 563       error("ERROR: First argument must be a cubic residue."),
 564    unless redei(a, b, p, om[1], p3) do b : b+1,
 565    c : mod(- f3power(a, b, p, (1 + p*(p + 1))/3), p),
 566    omega_1 : mod((-1 + om[1])*(?ash(p + 1, -1)), p),
 567    c_1 : mod(omega_1*c, p),
 568    c_2 : mod(- c - c_1, p),
 569    [c, c_1, c_2] )$
 570
 571 /* For each prime p in the form 4k+1 we find two positive integers  */
 572 /* a and b such that p=a^2+b^2. Useful to extend a Z-factorization  */
 573 /* into a Z[i]-factorization (order of Gaussian integers)           */
 574
 575 qsplit(p) := block([m, a, b:1, c,d,e,f, k,j],
 576    if p = 2 then return([1,1]),
 577    m : msqrt(-1, p),
 578    if m[1] < m[2] then a : m[1]
 579    else a : m[2],
 580    k : (a^2 + 1)/p,
 581    while k > 1 do (
 582       k : (a^2 + b^2)/p,
 583       c : mod(a, k),
 584       d : mod(b, k),
 585       j : (c^2 + d^2)/k,
 586       e : (a*c + b*d)/k,
 587       f : (b*c - a*d)/k,
 588       a : abs(e),
 589       b : abs(f),
 590       k : j ),
 591    [a, b] )$
 592
 593
 594 /* For each prime p in the form 3k+1 we find two positive integers     */
 595 /* a and b such that p=a^2+ab+b^2. Useful to extend a Z-factorization  */
 596 /* into a Z[omega]-factorization (order of Eisenstein integers)        */
 597 csplit(p):= block([m, a,b,c,d, t, x,y],
 598    if p = 3 then return([1,1]),
 599    m : msqrt(-3, p),
 600    if mod(m[1], 2) = 0 then b : m[2]
 601    else b : m[1],
 602    a : (b^2 + 3)/(4*p),
 603    c : p,
 604    x : 0,
 605    y : 1,
 606    while b # 1 do (
 607       if a > c then (
 608          d : a,
 609          a : c,
 610          c : d,
 611          b : - b,
 612          d : x,
 613          x : - y,
 614          y : d ),
 615       t : floor(b/?ash(a, 1)),
 616       b : b - ?ash(a*t, 1),
 617       c : c - t*(b + a*t),
 618       x : x + t*y ),
 619    if x < 0 then
 620       if y < - x then x : - (x + y)
 621       else (x : - x, y : y - x),
 622    if y < 0 then
 623       if x < - y then y : - (x + y)
 624       else (y : - y, x : x - y),
 625    [x, y] )$