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

Generated by: LCOV version 1.14