///////////////////////////////////////////////////// // Implementation in Magma of the algorithms in the // article: // Fast Jacobian group operations for C_{3,4} curves // over a large finite field // by Fatima K. Abu Salem and Kamal Khuri-Makdisi // preprint arxiv:math.NT/0610121 // to appear in LMS Journal of Computation and Mathematics // // General remark: the implementations here faithfully // follow the notations used in the above article. We // aimed primarily at having code that is easy for the // reader to follow side by side with the article. // This allows for easy verification of the formulas in // our article and for the number of multiplications // needed. We have not tried to optimize other aspects // of the code. For example, we frequently move data // around from one variable to another to bring out // how the computation unfolds in the different sections // of our article. This moving around of data can be // removed by renaming the variables, resulting in // greater runtime efficiency for the same number of // multiplications and inversions, but at the cost of a // loss of clarity in the source code. // // date: May 14, 2007 ///////////////////////////////////////////////////// // set up the finite field k p := NextPrime(10^8); k := GF(p); // Poor man's data structures: // 1) the pair of polynomials F,G representing a degree 3 divisor // is represented by a list [a,b,c,d,e,f,ainv]. newFG() returns a // blank vector whose contents we set to obtain a pair FG. geta := func; seta := procedure(~FG, val) FG[1]:=val; end procedure; getb := func; setb := procedure(~FG, val) FG[2]:=val; end procedure; getc := func; setc := procedure(~FG, val) FG[3]:=val; end procedure; getd := func; setd := procedure(~FG, val) FG[4]:=val; end procedure; gete := func; sete := procedure(~FG, val) FG[5]:=val; end procedure; getf := func; setf := procedure(~FG, val) FG[6]:=val; end procedure; getainv := func; setainv := procedure(~FG, val) FG[7]:=val; end procedure; newFG := func<| [k!0,k!0,k!0,k!0,k!0,k!0,k!0]>; // for debugging and testing. CAUTION: this only produces // a genuine divisor on the curve if the equation f of the // curve is contained in the ideal generated by F and // G in k[x,y]. What we did in practice was to find random // coefficients for F,G, F2, G2, and THEN to solve for // p0, p1, p2, q0, q1, q2 such that our polynomial f(x,y) // belonged to the intersection of with . randomFG := function() local answer; answer := newFG(); seta(~answer, Random(k)); // hope it's nonzero setb(~answer, Random(k)); setc(~answer, Random(k)); setd(~answer, Random(k)); sete(~answer, Random(k)); setf(~answer, Random(k)); setainv(~answer, 1/geta(answer)); return(answer); end function; // the above randomFG() function gave us these elements: FG := Vector(k, [31799142, 90609133, 49845484, 24754428, 79025059, 84333687, 80639977]); FG2 := Vector(k, [76509637, 62125676, 24876255, 68759361, 64391863, 63102834, 450027]); // and we fit our curve to pass through the points of D and D2. // This involved solving six linear equations for p0...q2 // yielding the answer cfs:=[64954369, 80199053, 62014001, 76871532, 68175610, 37961548]; q0 := k!cfs[1]; q1 := k!cfs[2]; q2 := k!cfs[3]; p0 := k!cfs[4]; p1 := k!cfs[5]; p2 := k!cfs[6]; /////////// // our original intention was to use the following to compare // the output of our algorithms with that of Magma's builtin // routines for function fields of algebraic curves: // // FF0 := FunctionField(k); // R := PolynomialRing(FF0); // f := y0^3 - x^4 + p2*x^2*y0 + p1*x*y0 + p0*y0 + q2*x^2 + q1*x + q0; // // FF1 := FunctionField(f); // // This was intended to allow us to test the results of // our algorithms to make sure that the implementation was // bug-free. // // // HOWEVER, we found it cleaner to do the testing and debugging // with ideals in the polynomial ring RR = k[x,y], and the Magma // builtins for computation with ideals using Groebner bases. // // To avoid name clashes with the originally intended approach // (and to also avoid name clashes between the variable f for which // G = xy + dy + ex + f on the one hand, and the completely different // from the equation f(x,y) of our curve on the other hand), we // called the equation of the curve ff(x,y) instead. /////////// RR := PolynomialRing(k,2); // we just copied the values of p0, p1, p2, q0, q1, q2 // from above into ff. Perhaps we should have written // this in terms of the variables p0...q2. // ff := RR!(y^3 - x^4 + 37961548*x^2*y + 68175610*x*y + 76871532*y + 62014001*x^2 + 80199053*x + 64954369); getF:=func; getG:=func; makeIdeal := func>; ////////////////////////////////////////////////////////// // 2) the vectors v1 = (alp, bet, gam, 1, 0) and // v2 = (del, eps, zet, 0, 1) // that give a basis for ker M' in sections 6 and 7 // are packaged into a single vector [alp, bet, gam, del, eps, zet] // and we provide a constructor newv1v2(). getalp := func; setalp := procedure(~v1v2, val) v1v2[1]:=val; end procedure; getbet := func; setbet := procedure(~v1v2, val) v1v2[2]:=val; end procedure; getgam := func; setgam := procedure(~v1v2, val) v1v2[3]:=val; end procedure; getdel := func; setdel := procedure(~v1v2, val) v1v2[4]:=val; end procedure; geteps := func; seteps := procedure(~v1v2, val) v1v2[5]:=val; end procedure; getzet := func; setzet := procedure(~v1v2, val) v1v2[6]:=val; end procedure; newv1v2 := func<| [k!0,k!0,k!0,k!0,k!0,k!0]>; // 3) We store the polynomials s,t corresponding to ker M': // s = x^3 + s1*y^2 + s2*xy + s3*x^2 + s4*y + s5*x + s6, // t = x^2y + t1*y^2 + t2*xy + t3*x^2 + t4*y + t5*x + t6 // as a vector [s1,s2,s3,s4,s5,s6,t1,t2,t3,t4,t5,t6], // and also include a constructor newst(). gets1 := func; sets1 := procedure(~st, val) st[1]:=val; end procedure; gets2 := func; sets2 := procedure(~st, val) st[2]:=val; end procedure; gets3 := func; sets3 := procedure(~st, val) st[3]:=val; end procedure; gets4 := func; sets4 := procedure(~st, val) st[4]:=val; end procedure; gets5 := func; sets5 := procedure(~st, val) st[5]:=val; end procedure; gets6 := func; sets6 := procedure(~st, val) st[6]:=val; end procedure; gett1 := func; sett1 := procedure(~st, val) st[7]:=val; end procedure; gett2 := func; sett2 := procedure(~st, val) st[8]:=val; end procedure; gett3 := func; sett3 := procedure(~st, val) st[9]:=val; end procedure; gett4 := func; sett4 := procedure(~st, val) st[10]:=val; end procedure; gett5 := func; sett5 := procedure(~st, val) st[11]:=val; end procedure; gett6 := func; sett6 := procedure(~st, val) st[12]:=val; end procedure; newst := func<| [k!0,k!0,k!0,k!0,k!0,k!0,k!0,k!0,k!0,k!0,k!0,k!0]>; ///////////////////////////////////////////////////////// // for testing and debugging: the following are similar // to getF() and getG() gets:=func; gett:=func; ///////////////////////////////////////////////////////// ///////////////////////////////////////////////////////// // Now we begin the actual algorithms from the article. ///////////////////////////////////////////////////////// // Functions to produce the matrices Tx and Ty for a given pair F,G. // These follow Proposition 3.2. // COST: 0M makeTx := func; // COST: 7M makeTy := function(FG) local a,ainv,b,c,d,e,f,g,h,i; a:=geta(FG); ainv:=getainv(FG); b:=getb(FG); c:=getc(FG); d:=getd(FG); e:=gete(FG); f:=getf(FG); g := ainv*(c + d*(d-b)) + e; h := ainv*(e*d - f); i := ainv*(e*c + f*(d-b)); return(Matrix(k,3,3, [k!0,-f,-i, k!0,-e,-h, k!1,-d,-g] ) ); end function; // COST: 6M // the "vector" below is actually a 3 x 1 matrix TxTimesVector := function(Tx,v) local v1, v2, v3; v1:=v[1,1]; v2:=v[2,1]; v3:=v[3,1]; return(Matrix(k,3,1, [ Tx[1,2]*v2 + Tx[1,3]*v3, v1 + Tx[2,2]*v2 + Tx[2,3]*v3, Tx[3,2]*v2 + Tx[3,3]*v3 ])); end function; // COST: 6M // the "vector" below is actually a 3 x 1 matrix TyTimesVector := function(Ty,v) local v1, v2, v3; v1:=v[1,1]; v2:=v[2,1]; v3:=v[3,1]; return(Matrix(k,3,1, [ Ty[1,2]*v2 + Ty[1,3]*v3, Ty[2,2]*v2 + Ty[2,3]*v3, v1 + Ty[3,2]*v2 + Ty[3,3]*v3 ])); end function; // COST: 11M // the "vector" below is actually a 3 x 1 matrix OneShotTyTimesVector := function(FG,v) local a,ainv,b,c,d,e,f,g,h,i, alp, bet, gam, del, eps, zet, eta, gamainv; a:=geta(FG); ainv:=getainv(FG); b:=getb(FG); c:=getc(FG); d:=getd(FG); e:=gete(FG); f:=getf(FG); alp:=v[1,1]; bet:=v[2,1]; gam:=v[3,1]; gamainv:=gam*ainv; del:=gamainv*f; eps:=alp-gamainv*c; zet:=gamainv*e; eta:=bet - gamainv*(b-d); return(Matrix(k,3,1, [ - zet*c - eta*f, del - zet*b - eta*e, eps - zet*a - eta*d ])); end function; // COST: 11M // the "two vectors" below are columns of a 3 x 2 matrix // Note that the matrix B is the BOTTOM 2 rows of v. // Similarly the matrix A consists of the RIGHTMOST 2 cols of Tx. // the reason for this is to minimize confusion in programming // Strassen's formulas. (N.B. "minimize" != "eliminate".) // Acknowledgement: Strassen's formulas and notation copied // off Wikipedia (I verified that they give the correct answer). TxTimesTwoVectors := function(Tx,v) local top1, top2, B11, B21, B12, B22, A11, A12, A21, A22, A31, A32, M1, M2, M3, M4, M5, M6, M7, C11, C12, C21, C22; top1:=v[1,1]; top2:=v[1,2]; B11:=v[2,1]; B12:=v[2,2]; B21:=v[3,1]; B22:=v[3,2]; A11:=Tx[1,2]; A12:=Tx[1,3]; A21:=Tx[2,2]; A22:=Tx[2,3]; A31:=Tx[3,2]; A32:=Tx[3,3]; M1 := (A11+A22)*(B11+B22); M2 := (A21+A22)*B11; M3 := A11*(B12-B22); M4 := A22*(B21-B11); M5 := (A11+A12)*B22; M6 := (A21-A11)*(B11+B12); M7 := (A12-A22)*(B21+B22); C11:= M1 + M4 - M5 + M7; C12:= M3 + M5; C21:= M2 + M4; C22:= M1 - M2 + M3 + M6; return(Matrix(k,3,2, [ C11, C12, top1 + C21, top2 + C22, A31*B11 + A32*B21, A31*B12 + A32*B22 ])); end function; ////////////////////////////////////////////////////////////// // Note: if we do not use Strassen, then we want to return // the following, for a total cost of 12M instead of 11M: // return(Matrix(k,3,2, // [ A11*B11 + A12*B21, A11*B12 + A12*B22, // top1 + A21*B11 + A22*B21, top2 + A21*B12 + A22*B22, // A31*B11 + A32*B21, A31*B12 + A32*B22 // ])); ////////////////////////////////////////////////////////////// // COST: 22M // This follows Proposition 4.1 EXCEPT that the roles of F,G and F',G' // are reversed. Thus the kernel of M gives the requisite combinations // of F,G that yield s,t. // Notation: BFBG is the 3x2 matrix whose columns are BF and BG (each // column represents the reduction of F or of G in R/). BxFBxG // similarly has columns BxF and BxG representing xF and xG. makeMprimeForAddition := function(FG, FGprime) local Tx, BFBG, BxFBxG, ByF, answer; BFBG := Matrix(k,3,2, [getc(FG) - getc(FGprime), getf(FG) - getf(FGprime), getb(FG) - getb(FGprime), gete(FG) - gete(FGprime), geta(FG) - geta(FGprime), getd(FG) - getd(FGprime) ]); Tx := makeTx(FGprime); BxFBxG := TxTimesTwoVectors(Tx, BFBG); // 11M so far ByF := OneShotTyTimesVector(FGprime, ColumnSubmatrix(BFBG,1,1)); // plus 11M // Now assemble into the answer, with columns // BF, BG, ByF-BxG, BxF, BxG in order: answer:= ZeroMatrix(k,3,5); InsertBlock(~answer, BFBG, 1,1); InsertBlock(~answer, ByF-ColumnSubmatrixRange(BxFBxG,2,2), 1,3); InsertBlock(~answer, BxFBxG, 1,4); return(answer); end function; ///////////////////////////////////////////////////////////// // We do inversion FIRST, before doubling, since there is // some repetition of formulas. Thus section 10 comes before // section 5. ///////////////////////////////////////////////////////////// // COST: 7M // Follows part 1 of Proposition 10.1 makeInverse := function(FG) local a, ainv, b, c, d, e, f, l, m, lainv, ab, d2, e2, f2, answer; a:=geta(FG); ainv:=getainv(FG); b:=getb(FG); c:=getc(FG); d:=getd(FG); e:=gete(FG); f:=getf(FG); m := e + a*(a + p2); l := c + (d - b)*d; lainv := l * ainv; ab := a * b; d2 := b - d; e2 := -(lainv + m); f2 := m*d + (lainv + e)*(d - b) + a*(ab - p1) - f; answer := newFG(); // has the same a, b, c, ainv as FG seta(~answer, a); setb(~answer, b); setc(~answer, c); setainv(~answer, ainv); setd(~answer, d2); sete(~answer, e2); setf(~answer, f2); return(answer); end function; // COST: 34M // This follows Proposition 5.4 including Proposition 10.1 part 2 makeMprimeForDoubling := function(FG) local a, ainv, b, c, d, e, f, l, m, lainv, ab, g, h, i, d2, e2, f2, u, Tx, Ty, BFBG, BxFBxG, ByF, answer; // CAUTION: the matrices BFBG, etc... // represent the differential forms dF,dG, etc.. // via polynomials G1 and -H1. See comment below. // step 1: find G1 and H1, costing 10M // here G1 = x*y + d2*y + e2*x + f2 // and H1 = -y^2 + a*x^2 + lainv*y - ab*x + u a:=geta(FG); ainv:=getainv(FG); b:=getb(FG); c:=getc(FG); d:=getd(FG); e:=gete(FG); f:=getf(FG); m := e + a*(a + p2); l := c + (d - b)*d; lainv := l * ainv; ab := a * b; d2 := b - d; e2 := -(lainv + m); f2 := m*d + (lainv + e)*(d - b) + a*(ab - p1) - f; u := (lainv + m)*e + a*(b*b -c - q2); // step 2: form Tx and find g,h,i to produce Ty. // Finding Ty costs only 5M instead of 7M, due // to partial results from step 1 Tx := makeTx(FG); g := lainv + e; h := ainv*(e*d - f); i := ainv*(e*c + f*(d-b)); Ty := Matrix(k,3,3, [k!0,-f,-i, k!0,-e,-h, k!1,-d,-g] ); // step 3: Form the matrix whose columns are BG1, BminusH1. // The polynomials G1 and -H1 play the roles that were // played by F and G when we did addition. Hence I will // ***still*** call the resulting matrix BFBG so as to reuse // the rest of the code from makeMprimeForAddition() below. // The punctilious reader of this code can think of BFBG as // really being BdFBdG, representing the differential forms // dF and dG via G1 and -H1, as in Lemma 5.3. // cost of this step: 2M. BFBG := Matrix(k,3,2, [f2 - f, -u - i + a*c, e2 - e, ab - h + ab, d2 - d, -lainv - g + a*a ]); // Step 4: continue just as in the algorithm for addition. // this part costs 11M+6M as indicated below. BxFBxG := TxTimesTwoVectors(Tx, BFBG); // costs 11M ByF := TyTimesVector(Ty, ColumnSubmatrix(BFBG,1,1)); // costs 6M // Now assemble into the answer, with columns // BF, BG, ByF-BxG, BxF, BxG in order: answer:= ZeroMatrix(k,3,5); InsertBlock(~answer, BFBG, 1,1); InsertBlock(~answer, ByF-ColumnSubmatrixRange(BxFBxG,2,2), 1,3); InsertBlock(~answer, BxFBxG, 1,4); return(answer); end function; // COST: 39M, 1I // find the kernel of M' as in Proposition 6.1. // returns the sequence of values // alpha, beta, gamma, delta, epsilon, zeta // that describe the basis elements v1' and v2' of the kernel. // This is returned as a data structure of type "v1v2". kernelOfMprime := function(Mprime) local A1, B1, C1, D1, E1, A2, B2, C2, D2, E2, A3, B3, C3, D3, E3, D, Del12, Del13, Del23, sig1, sig2, sig3, sig4, sig5, U, Q1, Q2, Q3, Q4, A1inv, Dinv, Uinv, alp, bet, gam, del, eps, zet, answer; // first set up variables as in the notation of // Proposition 6.1 (no cost involved) A1:=Mprime[1,1]; A2:=Mprime[2,1];A3:=Mprime[3,1]; B1:=Mprime[1,2]; B2:=Mprime[2,2];B3:=Mprime[3,2]; C1:=Mprime[1,3]; C2:=Mprime[2,3];C3:=Mprime[3,3]; D1:=Mprime[1,4]; D2:=Mprime[2,4];D3:=Mprime[3,4]; E1:=Mprime[1,5]; E2:=Mprime[2,5];E3:=Mprime[3,5]; // then find the echelon form (costs 21M) Del12 := A1*B2 - A2*B1; D := Del12; Del13 := A1*B3 - A3*B1; Del23 := A2*B3 - A3*B2; sig1 := A1*C2 - A2*C1; sig2 := A1*D2 - A2*D1; sig3 := A1*E2 - A2*E1; U := Del12*C3 - Del13*C2 + Del23*C1; sig4 := Del12*D3 - Del13*D2 + Del23*D1; sig5 := Del12*E3 - Del13*E2 + Del23*E1; // get three reciprocals using one inversion (6M,1I) Q1 := A1*D; Q2 := Q1*U; Q3 := 1/Q2; Uinv := Q1*Q3; Q4 := U*Q3; Dinv := A1*Q4; A1inv := D*Q4; // back substitution to find the basis elements (12M) gam := -Uinv*sig4; bet := -Dinv*(gam*sig1 + sig2); alp := -A1inv*(bet*B1 + gam*C1 + D1); zet := -Uinv*sig5; eps := -Dinv*(zet*sig1 + sig3); del := -A1inv*(eps*B1 + zet*C1 + E1); answer := newv1v2(); setalp(~answer, alp); setbet(~answer, bet); setgam(~answer, gam); setdel(~answer, del); seteps(~answer, eps); setzet(~answer, zet); return(answer); end function; // COST: 18M // Given the kernel of Mprime, find the polynomials s and t, // following Proposition 7.1. // HOWEVER, we do this by taking a combination of the // polynomials F and G and NOT of the polynomials F' and G'. // This is because we set up addition of divisors by reversing // the roles of {F,G} and {F',G'}. See the remarks above // the function makeMprimeForAddition(). makeCoefficientsOfst := function(FG,v1v2) local a, b, c, d, e, f, // we don't need ainv alp, bet, gam, del, eps, zet, gama, alpc, game, betf, zeta, delc, zete, epsf, s1, s2, s3, s4, s5, s6, t1, t2, t3, t4, t5, t6, answer; a := geta(FG); b := getb(FG); c := getc(FG); d := getd(FG); e := gete(FG); f := getf(FG); alp := getalp(v1v2); bet := getbet(v1v2); gam := getgam(v1v2); del := getdel(v1v2); eps := geteps(v1v2); zet := getzet(v1v2); // First half of computation: find // s = xF + (alp + gam*y)x^2 + (bet - gam*x)xy ...(I) // + gam*(b - d)xy + alp*b*x + bet*d*y ...(III)[sic] // + (alp+gam*y)(a*y+c) + (bet-gam*x)(ex+f)...(IV) // as in the proof of Proposition 7.1. // first initialize s = the term (I). Note the 'gam's cancel. s1 := 0; s2 := a + bet; s3 := b + alp; s4 := 0; s5 := c; s6 := 0; // now add (III) to the result, costing 3M: s2 +:= gam*(b - d); s4 +:= bet*d; s5 +:= alp*b; // add the first term of (IV), by Karatsuba, for 3M: gama := gam * a; alpc := alp * c; s1 +:= gama; s6 +:= alpc; s4 +:= (alp+gam)*(a+c) - gama - alpc; // similarly add the second term of (IV), for 3M: game := gam*e; betf := bet*f; s3 -:= game; s6 +:= betf; s5 +:= (bet-gam)*(e+f) + game - betf; // Second half of computation: find // t = xG + (alp + gam*y)x^2 + (bet - gam*x)xy ...(I) // + gam*(b - d)xy + alp*b*x + bet*d*y ...(III) // + (alp+gam*y)(a*y+c) + (bet-gam*x)(ex+f)...(IV) // note that xG has replaced xF in (I), and we have also // replaced alp -> del, bet -> eps, gam -> zet. // first initialize t = the term (I). Note the 'zet's cancel. t1 := 0; t2 := d + eps; t3 := e + del; t4 := 0; t5 := f; t6 := 0; // now add (III) to the result, costing 3M: t2 +:= zet*(b - d); t4 +:= eps*d; t5 +:= del*b; // add the first term of (IV), by Karatsuba, for 3M: zeta := zet * a; delc := del * c; t1 +:= zeta; t6 +:= delc; t4 +:= (del+zet)*(a+c) - zeta - delc; // similarly add the second term of (IV), for 3M: zete := zet*e; epsf := eps*f; t3 -:= zete; t6 +:= epsf; t5 +:= (eps-zet)*(e+f) + zete - epsf; // assemble together answer := newst(); sets1(~answer,s1); sets2(~answer,s2); sets3(~answer,s3); sets4(~answer,s4); sets5(~answer,s5); sets6(~answer,s6); sett1(~answer,t1); sett2(~answer,t2); sett3(~answer,t3); sett4(~answer,t4); sett5(~answer,t5); sett6(~answer,t6); return(answer); end function; // COST: 20M // Find the matrix M'', whose columns correspond to the // image of a basis for tW^7 inside W^17/(W^9 + sW^8). // This follows Lemma 9.2, which combines the formulas // of Lemmas 8.1 and 9.1 and arranges to save one // multiplication in the process. The bulk of the // formulas actually come from Lemma 9.1. // The result is M'' = [alp1, ..., alp5 // bet1, ..., bet5 // gam1, ..., gam5]. makeMdoublePrime := function(st) local s1, s2, s3, s4, s5, s6, t1, t2, t3, t4, t5, t6, s3p2, s1s2plusp2, s1s1, s1s3, s1s5plusq2, s1s4plusp1, m1, m2, m3, m4, m5, z1, z2, z3, z4, z5, z6, ell1, ell2, ell3, ell4, ell5, A11, A12, A21, A22, B11, B12, B21, B22, M1, M2, M3, M4, M5, M6, M7, C11, C12, C21, C22, alp1, alp2, alp3, alp4, alp5, bet1, bet2, bet3, bet4, bet5, gam1, gam2, gam3, gam4, gam5; // startup s1 := gets1(st); s2 := gets2(st); s3 := gets3(st); s4 := gets4(st); s5 := gets5(st); s6 := gets6(st); t1 := gett1(st); t2 := gett2(st); t3 := gett3(st); t4 := gett4(st); t5 := gett5(st); t6 := gett6(st); // some preliminary multiplications, namely those // for Lemmas 8.1 and 9.1 but without including // t3p2 and s1(s4+p1+s3p2), which are replaced // later by a single multiplication to find m1 // as in Lemma 9.2. // Cost: 6M for this first part. s3p2 := s3*p2; s1s2plusp2 := s1*(s2 + p2); s1s1 := s1*s1; s1s3 := s1*s3; s1s5plusq2 := s1*(s5 + q2); s1s4plusp1 := s1*(s4 + p1); // first column of M'', from the reduction // C1bar of C1. Free, of course. alp1 := k!1; bet1 := k!0; gam1 := k!0; // second column of M'', from C2bar, also free. alp2 := t2 - s3 + s1s2plusp2; bet2 := t1 - s2 + s1s1; gam2 := 0; // third column of M'', from C3bar. // Cost: 2M alp3 := t3 - t1*(s2 + p2); bet3 := t2 - t1*s1; gam3 := k!1; // preparation for column 4, as in Lemma 9.1. // nonzero entries m1...m5 of the column vector C4 - D10. // Note the special formula for m1, from Lemma 9.2. // Cost: 1M m1 := t4 - s5 + s1s4plusp1 + (t3+s1s3)*p2; m2 := - s4; // from p1 - (s4 + p1) m3 := t3 + s1s3; m4 := t2 - s3 + s1s2plusp2; m5 := t1 - s2 + s1s1; // p2 cancels similarly to m2 above // first half of preparation for column 5: // nonzero entries z1...z6 of the column vector C5 - D11. // Free. z1 := t5 + s1s5plusq2; z2 := t4 - s5 + s1s4plusp1; z3 := - s4; z4 := t3 + s1s3; z5 := t2 - s3 + s1s2plusp2; z6 := t1 - s2 + s1s1; // second half of preparation for column 5: // nonzero entries ell1...ell5 of C5 - D11 - z6*C9. // Cost: 4M ell1 := z1 - z6*(s4 + p1 + s3p2); ell2 := z2; ell3 := z3 - z6*s3; ell4 := z4 - z6*(s2 + p2); ell5 := z5 - z6*s1; // perform the 2x2 matrix multiplication of equation // (17) in Lemma 9.1, by reusing the code for Strassen // multiplication from TxTimesTwoVectors. // Cost: 7M A11 := s2 + p2; A12 := s3 - s1s2plusp2; A21 := s1; A22 := s2 - s1s1; B11 := m3; B12 := ell3; B21 := m4; B22 := ell4; M1 := (A11+A22)*(B11+B22); M2 := (A21+A22)*B11; M3 := A11*(B12-B22); M4 := A22*(B21-B11); M5 := (A11+A12)*B22; M6 := (A21-A11)*(B11+B12); M7 := (A12-A22)*(B21+B22); C11:= M1 + M4 - M5 + M7; C12:= M3 + M5; C21:= M2 + M4; C22:= M1 - M2 + M3 + M6; // finally use the above matrix product to get // columns 4 and 5 of M''. // Free. alp4 := m1 - C11; bet4 := m2 - C21; gam4 := m5; // alp5 := ell1 - C12; bet5 := ell2 - C22; gam5 := ell5; return(Matrix(k,3,5, [alp1, alp2, alp3, alp4, alp5, bet1, bet2, bet3, bet4, bet5, gam1, gam2, gam3, gam4, gam5] )); end function; // COST: 11M, 1I // This follows the second half of Proposition 9.3. // (The first half is the previous routine // makeMdoublePrime() above.) // // The routine below could have also been called // kernelOfMdoublePrime(). It calculates the // values a'', b'', c'', d'', e'', f'' (as well // as the inverse of a'') so that the kernel // of M'' is spanned by the transposes of // (c'',b'',a'',1,0) and (f'',e'',d'',0,1). // Then the desired divisor D'' such that // D + D' + D'' is principal is defined by // F'' = x^2 + a''*y + b''*x + c'', // G'' = x*y + d''*y + e''*x + f'' // However, for typographical convenience, // we shall just write a,b,c,d,e,f and drop // the double primes in the routine below. // Since M'' is in echelon form with the pivots // at M''[1,1] and M''[3,3] equal to 1, we // only need to do a simple back substitution // after we invert M''[2,2]. As explained // in the proof of Proposition 9.3, we combine // with another inversion to find ainv simultaneously. makeFGdoublePrime := function(MdoublePrime) local alp1, alp2, alp3, alp4, alp5, bet1, bet2, bet3, bet4, bet5, gam1, gam2, gam3, gam4, gam5, bet2gam4inv, gam4inv, bet2inv, a, b, c, d, e, f, ainv, answer; // initialize (although some entries // are zero or one, so superfluous; // see the proof of Proposition 9.3) alp1 := MdoublePrime[1,1]; // this is actually 1 bet1 := MdoublePrime[2,1]; // this is actually 0 gam1 := MdoublePrime[3,1]; // this is actually 0 alp2 := MdoublePrime[1,2]; bet2 := MdoublePrime[2,2]; gam2 := MdoublePrime[3,2]; // this is actually 0 alp3 := MdoublePrime[1,3]; bet3 := MdoublePrime[2,3]; gam3 := MdoublePrime[3,3]; // this is actually 1 alp4 := MdoublePrime[1,4]; bet4 := MdoublePrime[2,4]; gam4 := MdoublePrime[3,4]; alp5 := MdoublePrime[1,5]; bet5 := MdoublePrime[2,5]; gam5 := MdoublePrime[3,5]; // invert both bet2 and gam4. // cost: 3M, 1I bet2gam4inv := 1 / (bet2*gam4); gam4inv := bet2gam4inv * bet2; bet2inv := bet2gam4inv * gam4; // back substitution, and record 1/a. // cost: 8M a := -gam4; ainv := -gam4inv; b := -bet2inv*(bet3*a + bet4); c := -(alp2*b + alp3*a + alp4); d := -gam5; e := -bet2inv*(bet3*d + bet5); f := -(alp2*e + alp3*d + alp5); // put it all together answer := newFG(); seta(~answer,a); setainv(~answer,ainv); setb(~answer,b); setc(~answer,c); setd(~answer,d); sete(~answer,e); setf(~answer,f); return(answer); end function; // COST: 117M, 2I // Adds two distinct typical elements of the Jacobian. // Note that if we do not invert at the end, we can // save 7M from the above. This is reasonable if we // mainly find multiples of elements in the Jacobian, // using the "double-and-add" method. See the article // by Avanzi, Frey, Lange, and Oyono, "On using expansions // to the base of -2". // Caution: the roles of FG and FGprime have been reversed // throughout from the treatment of the article. This is // described in the comments above. makeAddition := function(FG, FGprime) local v1v2, MdoublePrime; // (22 + 39)M, 1I v1v2 := kernelOfMprime(makeMprimeForAddition(FG,FGprime)); // (18 + 20)M MdoublePrime := makeMdoublePrime(makeCoefficientsOfst(FG,v1v2)); // (11 + 7)M, 1I return(makeInverse(makeFGdoublePrime(MdoublePrime))); end function; // COST: 129M, 2I // The only change from the above is that we use // makeMprimeForDoubling() instead of makeMprimeForAddition() makeDoubling := function(FG) local v1v2, MdoublePrime; // (34 + 39)M, 1I v1v2 := kernelOfMprime(makeMprimeForDoubling(FG)); // (18 + 20)M MdoublePrime := makeMdoublePrime(makeCoefficientsOfst(FG,v1v2)); // (11 + 7)M, 1I return(makeInverse(makeFGdoublePrime(MdoublePrime))); end function;