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 - elltors.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21061-8feaff2) Lines: 392 403 97.3 %
Date: 2017-09-24 06:24:57 Functions: 22 24 91.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  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             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**         TORSION OF ELLIPTIC CURVES over NUMBER FIELDS          **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : static int
      22         336 : smaller_x(GEN p, GEN q)
      23             : {
      24         336 :   int s = abscmpii(denom(p), denom(q));
      25         336 :   return (s<0 || (s==0 && abscmpii(numer(p),numer(q)) < 0));
      26             : }
      27             : 
      28             : /* best generator in cycle of length k */
      29             : static GEN
      30         448 : best_in_cycle(GEN e, GEN p, long k)
      31             : {
      32         448 :   GEN p0 = p,q = p;
      33             :   long i;
      34             : 
      35         721 :   for (i=2; i+i<k; i++)
      36             :   {
      37         273 :     q = elladd(e,q,p0);
      38         273 :     if (ugcd(i,k)==1 && smaller_x(gel(q,1), gel(p,1))) p = q;
      39             :   }
      40         448 :   return (gsigne(ec_dmFdy_evalQ(e,p)) < 0)? ellneg(e,p): p;
      41             : }
      42             : 
      43             : /* <p,q> = E_tors, possibly NULL (= oo), p,q independent unless NULL
      44             :  * order p = k, order q = 2 unless NULL */
      45             : static GEN
      46         588 : tors(GEN e, long k, GEN p, GEN q, GEN v)
      47             : {
      48             :   GEN r;
      49         588 :   if (q)
      50             :   {
      51          98 :     long n = k>>1;
      52          98 :     GEN p1, best = q, np = ellmul(e,p,utoipos(n));
      53          98 :     if (n % 2 && smaller_x(gel(np,1), gel(best,1))) best = np;
      54          98 :     p1 = elladd(e,q,np);
      55          98 :     if (smaller_x(gel(p1,1), gel(best,1))) q = p1;
      56          77 :     else if (best == np) { p = elladd(e,p,q); q = np; }
      57          98 :     p = best_in_cycle(e,p,k);
      58          98 :     if (v)
      59             :     {
      60           0 :       p = ellchangepointinv(p,v);
      61           0 :       q = ellchangepointinv(q,v);
      62             :     }
      63          98 :     r = cgetg(4,t_VEC);
      64          98 :     gel(r,1) = utoipos(2*k);
      65          98 :     gel(r,2) = mkvec2(utoipos(k), gen_2);
      66          98 :     gel(r,3) = mkvec2copy(p, q);
      67             :   }
      68             :   else
      69             :   {
      70         490 :     if (p)
      71             :     {
      72         350 :       p = best_in_cycle(e,p,k);
      73         350 :       if (v) p = ellchangepointinv(p,v);
      74         350 :       r = cgetg(4,t_VEC);
      75         350 :       gel(r,1) = utoipos(k);
      76         350 :       gel(r,2) = mkvec( gel(r,1) );
      77         350 :       gel(r,3) = mkvec( gcopy(p) );
      78             :     }
      79             :     else
      80             :     {
      81         140 :       r = cgetg(4,t_VEC);
      82         140 :       gel(r,1) = gen_1;
      83         140 :       gel(r,2) = cgetg(1,t_VEC);
      84         140 :       gel(r,3) = cgetg(1,t_VEC);
      85             :     }
      86             :   }
      87         588 :   return r;
      88             : }
      89             : 
      90             : /* Finds a multiplicative upper bound for #E_tor; assume integral model */
      91             : static long
      92         588 : torsbound(GEN e)
      93             : {
      94         588 :   GEN D = ell_get_disc(e);
      95         588 :   pari_sp av = avma, av2;
      96             :   long m, b, bold, nb;
      97             :   forprime_t S;
      98         588 :   long CM = ellQ_get_CM(e);
      99         588 :   nb = expi(D) >> 3;
     100             :   /* nb = number of primes to try ~ 1 prime every 8 bits in D */
     101         588 :   b = bold = 5040; /* = 2^4 * 3^2 * 5 * 7 */
     102         588 :   m = 0;
     103             :   /* p > 2 has good reduction => E(Q) injects in E(Fp) */
     104         588 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
     105         588 :   av2 = avma;
     106        4179 :   while (m < nb || (b > 12 && b != 16))
     107             :   {
     108        3122 :     ulong p = u_forprime_next(&S);
     109        3122 :     if (!p) pari_err_BUG("torsbound [ran out of primes]");
     110        3122 :     if (!umodiu(D, p)) continue;
     111             : 
     112        2492 :     b = ugcd(b, p+1 - ellap_CM_fast(e,p,CM));
     113        2492 :     avma = av2;
     114        2492 :     if (b == 1) break;
     115        2373 :     if (b == bold) m++; else { bold = b; m = 0; }
     116             :   }
     117         588 :   avma = av; return b;
     118             : }
     119             : 
     120             : /* return a rational point of order pk = p^k on E, or NULL if E(Q)[k] = O.
     121             :  * *fk is either NULL (pk = 4 or prime) or elldivpol(p^(k-1)).
     122             :  * Set *fk to elldivpol(p^k) */
     123             : static GEN
     124         448 : tpoint(GEN E, long pk, GEN *fk)
     125             : {
     126         448 :   GEN f = elldivpol(E,pk,0), g = *fk, v;
     127             :   long i, l;
     128         448 :   *fk = f;
     129         448 :   if (g) f = RgX_div(f, g);
     130         448 :   v = nfrootsQ(f); l = lg(v);
     131         728 :   for (i = 1; i < l; i++)
     132             :   {
     133         630 :     GEN x = gel(v,i);
     134         630 :     GEN y = ellordinate(E,x,0);
     135         630 :     if (lg(y) != 1) return mkvec2(x,gel(y,1));
     136             :   }
     137          98 :   return NULL;
     138             : }
     139             : /* return E(Q)[2] */
     140             : static GEN
     141         322 : t2points(GEN E, GEN *f2)
     142             : {
     143             :   long i, l;
     144             :   GEN v;
     145         322 :   *f2 = ec_bmodel(E);
     146         322 :   v = nfrootsQ(*f2); l = lg(v);
     147         826 :   for (i = 1; i < l; i++)
     148             :   {
     149         504 :     GEN x = gel(v,i);
     150         504 :     GEN y = ellordinate(E,x,0);
     151         504 :     if (lg(y) != 1) gel(v,i) = mkvec2(x,gel(y,1));
     152             :   }
     153         322 :   return v;
     154             : }
     155             : 
     156             : static GEN
     157         588 : elltors_divpol(GEN E)
     158             : {
     159         588 :   GEN T2 = NULL, p, P, Q, v;
     160             :   long v2, r2, B;
     161             : 
     162         588 :   E = ellintegralmodel_i(E, &v);
     163         588 :   B = torsbound(E); /* #E_tor | B */
     164         588 :   if (B == 1) return tors(E,1,NULL,NULL, v);
     165         469 :   v2 = vals(B); /* bound for v_2(point order) */
     166         469 :   B >>= v2;
     167         469 :   p = const_vec(9, NULL);
     168         469 :   r2 = 0;
     169         469 :   if (v2) {
     170             :     GEN f;
     171         322 :     T2 = t2points(E, &f);
     172         322 :     switch(lg(T2)-1)
     173             :     {
     174          14 :       case 0:  v2 = 0; break;
     175         210 :       case 1:  r2 = 1; if (v2 == 4) v2 = 3; break;
     176          98 :       default: r2 = 2; v2--; break; /* 3 */
     177             :     }
     178         322 :     if (v2) gel(p,2) = gel(T2,1);
     179             :     /* f = f_2 */
     180         322 :     if (v2 > 1) { gel(p,4) = tpoint(E,4, &f); if (!gel(p,4)) v2 = 1; }
     181             :     /* if (v2>1) now f = f4 */
     182         322 :     if (v2 > 2) { gel(p,8) = tpoint(E,8, &f); if (!gel(p,8)) v2 = 2; }
     183             :   }
     184         469 :   B <<= v2;
     185         469 :   if (B % 3 == 0) {
     186         112 :     GEN f3 = NULL;
     187         112 :     gel(p,3) = tpoint(E,3,&f3);
     188         112 :     if (!gel(p,3)) B /= (B%9)? 3: 9;
     189         112 :     if (gel(p,3) && B % 9 == 0)
     190             :     {
     191           7 :       gel(p,9) = tpoint(E,9,&f3);
     192           7 :       if (!gel(p,9)) B /= 3;
     193             :     }
     194             :   }
     195         469 :   if (B % 5 == 0) {
     196          91 :     GEN junk = NULL;
     197          91 :     gel(p,5) = tpoint(E,5,&junk);
     198          91 :     if (!gel(p,5)) B /= 5;
     199             :   }
     200         469 :   if (B % 7 == 0) {
     201          21 :     GEN junk = NULL;
     202          21 :     gel(p,7) = tpoint(E,7,&junk);
     203          21 :     if (!gel(p,7)) B /= 7;
     204             :   }
     205             :   /* B is the exponent of E_tors(Q), r2 is the rank of its 2-Sylow,
     206             :    * for i > 1, p[i] is a point of order i if one exists and i is a prime power
     207             :    * and NULL otherwise */
     208         469 :   if (r2 == 2) /* 2 cyclic factors */
     209             :   { /* C2 x C2 */
     210          98 :     if (B == 2) return tors(E,2, gel(T2,1), gel(T2,2), v);
     211          42 :     else if (B == 6)
     212             :     { /* C2 x C6 */
     213          14 :       P = elladd(E, gel(p,3), gel(T2,1));
     214          14 :       Q = gel(T2,2);
     215             :     }
     216             :     else
     217             :     { /* C2 x C4 or C2 x C8 */
     218          28 :       P = gel(p, B);
     219          28 :       Q = gel(T2,2);
     220          28 :       if (gequal(Q, ellmul(E, P, utoipos(B>>1)))) Q = gel(T2,1);
     221             :     }
     222             :   }
     223             :   else /* cyclic */
     224             :   {
     225         371 :     Q = NULL;
     226         371 :     if (v2)
     227             :     {
     228         210 :       if (B>>v2 == 1)
     229         168 :         P = gel(p, B);
     230             :       else
     231          42 :         P = elladd(E, gel(p, B>>v2), gel(p,1<<v2));
     232             :     }
     233         161 :     else P = gel(p, B);
     234             :   }
     235         413 :   return tors(E,B, P, Q, v);
     236             : }
     237             : 
     238             : /* either return one prime of degree 1 above p or NULL (none or expensive) */
     239             : static GEN
     240        6125 : primedec_deg1(GEN K, GEN p)
     241             : {
     242        6125 :   GEN r, T, f = nf_get_index(K);
     243        6125 :   if (dvdii(f,p)) return NULL;
     244        6048 :   T = nf_get_pol(K);
     245        6048 :   r = FpX_oneroot(T, p); if (!r) return NULL;
     246        1106 :   r = deg1pol_shallow(gen_1, Fp_neg(r,p), varn(T));
     247        1106 :   return idealprimedec_kummer(K, r, 1, p);
     248             : }
     249             : 
     250             : /* Bound for the elementary divisors of the torsion group of elliptic curve
     251             :  * Reduce the curve modulo some small good primes */
     252             : static GEN
     253         322 : nftorsbound(GEN E)
     254             : {
     255             :   pari_sp av;
     256         322 :   long k = 0, g;
     257         322 :   GEN B1 = gen_0, B2 = gen_0, K = ellnf_get_nf(E);
     258         322 :   GEN D = ell_get_disc(E), ND = idealnorm(K,D);
     259             :   forprime_t S;
     260         322 :   if (typ(ND) == t_FRAC) ND = gel(ND,1);
     261         322 :   ND = mulii(ND, Q_denom(vecslice(E,1,5)));
     262         322 :   g = maxss(5, expi(ND) >> 3);
     263         322 :   if (g > 20) g = 20;
     264             :   /* P | p such that e(P/p) < p-1 => E(K) injects in E(k(P)) [otherwise
     265             :    * we may lose some p-torsion]*/
     266         322 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
     267         322 :   av = avma;
     268        9163 :   while (k < g) /* k = number of good primes already used */
     269             :   {
     270        8603 :     ulong p = u_forprime_next(&S);
     271             :     GEN P, gp;
     272             :     long j, l;
     273        8603 :     if (!umodiu(ND,p)) continue;
     274        7826 :     gp = utoipos(p);
     275             :     /* primes of degree 1 are easier and give smaller bounds */
     276        7826 :     if (typ(D) != t_POLMOD) /* E/Q */
     277             :     {
     278        5796 :       P = primedec_deg1(K, gp); /* single P|p has all the information */
     279        5796 :       if (!P) continue;
     280         973 :       P = mkvec(P);
     281             :     }
     282             :     else
     283        2030 :       P = idealprimedec_limit_f(K, utoipos(p), 1);
     284        3003 :     l = lg(P);
     285        5376 :     for (j = 1; j < l; j++,k++)
     286             :     {
     287        2457 :       GEN Q = gel(P,j), EQ, cyc;
     288             :       long n;
     289        2457 :       if ((ulong)pr_get_e(Q) >= p-1) continue;
     290        2457 :       EQ = ellinit(E,zkmodprinit(K,Q),0);
     291        2457 :       cyc = ellgroup(EQ, NULL);
     292        2457 :       n = lg(cyc)-1;
     293        2457 :       if (n == 0) return mkvec2(gen_1,gen_1);
     294        2450 :       B1 = gcdii(B1,gel(cyc,1));
     295        2450 :       B2 = (n == 1)? gen_1: gcdii(B2,gel(cyc,2));
     296        2450 :       obj_free(EQ);
     297             :       /* division by 2 is cheap when it fails, no need to have a sharp bound */
     298        2450 :       if (Z_ispow2(B1)) return mkvec2(B1,B2);
     299             :     }
     300        2919 :     if ((g & 15) == 0) gerepileall(av, 2, &B1, &B2);
     301             :   }
     302         238 :   if (abscmpiu(B2, 2) > 0)
     303             :   { /* if E(K) has full n-torsion then K contains the n-th roots of 1 */
     304          91 :     GEN n = gel(rootsof1(K), 1);
     305          91 :     B2 = gcdii(B2,n);
     306             :   }
     307         238 :   return mkvec2(B1,B2);
     308             : }
     309             : 
     310             : /* Checks whether the point P is divisible by n in E(K), where xn is
     311             :  * [phi_n, psi_n^2]
     312             :  * If true, returns a point Q such that nQ = P or -P. Else, returns NULL */
     313             : static GEN
     314         189 : ellnfis_divisible_by(GEN E, GEN K, GEN P, GEN xn)
     315             : {
     316         189 :   GEN r, x = gel(P,1);
     317             :   long i, l;
     318         189 :   if (ell_is_inf(P)) return P;
     319         189 :   r = nfroots(K, RgX_sub(RgX_Rg_mul(gel(xn,2), x), gel(xn,1)));
     320         189 :   l = lg(r);
     321         189 :   for(i=1; i<l; i++)
     322             :   {
     323         112 :     GEN a = gel(r,i), y = ellordinate(E,a,0);
     324         112 :     if (lg(y) != 1) return mkvec2(a, gel(y,1));
     325             :   }
     326          77 :   return NULL;
     327             : }
     328             : 
     329             : long
     330         140 : ellisdivisible(GEN E, GEN P, GEN n, GEN *pQ)
     331             : {
     332         140 :   pari_sp av = avma;
     333         140 :   GEN xP, R, K = NULL, N = NULL;
     334             :   long i, l;
     335         140 :   checkell(E);
     336         140 :   switch(ell_get_type(E))
     337             :   {
     338         140 :     case t_ELL_Q: break;
     339           0 :     case t_ELL_NF: K = ellnf_get_nf(E); break;
     340           0 :     default: pari_err_TYPE("ellisdivisible",E);
     341             :   }
     342         140 :   checkellpt(P);
     343         140 :   switch(typ(n))
     344             :   {
     345             :     case t_INT:
     346          63 :       N = n;
     347          63 :       if (!isprime(absi(n)))
     348             :       {
     349          42 :         GEN f = absZ_factor(n), LP = gel(f,1), LE = gel(f,2);
     350          42 :         l = lg(LP);
     351          84 :         for (i = 1; i < l; i++)
     352             :         {
     353          56 :           long j, e = itos(gel(LE,i));
     354          56 :           GEN xp = ellxn(E,itos(gel(LP,i)),0);
     355         112 :           for (j = 1; j <= e; j++)
     356          70 :             if (!ellisdivisible(E, P, xp, &P)) { avma = av; return 0; }
     357             :         }
     358          28 :         if (pQ)
     359             :         {
     360          28 :           if (signe(n) < 0) P = ellneg(E, P);
     361          28 :           *pQ = gerepilecopy(av,P);
     362             :         }
     363           0 :         else av = avma;
     364          28 :         return 1;
     365             :       }
     366          21 :       n = ellxn(E, itou(n), 0);
     367          21 :       break;
     368             :     case t_VEC:
     369          77 :       if (lg(n) == 3 && typ(gel(n,1)) == t_POL && typ(gel(n,2)) == t_POL) break;
     370             :     default:
     371           0 :       pari_err_TYPE("ellisdivisible",n);
     372           0 :       break;
     373             :   }
     374          98 :   if (ell_is_inf(P)) { if (pQ) *pQ = ellinf(); return 1; }
     375          91 :   if (!N) {
     376          70 :     long d, d2 = degpol(gel(n,1));
     377          70 :     if (d2 < 0)
     378           7 :       N = gen_0;
     379             :     else
     380             :     {
     381          63 :       if (!uissquareall(d2,(ulong*)&d)) pari_err_TYPE("ellisdivisible",n);
     382          63 :       N = stoi(d);
     383             :     }
     384             :   }
     385          91 :   if (!signe(N)) return 0;
     386          84 :   xP = gel(P,1);
     387          84 :   R = nfroots(K, RgX_sub(RgX_Rg_mul(gel(n,2), xP), gel(n,1)));
     388          84 :   l = lg(R);
     389         168 :   for(i=1; i<l; i++)
     390             :   {
     391          70 :     GEN Q,y, x = gel(R,i), a = ellordinate(E,x,0);
     392          70 :     if (lg(a) == 1) continue;
     393          70 :     y = gel(a,1);
     394          70 :     Q = mkvec2(x,y);
     395          70 :     if (!gequal(P,ellmul(E,Q,N))) Q = ellneg(E,Q); /* nQ = -P */
     396          70 :     if (pQ) *pQ = gerepilecopy(av,Q); else avma = av;
     397          70 :     return 1;
     398             :   }
     399          14 :   avma = av; return 0;
     400             : }
     401             : 
     402             : /* 2-torsion point of abscissa x */
     403             : static GEN
     404         126 : tor2(GEN E, GEN x) { return mkvec2(x, gmul2n(gneg(ec_h_evalx(E,x)), -1)); }
     405             : 
     406             : static GEN
     407          35 : ptor0(void)
     408          35 : { return mkvec2(mkvec(gen_1),cgetg(1,t_VEC)); }
     409             : static GEN
     410         126 : ptor1(long p, long n, GEN P)
     411         126 : { return mkvec2(mkvec(powuu(p,n)), mkvec(P)); }
     412             : static GEN
     413          98 : ptor2(long p, long n1, long n2, GEN P1, GEN P2)
     414          98 : { return mkvec2(mkvec2(powuu(p,n1), powuu(p,n2)), mkvec2(P1,P2)); }
     415             : 
     416             : /* Computes the p-primary torsion in E(K). Assume that p is small, should use
     417             :  * Weil pairing otherwise.
     418             :  * N1, N2 = upper bounds on the integers n1 >= n2 such that
     419             :  * E(K)[p^oo] = Z/p^n1 x Z/p^n2
     420             :  * Returns [cyc,gen], where E(K)[p^oo] = sum Z/cyc[i] gen[i] */
     421             : static GEN
     422         259 : ellnftorsprimary(GEN E, long p, long N1, long N2, long v)
     423             : {
     424         259 :   GEN X, P1, P2, Q1, Q2, xp, K = ellnf_get_nf(E);
     425             :   long n1, n2;
     426             : 
     427             :   /* compute E[p] = < P1 > or < P1, P2 > */
     428         259 :   P1 = P2 = ellinf();
     429         259 :   X = nfroots(K, elldivpol(E,p,v));
     430         259 :   if(lg(X) == 1) return ptor0();
     431         231 :   if (p==2)
     432             :   {
     433          77 :     P1 = tor2(E, gel(X,1));
     434          77 :     if (lg(X) > 2) P2 = tor2(E, gel(X,2)); /* E[2] = (Z/2Z)^2 */
     435             :   }
     436             :   else
     437             :   {
     438         154 :     long j, l = lg(X), nT, a;
     439         154 :     GEN T = vectrunc_init(l);
     440         539 :     for(j=1; j < l; j++)
     441             :     {
     442         385 :       GEN a = gel(X,j), Y = ellordinate(E,a,0);
     443         385 :       if (lg(Y) != 1) vectrunc_append(T, mkvec2(a,gel(Y,1)));
     444             :     }
     445         154 :     nT = lg(T)-1;
     446         154 :     if (!nT) return ptor0();
     447         147 :     P1 = gel(T,1);
     448         147 :     a = (p-1)/2;
     449         147 :     if (nT != a)
     450             :     { /* E[p] = (Z/pZ)^2 */
     451          49 :       GEN Z = cgetg(a+1,t_VEC), Q1 = P1;
     452             :       long k;
     453          49 :       gel(Z,1) = Q1;
     454          49 :       for (k=2; k <= a; k++) gel(Z,k) = elladd(E,Q1,P1);
     455          49 :       gen_sort_inplace(Z, (void*)&cmp_universal, &cmp_nodata, NULL);
     456          49 :       while (tablesearch(Z, gel(T,k), &cmp_universal)) k++;
     457          49 :       P2 = gel(T,k);
     458             :     }
     459             :   }
     460         224 :   xp = ellxn(E, p, v);
     461             : 
     462         224 :   if (ell_is_inf(P2))
     463             :   { /* E[p^oo] is cyclic, start from P1 and divide by p while possible */
     464         140 :     for (n1 = 1; n1 < N1; n1++)
     465             :     {
     466          14 :       GEN Q = ellnfis_divisible_by(E,K,P1,xp);
     467          14 :       if (!Q) break;
     468          14 :       P1 = Q;
     469             :     }
     470         126 :     return ptor1(p, n1, P1);
     471             :   }
     472             : 
     473             :   /* E[p] = (Z/pZ)^2, compute n2 and E[p^n2] */
     474          98 :   Q1 = NULL;
     475         119 :   for (n2 = 1; n2 < N2; n2++)
     476             :   {
     477          21 :     Q1 = ellnfis_divisible_by(E,K,P1,xp);
     478          21 :     Q2 = ellnfis_divisible_by(E,K,P2,xp);
     479          21 :     if (!Q1 || !Q2) break;
     480          21 :     P1 = Q1;
     481          21 :     P2 = Q2;
     482             :   }
     483             : 
     484             :   /* compute E[p^\infty] = < P1, P2 > */
     485          98 :   n1 = n2;
     486          98 :   if (n2 == N2)
     487             :   {
     488          98 :     if (N1 == N2) return ptor2(p, n2,n2, P1,P2);
     489          42 :     Q1 = ellnfis_divisible_by(E,K,P1,xp);
     490             :   }
     491          42 :   if (Q1) { P1 = Q1; n1++; }
     492             :   else
     493             :   {
     494          28 :     Q2 = ellnfis_divisible_by(E,K,P2,xp);
     495          28 :     if (Q2) { P2 = P1; P1 = Q2; n1++; }
     496             :     else
     497             :     {
     498             :       long k;
     499          35 :       for (k = 1; k < p; k++)
     500             :       {
     501          28 :         P1 = elladd(E,P1,P2);
     502          28 :         Q1 = ellnfis_divisible_by(E,K,P1,xp);
     503          28 :         if (Q1) { P1 = Q1; n1++; break; }
     504             :       }
     505          28 :       if (k == p) return ptor2(p, n2,n2, P1,P2);
     506             :     }
     507             :   }
     508             :   /* P1,P2 of order p^n1,p^n2 with n1=n2+1.
     509             :    * Keep trying to divide P1 + k P2 with 0 <= k < p by p */
     510          91 :   while (n1 < N1)
     511             :   {
     512          21 :     Q1 = ellnfis_divisible_by(E,K,P1,xp);
     513          21 :     if (Q1) { P1 = Q1; n1++; }
     514             :     else
     515             :     {
     516             :       long k;
     517          14 :       for (k = 1; k < p; k++)
     518             :       {
     519          14 :         P1 = elladd(E,P1,P2);
     520          14 :         Q1 = ellnfis_divisible_by(E,K,P1,xp);
     521          14 :         if (Q1) { P1 = Q1; n1++; break; }
     522             :       }
     523          14 :       if (k == p) break;
     524             :     }
     525             :   }
     526          35 :   return ptor2(p, n1,n2, P1,P2);
     527             : }
     528             : 
     529             : /* P affine point */
     530             : static GEN
     531         259 : nfpt(GEN e, GEN P)
     532             : {
     533         259 :   GEN T = nf_get_pol(ellnf_get_nf(e));
     534         259 :   GEN x = gel(P,1), y = gel(P,2);
     535         259 :   long tx = typ(x), ty = typ(y);
     536         259 :   if (tx == ty) return P;
     537          98 :   if (tx != t_POLMOD) x = mkpolmod(x,T); else y = mkpolmod(y,T);
     538          98 :   return mkvec2(x,y);
     539             : }
     540             : /* Computes the torsion subgroup of E(K), as [order, cyc, gen] */
     541             : static GEN
     542         189 : ellnftors(GEN e)
     543             : {
     544         189 :   GEN B = nftorsbound(e), B1 = gel(B,1), B2 = gel(B,2), d1,d2, P1,P2;
     545         189 :   GEN f = Z_factor(B1), P = gel(f,1), E = gel(f,2);
     546         189 :   long i, l = lg(P), v = fetch_var_higher();
     547             : 
     548         189 :   d1 = d2 = gen_1; P1 = P2 = ellinf();
     549         448 :   for (i=1; i<l; i++)
     550             :   {
     551         259 :     long p = itos(gel(P,i)); /* Compute p-primary torsion */
     552         259 :     long N1 = itos(gel(E,i)); /* >= n1 */
     553         259 :     long N2 = Z_lval(B2,p); /* >= n2 */
     554         259 :     GEN T = ellnftorsprimary(e, p, N1, N2, v), cyc = gel(T,1), gen = gel(T,2);
     555         259 :     if (is_pm1(gel(cyc,1))) continue;
     556             :     /* update generators P1,P2 and their respective orders d1,d2 */
     557         224 :     P1 = elladd(e, P1, gel(gen,1)); d1 = mulii(d1, gel(cyc,1));
     558         224 :     if (lg(cyc) > 2)
     559          98 :     { P2 = elladd(e, P2, gel(gen,2)); d2 = mulii(d2, gel(cyc,2)); }
     560             :   }
     561         189 :   (void)delete_var();
     562         189 :   if (is_pm1(d1)) return mkvec3(gen_1,cgetg(1,t_VEC),cgetg(1,t_VEC));
     563         168 :   if (is_pm1(d2)) return mkvec3(d1, mkvec(d1), mkvec(nfpt(e,P1)));
     564          91 :   return mkvec3(mulii(d1,d2), mkvec2(d1,d2), mkvec2(nfpt(e,P1),nfpt(e,P2)));
     565             : }
     566             : 
     567             : GEN
     568         777 : elltors(GEN e)
     569             : {
     570         777 :   pari_sp av = avma;
     571         777 :   GEN t = NULL;
     572         777 :   checkell(e);
     573         777 :   switch(ell_get_type(e))
     574             :   {
     575         588 :     case t_ELL_Q:  t = elltors_divpol(e); break;
     576         189 :     case t_ELL_NF: t = ellnftors(e); break;
     577             :     case t_ELL_Fp:
     578           0 :     case t_ELL_Fq: return ellgroup0(e,NULL,1);
     579           0 :     default: pari_err_TYPE("elltors",e);
     580             :   }
     581         777 :   return gerepilecopy(av, t);
     582             : }
     583             : 
     584             : GEN
     585           0 : elltors0(GEN e, long flag) { (void)flag; return elltors(e); }
     586             : 
     587             : /********************************************************************/
     588             : /**                                                                **/
     589             : /**                ORDER OF POINTS over NUMBER FIELDS              **/
     590             : /**                                                                **/
     591             : /********************************************************************/
     592             : /* E a t_ELL_Q (use Mazur's theorem) */
     593             : long
     594         301 : ellorder_Q(GEN E, GEN P)
     595             : {
     596         301 :   pari_sp av = avma;
     597             :   GEN dx, dy, d4, d6, D, Pp, Q;
     598             :   forprime_t S;
     599             :   ulong a4, p;
     600             :   long k;
     601         301 :   if (ell_is_inf(P)) return 1;
     602         301 :   if (gequal(P, ellneg(E,P))) return 2;
     603             : 
     604         280 :   dx = Q_denom(gel(P,1));
     605         280 :   dy = Q_denom(gel(P,2));
     606         280 :   if (ell_is_integral(E)) /* integral model, try Nagell Lutz */
     607         231 :     if (abscmpiu(dx, 4) > 0 || abscmpiu(dy, 8) > 0) return 0;
     608             : 
     609         238 :   d4 = Q_denom(ell_get_c4(E));
     610         238 :   d6 = Q_denom(ell_get_c6(E));
     611         238 :   D = ell_get_disc (E);
     612             :   /* choose not too small prime p dividing neither a coefficient of the
     613             :      short Weierstrass form nor of P and leading to good reduction */
     614         238 :   u_forprime_init(&S, 100003, ULONG_MAX);
     615         238 :   while ( (p = u_forprime_next(&S)) )
     616         238 :     if (umodiu(d4, p) && umodiu(d6, p) && Rg_to_Fl(D, p)
     617         238 :      && umodiu(dx, p) && umodiu(dy, p)) break;
     618             : 
     619             :   /* transform E into short Weierstrass form Ep modulo p and P to Pp on Ep,
     620             :    * check whether the order of Pp on Ep is <= 12 */
     621         238 :   Pp = point_to_a4a6_Fl(E, P, p, &a4);
     622        2646 :   for (Q = Fle_dbl(Pp, a4, p), k = 2;
     623        4753 :        !ell_is_inf(Q) && k <= 12;
     624        2170 :        Q = Fle_add(Q, Pp, a4, p), k++) /* empty */;
     625             : 
     626         238 :   if (k == 13) k = 0;
     627             :   else
     628             :   { /* check whether [k]P = O over Q. Save potentially costly last elladd */
     629             :     GEN R;
     630          63 :     Q = ellmul(E, P, utoipos(k>>1));
     631          63 :     R = odd(k)? elladd(E, P,Q): Q;
     632          63 :     if (!gequal(Q, ellneg(E,R))) k = 0;
     633             :   }
     634         238 :   avma = av; return k;
     635             : }
     636             : /* E a t_ELL_NF */
     637             : static GEN
     638         161 : ellorder_nf(GEN E, GEN P)
     639             : {
     640         161 :   GEN K = ellnf_get_nf(E), B;
     641         161 :   pari_sp av = avma;
     642             :   GEN dx, dy, d4, d6, D, ND, Ep, Pp, Q, gp, modpr, pr, T, k;
     643             :   forprime_t S;
     644             :   ulong a4, p;
     645         161 :   if (ell_is_inf(P)) return gen_1;
     646         161 :   if (gequal(P, ellneg(E,P))) return gen_2;
     647             : 
     648         133 :   B = gel(nftorsbound(E), 1);
     649         133 :   dx = Q_denom(gel(P,1));
     650         133 :   dy = Q_denom(gel(P,2));
     651         133 :   d4 = Q_denom(ell_get_c4(E));
     652         133 :   d6 = Q_denom(ell_get_c6(E));
     653         133 :   D = ell_get_disc(E);
     654         133 :   ND = idealnorm(K,D);
     655         133 :   if (typ(ND) == t_FRAC) ND = gel(ND,1);
     656             : 
     657             :   /* choose not too small prime p of degree 1 dividing neither a coefficient of
     658             :    * the short Weierstrass form nor of P and leading to good reduction */
     659         133 :   u_forprime_init(&S, 100003, ULONG_MAX);
     660         133 :   while ( (p = u_forprime_next(&S)) )
     661             :   {
     662         329 :     if (!umodiu(d4, p) || !umodiu(d6, p) || !umodiu(ND, p)
     663         329 :      || !umodiu(dx, p) || !umodiu(dy, p)) continue;
     664         329 :     gp = utoipos(p);
     665         329 :     pr = primedec_deg1(K, gp);
     666         329 :     if (pr) break;
     667             :   }
     668             : 
     669         133 :   modpr = nf_to_Fq_init(K, &pr,&T,&gp);
     670         133 :   Ep = ellinit(E, pr, 0);
     671         133 :   Pp = nfV_to_FqV(P, K, modpr);
     672             : 
     673             :   /* transform E into short Weierstrass form Ep modulo p and P to Pp on Ep,
     674             :    * check whether the order of Pp on Ep divides B */
     675         133 :   Pp = point_to_a4a6_Fl(Ep, Pp, p, &a4);
     676         133 :   if (!ell_is_inf(Fle_mul(Pp, B, a4, p))) { avma = av; return gen_0; }
     677         133 :   k = Fle_order(Pp, B, a4, p);
     678             :   { /* check whether [k]P = O over K. Save potentially costly last elladd */
     679             :     GEN R;
     680         133 :     Q = ellmul(E, P, shifti(k,-1));
     681         133 :     R = mod2(k)? elladd(E, P,Q): Q;
     682         133 :     if (!gequal(Q, ellneg(E,R))) k = gen_0;
     683             :   }
     684         133 :   return gerepileuptoint(av, k);
     685             : }
     686             : 
     687             : GEN
     688        2765 : ellorder(GEN E, GEN P, GEN o)
     689             : {
     690        2765 :   pari_sp av = avma;
     691        2765 :   GEN fg, r, E0 = E;
     692        2765 :   checkell(E); checkellpt(P);
     693        2765 :   if (ell_is_inf(P)) return gen_1;
     694        2751 :   if (ell_get_type(E)==t_ELL_Q)
     695             :   {
     696        1869 :     GEN p = NULL;
     697        1869 :     if (is_rational_t(typ(gel(P,1))) && is_rational_t(typ(gel(P,2))))
     698          91 :       return utoi( ellorder_Q(E, P) );
     699        1778 :     if (RgV_is_FpV(P,&p) && p)
     700             :     {
     701        1778 :       E = ellinit(E,p,0);
     702        1778 :       if (lg(E)==1) pari_err_IMPL("ellorder for curve with singular reduction");
     703             :     }
     704             :   }
     705        2653 :   if (ell_get_type(E)==t_ELL_NF) return ellorder_nf(E, P);
     706        2492 :   checkell_Fq(E);
     707        2492 :   fg = ellff_get_field(E);
     708        2492 :   if (!o) o = ellff_get_o(E);
     709        2492 :   if (typ(fg)==t_FFELT)
     710        1750 :     r = FF_ellorder(E, P, o);
     711             :   else
     712             :   {
     713         742 :     GEN p = fg, e = ellff_get_a4a6(E);
     714         742 :     GEN Pp = FpE_changepointinv(RgE_to_FpE(P,p), gel(e,3), p);
     715         742 :     r = FpE_order(Pp, o, gel(e,1), p);
     716             :   }
     717        2492 :   if (E != E0) obj_free(E);
     718        2492 :   return gerepileuptoint(av, r);
     719             : }
     720             : 
     721             : GEN
     722           0 : orderell(GEN e, GEN z) { return ellorder(e,z,NULL); }

Generated by: LCOV version 1.11