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 to exceed 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.16.2 lcov report (development 29115-f22e516b23) Lines: 427 456 93.6 %
Date: 2024-03-19 08:06:49 Functions: 57 64 89.1 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : /* Not so fast arithmetic with points over elliptic curves over Fl */
      19             : 
      20             : /***********************************************************************/
      21             : /**                              Flj                                  **/
      22             : /***********************************************************************/
      23             : /* Jacobian coordinates: we represent a projective point (x:y:z) on E by
      24             :  * [z*x, z^2*y, z]. Not the fastest representation available for the problem,
      25             :  * but easy to implement and up to 60% faster than the school-book method. */
      26             : 
      27             : /* Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
      28             :  * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl */
      29             : INLINE void
      30    25129605 : Flj_dbl_indir_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
      31             : {
      32             :   ulong X1, Y1, Z1;
      33             :   ulong XX, YY, YYYY, ZZ, S, M, T;
      34             : 
      35    25129605 :   X1 = P[1]; Y1 = P[2]; Z1 = P[3];
      36             : 
      37    25129605 :   if (Z1 == 0) { Q[1] = X1; Q[2] = Y1; Q[3] = Z1; return; }
      38             : 
      39    24947649 :   XX = Fl_sqr_pre(X1, p, pi);
      40    24949765 :   YY = Fl_sqr_pre(Y1, p, pi);
      41    24933102 :   YYYY = Fl_sqr_pre(YY, p, pi);
      42    24929240 :   ZZ = Fl_sqr_pre(Z1, p, pi);
      43    24928681 :   S = Fl_double(Fl_sub(Fl_sqr_pre(Fl_add(X1, YY, p), p, pi),
      44             :                        Fl_add(XX, YYYY, p), p), p);
      45    24938016 :   M = Fl_addmul_pre(Fl_triple(XX, p), a4, Fl_sqr_pre(ZZ, p, pi), p, pi);
      46    24959340 :   T = Fl_sub(Fl_sqr_pre(M, p, pi), Fl_double(S, p), p);
      47    24945893 :   Q[1] = T;
      48    24945893 :   Q[2] = Fl_sub(Fl_mul_pre(M, Fl_sub(S, T, p), p, pi),
      49             :                 Fl_double(Fl_double(Fl_double(YYYY, p), p), p), p);
      50    24939111 :   Q[3] = Fl_sub(Fl_sqr_pre(Fl_add(Y1, Z1, p), p, pi),
      51             :                 Fl_add(YY, ZZ, p), p);
      52             : }
      53             : 
      54             : INLINE void
      55    21215033 : Flj_dbl_pre_inplace(GEN P, ulong a4, ulong p, ulong pi)
      56             : {
      57    21215033 :   Flj_dbl_indir_pre(P, P, a4, p, pi);
      58    21214421 : }
      59             : 
      60             : GEN
      61     3919444 : Flj_dbl_pre(GEN P, ulong a4, ulong p, ulong pi)
      62             : {
      63     3919444 :   GEN Q = cgetg(4, t_VECSMALL);
      64     3919241 :   Flj_dbl_indir_pre(P, Q, a4, p, pi);
      65     3918636 :   return Q;
      66             : }
      67             : 
      68             : /* Cost: 11M + 5S + 9add + 4*2.
      69             :  * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl */
      70             : INLINE void
      71     7606464 : Flj_add_indir_pre(GEN P, GEN Q, GEN R, ulong a4, ulong p, ulong pi)
      72             : {
      73             :   ulong X1, Y1, Z1, X2, Y2, Z2;
      74             :   ulong Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W;
      75     7606464 :   X1 = P[1]; Y1 = P[2]; Z1 = P[3];
      76     7606464 :   X2 = Q[1]; Y2 = Q[2]; Z2 = Q[3];
      77             : 
      78     7606464 :   if (Z2 == 0) { R[1] = X1; R[2] = Y1; R[3] = Z1; return; }
      79     7606153 :   if (Z1 == 0) { R[1] = X2; R[2] = Y2; R[3] = Z2; return; }
      80             : 
      81     7589815 :   Z1Z1 = Fl_sqr_pre(Z1, p, pi);
      82     7590198 :   Z2Z2 = Fl_sqr_pre(Z2, p, pi);
      83     7588644 :   U1 = Fl_mul_pre(X1, Z2Z2, p, pi);
      84     7588745 :   U2 = Fl_mul_pre(X2, Z1Z1, p, pi);
      85     7588763 :   S1 = Fl_mul_pre(Y1, Fl_mul_pre(Z2, Z2Z2, p, pi), p, pi);
      86     7589214 :   S2 = Fl_mul_pre(Y2, Fl_mul_pre(Z1, Z1Z1, p, pi), p, pi);
      87     7589047 :   H = Fl_sub(U2, U1, p);
      88     7589460 :   r = Fl_double(Fl_sub(S2, S1, p), p);
      89             : 
      90     7590074 :   if (H == 0) {
      91      598385 :     if (r == 0) Flj_dbl_indir_pre(P, R, a4, p, pi); /* P = Q */
      92      593518 :     else { R[1] = R[2] = 1; R[3] = 0; } /* P = -Q */
      93      598385 :     return;
      94             :   }
      95     6991689 :   I = Fl_sqr_pre(Fl_double(H, p), p, pi);
      96     6992026 :   J = Fl_mul_pre(H, I, p, pi);
      97     6991375 :   V = Fl_mul_pre(U1, I, p, pi);
      98     6991395 :   W = Fl_sub(Fl_sqr_pre(r, p, pi), Fl_add(J, Fl_double(V, p), p), p);
      99     6992118 :   R[1] = W;
     100     6992118 :   R[2] = Fl_sub(Fl_mul_pre(r, Fl_sub(V, W, p), p, pi),
     101             :                 Fl_double(Fl_mul_pre(S1, J, p, pi), p), p);
     102     6991592 :   R[3] = Fl_mul_pre(Fl_sub(Fl_sqr_pre(Fl_add(Z1, Z2, p), p, pi),
     103             :                            Fl_add(Z1Z1, Z2Z2, p), p), H, p, pi);
     104             : }
     105             : 
     106             : INLINE void
     107     7514435 : Flj_add_pre_inplace(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
     108     7514435 : { Flj_add_indir_pre(P, Q, P, a4, p, pi); }
     109             : 
     110             : GEN
     111       91290 : Flj_add_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
     112             : {
     113       91290 :   GEN R = cgetg(4, t_VECSMALL);
     114       91290 :   Flj_add_indir_pre(P, Q, R, a4, p, pi);
     115       91290 :   return R;
     116             : }
     117             : 
     118             : GEN
     119     2446207 : Flj_neg(GEN Q, ulong p)
     120     2446207 : { return mkvecsmall3(Q[1], Fl_neg(Q[2], p), Q[3]); }
     121             : 
     122             : typedef struct {
     123             :   ulong pbits, nbits;  /* Positive bits and negative bits */
     124             :   ulong lnzb; /* Leading nonzero bit */
     125             : } naf_t;
     126             : 
     127             : /* Return the signed binary representation (i.e. the Non-Adjacent Form
     128             :  * in base 2) of a; a = x.pbits - x.nbits (+ 2^BILif < 0; this
     129             :  * exceptional case can happen if a > 2^(BIL-1)) */
     130             : static void
     131     4457373 : naf_repr(naf_t *x, ulong a)
     132             : {
     133     4457373 :   ulong pbits = 0, nbits = 0, c0 = 0, c1, a0;
     134             :   long t, i;
     135             : 
     136    44843619 :   for (i = 0; a; a >>= 1, ++i) {
     137    40386246 :     a0 = a & 1;
     138    40386246 :     c1 = (c0 + a0 + ((a & 2) >> 1)) >> 1;
     139    40386246 :     t = c0 + a0 - (c1 << 1);
     140    40386246 :     if (t < 0)
     141     6264968 :       nbits |= (1UL << i);
     142    34121278 :     else if (t > 0)
     143     7133252 :       pbits |= (1UL << i);
     144    40386246 :     c0 = c1;
     145             :   }
     146     4457373 :   c1 = c0 >> 1;
     147     4457373 :   t = c0 - (c1 << 1);
     148             :   /* since a >= 0, we have t >= 0; if i = BIL, pbits (virtually) overflows;
     149             :    * that leading overflowed bit is implied and not recorded in pbits */
     150     4457373 :   if (t > 0 && i != BITS_IN_LONG) pbits |= (1UL << i);
     151     4457373 :   x->pbits = pbits;
     152     4457373 :   x->nbits = nbits;
     153     4457373 :   x->lnzb = t? i-2: i-3;
     154     4457373 : }
     155             : 
     156             : /* Standard left-to-right signed double-and-add to compute [n]P. */
     157             : static GEN
     158     3981219 : Flj_mulu_pre_naf(GEN P, ulong n, ulong a4, ulong p, ulong pi, const naf_t *x)
     159             : {
     160             :   GEN R, Pi;
     161             :   ulong pbits, nbits;
     162             :   ulong m;
     163             : 
     164     3981219 :   if (n == 0) return mkvecsmall3(1, 1, 0);
     165     3981219 :   if (n == 1) return Flv_copy(P);
     166             : 
     167     3919677 :   R = Flj_dbl_pre(P, a4, p, pi);
     168     3918558 :   if (n == 2) return R;
     169             : 
     170     3082144 :   pbits = x->pbits;
     171     3082144 :   nbits = x->nbits;
     172     3082144 :   Pi = nbits? Flj_neg(P, p): NULL;
     173     3083043 :   m = (1UL << x->lnzb);
     174    24297881 :   for ( ; m; m >>= 1) {
     175    21214373 :     Flj_dbl_pre_inplace(R, a4, p, pi);
     176    21214068 :     if (m & pbits)
     177     3295220 :       Flj_add_pre_inplace(R, P, a4, p, pi);
     178    17918848 :     else if (m & nbits)
     179     4222087 :       Flj_add_pre_inplace(R, Pi, a4, p, pi);
     180             :   }
     181     3083508 :   return gc_const((pari_sp)R, R);
     182             : }
     183             : 
     184             : GEN
     185     2749798 : Flj_mulu_pre(GEN P, ulong n, ulong a4, ulong p, ulong pi)
     186             : {
     187     2749798 :   naf_t x; naf_repr(&x, n);
     188     2750117 :   return Flj_mulu_pre_naf(P, n, a4, p, pi, &x);
     189             : }
     190             : 
     191             : struct _Flj { ulong a4, p, pi; };
     192             : 
     193             : static GEN
     194       91290 : _Flj_add(void *E, GEN P, GEN Q)
     195             : {
     196       91290 :   struct _Flj *ell=(struct _Flj *) E;
     197       91290 :   return Flj_add_pre(P, Q, ell->a4, ell->p, ell->pi);
     198             : }
     199             : 
     200             : static GEN
     201      106612 : _Flj_mul(void *E, GEN P, GEN n)
     202             : {
     203      106612 :   struct _Flj *ell = (struct _Flj *) E;
     204      106612 :   long s = signe(n);
     205             :   GEN Q;
     206      106612 :   if (s==0) return mkvecsmall3(1, 1, 0);
     207      106612 :   Q = Flj_mulu_pre(P, itou(n), ell->a4, ell->p, ell->pi);
     208      106612 :   return s>0 ? Q : Flj_neg(Q, ell->p);
     209             : }
     210             : static GEN
     211           0 : _Flj_one(void *E)
     212           0 : { (void) E; return mkvecsmall3(1, 1, 0); }
     213             : 
     214             : GEN
     215       15322 : FljV_factorback_pre(GEN P, GEN L, ulong a4, ulong p, ulong pi)
     216             : {
     217             :   struct _Flj E;
     218       15322 :   E.a4 = a4; E.p = p; E.pi = pi;
     219       15322 :   return gen_factorback(P, L, (void*)&E, &_Flj_add, &_Flj_mul, &_Flj_one);
     220             : }
     221             : 
     222             : ulong
     223      291224 : Flj_order_ufact(GEN P, ulong n, GEN fa, ulong a4, ulong p, ulong pi)
     224             : {
     225      291224 :   pari_sp av = avma;
     226      291224 :   GEN T = gel(fa,1), E = gel(fa,2);
     227      291224 :   long i, l = lg(T);
     228      291224 :   ulong res = 1;
     229             : 
     230      821763 :   for (i = 1; i < l; i++, set_avma(av))
     231             :   {
     232      656838 :     ulong j, t = T[i], e = E[i];
     233      656838 :     GEN b = P;
     234      656838 :     naf_t x; naf_repr(&x, t);
     235      656954 :     if (l != 2) b = Flj_mulu_pre(b, n / upowuu(t,e), a4, p, pi);
     236     1888215 :     for (j = 0; j < e && b[3]; j++) b = Flj_mulu_pre_naf(b, t, a4, p, pi, &x);
     237      656831 :     if (b[3]) return 0;
     238      530524 :     res *= upowuu(t, j);
     239             :   }
     240      164925 :   return res;
     241             : }
     242             : 
     243             : GEN
     244     1577405 : Fle_to_Flj(GEN P)
     245     3154744 : { return ell_is_inf(P) ? mkvecsmall3(1UL, 1UL, 0UL):
     246     1577380 :                          mkvecsmall3(P[1], P[2], 1UL); }
     247             : 
     248             : GEN
     249       68257 : Flj_to_Fle(GEN P, ulong p)
     250             : {
     251       68257 :   if (P[3] == 0) return ellinf();
     252             :   else
     253             :   {
     254       67557 :     ulong Z = Fl_inv(P[3], p);
     255       67557 :     ulong Z2 = Fl_sqr(Z, p);
     256       67557 :     ulong X3 = Fl_mul(P[1], Z2, p);
     257       67557 :     ulong Y3 = Fl_mul(P[2], Fl_mul(Z, Z2, p), p);
     258       67557 :     return mkvecsmall2(X3, Y3);
     259             :   }
     260             : }
     261             : 
     262             : GEN
     263     1592597 : Flj_to_Fle_pre(GEN P, ulong p, ulong pi)
     264             : {
     265     1592597 :   if (P[3] == 0) return ellinf();
     266             :   else
     267             :   {
     268     1406087 :     ulong Z = Fl_inv(P[3], p);
     269     1406337 :     ulong Z2 = Fl_sqr_pre(Z, p, pi);
     270     1406218 :     ulong X3 = Fl_mul_pre(P[1], Z2, p, pi);
     271     1406155 :     ulong Y3 = Fl_mul_pre(P[2], Fl_mul_pre(Z, Z2, p, pi), p, pi);
     272     1406160 :     return mkvecsmall2(X3, Y3);
     273             :   }
     274             : }
     275             : 
     276             : INLINE void
     277     8263398 : random_Fle_pre_indir(ulong a4, ulong a6, ulong p, ulong pi,
     278             :                      ulong *pt_x, ulong *pt_y)
     279             : {
     280             :   ulong x, x2, y, rhs;
     281             :   do
     282             :   {
     283     8263398 :     x   = random_Fl(p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     284     8273036 :     x2  = Fl_sqr_pre(x, p, pi);
     285     8253938 :     rhs = Fl_addmul_pre(a6, x, Fl_add(x2, a4, p), p, pi);
     286     8252727 :   } while ((!rhs && !Fl_add(Fl_triple(x2,p),a4,p)) || krouu(rhs, p) < 0);
     287     4139811 :   y = Fl_sqrt_pre(rhs, p, pi);
     288     4140795 :   *pt_x = x; *pt_y = y;
     289     4140795 : }
     290             : 
     291             : GEN
     292      490470 : random_Flj_pre(ulong a4, ulong a6, ulong p, ulong pi)
     293             : {
     294             :   ulong x, y;
     295      490470 :   random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
     296      490468 :   return mkvecsmall3(x, y, 1);
     297             : }
     298             : 
     299             : GEN
     300      153252 : Flj_changepointinv_pre(GEN P, GEN ch, ulong p, ulong pi)
     301             : {
     302             :   ulong c, u, r, s, t, u2, u3, z2;
     303      153252 :   ulong x  = uel(P,1), y = uel(P,2), z = uel(P,3);
     304             :   GEN w;
     305      153252 :   if (z == 0) return Flv_copy(P);
     306      153224 :   u = ch[1]; r = ch[2];
     307      153224 :   s = ch[3]; t = ch[4];
     308      153224 :   u2 = Fl_sqr_pre(u, p, pi); u3 = Fl_mul_pre(u, u2, p, pi);
     309      153224 :   c = Fl_mul_pre(u2, x, p, pi);
     310      153224 :   z2 = Fl_sqr_pre(z, p, pi);
     311      153224 :   w = cgetg(4, t_VECSMALL);
     312      153224 :   uel(w,1) = Fl_add(c, Fl_mul_pre(r, z2, p, pi), p);
     313      153224 :   uel(w,2) = Fl_add(Fl_mul_pre(u3 ,y, p, pi),
     314             :                     Fl_mul_pre(z, Fl_add(Fl_mul_pre(s,c,p,pi),
     315             :                                          Fl_mul_pre(z2,t,p,pi), p), p, pi), p);
     316      153224 :   uel(w,3) = z;
     317      153224 :   return w;
     318             : }
     319             : 
     320             : /***********************************************************************/
     321             : /**                              Fle                                  **/
     322             : /***********************************************************************/
     323             : GEN
     324       16041 : Fle_changepoint(GEN P, GEN ch, ulong p)
     325             : {
     326             :   ulong c, u, r, s, t, v, v2, v3;
     327             :   GEN z;
     328       16041 :   if (ell_is_inf(P)) return ellinf();
     329       16041 :   u = ch[1]; r = ch[2];
     330       16041 :   s = ch[3]; t = ch[4];
     331       16041 :   v = Fl_inv(u, p); v2 = Fl_sqr(v,p); v3 = Fl_mul(v,v2,p);
     332       16041 :   c = Fl_sub(uel(P,1),r,p);
     333       16041 :   z = cgetg(3,t_VECSMALL);
     334       16041 :   z[1] = Fl_mul(v2, c, p);
     335       16041 :   z[2] = Fl_mul(v3, Fl_sub(uel(P,2), Fl_add(Fl_mul(s,c, p),t, p),p),p);
     336       16041 :   return z;
     337             : }
     338             : 
     339             : GEN
     340      134231 : Fle_changepointinv(GEN P, GEN ch, ulong p)
     341             : {
     342             :   ulong c, u, r, s, t, u2, u3;
     343             :   GEN z;
     344      134231 :   if (ell_is_inf(P)) return ellinf();
     345      133531 :   u = ch[1]; r = ch[2];
     346      133531 :   s = ch[3]; t = ch[4];
     347      133531 :   u2 = Fl_sqr(u, p); u3 = Fl_mul(u,u2,p);
     348      133531 :   c = Fl_mul(u2,uel(P,1), p);
     349      133531 :   z = cgetg(3, t_VECSMALL);
     350      133531 :   z[1] = Fl_add(c,r,p);
     351      133531 :   z[2] = Fl_add(Fl_mul(u3,uel(P,2),p), Fl_add(Fl_mul(s,c,p), t, p), p);
     352      133531 :   return z;
     353             : }
     354             : static GEN
     355      715240 : Fle_dbl_slope(GEN P, ulong a4, ulong p, ulong *slope)
     356             : {
     357             :   ulong x, y, Qx, Qy;
     358      715240 :   if (ell_is_inf(P) || !P[2]) return ellinf();
     359      645006 :   x = P[1]; y = P[2];
     360      645006 :   *slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
     361             :                   Fl_double(y, p), p);
     362      645008 :   Qx = Fl_sub(Fl_sqr(*slope, p), Fl_double(x, p), p);
     363      645001 :   Qy = Fl_sub(Fl_mul(*slope, Fl_sub(x, Qx, p), p), y, p);
     364      644996 :   return mkvecsmall2(Qx, Qy);
     365             : }
     366             : 
     367             : GEN
     368      466562 : Fle_dbl(GEN P, ulong a4, ulong p)
     369             : {
     370             :   ulong slope;
     371      466562 :   return Fle_dbl_slope(P,a4,p,&slope);
     372             : }
     373             : 
     374             : static GEN
     375     1514073 : Fle_add_slope(GEN P, GEN Q, ulong a4, ulong p, ulong *slope)
     376             : {
     377             :   ulong Px, Py, Qx, Qy, Rx, Ry;
     378     1514073 :   if (ell_is_inf(P)) return Q;
     379     1514079 :   if (ell_is_inf(Q)) return P;
     380     1514083 :   Px = P[1]; Py = P[2];
     381     1514083 :   Qx = Q[1]; Qy = Q[2];
     382     1514083 :   if (Px==Qx) return Py==Qy ? Fle_dbl_slope(P, a4, p, slope): ellinf();
     383     1382891 :   *slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
     384     1382996 :   Rx = Fl_sub(Fl_sub(Fl_sqr(*slope, p), Px, p), Qx, p);
     385     1382963 :   Ry = Fl_sub(Fl_mul(*slope, Fl_sub(Px, Rx, p), p), Py, p);
     386     1382963 :   return mkvecsmall2(Rx, Ry);
     387             : }
     388             : 
     389             : GEN
     390     1498472 : Fle_add(GEN P, GEN Q, ulong a4, ulong p)
     391             : {
     392             :   ulong slope;
     393     1498472 :   return Fle_add_slope(P,Q,a4,p,&slope);
     394             : }
     395             : 
     396             : static GEN
     397      131208 : Fle_neg(GEN P, ulong p)
     398             : {
     399      131208 :   if (ell_is_inf(P)) return P;
     400      131208 :   return mkvecsmall2(P[1], Fl_neg(P[2], p));
     401             : }
     402             : 
     403             : GEN
     404           0 : Fle_sub(GEN P, GEN Q, ulong a4, ulong p)
     405             : {
     406           0 :   pari_sp av = avma;
     407             :   ulong slope;
     408           0 :   return gerepileupto(av, Fle_add_slope(P, Fle_neg(Q, p), a4, p, &slope));
     409             : }
     410             : 
     411             : struct _Fle { ulong a4, a6, p; };
     412             : 
     413             : static GEN
     414           0 : _Fle_dbl(void *E, GEN P)
     415             : {
     416           0 :   struct _Fle *ell = (struct _Fle *) E;
     417           0 :   return Fle_dbl(P, ell->a4, ell->p);
     418             : }
     419             : 
     420             : static GEN
     421      365533 : _Fle_add(void *E, GEN P, GEN Q)
     422             : {
     423      365533 :   struct _Fle *ell=(struct _Fle *) E;
     424      365533 :   return Fle_add(P, Q, ell->a4, ell->p);
     425             : }
     426             : 
     427             : GEN
     428     2021210 : Fle_mulu(GEN P, ulong n, ulong a4, ulong p)
     429             : {
     430             :   ulong pi;
     431     2021210 :   if (!n || ell_is_inf(P)) return ellinf();
     432     2021202 :   if (n==1) return zv_copy(P);
     433     2008476 :   if (n==2) return Fle_dbl(P, a4, p);
     434     1577197 :   pi = get_Fl_red(p);
     435     1577404 :   return Flj_to_Fle_pre(Flj_mulu_pre(Fle_to_Flj(P), n, a4, p, pi), p, pi);
     436             : }
     437             : 
     438             : static GEN
     439     1512573 : _Fle_mul(void *E, GEN P, GEN n)
     440             : {
     441     1512573 :   pari_sp av = avma;
     442     1512573 :   struct _Fle *e=(struct _Fle *) E;
     443     1512573 :   long s = signe(n);
     444             :   GEN Q;
     445     1512573 :   if (!s || ell_is_inf(P)) return ellinf();
     446     1493719 :   if (s < 0) P = Fle_neg(P, e->p);
     447     1493719 :   if (is_pm1(n)) return s > 0? zv_copy(P): P;
     448     1166846 :   Q = (lgefint(n)==3) ? Fle_mulu(P, uel(n,2), e->a4, e->p):
     449           0 :                         gen_pow(P, n, (void*)e, &_Fle_dbl, &_Fle_add);
     450     1166926 :   return s > 0? Q: gerepileuptoleaf(av, Q);
     451             : }
     452             : 
     453             : GEN
     454       28959 : Fle_mul(GEN P, GEN n, ulong a4, ulong p)
     455             : {
     456             :   struct _Fle E;
     457       28959 :   E.a4 = a4; E.p = p;
     458       28959 :   return _Fle_mul(&E, P, n);
     459             : }
     460             : 
     461             : /* Finds a random nonsingular point on E */
     462             : GEN
     463     3648908 : random_Fle_pre(ulong a4, ulong a6, ulong p, ulong pi)
     464             : {
     465             :   ulong x, y;
     466     3648908 :   random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
     467     3650720 :   return mkvecsmall2(x, y);
     468             : }
     469             : 
     470             : GEN
     471       26040 : random_Fle(ulong a4, ulong a6, ulong p)
     472       26040 : { return random_Fle_pre(a4, a6, p, get_Fl_red(p)); }
     473             : 
     474             : static GEN
     475           0 : _Fle_rand(void *E)
     476             : {
     477           0 :   struct _Fle *e=(struct _Fle *) E;
     478           0 :   return random_Fle(e->a4, e->a6, e->p);
     479             : }
     480             : 
     481             : static const struct bb_group Fle_group={_Fle_add,_Fle_mul,_Fle_rand,hash_GEN,zv_equal,ell_is_inf,NULL};
     482             : 
     483             : GEN
     484      288787 : Fle_order(GEN z, GEN o, ulong a4, ulong p)
     485             : {
     486      288787 :   pari_sp av = avma;
     487             :   struct _Fle e;
     488      288787 :   e.a4=a4;
     489      288787 :   e.p=p;
     490      288787 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &Fle_group));
     491             : }
     492             : 
     493             : GEN
     494       54327 : Fle_log(GEN a, GEN b, GEN o, ulong a4, ulong p)
     495             : {
     496       54327 :   pari_sp av = avma;
     497             :   struct _Fle e;
     498       54327 :   e.a4=a4;
     499       54327 :   e.p=p;
     500       54327 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &Fle_group));
     501             : }
     502             : 
     503             : ulong
     504           0 : Fl_ellj(ulong a4, ulong a6, ulong p)
     505             : {
     506           0 :   if (SMALL_ULONG(p))
     507             :   { /* a43 = 4 a4^3 */
     508           0 :     ulong a43 = Fl_double(Fl_double(Fl_mul(a4, Fl_sqr(a4, p), p), p), p);
     509             :     /* a62 = 27 a6^2 */
     510           0 :     ulong a62 = Fl_mul(Fl_sqr(a6, p), 27 % p, p);
     511           0 :     ulong z1 = Fl_mul(a43, 1728 % p, p);
     512           0 :     ulong z2 = Fl_add(a43, a62, p);
     513           0 :     return Fl_div(z1, z2, p);
     514             :   }
     515           0 :   return Fl_ellj_pre(a4, a6, p, get_Fl_red(p));
     516             : }
     517             : 
     518             : void
     519      174495 : Fl_ellj_to_a4a6(ulong j, ulong p, ulong *pt_a4, ulong *pt_a6)
     520             : {
     521      174495 :   ulong zagier = 1728 % p;
     522      174495 :   if (j == 0)           { *pt_a4 = 0; *pt_a6 =1; }
     523      174481 :   else if (j == zagier) { *pt_a4 = 1; *pt_a6 =0; }
     524             :   else
     525             :   {
     526      174467 :     ulong k = Fl_sub(zagier, j, p);
     527      174457 :     ulong kj = Fl_mul(k, j, p);
     528      174451 :     ulong k2j = Fl_mul(kj, k, p);
     529      174441 :     *pt_a4 = Fl_triple(kj, p);
     530      174443 :     *pt_a6 = Fl_double(k2j, p);
     531             :   }
     532      174467 : }
     533             : 
     534             : ulong
     535     2116438 : Fl_elldisc_pre(ulong a4, ulong a6, ulong p, ulong pi)
     536             : { /* D = -(4A^3 + 27B^2) */
     537             :   ulong t1, t2;
     538     2116438 :   t1 = Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi);
     539     2114959 :   t1 = Fl_double(Fl_double(t1, p), p);
     540     2114551 :   t2 = Fl_mul_pre(27 % p, Fl_sqr_pre(a6, p, pi), p, pi);
     541     2114521 :   return Fl_neg(Fl_add(t1, t2, p), p);
     542             : }
     543             : 
     544             : ulong
     545           0 : Fl_elldisc(ulong a4, ulong a6, ulong p)
     546             : {
     547           0 :   if (SMALL_ULONG(p))
     548             :   { /* D = -(4A^3 + 27B^2) */
     549             :     ulong t1, t2;
     550           0 :     t1 = Fl_mul(a4, Fl_sqr(a4, p), p);
     551           0 :     t1 = Fl_double(Fl_double(t1, p), p);
     552           0 :     t2 = Fl_mul(27 % p, Fl_sqr(a6, p), p);
     553           0 :     return Fl_neg(Fl_add(t1, t2, p), p);
     554             :   }
     555           0 :   return Fl_elldisc_pre(a4, a6, p, get_Fl_red(p));
     556             : }
     557             : 
     558             : void
     559      105135 : Fl_elltwist_disc(ulong a4, ulong a6, ulong D, ulong p, ulong *pa4, ulong *pa6)
     560             : {
     561      105135 :   ulong D2 = Fl_sqr(D, p);
     562      105135 :   *pa4 = Fl_mul(a4, D2, p);
     563      105135 :   *pa6 = Fl_mul(a6, Fl_mul(D, D2, p), p);
     564      105134 : }
     565             : 
     566             : void
     567           0 : Fl_elltwist(ulong a4, ulong a6, ulong p, ulong *pt_a4, ulong *pt_a6)
     568           0 : { Fl_elltwist_disc(a4, a6, nonsquare_Fl(p), p, pt_a4, pt_a6); }
     569             : 
     570             : static void
     571    52402051 : Fle_dbl_sinv_pre_inplace(GEN P, ulong a4, ulong sinv, ulong p, ulong pi)
     572             : {
     573             :   ulong x, y, slope;
     574    52402051 :   if (uel(P,1)==p) return;
     575    52159193 :   if (!P[2]) { P[1] = p; return; }
     576    52011157 :   x = P[1]; y = P[2];
     577    52011157 :   slope = Fl_mul_pre(Fl_add(Fl_triple(Fl_sqr_pre(x, p, pi), p), a4, p),
     578             :                 sinv, p, pi);
     579    52015553 :   P[1] = Fl_sub(Fl_sqr_pre(slope, p, pi), Fl_double(x, p), p);
     580    51913760 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(x, P[1], p), p, pi), y, p);
     581             : }
     582             : 
     583             : static void
     584     7217425 : Fle_add_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
     585             : {
     586             :   ulong Px, Py, Qx, Qy, slope;
     587     7217425 :   if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Q[2]; }
     588     7217425 :   if (ell_is_inf(Q)) return;
     589     7216846 :   Px = P[1]; Py = P[2];
     590     7216846 :   Qx = Q[1]; Qy = Q[2];
     591     7216846 :   if (Px==Qx)
     592             :   {
     593       29243 :     if (Py==Qy) Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
     594       14027 :     else P[1] = p;
     595       29243 :     return;
     596             :   }
     597     7187603 :   slope = Fl_mul_pre(Fl_sub(Py, Qy, p), sinv, p, pi);
     598     7186954 :   P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
     599     7177339 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
     600             : }
     601             : 
     602             : static void
     603     8052669 : Fle_sub_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
     604             : {
     605             :   ulong Px, Py, Qx, Qy, slope;
     606     8052669 :   if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Fl_neg(Q[2], p); }
     607     8052669 :   if (ell_is_inf(Q)) return;
     608     8054837 :   Px = P[1]; Py = P[2];
     609     8054837 :   Qx = Q[1]; Qy = Q[2];
     610     8054837 :   if (Px==Qx)
     611             :   {
     612       34939 :     if (Py==Qy) P[1] = p;
     613             :     else
     614       13668 :       Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
     615       34939 :     return;
     616             :   }
     617     8019898 :   slope = Fl_mul_pre(Fl_add(Py, Qy, p), sinv, p, pi);
     618     8012571 :   P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
     619     7996364 :   P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
     620             : }
     621             : 
     622             : static long
     623    67909866 : skipzero(long n) { return n ? n:1; }
     624             : 
     625             : void
     626     1797364 : FleV_add_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
     627             : {
     628     1797364 :   long i, l=lg(a4);
     629     1797364 :   GEN sinv = cgetg(l, t_VECSMALL);
     630     9037384 :   for(i=1; i<l; i++)
     631     7241075 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
     632     1796309 :   Flv_inv_pre_inplace(sinv, p, pi);
     633     9013094 :   for (i=1; i<l; i++)
     634     7218469 :     Fle_add_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
     635     1794625 : }
     636             : 
     637             : void
     638     2065919 : FleV_sub_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
     639             : {
     640     2065919 :   long i, l=lg(a4);
     641     2065919 :   GEN sinv = cgetg(l, t_VECSMALL);
     642    10137328 :   for(i=1; i<l; i++)
     643     8071658 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
     644     2065670 :   Flv_inv_pre_inplace(sinv, p, pi);
     645    10116007 :   for (i=1; i<l; i++)
     646     8053108 :     Fle_sub_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
     647     2062899 : }
     648             : 
     649             : void
     650    13540353 : FleV_dbl_pre_inplace(GEN P, GEN a4, ulong p, ulong pi)
     651             : {
     652    13540353 :   long i, l=lg(a4);
     653    13540353 :   GEN sinv = cgetg(l, t_VECSMALL);
     654    66874125 :   for(i=1; i<l; i++)
     655    53355942 :     uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_double(umael(P,i,2), p));
     656    13518183 :   Flv_inv_pre_inplace(sinv, p, pi);
     657    65912975 :   for(i=1; i<l; i++)
     658    52428562 :     Fle_dbl_sinv_pre_inplace(gel(P,i), uel(a4,i), uel(sinv,i), p, pi);
     659    13484413 : }
     660             : 
     661             : static void
     662     1051867 : FleV_mulu_pre_naf_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi, const naf_t *x)
     663             : {
     664     1051867 :   pari_sp av = avma;
     665             :   ulong pbits, nbits, m;
     666             :   GEN R;
     667     1051867 :   if (n == 1) return;
     668             : 
     669     1051867 :   R = P; P = gcopy(P);
     670     1051597 :   FleV_dbl_pre_inplace(R, a4, p, pi);
     671     1050290 :   if (n == 2) return;
     672             : 
     673     1050196 :   pbits = x->pbits;
     674     1050196 :   nbits = x->nbits;
     675     1050196 :   m = (1UL << x->lnzb);
     676    13563180 :   for ( ; m; m >>= 1) {
     677    12512714 :     FleV_dbl_pre_inplace(R, a4, p, pi);
     678    12505353 :     if (m & pbits)
     679     1797494 :       FleV_add_pre_inplace(R, P, a4, p, pi);
     680    10707859 :     else if (m & nbits)
     681     2066219 :       FleV_sub_pre_inplace(R, P, a4, p, pi);
     682             :   }
     683     1050466 :   set_avma(av);
     684             : }
     685             : 
     686             : void
     687     1051348 : FleV_mulu_pre_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi)
     688             : {
     689     1051348 :   naf_t x; naf_repr(&x, n);
     690     1051880 :   FleV_mulu_pre_naf_inplace(P, n, a4, p, pi, &x);
     691     1050120 : }
     692             : 
     693             : /***********************************************************************/
     694             : /**                            Pairings                               **/
     695             : /**               Derived from APIP by Jerome Milan, 2012             **/
     696             : /***********************************************************************/
     697             : static ulong
     698      419084 : Fle_vert(GEN P, GEN Q, ulong a4, ulong p)
     699             : {
     700      419084 :   if (ell_is_inf(P))
     701      154989 :     return 1;
     702      264095 :   if (uel(Q, 1) != uel(P, 1))
     703      244720 :     return Fl_sub(uel(Q, 1), uel(P, 1), p);
     704       19375 :   if (uel(P,2) != 0 ) return 1;
     705       14435 :   return Fl_inv(Fl_add(Fl_triple(Fl_sqr(uel(P,1),p), p), a4, p), p);
     706             : }
     707             : 
     708             : static ulong
     709      138991 : Fle_Miller_line(GEN R, GEN Q, ulong slope, ulong a4, ulong p)
     710             : {
     711      138991 :   ulong x = uel(Q, 1), y = uel(Q, 2);
     712      138991 :   ulong tmp1 = Fl_sub(x, uel(R, 1), p);
     713      138991 :   ulong tmp2 = Fl_add(Fl_mul(tmp1, slope, p), uel(R,2), p);
     714      138991 :   if (y != tmp2)
     715      128279 :     return Fl_sub(y, tmp2, p);
     716       10712 :   if (y == 0)
     717        6934 :     return 1;
     718             :   else
     719             :   {
     720             :     ulong s1, s2;
     721        3778 :     ulong y2i = Fl_inv(Fl_double(y, p), p);
     722        3778 :     s1 = Fl_mul(Fl_add(Fl_triple(Fl_sqr(x, p), p), a4, p), y2i, p);
     723        3778 :     if (s1 != slope)
     724        2119 :       return Fl_sub(s1, slope, p);
     725        1659 :     s2 = Fl_mul(Fl_sub(Fl_triple(x, p), Fl_sqr(s1, p), p), y2i, p);
     726        1659 :     return s2 ? s2: y2i;
     727             :   }
     728             : }
     729             : 
     730             : /* Computes the equation of the line tangent to R and returns its
     731             :  * evaluation at the point Q. Also doubles the point R. */
     732             : static ulong
     733      258513 : Fle_tangent_update(GEN R, GEN Q, ulong a4, ulong p, GEN *pt_R)
     734             : {
     735      258513 :   if (ell_is_inf(R)) { *pt_R = ellinf(); return 1; }
     736      226122 :   else if (uel(R,2) == 0) { *pt_R = ellinf(); return Fle_vert(R, Q, a4, p); }
     737             :   else
     738             :   {
     739             :     ulong slope;
     740      123397 :     *pt_R = Fle_dbl_slope(R, a4, p, &slope);
     741      123397 :     return Fle_Miller_line(R, Q, slope, a4, p);
     742             :   }
     743             : }
     744             : 
     745             : /* Computes the equation of the line through R and P, and returns its
     746             :  * evaluation at the point Q. Also adds P to the point R. */
     747             : static ulong
     748       36720 : Fle_chord_update(GEN R, GEN P, GEN Q, ulong a4, ulong p, GEN *pt_R)
     749             : {
     750       36720 :   if (ell_is_inf(R)) { *pt_R = P; return Fle_vert(P, Q, a4, p); }
     751       35467 :   else if (ell_is_inf(P)) { *pt_R = R; return Fle_vert(R, Q, a4, p); }
     752       35467 :   else if (uel(P, 1) == uel(R, 1))
     753             :   {
     754       19873 :     if (uel(P, 2) == uel(R, 2)) return Fle_tangent_update(R, Q, a4, p, pt_R);
     755       19873 :     else { *pt_R = ellinf(); return Fle_vert(R, Q, a4, p); }
     756             :   }
     757             :   else
     758             :   {
     759             :     ulong slope;
     760       15594 :     *pt_R = Fle_add_slope(P, R, a4, p, &slope);
     761       15594 :     return Fle_Miller_line(R, Q, slope, a4, p);
     762             :   }
     763             : }
     764             : 
     765             : struct _Fle_miller { ulong p, a4; GEN P; };
     766             : static GEN
     767      258513 : Fle_Miller_dbl(void* E, GEN d)
     768             : {
     769      258513 :   struct _Fle_miller *m = (struct _Fle_miller *)E;
     770      258513 :   ulong p = m->p, a4 = m->a4;
     771      258513 :   GEN P = m->P;
     772             :   ulong v, line;
     773      258513 :   ulong N = Fl_sqr(umael(d,1,1), p);
     774      258513 :   ulong D = Fl_sqr(umael(d,1,2), p);
     775      258513 :   GEN point = gel(d,2);
     776      258513 :   line = Fle_tangent_update(point, P, a4, p, &point);
     777      258513 :   N  = Fl_mul(N, line, p);
     778      258513 :   v = Fle_vert(point, P, a4, p);
     779      258513 :   D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
     780             : }
     781             : static GEN
     782       36720 : Fle_Miller_add(void* E, GEN va, GEN vb)
     783             : {
     784       36720 :   struct _Fle_miller *m = (struct _Fle_miller *)E;
     785       36720 :   ulong p = m->p, a4= m->a4;
     786       36720 :   GEN P = m->P, point;
     787             :   ulong v, line;
     788       36720 :   ulong na = umael(va,1,1), da = umael(va,1,2);
     789       36720 :   ulong nb = umael(vb,1,1), db = umael(vb,1,2);
     790       36720 :   GEN pa = gel(va,2), pb = gel(vb,2);
     791       36720 :   ulong N = Fl_mul(na, nb, p);
     792       36720 :   ulong D = Fl_mul(da, db, p);
     793       36720 :   line = Fle_chord_update(pa, pb, P, a4, p, &point);
     794       36720 :   N = Fl_mul(N, line, p);
     795       36720 :   v = Fle_vert(point, P, a4, p);
     796       36720 :   D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
     797             : }
     798             : 
     799             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     800             :  * the standard Miller algorithm. */
     801             : static ulong
     802      121345 : Fle_Miller(GEN Q, GEN P, ulong m, ulong a4, ulong p)
     803             : {
     804      121345 :   pari_sp av = avma;
     805             :   struct _Fle_miller d;
     806             :   GEN v;
     807             :   ulong N, D;
     808             : 
     809      121345 :   d.a4 = a4; d.p = p; d.P = P;
     810      121345 :   v = gen_powu_i(mkvec2(mkvecsmall2(1,1), Q), m, (void*)&d,
     811             :                 Fle_Miller_dbl, Fle_Miller_add);
     812      121345 :   N = umael(v,1,1); D = umael(v,1,2);
     813      121345 :   return gc_ulong(av, Fl_div(N, D, p));
     814             : }
     815             : 
     816             : ulong
     817       53361 : Fle_weilpairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
     818             : {
     819       53361 :   pari_sp ltop = avma;
     820             :   ulong N, D, w;
     821       53361 :   if (ell_is_inf(P) || ell_is_inf(Q) || zv_equal(P,Q)) return 1;
     822       52234 :   N = Fle_Miller(P, Q, m, a4, p);
     823       52234 :   D = Fle_Miller(Q, P, m, a4, p);
     824       52234 :   w = Fl_div(N, D, p);
     825       52234 :   if (odd(m)) w  = Fl_neg(w, p);
     826       52234 :   return gc_ulong(ltop, w);
     827             : }
     828             : 
     829             : ulong
     830       16877 : Fle_tatepairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
     831             : {
     832       16877 :   if (ell_is_inf(P) || ell_is_inf(Q)) return 1;
     833       16877 :   return Fle_Miller(P, Q, m, a4, p);
     834             : }
     835             : 
     836             : GEN
     837       11438 : Fl_ellptors(ulong l, ulong N, ulong a4, ulong a6, ulong p)
     838             : {
     839       11438 :   long v = z_lval(N, l);
     840             :   ulong N0, N1;
     841             :   GEN F;
     842       11438 :   if (v==0) return cgetg(1,t_VEC);
     843       11438 :   N0 = upowuu(l, v); N1 = N/N0;
     844       11438 :   F = mkmat2(mkcols(l), mkcols(v));
     845             :   while(1)
     846        1582 :   {
     847       13020 :     pari_sp av2 = avma;
     848             :     GEN P, Q;
     849             :     ulong d, s, t;
     850             : 
     851       13020 :     P = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
     852       13020 :     Q = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
     853       13020 :     s = itou(Fle_order(P, F, a4, p));
     854       13020 :     t = itou(Fle_order(Q, F, a4, p));
     855       13020 :     if (s < t) { swap(P,Q); lswap(s,t); }
     856       13020 :     if (s == N0) retmkvec(Fle_mulu(P, s/l, a4, p));
     857        2394 :     d = Fl_order(Fle_weilpairing(P, Q, s, a4, p), s, p);
     858        2394 :     if (s*d == N0)
     859         812 :       retmkvec2(Fle_mulu(P, s/l, a4, p), Fle_mulu(Q, t/l, a4, p));
     860        1582 :     set_avma(av2);
     861             :   }
     862             : }

Generated by: LCOV version 1.14