Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - Fle.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.11.0 lcov report (development 22860-5579deb0b) Lines: 296 341 86.8 %
Date: 2018-07-19 05:36:41 Functions: 39 49 79.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /* Not so fast arithmetic with points over elliptic curves over Fl */
      18             : 
      19             : /***********************************************************************/
      20             : /**                                                                   **/
      21             : /**                              Flj                                  **/
      22             : /**                                                                   **/
      23             : /***********************************************************************/
      24             : 
      25             : /* Arithmetic is implemented using Jacobian coordinates, representing
      26             :  * a projective point (x : y : z) on E by [z*x , z^2*y , z].  This is
      27             :  * probably not the fastest representation available for the given
      28             :  * problem, but they're easy to implement and up to 60% faster than
      29             :  * the school-book method used in Fle_mulu(). */
      30             : 
      31             : /* Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
      32             :  * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl */
      33             : INLINE void
      34    20682142 : Flj_dbl_indir_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
      35             : {
      36             :   ulong X1, Y1, Z1;
      37             :   ulong XX, YY, YYYY, ZZ, S, M, T;
      38             : 
      39    20682142 :   X1 = P[1]; Y1 = P[2]; Z1 = P[3];
      40             : 
      41    20682142 :   if (Z1 == 0) { Q[1] = X1; Q[2] = Y1; Q[3] = Z1; return; }
      42             : 
      43    20542717 :   XX = Fl_sqr_pre(X1, p, pi);
      44    20542275 :   YY = Fl_sqr_pre(Y1, p, pi);
      45    20539402 :   YYYY = Fl_sqr_pre(YY, p, pi);
      46    20538189 :   ZZ = Fl_sqr_pre(Z1, p, pi);
      47    20537749 :   S = Fl_double(Fl_sub(Fl_sqr_pre(Fl_add(X1, YY, p), p, pi),
      48             :                        Fl_add(XX, YYYY, p), p), p);
      49    20538150 :   M = Fl_addmul_pre(Fl_triple(XX, p), a4, Fl_sqr_pre(ZZ, p, pi), p, pi);
      50    20541366 :   T = Fl_sub(Fl_sqr_pre(M, p, pi), Fl_double(S, p), p);
      51    20539895 :   Q[1] = T;
      52    20539895 :   Q[2] = Fl_sub(Fl_mul_pre(M, Fl_sub(S, T, p), p, pi),
      53             :                 Fl_double(Fl_double(Fl_double(YYYY, p), p), p), p);
      54    20538016 :   Q[3] = Fl_sub(Fl_sqr_pre(Fl_add(Y1, Z1, p), p, pi),
      55             :                 Fl_add(YY, ZZ, p), p);
      56             : }
      57             : 
      58             : INLINE void
      59    17539148 : Flj_dbl_pre_inplace(GEN P, ulong a4, ulong p, ulong pi)
      60             : {
      61    17539148 :   Flj_dbl_indir_pre(P, P, a4, p, pi);
      62    17537492 : }
      63             : 
      64             : GEN
      65     3139406 : Flj_dbl_pre(GEN P, ulong a4, ulong p, ulong pi)
      66             : {
      67     3139406 :   GEN Q = cgetg(4, t_VECSMALL);
      68     3139377 :   Flj_dbl_indir_pre(P, Q, a4, p, pi);
      69     3139314 :   return Q;
      70             : }
      71             : 
      72             : /* Cost: 11M + 5S + 9add + 4*2.
      73             :  * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl */
      74             : INLINE void
      75     6201563 : Flj_add_indir_pre(GEN P, GEN Q, GEN R, ulong a4, ulong p, ulong pi)
      76             : {
      77             :   ulong X1, Y1, Z1, X2, Y2, Z2;
      78             :   ulong Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W;
      79     6201563 :   X1 = P[1]; Y1 = P[2]; Z1 = P[3];
      80     6201563 :   X2 = Q[1]; Y2 = Q[2]; Z2 = Q[3];
      81             : 
      82     6201563 :   if (Z2 == 0) { R[1] = X1; R[2] = Y1; R[3] = Z1; return; }
      83     6201454 :   if (Z1 == 0) { R[1] = X2; R[2] = Y2; R[3] = Z2; return; }
      84             : 
      85     6191402 :   Z1Z1 = Fl_sqr_pre(Z1, p, pi);
      86     6191362 :   Z2Z2 = Fl_sqr_pre(Z2, p, pi);
      87     6191246 :   U1 = Fl_mul_pre(X1, Z2Z2, p, pi);
      88     6191310 :   U2 = Fl_mul_pre(X2, Z1Z1, p, pi);
      89     6191250 :   S1 = Fl_mul_pre(Y1, Fl_mul_pre(Z2, Z2Z2, p, pi), p, pi);
      90     6191329 :   S2 = Fl_mul_pre(Y2, Fl_mul_pre(Z1, Z1Z1, p, pi), p, pi);
      91     6191423 :   H = Fl_sub(U2, U1, p);
      92     6191380 :   r = Fl_double(Fl_sub(S2, S1, p), p);
      93             : 
      94     6191451 :   if (H == 0) {
      95      469571 :     if (r == 0) Flj_dbl_indir_pre(P, R, a4, p, pi); /*Points are equal, double*/
      96      467051 :     else { R[1] = R[2] = 1; R[3] = 0; } /* Points are opposite, return zero */
      97      469571 :     return;
      98             :   }
      99     5721880 :   I = Fl_sqr_pre(Fl_double(H, p), p, pi);
     100     5722070 :   J = Fl_mul_pre(H, I, p, pi);
     101     5721942 :   V = Fl_mul_pre(U1, I, p, pi);
     102     5721967 :   W = Fl_sub(Fl_sqr_pre(r, p, pi), Fl_add(J, Fl_double(V, p), p), p);
     103     5722051 :   R[1] = W;
     104     5722051 :   R[2] = Fl_sub(Fl_mul_pre(r, Fl_sub(V, W, p), p, pi),
     105             :                 Fl_double(Fl_mul_pre(S1, J, p, pi), p), p);
     106     5722088 :   R[3] = Fl_mul_pre(Fl_sub(Fl_sqr_pre(Fl_add(Z1, Z2, p), p, pi),
     107             :                            Fl_add(Z1Z1, Z2Z2, p), p), H, p, pi);
     108             : }
     109             : 
     110             : INLINE void
     111     6201312 : Flj_add_pre_inplace(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
     112     6201312 : { Flj_add_indir_pre(P, Q, P, a4, p, pi); }
     113             : 
     114             : GEN
     115           0 : Flj_add_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
     116             : {
     117           0 :   GEN R = cgetg(4, t_VECSMALL);
     118           0 :   Flj_add_indir_pre(P, Q, R, a4, p, pi);
     119           0 :   return R;
     120             : }
     121             : 
     122             : GEN
     123     2435759 : Flj_neg(GEN Q, ulong p)
     124     2435759 : { return mkvecsmall3(Q[1], Fl_neg(Q[2], p), Q[3]); }
     125             : 
     126             : typedef struct {
     127             :   ulong n; /* The number being represented */
     128             :   ulong pbits, nbits;  /* Positive bits and negative bits */
     129             :   ulong lnzb; /* Leading non-zero bit */
     130             : } naf_t;
     131             : 
     132             : /* Return the signed binary representation (i.e. the Non-Adjacent Form
     133             :  * in base 2) of 0 <= a < 2^63 */
     134             : static void
     135     3504833 : naf_repr(naf_t *x, ulong a)
     136             : {
     137             :   long t, i;
     138             :   ulong pbits, nbits;
     139     3504833 :   ulong c0 = 0, c1, a0;
     140             : 
     141     3504833 :   x->n = a;
     142     3504833 :   pbits = nbits = 0;
     143    32930319 :   for (i = 0; a; a >>= 1, ++i) {
     144    29425486 :     a0 = a & 1;
     145    29425486 :     c1 = (c0 + a0 + ((a & 2) >> 1)) >> 1;
     146    29425486 :     t = c0 + a0 - (c1 << 1);
     147    29425486 :     if (t < 0)
     148     4532261 :       nbits |= (1UL << i);
     149    24893225 :     else if (t > 0)
     150     5378221 :       pbits |= (1UL << i);
     151    29425486 :     c0 = c1;
     152             :   }
     153     3504833 :   c1 = c0 >> 1;
     154     3504833 :   t = c0 - (c1 << 1);
     155             :   /* Note that we don't need to check whether t < 0, since a >= 0 implies
     156             :    * that this most significant signed bit must be non-negative. */
     157     3504833 :   if (t > 0) pbits |= (1UL << i);
     158             : 
     159     3504833 :   x->pbits = pbits;
     160     3504833 :   x->nbits = nbits;
     161             :   /* Note that expu returns the least nonzero bit in the argument,
     162             :    * like the bit-scan-rev instruction on Intel architectures. */
     163             :   /* Using pbits here is justified by the fact that a >= 0, so the
     164             :    * most significant bit must be positive. */
     165     3504833 :   x->lnzb = expu(pbits) - 2;
     166     3504820 : }
     167             : 
     168             : /* Standard left-to-right signed double-and-add to compute [n]P. */
     169             : static GEN
     170     3386235 : Flj_mulu_pre_naf(GEN P, ulong n, ulong a4, ulong p, ulong pi, const naf_t *x)
     171             : {
     172             :   GEN R, Pinv;
     173             :   ulong pbits, nbits, lnzb;
     174             :   ulong m;
     175             : 
     176     3386235 :   if (n == 0) return mkvecsmall3(1, 1, 0);
     177     3386203 :   if (n == 1) return Flv_copy(P);
     178             : 
     179     3139415 :   R = Flj_dbl_pre(P, a4, p, pi);
     180     3139303 :   if (n == 2) return R;
     181             : 
     182     2435766 :   pbits = x->pbits;
     183     2435766 :   nbits = x->nbits;
     184     2435766 :   lnzb = x->lnzb;
     185             : 
     186     2435766 :   Pinv = Flj_neg(P, p);
     187     2435924 :   m = (1UL << lnzb);
     188    19975227 :   for ( ; m; m >>= 1) {
     189    17539404 :     Flj_dbl_pre_inplace(R, a4, p, pi);
     190    17539814 :     if (m & pbits)
     191     2725053 :       Flj_add_pre_inplace(R, P, a4, p, pi);
     192    14814761 :     else if (m & nbits)
     193     3476928 :       Flj_add_pre_inplace(R, Pinv, a4, p, pi);
     194             :   }
     195     2435823 :   avma = (pari_sp)R; return R;
     196             : }
     197             : 
     198             : GEN
     199     2323758 : Flj_mulu_pre(GEN P, ulong n, ulong a4, ulong p, ulong pi)
     200             : {
     201             :   naf_t x;
     202     2323758 :   naf_repr(&x, n);
     203     2323824 :   return Flj_mulu_pre_naf(P, n, a4, p, pi, &x);
     204             : }
     205             : 
     206             : ulong
     207      245053 : Flj_order_ufact(GEN P, ulong n, GEN F, ulong a4, ulong p, ulong pi)
     208             : {
     209      245053 :   pari_sp av = avma;
     210      245053 :   ulong res = 1;
     211             :   long i, nfactors;
     212             :   GEN primes, exps;
     213             : 
     214      245053 :   primes = gel(F, 1);
     215      245053 :   nfactors = lg(primes);
     216      245053 :   exps = gel(F, 2);
     217             : 
     218      699704 :   for (i = 1; i < nfactors; ++i) {
     219      560259 :     ulong q, pp = primes[i];
     220      560259 :     long j, k, ei = exps[i];
     221             :     naf_t x;
     222             :     GEN b;
     223             : 
     224      560259 :     for (q = pp, j = 1; j < ei; ++j) q *= pp;
     225      560259 :     b = Flj_mulu_pre(P, n / q, a4, p, pi);
     226             : 
     227      560259 :     naf_repr(&x, pp);
     228     1622670 :     for (j = 0; j < ei && b[3] != 0; ++j)
     229     1062411 :       b = Flj_mulu_pre_naf(b, pp, a4, p, pi, &x);
     230      560259 :     if (b[3] != 0) return 0;
     231      454651 :     for (k = 0; k < j; ++k) res *= pp;
     232      454651 :     avma = av;
     233             :   }
     234      139445 :   return res;
     235             : }
     236             : 
     237             : GEN
     238     1164732 : Fle_to_Flj(GEN P)
     239     2329464 : { return ell_is_inf(P) ? mkvecsmall3(1UL, 1UL, 0UL):
     240     1164732 :                          mkvecsmall3(P[1], P[2], 1UL);
     241             : }
     242             : 
     243             : GEN
     244     1164732 : Flj_to_Fle_pre(GEN P, ulong p, ulong pi)
     245             : {
     246     1164732 :   if (P[3] == 0) return ellinf();
     247             :   else
     248             :   {
     249     1055642 :     ulong Z = Fl_inv(P[3], p);
     250     1055642 :     ulong Z2 = Fl_sqr_pre(Z, p, pi);
     251     1055642 :     ulong X3 = Fl_mul_pre(P[1], Z2, p, pi);
     252     1055642 :     ulong Y3 = Fl_mul_pre(P[2], Fl_mul_pre(Z, Z2, p, pi), p, pi);
     253     1055642 :     return mkvecsmall2(X3, Y3);
     254             :   }
     255             : }
     256             : 
     257             : INLINE void
     258     6469868 : random_Fle_pre_indir(ulong a4, ulong a6, ulong p, ulong pi,
     259             :                      ulong *pt_x, ulong *pt_y)
     260             : {
     261             :   ulong x, x2, y, rhs;
     262             :   do
     263             :   {
     264     6469868 :     x   = random_Fl(p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     265     6469909 :     x2  = Fl_sqr_pre(x, p, pi);
     266     6469882 :     rhs = Fl_addmul_pre(a6, x, Fl_add(x2, a4, p), p, pi);
     267     6469888 :   } while ((!rhs && !Fl_add(Fl_triple(x2,p),a4,p)) || krouu(rhs, p) < 0);
     268     3241794 :   y = Fl_sqrt_pre(rhs, p, pi);
     269     3241789 :   *pt_x = x; *pt_y = y;
     270     3241789 : }
     271             : 
     272             : GEN
     273      415600 : random_Flj_pre(ulong a4, ulong a6, ulong p, ulong pi)
     274             : {
     275             :   ulong x, y;
     276      415600 :   random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
     277      415602 :   return mkvecsmall3(x, y, 1);
     278             : }
     279             : 
     280             : /***********************************************************************/
     281             : /**                                                                   **/
     282             : /**                              Fle                                  **/
     283             : /**                                                                   **/
     284             : /***********************************************************************/
     285             : GEN
     286           0 : Fle_changepoint(GEN x, GEN ch, ulong p)
     287             : {
     288             :   ulong p1,u,r,s,t,v,v2,v3;
     289             :   GEN z;
     290           0 :   if (ell_is_inf(x)) return x;
     291           0 :   u = ch[1]; r = ch[2];
     292           0 :   s = ch[3]; t = ch[4];
     293           0 :   v = Fl_inv(u, p); v2 = Fl_sqr(v,p); v3 = Fl_mul(v,v2,p);
     294           0 :   p1 = Fl_sub(x[1],r,p);
     295           0 :   z = cgetg(3,t_VECSMALL);
     296           0 :   z[1] = Fl_mul(v2, p1, p);
     297           0 :   z[2] = Fl_mul(v3, Fl_sub(x[2], Fl_add(Fl_mul(s,p1, p),t, p),p),p);
     298           0 :   return z;
     299             : }
     300             : 
     301             : GEN
     302         294 : Fle_changepointinv(GEN x, GEN ch, ulong p)
     303             : {
     304             :   ulong u, r, s, t, X, Y, u2, u3, u2X;
     305             :   GEN z;
     306         294 :   if (ell_is_inf(x)) return x;
     307         294 :   X = x[1]; Y = x[2];
     308         294 :   u = ch[1]; r = ch[2];
     309         294 :   s = ch[3]; t = ch[4];
     310         294 :   u2 = Fl_sqr(u, p); u3 = Fl_mul(u,u2,p);
     311         294 :   u2X = Fl_mul(u2,X, p);
     312         294 :   z = cgetg(3, t_VECSMALL);
     313         294 :   z[1] = Fl_add(u2X,r,p);
     314         294 :   z[2] = Fl_add(Fl_mul(u3,Y,p), Fl_add(Fl_mul(s,u2X,p), t, p), p);
     315         294 :   return z;
     316             : }
     317             : static GEN
     318      420774 : Fle_dbl_slope(GEN P, ulong a4, ulong p, ulong *slope)
     319             : {
     320             :   ulong x, y, Qx, Qy;
     321      420774 :   if (ell_is_inf(P) || !P[2]) return ellinf();
     322      368207 :   x = P[1]; y = P[2];
     323      368207 :   *slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
     324             :                   Fl_double(y, p), p);
     325      368211 :   Qx = Fl_sub(Fl_sqr(*slope, p), Fl_double(x, p), p);
     326      368209 :   Qy = Fl_sub(Fl_mul(*slope, Fl_sub(x, Qx, p), p), y, p);
     327      368207 :   return mkvecsmall2(Qx, Qy);
     328             : }
     329             : 
     330             : GEN
     331      343543 : Fle_dbl(GEN P, ulong a4, ulong p)
     332             : {
     333             :   ulong slope;
     334      343543 :   return Fle_dbl_slope(P,a4,p,&slope);
     335             : }
     336             : 
     337             : static GEN
     338      672567 : Fle_add_slope(GEN P, GEN Q, ulong a4, ulong p, ulong *slope)
     339             : {
     340             :   ulong Px, Py, Qx, Qy, Rx, Ry;
     341      672567 :   if (ell_is_inf(P)) return Q;
     342      672581 :   if (ell_is_inf(Q)) return P;
     343      672628 :   Px = P[1]; Py = P[2];
     344      672628 :   Qx = Q[1]; Qy = Q[2];
     345      672628 :   if (Px==Qx) return Py==Qy ? Fle_dbl_slope(P, a4, p, slope): ellinf();
     346      595325 :   *slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
     347      595389 :   Rx = Fl_sub(Fl_sub(Fl_sqr(*slope, p), Px, p), Qx, p);
     348      595343 :   Ry = Fl_sub(Fl_mul(*slope, Fl_sub(Px, Rx, p), p), Py, p);
     349      595341 :   return mkvecsmall2(Rx, Ry);
     350             : }
     351             : 
     352             : GEN
     353      672551 : Fle_add(GEN P, GEN Q, ulong a4, ulong p)
     354             : {
     355             :   ulong slope;
     356      672551 :   return Fle_add_slope(P,Q,a4,p,&slope);
     357             : }
     358             : 
     359             : static GEN
     360          98 : Fle_neg(GEN P, ulong p)
     361             : {
     362          98 :   if (ell_is_inf(P)) return P;
     363          98 :   return mkvecsmall2(P[1], Fl_neg(P[2], p));
     364             : }
     365             : 
     366             : GEN
     367           0 : Fle_sub(GEN P, GEN Q, ulong a4, ulong p)
     368             : {
     369           0 :   pari_sp av = avma;
     370             :   ulong slope;
     371           0 :   return gerepileupto(av, Fle_add_slope(P, Fle_neg(Q, p), a4, p, &slope));
     372             : }
     373             : 
     374             : struct _Fle { ulong a4, a6, p; };
     375             : 
     376             : static GEN
     377           0 : _Fle_dbl(void *E, GEN P)
     378             : {
     379           0 :   struct _Fle *ell = (struct _Fle *) E;
     380           0 :   return Fle_dbl(P, ell->a4, ell->p);
     381             : }
     382             : 
     383             : static GEN
     384        7777 : _Fle_add(void *E, GEN P, GEN Q)
     385             : {
     386        7777 :   struct _Fle *ell=(struct _Fle *) E;
     387        7777 :   return Fle_add(P, Q, ell->a4, ell->p);
     388             : }
     389             : 
     390             : GEN
     391     1509311 : Fle_mulu(GEN P, ulong n, ulong a4, ulong p)
     392             : {
     393             :   ulong pi;
     394     1509311 :   if (!n || ell_is_inf(P)) return ellinf();
     395     1509311 :   if (n==1) return zv_copy(P);
     396     1508114 :   if (n==2) return Fle_dbl(P, a4, p);
     397     1164732 :   pi = get_Fl_red(p);
     398     1164732 :   return Flj_to_Fle_pre(Flj_mulu_pre(Fle_to_Flj(P), n, a4, p, pi), p, pi);
     399             : }
     400             : 
     401             : static GEN
     402      872906 : _Fle_mul(void *E, GEN P, GEN n)
     403             : {
     404      872906 :   pari_sp av = avma;
     405      872906 :   struct _Fle *e=(struct _Fle *) E;
     406      872906 :   long s = signe(n);
     407             :   GEN Q;
     408      872906 :   if (!s || ell_is_inf(P)) return ellinf();
     409      872892 :   if (s < 0) P = Fle_neg(P, e->p);
     410      872892 :   if (is_pm1(n)) return s > 0? zv_copy(P): P;
     411      872402 :   Q = (lgefint(n)==3) ? Fle_mulu(P, uel(n,2), e->a4, e->p):
     412             :                         gen_pow(P, n, (void*)e, &_Fle_dbl, &_Fle_add);
     413      872402 :   return s > 0? Q: gerepileuptoleaf(av, Q);
     414             : }
     415             : 
     416             : GEN
     417         133 : Fle_mul(GEN P, GEN n, ulong a4, ulong p)
     418             : {
     419             :   struct _Fle E;
     420         133 :   E.a4 = a4; E.p = p;
     421         133 :   return _Fle_mul(&E, P, n);
     422             : }
     423             : 
     424             : /* Finds a random non-singular point on E */
     425             : GEN
     426     2826186 : random_Fle_pre(ulong a4, ulong a6, ulong p, ulong pi)
     427             : {
     428             :   ulong x, y;
     429     2826186 :   random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
     430     2826186 :   return mkvecsmall2(x, y);
     431             : }
     432             : 
     433             : GEN
     434           0 : random_Fle(ulong a4, ulong a6, ulong p)
     435           0 : { return random_Fle_pre(a4, a6, p, get_Fl_red(p)); }
     436             : 
     437             : static GEN
     438           0 : _Fle_rand(void *E)
     439             : {
     440           0 :   struct _Fle *e=(struct _Fle *) E;
     441           0 :   return random_Fle(e->a4, e->a6, e->p);
     442             : }
     443             : 
     444             : static const struct bb_group Fle_group={_Fle_add,_Fle_mul,_Fle_rand,hash_GEN,zv_equal,ell_is_inf,NULL};
     445             : 
     446             : GEN
     447      219106 : Fle_order(GEN z, GEN o, ulong a4, ulong p)
     448             : {
     449      219106 :   pari_sp av = avma;
     450             :   struct _Fle e;
     451      219106 :   e.a4=a4;
     452      219106 :   e.p=p;
     453      219106 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &Fle_group));
     454             : }
     455             : 
     456             : GEN
     457          42 : Fle_log(GEN a, GEN b, GEN o, ulong a4, ulong p)
     458             : {
     459          42 :   pari_sp av = avma;
     460             :   struct _Fle e;
     461          42 :   e.a4=a4;
     462          42 :   e.p=p;
     463          42 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &Fle_group));
     464             : }
     465             : 
     466             : ulong
     467           0 : Fl_ellj(ulong a4, ulong a6, ulong p)
     468             : {
     469           0 :   if (SMALL_ULONG(p))
     470             :   { /* a43 = 4 a4^3 */
     471           0 :     ulong a43 = Fl_double(Fl_double(Fl_mul(a4, Fl_sqr(a4, p), p), p), p);
     472             :     /* a62 = 27 a6^2 */
     473           0 :     ulong a62 = Fl_mul(Fl_sqr(a6, p), 27 % p, p);
     474           0 :     ulong z1 = Fl_mul(a43, 1728 % p, p);
     475           0 :     ulong z2 = Fl_add(a43, a62, p);
     476           0 :     return Fl_div(z1, z2, p);
     477             :   }
     478           0 :   return Fl_ellj_pre(a4, a6, p, get_Fl_red(p));
     479             : }
     480             : 
     481             : void
     482      130181 : Fl_ellj_to_a4a6(ulong j, ulong p, ulong *pt_a4, ulong *pt_a6)
     483             : {
     484      130181 :   ulong zagier = 1728 % p;
     485      130181 :   if (j == 0)           { *pt_a4 = 0; *pt_a6 =1; }
     486      130181 :   else if (j == zagier) { *pt_a4 = 1; *pt_a6 =0; }
     487             :   else
     488             :   {
     489      130181 :     ulong k = Fl_sub(zagier, j, p);
     490      130180 :     ulong kj = Fl_mul(k, j, p);
     491      130185 :     ulong k2j = Fl_mul(kj, k, p);
     492      130182 :     *pt_a4 = Fl_triple(kj, p);
     493      130180 :     *pt_a6 = Fl_double(k2j, p);
     494             :   }
     495      130183 : }
     496             : 
     497             : ulong
     498     2439464 : Fl_elldisc_pre(ulong a4, ulong a6, ulong p, ulong pi)
     499             : { /* D = -(4A^3 + 27B^2) */
     500             :   ulong t1, t2;
     501     2439464 :   t1 = Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi);
     502     2439464 :   t1 = Fl_double(Fl_double(t1, p), p);
     503     2439464 :   t2 = Fl_mul_pre(27 % p, Fl_sqr_pre(a6, p, pi), p, pi);
     504     2439464 :   return Fl_neg(Fl_add(t1, t2, p), p);
     505             : }
     506             : 
     507             : ulong
     508           0 : Fl_elldisc(ulong a4, ulong a6, ulong p)
     509             : {
     510           0 :   if (SMALL_ULONG(p))
     511             :   { /* D = -(4A^3 + 27B^2) */
     512             :     ulong t1, t2;
     513           0 :     t1 = Fl_mul(a4, Fl_sqr(a4, p), p);
     514           0 :     t1 = Fl_double(Fl_double(t1, p), p);
     515           0 :     t2 = Fl_mul(27 % p, Fl_sqr(a6, p), p);
     516           0 :     return Fl_neg(Fl_add(t1, t2, p), p);
     517             :   }
     518           0 :   return Fl_elldisc_pre(a4, a6, p, get_Fl_red(p));
     519             : }
     520             : 
     521             : static ulong
     522           0 : nonsquare_Fl(ulong p)
     523             : {
     524             :   ulong a;
     525           0 :   do a = random_Fl(p); while (krouu(a, p) >= 0);
     526           0 :   return a;
     527             : }
     528             : 
     529             : void
     530      210146 : Fl_elltwist_disc(ulong a4, ulong a6, ulong D, ulong p, ulong *pa4, ulong *pa6)
     531             : {
     532      210146 :   ulong D2 = Fl_sqr(D, p);
     533      210146 :   *pa4 = Fl_mul(a4, D2, p);
     534      210148 :   *pa6 = Fl_mul(a6, Fl_mul(D, D2, p), p);
     535      210143 : }
     536             : 
     537             : void
     538           0 : Fl_elltwist(ulong a4, ulong a6, ulong p, ulong *pt_a4, ulong *pt_a6)
     539           0 : { Fl_elltwist_disc(a4, a6, nonsquare_Fl(p), p, pt_a4, pt_a6); }
     540             : 
     541             : static void
     542    41560911 : Fle_dbl_sinv_pre_inplace(GEN P, ulong a4, ulong sinv, ulong p, ulong pi)
     543             : {
     544             :   ulong x, y, slope;
     545    41560911 :   if (uel(P,1)==p) return;
     546    41356941 :   if (!P[2]) { P[1] = p; return; }
     547    41234517 :   x = P[1]; y = P[2];
     548    41234517 :   slope = Fl_mul_pre(Fl_add(Fl_triple(Fl_sqr_pre(x, p, pi), p), a4, p),
     549             :                 sinv, p, pi);
     550    41234517 :   P[1] = Fl_sub(Fl_sqr_pre(slope, p, pi), Fl_double(x, p), p);
     551    41234517 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(x, P[1], p), p, pi), y, p);
     552             : }
     553             : 
     554             : static void
     555     5494593 : Fle_add_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
     556             : {
     557             :   ulong Px, Py, Qx, Qy, slope;
     558     5494593 :   if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Q[2]; }
     559     5494593 :   if (ell_is_inf(Q)) return;
     560     5494593 :   Px = P[1]; Py = P[2];
     561     5494593 :   Qx = Q[1]; Qy = Q[2];
     562     5494593 :   if (Px==Qx)
     563             :   {
     564       22678 :     if (Py==Qy) Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
     565       11112 :     else P[1] = p;
     566       22678 :     return;
     567             :   }
     568     5471915 :   slope = Fl_mul_pre(Fl_sub(Py, Qy, p), sinv, p, pi);
     569     5471915 :   P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
     570     5471915 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
     571             : }
     572             : 
     573             : static void
     574     6159246 : Fle_sub_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
     575             : {
     576             :   ulong Px, Py, Qx, Qy, slope;
     577     6159246 :   if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Fl_neg(Q[2], p); }
     578     6159246 :   if (ell_is_inf(Q)) return;
     579     6159246 :   Px = P[1]; Py = P[2];
     580     6159246 :   Qx = Q[1]; Qy = Q[2];
     581     6159246 :   if (Px==Qx)
     582             :   {
     583       27223 :     if (Py==Qy) P[1] = p;
     584             :     else
     585       10325 :       Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
     586       27223 :     return;
     587             :   }
     588     6132023 :   slope = Fl_mul_pre(Fl_add(Py, Qy, p), sinv, p, pi);
     589     6132023 :   P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
     590     6132023 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
     591             : }
     592             : 
     593             : static long
     594    52968113 : skipzero(long n) { return n ? n:1; }
     595             : 
     596             : void
     597      951165 : FleV_add_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
     598             : {
     599      951165 :   long i, l=lg(a4);
     600      951165 :   GEN sinv = cgetg(l, t_VECSMALL);
     601     6445758 :   for(i=1; i<l; i++)
     602     5494593 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
     603      951165 :   Flv_inv_pre_inplace(sinv, p, pi);
     604     6445758 :   for (i=1; i<l; i++)
     605     5494593 :     Fle_add_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
     606      951165 : }
     607             : 
     608             : void
     609     1087892 : FleV_sub_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
     610             : {
     611     1087892 :   long i, l=lg(a4);
     612     1087892 :   GEN sinv = cgetg(l, t_VECSMALL);
     613     7247138 :   for(i=1; i<l; i++)
     614     6159246 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
     615     1087892 :   Flv_inv_pre_inplace(sinv, p, pi);
     616     7247138 :   for (i=1; i<l; i++)
     617     6159246 :     Fle_sub_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
     618     1087892 : }
     619             : 
     620             : void
     621     7565264 : FleV_dbl_pre_inplace(GEN P, GEN a4, ulong p, ulong pi)
     622             : {
     623     7565264 :   long i, l=lg(a4);
     624     7565264 :   GEN sinv = cgetg(l, t_VECSMALL);
     625    49104284 :   for(i=1; i<l; i++)
     626    41539020 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_double(umael(P,i,2), p));
     627     7565264 :   Flv_inv_pre_inplace(sinv, p, pi);
     628    49104284 :   for(i=1; i<l; i++)
     629    41539020 :     Fle_dbl_sinv_pre_inplace(gel(P,i), uel(a4,i), uel(sinv,i), p, pi);
     630     7565264 : }
     631             : 
     632             : static void
     633      620733 : FleV_mulu_pre_naf_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi, const naf_t *x)
     634             : {
     635      620733 :   pari_sp av = avma;
     636             :   ulong pbits, nbits, lnzb, m;
     637             :   GEN R;
     638      620733 :   if (n == 1) return;
     639             : 
     640      620733 :   R = P; P = gcopy(P);
     641      620733 :   FleV_dbl_pre_inplace(R, a4, p, pi);
     642      620733 :   if (n == 2) return;
     643             : 
     644      620629 :   pbits = x->pbits;
     645      620629 :   nbits = x->nbits;
     646      620629 :   lnzb = x->lnzb;
     647             : 
     648      620629 :   m = (1UL << lnzb);
     649     7565160 :   for ( ; m; m >>= 1) {
     650     6944531 :     FleV_dbl_pre_inplace(R, a4, p, pi);
     651     6944531 :     if (m & pbits)
     652      951165 :       FleV_add_pre_inplace(R, P, a4, p, pi);
     653     5993366 :     else if (m & nbits)
     654     1087892 :       FleV_sub_pre_inplace(R, P, a4, p, pi);
     655             :   }
     656      620629 :   avma = av;
     657             : }
     658             : 
     659             : void
     660      620733 : FleV_mulu_pre_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi)
     661             : {
     662             :   naf_t x;
     663      620733 :   naf_repr(&x, n);
     664      620733 :   FleV_mulu_pre_naf_inplace(P, n, a4, p, pi, &x);
     665      620733 : }

Generated by: LCOV version 1.13