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 - ellanal.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 726 788 92.1 %
Date: 2022-07-03 07:33:15 Functions: 67 71 94.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2010  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             : /**                  L functions of elliptic curves                **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : #define DEBUGLEVEL DEBUGLEVEL_ellanal
      24             : 
      25             : struct baby_giant
      26             : {
      27             :   GEN baby, giant, sum;
      28             :   GEN bnd, rbnd;
      29             : };
      30             : 
      31             : /* Generic Buhler-Gross algorithm */
      32             : 
      33             : struct bg_data
      34             : {
      35             :   GEN E, N; /* ell, conductor */
      36             :   GEN bnd; /* t_INT; will need all an for n <= bnd */
      37             :   ulong rootbnd; /* sqrt(bnd) */
      38             :   GEN an; /* t_VECSMALL: cache of an, n <= rootbnd */
      39             :   GEN p;  /* t_VECSMALL: primes <= rootbnd */
      40             : };
      41             : 
      42             : typedef void bg_fun(void*el, GEN n, GEN a);
      43             : 
      44             : /* a = a_n, where p = bg->pp[i] divides n, and lasta = a_{n/p}.
      45             :  * Call fun(E, N, a_N), for all N, n | N, P^+(N) <= p, a_N != 0,
      46             :  * i.e. assumes that fun accumulates a_N * w(N) */
      47             : 
      48             : static void
      49     1528191 : gen_BG_add(void *E, bg_fun *fun, struct bg_data *bg, GEN n, long i, GEN a, GEN lasta)
      50             : {
      51     1528191 :   pari_sp av = avma;
      52             :   long j;
      53     1528191 :   ulong nn = itou_or_0(n);
      54     1528191 :   if (nn && nn <= bg->rootbnd) bg->an[nn] = itos(a);
      55             : 
      56     1528191 :   if (signe(a))
      57             :   {
      58      367465 :     fun(E, n, a);
      59      367465 :     j = 1;
      60             :   }
      61             :   else
      62     1160726 :     j = i;
      63     3051076 :   for(; j <= i; j++)
      64             :   {
      65     2420509 :     ulong p = bg->p[j];
      66     2420509 :     GEN nexta, pn = mului(p, n);
      67     2420509 :     if (cmpii(pn, bg->bnd) > 0) return;
      68     1522885 :     nexta = mulis(a, bg->an[p]);
      69     1522885 :     if (i == j && umodiu(bg->N, p)) nexta = subii(nexta, mului(p, lasta));
      70     1522885 :     gen_BG_add(E, fun, bg, pn, j, nexta, a);
      71     1522885 :     set_avma(av);
      72             :   }
      73             : }
      74             : 
      75             : static void
      76          77 : gen_BG_init(struct bg_data *bg, GEN E, GEN N, GEN bnd)
      77             : {
      78          77 :   bg->E = E;
      79          77 :   bg->N = N;
      80          77 :   bg->bnd = bnd;
      81          77 :   bg->rootbnd = itou(sqrtint(bnd));
      82          77 :   bg->p = primes_upto_zv(bg->rootbnd);
      83          77 :   bg->an = ellanQ_zv(E, bg->rootbnd);
      84          77 : }
      85             : 
      86             : static void
      87          77 : gen_BG_rec(void *E, bg_fun *fun, struct bg_data *bg)
      88             : {
      89          77 :   long i, j, lp = lg(bg->p)-1;
      90          77 :   GEN bndov2 = shifti(bg->bnd, -1);
      91          77 :   pari_sp av = avma, av2;
      92             :   GEN p;
      93             :   forprime_t S;
      94          77 :   (void)forprime_init(&S, utoipos(bg->p[lp]+1), bg->bnd);
      95          77 :   av2 = avma;
      96          77 :   if (DEBUGLEVEL)
      97           0 :     err_printf("1st stage, using recursion for p <= %ld\n", bg->p[lp]);
      98        5383 :   for (i = 1; i <= lp; i++)
      99             :   {
     100        5306 :     ulong pp = bg->p[i];
     101        5306 :     long ap = bg->an[pp];
     102        5306 :     gen_BG_add(E, fun, bg, utoipos(pp), i, stoi(ap), gen_1);
     103        5306 :     set_avma(av2);
     104             :   }
     105          77 :   if (DEBUGLEVEL) err_printf("2nd stage, looping for p <= %Ps\n", bndov2);
     106     1113672 :   while ( (p = forprime_next(&S)) )
     107             :   {
     108             :     long jmax;
     109     1113672 :     GEN ap = ellap(bg->E, p);
     110     1113672 :     pari_sp av3 = avma;
     111     1113672 :     if (!signe(ap)) continue;
     112             : 
     113      556178 :     jmax = itou( divii(bg->bnd, p) ); /* 2 <= jmax <= el->rootbound */
     114      556178 :     fun(E, p, ap);
     115     8337322 :     for (j = 2;  j <= jmax; j++)
     116             :     {
     117     7781144 :       long aj = bg->an[j];
     118             :       GEN a, n;
     119     7781144 :       if (!aj) continue;
     120     1136044 :       a = mulis(ap, aj);
     121     1136044 :       n = muliu(p, j);
     122     1136044 :       fun(E, n, a);
     123     1136044 :       set_avma(av3);
     124             :     }
     125      556178 :     set_avma(av2);
     126      556178 :     if (abscmpii(p, bndov2) >= 0) break;
     127             :   }
     128          77 :   if (DEBUGLEVEL) err_printf("3nd stage, looping for p <= %Ps\n", bg->bnd);
     129      997157 :   while ( (p = forprime_next(&S)) )
     130             :   {
     131      997080 :     GEN ap = ellap(bg->E, p);
     132      997080 :     if (!signe(ap)) continue;
     133      497896 :     fun(E, p, ap);
     134      497896 :     set_avma(av2);
     135             :   }
     136          77 :   set_avma(av);
     137          77 : }
     138             : 
     139             : /******************************************************************
     140             :  *
     141             :  * L functions of elliptic curves
     142             :  * Pascal Molin (molin.maths@gmail.com) 2014
     143             :  *
     144             :  ******************************************************************/
     145             : 
     146             : struct lcritical
     147             : {
     148             :   GEN h;    /* real */
     149             :   long cprec; /* computation prec */
     150             :   long L; /* number of points */
     151             :   GEN  K; /* length of series */
     152             :   long real;
     153             : };
     154             : 
     155             : static double
     156         245 : logboundG0(long e, double aY)
     157             : {
     158             :   double cla, loggam;
     159         245 :   cla = 1 + 1/sqrt(aY);
     160         245 :   if (e) cla = ( cla + 1/(2*aY) ) / (2*sqrt(aY));
     161         245 :   loggam = (e) ? M_LN2-aY : -aY + log( log( 1+1/aY) );
     162         245 :   return log(cla) + loggam;
     163             : }
     164             : 
     165             : static void
     166         245 : param_points(GEN N, double Y, double tmax, long bprec, long *cprec, long *L,
     167             :              GEN *K, double *h)
     168             : {
     169             :   double D, a, aY, X, logM;
     170         245 :   long d = 2, w = 1;
     171         245 :   tmax *= d;
     172         245 :   D = bprec * M_LN2 + M_PI/4*tmax + 2;
     173         245 :   *cprec = nbits2prec(ceil(D / M_LN2) + 5);
     174         245 :   a = 2 * M_PI / sqrt(gtodouble(N));
     175         245 :   aY = a * cos(M_PI/2*Y);
     176         245 :   logM = 2*M_LN2 + logboundG0(w+1, aY) + tmax * Y * M_PI/2;
     177         245 :   *h = ( 2 * M_PI * M_PI / 2 * Y ) / ( D + logM );
     178         245 :   X = log( D / a);
     179         245 :   *L = ceil( X / *h);
     180         245 :   *K = ceil_safe(dbltor( D / a ));
     181         245 : }
     182             : 
     183             : static GEN
     184         210 : vecF2_lk(GEN E, GEN K, GEN rbnd, GEN Q, GEN sleh, long prec)
     185             : {
     186             :   pari_sp av;
     187         210 :   long l, L  = lg(K)-1;
     188         210 :   GEN a = ellanQ_zv(E, itos(gel(K,1)));
     189         210 :   GEN S = cgetg(L+1, t_VEC);
     190             : 
     191       12097 :   for (l = 1; l <= L; l++) gel(S,l) = cgetr(prec);
     192         210 :   av = avma;
     193       12097 :   for (l = 1; l <= L; l++)
     194             :   {
     195             :     GEN e1, Sl, z, zB;
     196       11887 :     long aB, b, A, B, Kl = itou(gel(K,l));
     197             :     pari_sp av2;
     198             :     /* FIXME: could reduce prec here (useful for large prec) */
     199       11887 :     e1 = gel(Q, l);
     200       11887 :     Sl = real_0(prec);;
     201             :     /* baby-step giant step */
     202       11887 :     B = A = rbnd[l];
     203       11887 :     z = powersr(e1, B); zB = gel(z, B+1);
     204       11887 :     av2 = avma;
     205      223054 :     for (aB = A*B; aB >= 0; aB -= B)
     206             :     {
     207      211167 :       GEN s = real_0(prec); /* could change also prec here */
     208    14832771 :       for (b = B; b > 0; --b)
     209             :       {
     210    14621604 :         long k = aB+b;
     211    14621604 :         if (k <= Kl && a[k]) s = addrr(s, mulsr(a[k], gel(z, b+1)));
     212    14621604 :         if (gc_needed(av2, 1)) gerepileall(av2, 2, &s, &Sl);
     213             :       }
     214      211167 :       Sl = addrr(mulrr(Sl, zB), s);
     215             :     }
     216       11887 :     affrr(mulrr(Sl, gel(sleh,l)), gel(S, l)); /* to avoid copying all S */
     217       11887 :     set_avma(av);
     218             :   }
     219         210 :   return S;
     220             : }
     221             : 
     222             : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
     223             : static void
     224          35 : baby_init(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
     225             : {
     226          35 :   long i, j, l = lg(Q);
     227             :   GEN R, C, r0;
     228          35 :   C = cgetg(l,t_VEC);
     229        1218 :   for (i = 1; i < l; ++i)
     230        1183 :     gel(C, i) = powersr(gel(Q, i), rbnd[i]);
     231          35 :   R = cgetg(l,t_VEC);
     232          35 :   r0 = real_0(prec);
     233        1218 :   for (i = 1; i < l; ++i)
     234             :   {
     235        1183 :     gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
     236        1183 :     gmael(R, i, 1) = cgetr(prec);
     237        1183 :     affrr(gmael(C, i, 2),gmael(R, i, 1));
     238       80234 :     for (j = 2; j <= rbnd[i]; j++)
     239             :     {
     240       79051 :       gmael(R, i, j) = cgetr(prec);
     241       79051 :       affrr(r0, gmael(R, i, j));
     242             :     }
     243             :   }
     244          35 :   bb->baby = C; bb->giant = R;
     245          35 :   bb->bnd = bnd; bb->rbnd = rbnd;
     246          35 : }
     247             : 
     248             : static long
     249         245 : baby_size(GEN rbnd, long Ks, long prec)
     250             : {
     251         245 :   long i, s, m, l = lg(rbnd);
     252       13315 :   for (s = 0, i = 1; i < l; ++i)
     253       13070 :     s += rbnd[i];
     254         245 :   m = 2*s*prec + 3*l + s;
     255         245 :   if (DEBUGLEVEL)
     256           0 :     err_printf("ellL1: BG_add: %ld words, ellan: %ld words\n", m, Ks);
     257         245 :   return m;
     258             : }
     259             : 
     260             : static void
     261      454972 : ellL1_add(void *E, GEN n, GEN a)
     262             : {
     263      454972 :   pari_sp av = avma;
     264      454972 :   struct baby_giant *bb = (struct baby_giant*) E;
     265      454972 :   long j, l = lg(bb->giant);
     266     2292864 :   for (j = 1; j < l; j++)
     267     2292864 :     if (cmpii(n, gel(bb->bnd,j)) <= 0)
     268             :     {
     269     1837892 :       ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
     270     1837892 :       GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
     271     1837892 :       affrr(addrr(gel(giant, q+1), mulri(gel(baby, r+1), a)), gel(giant, q+1));
     272     1837892 :       set_avma(av);
     273      454972 :     } else break;
     274      454972 : }
     275             : 
     276             : static GEN
     277          35 : vecF2_lk_bsgs(GEN E, GEN bnd, GEN rbnd, GEN Q, GEN sleh, GEN N, long prec)
     278             : {
     279             :   struct bg_data bg;
     280             :   struct baby_giant bb;
     281          35 :   long k, L = lg(bnd)-1;
     282             :   GEN S;
     283          35 :   baby_init(&bb, Q, bnd, rbnd, prec);
     284          35 :   gen_BG_init(&bg, E, N, gel(bnd,1));
     285          35 :   gen_BG_rec((void*) &bb, ellL1_add, &bg);
     286          35 :   S = cgetg(L+1, t_VEC);
     287        1218 :   for (k = 1; k <= L; ++k)
     288             :   {
     289        1183 :     pari_sp av = avma;
     290        1183 :     long j, g = rbnd[k];
     291        1183 :     GEN giant = gmael(bb.baby, k, g+1), Sl = gmael(bb.giant, k, g);
     292       80234 :     for (j = g-1; j >=1; j--) Sl = addrr(mulrr(Sl, giant), gmael(bb.giant,k,j));
     293        1183 :     gel(S, k) = gerepileuptoleaf(av, mulrr(gel(sleh,k), Sl));
     294             :   }
     295          35 :   return S;
     296             : }
     297             : 
     298             : static long
     299       13070 : _sqrt(GEN x) { pari_sp av = avma; return gc_long(av, itou(sqrtint(x))); }
     300             : 
     301             : static GEN
     302         245 : vecF(struct lcritical *C, GEN E)
     303             : {
     304         245 :   pari_sp av = avma;
     305         245 :   long prec = C->cprec, Ks = itos_or_0(C->K), L = C->L, l;
     306         245 :   GEN N = ellQ_get_N(E), PiN;
     307         245 :   GEN e = mpexp(C->h), elh = powersr(e, L-1), Q, bnd, rbnd, vec;
     308             : 
     309         245 :   PiN = divrr(Pi2n(1,prec), sqrtr_abs(itor(N, prec)));
     310         245 :   setsigne(PiN, -1); /* - 2Pi/sqrt(N) */
     311         245 :   bnd = gpowers0(invr(e), L-1, C->K); /* bnd[i] = K exp(-(i-1)h) */
     312         245 :   rbnd = cgetg(L+1, t_VECSMALL);
     313         245 :   Q  = cgetg(L+1, t_VEC);
     314       13315 :   for (l = 1; l <= L; l++)
     315             :   {
     316       13070 :     gel(bnd,l) = ceil_safe(gel(bnd,l));
     317       13070 :     rbnd[l] = _sqrt(gel(bnd,l)) + 1;
     318       13070 :     gel(Q, l) = mpexp(mulrr(PiN, gel(elh, l)));
     319             :   }
     320         245 :   if (Ks && baby_size(rbnd, Ks, prec) > (Ks>>1))
     321         210 :     vec = vecF2_lk(E, bnd, rbnd, Q, elh, prec);
     322             :   else
     323          35 :     vec = vecF2_lk_bsgs(E, bnd, rbnd, Q, elh, N, prec);
     324         245 :   return gerepileupto(av, vec);
     325             : }
     326             : 
     327             : /* Lambda function by Fourier inversion. vec is a grid, t a scalar or t_SER */
     328             : static GEN
     329         273 : glambda(GEN t, GEN vec, GEN h, long real, long prec)
     330             : {
     331         273 :   GEN z, r, e = gexp(gmul(mkcomplex(gen_0,h), t), prec);
     332         273 :   long n = lg(vec)-1, i;
     333             : 
     334         273 :   r = real == 1? gmul2n(real_i(gel(vec, 1)), -1): gen_0;
     335         273 :   z = real == 1? e: gmul(powIs(3), e);
     336             :   /* FIXME: summing backward may be more stable */
     337       15807 :   for (i = 2; i <= n; i++)
     338             :   {
     339       15534 :     r = gadd(r, real_i(gmul(gel(vec,i), z)));
     340       15534 :     if (i < n) z = gmul(z, e);
     341             :   }
     342         273 :   return gmul(mulsr(4, h), r);
     343             : }
     344             : 
     345             : static GEN
     346         245 : Lpoints(struct lcritical *C, GEN e, double tmax, long bprec)
     347             : {
     348         245 :   double h = 0, Y = .97;
     349         245 :   GEN N = ellQ_get_N(e);
     350         245 :   param_points(N, Y, tmax, bprec, &C->cprec, &C->L, &C->K, &h);
     351         245 :   C->real = ellrootno_global(e);
     352         245 :   C->h = rtor(dbltor(h), C->cprec);
     353         245 :   return vecF(C, e);
     354             : }
     355             : 
     356             : static GEN
     357         273 : Llambda(GEN vec, struct lcritical *C, GEN t, long prec)
     358             : {
     359         273 :   GEN lambda = glambda(gprec_w(t, C->cprec), vec, C->h, C->real, C->cprec);
     360         273 :   return gprec_w(lambda, prec);
     361             : }
     362             : 
     363             : /* 2*(2*Pi)^(-s)*gamma(s)*N^(s/2); */
     364             : static GEN
     365         273 : ellgammafactor(GEN N, GEN s, long prec)
     366             : {
     367         273 :   GEN c = gpow(divrr(gsqrt(N,prec), Pi2n(1,prec)), s, prec);
     368         273 :   return gmul(gmul2n(c,1), ggamma(s, prec));
     369             : }
     370             : 
     371             : static GEN
     372         273 : ellL1_eval(GEN e, GEN vec, struct lcritical *C, GEN t, long prec)
     373             : {
     374         273 :   GEN g = ellgammafactor(ellQ_get_N(e), gaddgs(gmul(gen_I(),t), 1), prec);
     375         273 :   return gdiv(Llambda(vec, C, t, prec), g);
     376             : }
     377             : 
     378             : static GEN
     379         273 : ellL1_der(GEN e, GEN vec, struct lcritical *C, GEN t, long der, long prec)
     380             : {
     381         273 :   GEN r = polcoef_i(ellL1_eval(e, vec, C, t, prec), der, 0);
     382         273 :   r = gmul(r,powIs(C->real == 1 ? -der: 1-der));
     383         273 :   return gmul(real_i(r), mpfact(der));
     384             : }
     385             : 
     386             : GEN
     387         231 : ellL1_bitprec(GEN E, long r, long bitprec)
     388             : {
     389         231 :   pari_sp av = avma;
     390             :   struct lcritical C;
     391         231 :   long prec = nbits2prec(bitprec);
     392             :   GEN e, vec, t;
     393         231 :   if (r < 0)
     394           7 :     pari_err_DOMAIN("ellL1", "derivative order", "<", gen_0, stoi(r));
     395         224 :   e = ellanal_globalred(E, NULL);
     396         224 :   if (r == 0 && ellrootno_global(e) < 0) { set_avma(av); return gen_0; }
     397         210 :   vec = Lpoints(&C, e, 0., bitprec);
     398         210 :   t = r ? scalarser(gen_1, 0, r):  zeroser(0, 0);
     399         210 :   setvalp(t, 1);
     400         210 :   return gerepileupto(av, ellL1_der(e, vec, &C, t, r, prec));
     401             : }
     402             : 
     403             : GEN
     404           0 : ellL1(GEN E, long r, long prec) { return ellL1_bitprec(E, r, prec2nbits(prec)); }
     405             : 
     406             : GEN
     407          35 : ellanalyticrank_bitprec(GEN E, GEN eps, long bitprec)
     408             : {
     409          35 :   pari_sp av = avma, av2;
     410          35 :   long prec = nbits2prec(bitprec);
     411             :   struct lcritical C;
     412             :   pari_timer ti;
     413             :   GEN e, vec;
     414             :   long rk;
     415          35 :   if (DEBUGLEVEL) timer_start(&ti);
     416          35 :   if (!eps)
     417          35 :     eps = real2n(-bitprec/2+1, DEFAULTPREC);
     418             :   else
     419           0 :     if (typ(eps) != t_REAL) {
     420           0 :       eps = gtofp(eps, DEFAULTPREC);
     421           0 :       if (typ(eps) != t_REAL) pari_err_TYPE("ellanalyticrank", eps);
     422             :     }
     423          35 :   e = ellanal_globalred(E, NULL);
     424          35 :   vec = Lpoints(&C, e, 0., bitprec);
     425          35 :   if (DEBUGLEVEL) timer_printf(&ti, "init L");
     426          35 :   av2 = avma;
     427          35 :   for (rk = C.real>0 ? 0: 1;  ; rk += 2)
     428          28 :   {
     429             :     GEN Lrk;
     430          63 :     GEN t = rk ? scalarser(gen_1, 0, rk):  zeroser(0, 0);
     431          63 :     setvalp(t, 1);
     432          63 :     Lrk = ellL1_der(e, vec, &C, t, rk, prec);
     433          63 :     if (DEBUGLEVEL) timer_printf(&ti, "L^(%ld)=%Ps", rk, Lrk);
     434          63 :     if (abscmprr(Lrk, eps) > 0)
     435          35 :       return gerepilecopy(av, mkvec2(stoi(rk), Lrk));
     436          28 :     set_avma(av2);
     437             :   }
     438             : }
     439             : 
     440             : GEN
     441           0 : ellanalyticrank(GEN E, GEN eps, long prec)
     442             : {
     443           0 :   return ellanalyticrank_bitprec(E, eps, prec2nbits(prec));
     444             : }
     445             : 
     446             : /*        Heegner point computation
     447             : 
     448             :    This section is a C version by Bill Allombert of a GP script by
     449             :    Christophe Delaunay which was based on a GP script by John Cremona.
     450             :    Reference: Henri Cohen's book GTM 239.
     451             : */
     452             : 
     453             : static void
     454           0 : heegner_L1_bg(void*E, GEN n, GEN a)
     455             : {
     456           0 :   struct baby_giant *bb = (struct baby_giant*) E;
     457           0 :   long j, l = lg(bb->giant);
     458           0 :   for (j = 1; j < l; j++)
     459           0 :     if (cmpii(n, gel(bb->bnd,j)) <= 0)
     460             :     {
     461           0 :       ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
     462           0 :       GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
     463           0 :       gaffect(gadd(gel(giant, q+1), gdiv(gmul(gel(baby, r+1), a), n)), gel(giant, q+1));
     464             :     }
     465           0 : }
     466             : 
     467             : static void
     468     2102611 : heegner_L1(void*E, GEN n, GEN a)
     469             : {
     470     2102611 :   struct baby_giant *bb = (struct baby_giant*) E;
     471     2102611 :   long j, l = lg(bb->giant);
     472    12033980 :   for (j = 1; j < l; j++)
     473     9931369 :     if (cmpii(n, gel(bb->bnd,j)) <= 0)
     474             :     {
     475     8349474 :       ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
     476     8349474 :       GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
     477     8349474 :       GEN ex = mulreal(gel(baby, r+1), gel(giant, q+1));
     478     8349474 :       affrr(addrr(gel(bb->sum, j), divri(mulri(ex, a), n)), gel(bb->sum, j));
     479             :     }
     480     2102611 : }
     481             : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
     482             : static void
     483           0 : baby_init2(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
     484             : {
     485           0 :   long i, j, l = lg(Q);
     486             :   GEN R, C, r0;
     487           0 :   C = cgetg(l,t_VEC);
     488           0 :   for (i = 1; i < l; ++i)
     489           0 :     gel(C, i) = gpowers(gel(Q, i), rbnd[i]);
     490           0 :   R = cgetg(l,t_VEC);
     491           0 :   r0 = mkcomplex(real_0(prec),real_0(prec));
     492           0 :   for (i = 1; i < l; ++i)
     493             :   {
     494           0 :     gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
     495           0 :     gmael(R, i, 1) = cgetc(prec);
     496           0 :     gaffect(gmael(C, i, 2),gmael(R, i, 1));
     497           0 :     for (j = 2; j <= rbnd[i]; j++)
     498             :     {
     499           0 :       gmael(R, i, j) = cgetc(prec);
     500           0 :       gaffect(r0, gmael(R, i, j));
     501             :     }
     502             :   }
     503           0 :   bb->baby = C; bb->giant = R;
     504           0 :   bb->bnd = bnd; bb->rbnd = rbnd;
     505           0 : }
     506             : 
     507             : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
     508             : static void
     509          42 : baby_init3(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
     510             : {
     511          42 :   long i, l = lg(Q);
     512             :   GEN R, C, S;
     513          42 :   C = cgetg(l,t_VEC);
     514         189 :   for (i = 1; i < l; ++i)
     515         147 :     gel(C, i) = gpowers(gel(Q, i), rbnd[i]);
     516          42 :   R = cgetg(l,t_VEC);
     517         189 :   for (i = 1; i < l; ++i)
     518         147 :     gel(R, i) = gpowers(gmael(C, i, 1+rbnd[i]), rbnd[i]);
     519          42 :   S = cgetg(l,t_VEC);
     520         189 :   for (i = 1; i < l; ++i)
     521             :   {
     522         147 :     gel(S, i) = cgetr(prec);
     523         147 :     affrr(real_i(gmael(C, i, 2)), gel(S, i));
     524             :   }
     525          42 :   bb->baby = C; bb->giant = R; bb->sum = S;
     526          42 :   bb->bnd = bnd; bb->rbnd = rbnd;
     527          42 : }
     528             : 
     529             : /* ymin a t_REAL */
     530             : static GEN
     531          42 : heegner_psi(GEN E, GEN N, GEN points, long bitprec)
     532             : {
     533          42 :   pari_sp av = avma, av2;
     534             :   struct baby_giant bb;
     535             :   struct bg_data bg;
     536          42 :   long l, k, L = lg(points)-1, prec = nbits2prec(bitprec)+EXTRAPRECWORD;
     537          42 :   GEN Q, pi2 = Pi2n(1, prec), bnd, rbnd, bndmax;
     538          42 :   GEN B = divrr(mulur(bitprec,mplog2(DEFAULTPREC)), pi2);
     539             : 
     540          42 :   rbnd = cgetg(L+1, t_VECSMALL); av2 = avma;
     541          42 :   bnd = cgetg(L+1, t_VEC);
     542          42 :   Q  = cgetg(L+1, t_VEC);
     543         189 :   for (l = 1; l <= L; ++l)
     544             :   {
     545         147 :     gel(bnd,l) = ceil_safe(divrr(B,imag_i(gel(points, l))));
     546         147 :     rbnd[l] = itou(sqrtint(gel(bnd,l)))+1;
     547         147 :     gel(Q, l) = expIxy(pi2, gel(points, l), prec);
     548             :   }
     549          42 :   gerepileall(av2, 2, &bnd, &Q);
     550          42 :   bndmax = gel(bnd,vecindexmax(bnd));
     551          42 :   gen_BG_init(&bg, E, N, bndmax);
     552          42 :   if (bitprec >= 1900)
     553             :   {
     554           0 :     GEN S = cgetg(L+1, t_VEC);
     555           0 :     baby_init2(&bb, Q, bnd, rbnd, prec);
     556           0 :     gen_BG_rec((void*)&bb, heegner_L1_bg, &bg);
     557           0 :     for (k = 1; k <= L; ++k)
     558             :     {
     559           0 :       pari_sp av2 = avma;
     560           0 :       long j, g = rbnd[k];
     561           0 :       GEN giant = gmael(bb.baby, k, g+1), Sl = real_0(prec);
     562           0 :       for (j = g; j >=1; j--) Sl = gadd(gmul(Sl, giant), gmael(bb.giant,k,j));
     563           0 :       gel(S, k) = gerepileupto(av2, real_i(Sl));
     564             :     }
     565           0 :     return gerepileupto(av, S);
     566             :   }
     567             :   else
     568             :   {
     569          42 :     baby_init3(&bb, Q, bnd, rbnd, prec);
     570          42 :     gen_BG_rec((void*)&bb, heegner_L1, &bg);
     571          42 :     return gerepilecopy(av, bb.sum);
     572             :   }
     573             : }
     574             : 
     575             : /*Returns lambda_bad list for one prime p, nv = localred(E, p) */
     576             : static GEN
     577          91 : lambda1(GEN E, GEN nv, GEN p, long prec)
     578             : {
     579             :   pari_sp av;
     580             :   GEN res, lp;
     581          91 :   long kod = itos(gel(nv, 2));
     582          91 :   if (kod==2 || kod ==-2) return cgetg(1,t_VEC);
     583          91 :   av = avma; lp = glog(p, prec);
     584          91 :   if (kod > 4)
     585             :   {
     586          14 :     long n = Z_pval(ell_get_disc(E), p);
     587          14 :     long j, m = kod - 4, nl = 1 + (m >> 1L);
     588          14 :     res = cgetg(nl, t_VEC);
     589          35 :     for (j = 1; j < nl; j++)
     590          21 :       gel(res, j) = gmul(lp, gsubgs(gdivgu(sqru(j), n), j)); /* j^2/n - j */
     591             :   }
     592          77 :   else if (kod < -4)
     593          14 :     res = mkvec2(negr(lp), shiftr(mulrs(lp, kod), -2));
     594             :   else
     595             :   {
     596          63 :     const long lam[] = {8,9,0,6,0,0,0,3,4};
     597          63 :     long m = -lam[kod+4];
     598          63 :     res = mkvec(divru(mulrs(lp, m), 6));
     599             :   }
     600          91 :   return gerepilecopy(av, res);
     601             : }
     602             : 
     603             : static GEN
     604          42 : lambdalist(GEN E, long prec)
     605             : {
     606          42 :   pari_sp ltop = avma;
     607          42 :   GEN glob = ellglobalred(E), plist = gmael(glob,4,1), L = gel(glob,5);
     608          42 :   GEN res, v, D = ell_get_disc(E);
     609          42 :   long i, j, k, l, m, n, np = lg(plist), lr = 1;
     610          42 :   v = cgetg(np, t_VEC);
     611         147 :   for (j = 1, i = 1 ; j < np; ++j)
     612             :   {
     613         105 :     GEN p = gel(plist, j);
     614         105 :     if (dvdii(D, sqri(p)))
     615             :     {
     616          91 :       GEN la = lambda1(E, gel(L,j), p, prec);
     617          91 :       gel(v, i++) = la;
     618          91 :       lr *= lg(la);
     619             :     }
     620             :   }
     621          42 :   np = i;
     622          42 :   res = cgetg(lr+1, t_VEC);
     623          42 :   gel(res, 1) = gen_0; n = 1; m = 1;
     624         133 :   for (j = 1; j < np; ++j)
     625             :   {
     626          91 :     GEN w = gel(v, j);
     627          91 :     long lw = lg(w);
     628         308 :     for (k = 1; k <= n; k++)
     629             :     {
     630         217 :       GEN t = gel(res, k);
     631         455 :       for (l = 1, m = n; l < lw; l++, m+=n)
     632         238 :         gel(res, k + m) = mpadd(t, gel(w, l));
     633             :     }
     634          91 :     n = m;
     635             :   }
     636          42 :   return gerepilecopy(ltop, res);
     637             : }
     638             : 
     639             : /* P a t_INT or t_FRAC, return its logarithmic height */
     640             : static GEN
     641          98 : heightQ(GEN P, long prec)
     642             : {
     643             :   long s;
     644          98 :   if (typ(P) == t_FRAC)
     645             :   {
     646          56 :     GEN a = gel(P,1), b = gel(P,2);
     647          56 :     P = abscmpii(a,b) > 0 ? a: b;
     648             :   }
     649          98 :   s = signe(P);
     650          98 :   if (!s) return real_0(prec);
     651          84 :   if (s < 0) P = negi(P);
     652          84 :   return glog(P, prec);
     653             : }
     654             : 
     655             : /* t a t_INT or t_FRAC, returns max(1, log |t|), returns a t_REAL */
     656             : static GEN
     657         147 : logplusQ(GEN t, long prec)
     658             : {
     659         147 :   if (typ(t) == t_INT)
     660             :   {
     661          42 :     if (!signe(t)) return real_1(prec);
     662          28 :     if (signe(t) < 0) t = negi(t);
     663             :   }
     664             :   else
     665             :   {
     666         105 :     GEN a = gel(t,1), b = gel(t,2);
     667         105 :     if (abscmpii(a, b) < 0) return real_1(prec);
     668          56 :     if (signe(a) < 0) t = gneg(t);
     669             :   }
     670          84 :   return glog(t, prec);
     671             : }
     672             : 
     673             : /* See GTM239, p532, Th 8.1.18
     674             :  * Return M such that h_naive <= M */
     675             : GEN
     676          98 : hnaive_max(GEN ell, GEN ht)
     677             : {
     678          98 :   const long prec = LOWDEFAULTPREC; /* minimal accuracy */
     679          98 :   GEN b2     = ell_get_b2(ell), j = ell_get_j(ell);
     680          98 :   GEN logd   = glog(absi_shallow(ell_get_disc(ell)), prec);
     681          98 :   GEN logj   = logplusQ(j, prec);
     682          98 :   GEN hj     = heightQ(j, prec);
     683          49 :   GEN logb2p = signe(b2)? addrr(logplusQ(gdivgu(b2, 12),prec), mplog2(prec))
     684          98 :                         : real_1(prec);
     685          98 :   GEN mu     = addrr(divru(addrr(logd, logj),6), logb2p);
     686          98 :   return addrs(addrr(addrr(ht, divru(hj,12)), mu), 2);
     687             : }
     688             : 
     689             : static GEN
     690         147 : qfb_root(GEN Q, GEN vDi)
     691             : {
     692         147 :   GEN a2 = shifti(gel(Q, 1),1), b = gel(Q, 2);
     693         147 :   return mkcomplex(gdiv(negi(b),a2),divri(vDi,a2));
     694             : }
     695             : 
     696             : static GEN
     697       24668 : qimag2(GEN Q)
     698             : {
     699       24668 :   pari_sp av = avma;
     700       24668 :   GEN z = gdiv(negi(qfb_disc(Q)), shifti(sqri(gel(Q, 1)),2));
     701       24668 :   return gerepileupto(av, z);
     702             : }
     703             : 
     704             : /***************************************************/
     705             : /*Routines for increasing the imaginary parts using*/
     706             : /*Atkin-Lehner operators                           */
     707             : /***************************************************/
     708             : 
     709             : static GEN
     710       24668 : qfb_mult(GEN Q, GEN a, GEN b, GEN c, GEN d)
     711             : {
     712       24668 :   GEN A = gel(Q, 1) , B = gel(Q, 2), C = gel(Q, 3), D = qfb_disc(Q);
     713       24668 :   GEN a2 = sqri(a), b2 = sqri(b), c2 = sqri(c), d2 = sqri(d);
     714       24668 :   GEN ad = mulii(d, a), bc = mulii(b, c), e = subii(ad, bc);
     715       24668 :   GEN W1 = addii(addii(mulii(a2, A), mulii(mulii(c, a), B)), mulii(c2, C));
     716       24668 :   GEN W3 = addii(addii(mulii(b2, A), mulii(mulii(d, b), B)), mulii(d2, C));
     717       24668 :   GEN W2 = addii(addii(mulii(mulii(shifti(b,1), a), A),
     718             :                        mulii(addii(ad, bc), B)),
     719             :                  mulii(mulii(shifti(d,1), c), C));
     720       24668 :   if (!equali1(e)) {
     721       22190 :     W1 = diviiexact(W1,e);
     722       22190 :     W2 = diviiexact(W2,e);
     723       22190 :     W3 = diviiexact(W3,e);
     724             :   }
     725       24668 :   return mkqfb(W1, W2, W3, D);
     726             : }
     727             : 
     728             : #ifdef DEBUG
     729             : static void
     730             : best_point_old(GEN Q, GEN NQ, GEN f, GEN *u, GEN *v)
     731             : {
     732             :   long n, k;
     733             :   GEN U, c, d, A = gel(f,1), B = gel(f,2), C = gel(f,3), D = qfb_disc(f);
     734             :   GEN q = mkqfb(mulii(NQ, C), negi(B), diviiexact(A, NQ), D);
     735             :   redimagsl2(q, &U);
     736             :   *u = c = gcoeff(U, 1, 1);
     737             :   *v = d = gcoeff(U, 2, 1);
     738             :   if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
     739             :   for (n = 1;; n++)
     740             :   {
     741             :     for (k = -n; k <= n; k++)
     742             :     {
     743             :       *u = addis(c, k); *v = addiu(d, n);
     744             :       if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
     745             :       *v = subiu(d, n);
     746             :       if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
     747             :       *u = addiu(c, n); *v = addis(d, k);
     748             :       if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
     749             :       *u = subiu(c, n);
     750             :       if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
     751             :     }
     752             :   }
     753             : }
     754             : /* q(x,y) = ax^2 + bxy + cy^2 */
     755             : static GEN
     756             : qfb_eval(GEN q, GEN x, GEN y)
     757             : {
     758             :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     759             :   GEN x2 = sqri(x), y2 = sqri(y), xy = mulii(x,y);
     760             :   return addii(addii(mulii(a, x2), mulii(b,xy)), mulii(c, y2));
     761             : }
     762             : #endif
     763             : 
     764             : static long
     765        6580 : nexti(long i) { return i>0 ? -i : 1-i; }
     766             : 
     767             : /* q0 + i q1 + i^2 q2 */
     768             : static GEN
     769       12313 : qfmin_eval(GEN q0, GEN q1, GEN q2, long i)
     770       12313 : { return addii(mulis(addii(mulis(q2, i), q1), i), q0); }
     771             : 
     772             : /* assume a > 0, return gcd(a,b,c) */
     773             : static ulong
     774       16436 : gcduii(ulong a, GEN b, GEN c)
     775             : {
     776       16436 :   a = ugcdiu(b, a);
     777       16436 :   return a == 1? 1: ugcdiu(c, a);
     778             : }
     779             : 
     780             : static void
     781       24668 : best_point(GEN Q, GEN NQ, GEN f, GEN *pu, GEN *pv)
     782             : {
     783       24668 :   GEN a = mulii(NQ, gel(f,3)), b = negi(gel(f,2)), c = diviiexact(gel(f,1), NQ);
     784       24668 :   GEN D = qfb_disc(f);
     785       24668 :   GEN U, qr = redimagsl2(mkqfb(a, b, c, D), &U);
     786       24668 :   GEN A = gel(qr,1), B = gel(qr,2), A2 = shifti(A,1), AA4 = sqri(A2);
     787             :   GEN V, best;
     788             :   long y;
     789             : 
     790       24668 :   D = absi_shallow(D);
     791             :   /* 4A qr(x,y) = (2A x + By)^2 + D y^2
     792             :    * Write x = x0(y) + i, where x0 is an integer minimum
     793             :    * (the smallest in case of tie) of x-> qr(x,y), for given y.
     794             :    * 4A qr(x,y) = ((2A x0 + By)^2 + Dy^2) + 4A i (2A x0 + By) + 4A^2 i^2
     795             :    *            = q0(y) + q1(y) i + q2 i^2
     796             :    * Loop through (x,y), y>0 by (roughly) increasing values of qr(x,y) */
     797             : 
     798             :   /* We must test whether [X,Y]~ := U * [x,y]~ satisfy (X NQ, Y Q) = 1
     799             :    * This is equivalent to (X,Y) = 1 (note that (X,Y) = (x,y)), and
     800             :    * (X, Q) = (Y, NQ) = 1.
     801             :    * We have U * [x0+i, y]~ = U * [x0,y]~ + i U[,1] =: V0 + i U[,1] */
     802             : 
     803             :   /* try [1,0]~ = first minimum */
     804       24668 :   V = gel(U,1); /* U *[1,0]~ */
     805       24668 :   *pu = gel(V,1);
     806       24668 :   *pv = gel(V,2);
     807       30947 :   if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
     808             : 
     809             :   /* try [0,1]~ = second minimum */
     810       11935 :   V = gel(U,2); /* U *[0,1]~ */
     811       11935 :   *pu = gel(V,1);
     812       11935 :   *pv = gel(V,2);
     813       11935 :   if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
     814             : 
     815             :   /* (X,Y) = (1, \pm1) always works. Try to do better now */
     816        5656 :   best = subii(addii(a, c), absi_shallow(b));
     817        5656 :   *pu = gen_1;
     818        5656 :   *pv = signe(b) < 0? gen_1: gen_m1;
     819             : 
     820        5656 :   for (y = 1;; y++)
     821        9128 :   {
     822             :     GEN Dy2, r, By, x0, q0, q1, V0;
     823             :     long i;
     824       14784 :     if (y > 1)
     825             :     {
     826       10703 :       if (gcduii(y, gcoeff(U,1,1),  Q) != 1) continue;
     827        7308 :       if (gcduii(y, gcoeff(U,2,1), NQ) != 1) continue;
     828             :     }
     829       11396 :     Dy2 = mulii(D, sqru(y));
     830       11396 :     if (cmpii(Dy2, best) >= 0) break; /* we won't improve. STOP */
     831        5740 :     By = muliu(B,y), x0 = truedvmdii(negi(By), A2, &r);
     832        5740 :     if (cmpii(r, A) >= 0) { x0 = subiu(x0,1); r = subii(r, A2); }
     833             :     /* (2A x + By)^2 + Dy^2, minimal at x = x0. Assume A2 > 0 */
     834             :     /* r = 2A x0 + By */
     835        5740 :     q0 = addii(sqri(r), Dy2); /* minimal value for this y, at x0 */
     836        5740 :     if (cmpii(q0, best) >= 0) continue; /* we won't improve for this y */
     837        5733 :     q1 = shifti(mulii(A2, r), 1);
     838             : 
     839        5733 :     V0 = ZM_ZC_mul(U, mkcol2(x0, utoipos(y)));
     840       12313 :     for (i = 0;; i = nexti(i))
     841        6580 :     {
     842       12313 :       pari_sp av2 = avma;
     843       12313 :       GEN x, N = qfmin_eval(q0, q1, AA4, i);
     844       12313 :       if (cmpii(N , best) >= 0) break;
     845       12271 :       x = addis(x0, i);
     846       12271 :       if (ugcdiu(x, y) == 1)
     847             :       {
     848             :         GEN u, v;
     849       12229 :         V = ZC_add(V0, ZC_z_mul(gel(U,1), i)); /* [X, Y] */
     850       12229 :         u = gel(V,1);
     851       12229 :         v = gel(V,2);
     852       12229 :         if (is_pm1(gcdii(u, Q)) && is_pm1(gcdii(v, NQ)))
     853             :         {
     854        5691 :           *pu = u;
     855        5691 :           *pv = v;
     856        5691 :           best = N; break;
     857             :         }
     858             :       }
     859        6580 :       set_avma(av2);
     860             :     }
     861             :   }
     862             : #ifdef DEBUG
     863             :   {
     864             :     GEN oldu, oldv, F = mkqfb(a, b, c, qfb_disc(f));
     865             :     best_point_old(Q, NQ, f, &oldu, &oldv);
     866             :     if (!equalii(oldu, *pu) || !equalii(oldv, *pv))
     867             :     {
     868             :       if (!equali1(gcdii(mulii(*pu, NQ), mulii(*pv, Q))))
     869             :         pari_err_BUG("best_point (gcd)");
     870             :       if (cmpii(qfb_eval(F, *pu,*pv), qfb_eval(F, oldu, oldv)) > 0)
     871             :       {
     872             :         pari_warn(warner, "%Ps,%Ps,%Ps, %Ps > %Ps",
     873             :                           Q,NQ,f, mkvec2(*pu,*pv), mkvec2(oldu,oldv));
     874             :         pari_err_BUG("best_point (too large)");
     875             :       }
     876             :     }
     877             :   }
     878             : #endif
     879             : }
     880             : 
     881             : static GEN
     882       24668 : best_lift(GEN Q, GEN NQ, GEN f)
     883             : {
     884             :   GEN a, b, c, d, dQ, cNQ;
     885       24668 :   best_point(Q, NQ, f, &c, &d);
     886       24668 :   dQ = mulii(d, Q); cNQ = mulii(NQ, c);
     887       24668 :   (void)bezout(dQ, cNQ, &a, &b);
     888       24668 :   return qfb_mult(f, dQ, b, mulii(negi(Q),cNQ), mulii(a,Q));
     889             : }
     890             : 
     891             : static GEN
     892        2478 : lift_points(GEN listQ, GEN f, GEN *pt, GEN *pQ)
     893             : {
     894        2478 :   pari_sp av = avma;
     895        2478 :   GEN yf = gen_0, tf = NULL, Qf = NULL;
     896        2478 :   long k, l = lg(listQ);
     897       27146 :   for (k = 1; k < l; ++k)
     898             :   {
     899       24668 :     GEN c = gel(listQ, k), Q = gel(c,1), NQ = gel(c,2);
     900       24668 :     GEN t = best_lift(Q, NQ, f), y = qimag2(t);
     901       24668 :     if (gcmp(y, yf) > 0) { yf = y; Qf = Q; tf = t; }
     902             :   }
     903        2478 :   *pt = tf; *pQ = Qf; return gc_all(av, 3, &yf, pt, pQ);
     904             : }
     905             : 
     906             : /***************************/
     907             : /*         Twists          */
     908             : /***************************/
     909             : 
     910             : static GEN
     911          56 : ltwist1(GEN E, GEN d, long bitprec)
     912             : {
     913          56 :   pari_sp av = avma;
     914          56 :   GEN Ed = elltwist(E, d), z = ellL1_bitprec(Ed, 0, bitprec);
     915          56 :   obj_free(Ed); return gerepileuptoleaf(av, z);
     916             : }
     917             : 
     918             : /* Return O_re*c(E)/(4*O_vol*|E_t|^2) */
     919             : 
     920             : static GEN
     921          42 : heegner_indexmult(GEN om, long t, GEN tam, long prec)
     922             : {
     923          42 :   pari_sp av = avma;
     924          42 :   GEN Ovr = gabs(imag_i(gel(om, 2)), prec); /* O_vol/O_re, t_REAL */
     925          42 :   return gerepileupto(av, divru(divir(tam, Ovr), 4*t*t));
     926             : }
     927             : 
     928             : /* omega(gcd(D, N)), given faN = factor(N) */
     929             : static long
     930          56 : omega_N_D(GEN faN, ulong D)
     931             : {
     932          56 :   GEN P = gel(faN, 1);
     933          56 :   long i, l = lg(P), w = 0;
     934         196 :   for (i = 1; i < l; i++)
     935         140 :     if (dvdui(D, gel(P,i))) w++;
     936          56 :   return w;
     937             : }
     938             : 
     939             : static GEN
     940          56 : heegner_indexmultD(GEN faN, GEN a, long D, GEN sqrtD)
     941             : {
     942          56 :   pari_sp av = avma;
     943             :   GEN c;
     944             :   long w;
     945          56 :   switch(D)
     946             :   {
     947           0 :     case -3: w = 9; break;
     948           0 :     case -4: w = 4; break;
     949          56 :     default: w = 1;
     950             :   }
     951          56 :   c = shifti(stoi(w), omega_N_D(faN,-D)); /* (w(D)/2)^2 * 2^omega(gcd(D,N)) */
     952          56 :   return gerepileupto(av, mulri(mulrr(a, sqrtD), c));
     953             : }
     954             : 
     955             : static GEN
     956         399 : nf_to_basis(GEN nf, GEN x)
     957             : {
     958         399 :   x = nf_to_scalar_or_basis(nf, x);
     959         399 :   if (typ(x)!=t_COL)
     960         287 :     x = scalarcol(x, nf_get_degree(nf));
     961         399 :   return x;
     962             : }
     963             : 
     964             : static GEN
     965         196 : etnf_to_basis(GEN et, GEN x)
     966             : {
     967         196 :   long i, l = lg(et);
     968         196 :   GEN V = cgetg(l, t_VEC);
     969         595 :   for (i = 1; i < l; i++)
     970         399 :     gel(V,i) = nf_to_basis(gel(et,i), x);
     971         196 :   return shallowconcat1(V);
     972             : }
     973             : 
     974             : static GEN
     975         140 : etnf_get_M(GEN et)
     976             : {
     977         140 :   long i, l = lg(et);
     978         140 :   GEN V = cgetg(l, t_VEC);
     979         448 :   for (i = 1; i < l; i++)
     980         308 :     gel(V,i)=nf_get_M(gel(et,i));
     981         140 :   return shallowmatconcat(diagonal(V));
     982             : }
     983             : 
     984             : static long
     985          49 : etnf_get_varn(GEN et)
     986             : {
     987          49 :   return nf_get_varn(gel(et,1));
     988             : }
     989             : 
     990             : static GEN
     991          98 : heegner_descent_try_point(GEN nfA, GEN z, GEN den, long prec)
     992             : {
     993          98 :   pari_sp av = avma;
     994          98 :   GEN etal = gel(nfA,1), A = gel(nfA,2), cb = gel(nfA,3);
     995          98 :   GEN al = gel(nfA,4), th = gel(nfA, 5);
     996          98 :   GEN et = gel(etal,1), zk = gel(etal, 2), T = gel(etal,3);
     997          98 :   GEN M = etnf_get_M(et);
     998          98 :   long i, j, n = lg(th)-1, l = lg(al);
     999          98 :   GEN u2 = gsqr(gel(cb,1)), r = gel(cb,2);
    1000          98 :   GEN zz = gdiv(gsub(z,r), u2);
    1001          98 :   GEN be = cgetg(n+1, t_COL);
    1002         161 :   for (j = 1; j < l; j++)
    1003             :   {
    1004         105 :     GEN aj = gel(al, j), Aj = gel(A,j);
    1005         392 :     for (i = 1; i <= n; i++)
    1006         287 :       gel(be,i) = gsqrt(gmul(gsub(zz, gel(th,i)), gel(aj,i)), prec);
    1007         357 :     for (i = 0; i <= (1<<(n-1))-1; i++)
    1008             :     {
    1009             :       long eps;
    1010         294 :       GEN s = gmul(den, RgM_solve_realimag(M, be));
    1011         294 :       GEN S = grndtoi(s, &eps), V, S2;
    1012         294 :       gel(be,1+odd(i)) = gneg(gel(be,1+odd(i)));
    1013         294 :       if (eps > -7) continue;
    1014          42 :       S2 = QXQ_sqr(RgV_RgC_mul(zk, S), T);
    1015          42 :       V = gdiv(QXQ_mul(S2, Aj, T), sqri(den));
    1016          42 :       if (typ(V) != t_POL || degpol(V) != 1) continue;
    1017          42 :       if (gequalm1(gel(V,3)))
    1018          42 :         return gerepileupto(av,gadd(gmul(gel(V,2),u2),r));
    1019             :     }
    1020             :   }
    1021          56 :   return gc_NULL(av);
    1022             : }
    1023             : 
    1024             : static GEN
    1025        1785 : heegner_try_point(GEN E, GEN nfA, GEN lambdas, GEN ht, GEN z, long prec)
    1026             : {
    1027        1785 :   long l = lg(lambdas);
    1028             :   long i, eps;
    1029        1785 :   GEN P = real_i(pointell(E, z, prec)), x = gel(P,1);
    1030        1785 :   GEN rh = subrr(ht, shiftr(ellheightoo(E, P, prec),1));
    1031       26572 :   for (i = 1; i < l; ++i)
    1032             :   {
    1033       24829 :     GEN logd = shiftr(gsub(rh, gel(lambdas, i)), -1);
    1034       24829 :     GEN d, approxd = gexp(logd, prec);
    1035       24829 :     d = grndtoi(approxd, &eps);
    1036       24829 :     if (signe(d) > 0 && eps<-10)
    1037             :     {
    1038             :       GEN X, ylist;
    1039          98 :       if (DEBUGLEVEL > 2)
    1040           0 :         err_printf("\nTrying lambda number %ld, logd=%Ps, approxd=%Ps\n", i, logd, approxd);
    1041          98 :       X = heegner_descent_try_point(nfA, x, d, prec);
    1042          98 :       if (X)
    1043             :       {
    1044          42 :         ylist = ellordinate(E, X, prec);
    1045          42 :         if (lg(ylist) > 1)
    1046             :         {
    1047          42 :           GEN P = mkvec2(X, gel(ylist, 1));
    1048          42 :           GEN hp = ellheight(E,P,prec);
    1049          42 :           if (signe(hp) && cmprr(hp, shiftr(ht,1)) < 0 && cmprr(hp, shiftr(ht,-1)) > 0)
    1050          42 :             return P;
    1051           0 :           if (DEBUGLEVEL)
    1052           0 :             err_printf("found non-Heegner point %Ps\n", P);
    1053             :         }
    1054             :       }
    1055             :     }
    1056             :   }
    1057        1743 :   return NULL;
    1058             : }
    1059             : 
    1060             : static GEN
    1061          42 : heegner_find_point(GEN e, GEN nfA, GEN om, GEN ht, GEN z1, long k, long prec)
    1062             : {
    1063          42 :   GEN lambdas = lambdalist(e, prec);
    1064          42 :   pari_sp av = avma;
    1065             :   long m;
    1066          42 :   GEN Ore = gel(om, 1), Oim = gel(om, 2);
    1067          42 :   if (DEBUGLEVEL)
    1068           0 :     err_printf("%ld*%ld multipliers to test: ",k,lg(lambdas)-1);
    1069         966 :   for (m = 0; m < k; m++)
    1070             :   {
    1071         966 :     GEN P, z2 = divru(addrr(z1, mulsr(m, Ore)), k);
    1072         966 :     if (DEBUGLEVEL > 2)
    1073           0 :       err_printf("%ld ",m);
    1074         966 :     P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
    1075         966 :     if (P) return P;
    1076         931 :     if (signe(ell_get_disc(e)) > 0)
    1077             :     {
    1078         819 :       z2 = gadd(z2, gmul2n(Oim, -1));
    1079         819 :       P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
    1080         819 :       if (P) return P;
    1081             :     }
    1082         924 :     set_avma(av);
    1083             :   }
    1084           0 :   pari_err_BUG("ellheegner, point not found");
    1085             :   return NULL; /* LCOV_EXCL_LINE */
    1086             : }
    1087             : 
    1088             : /* N > 1, fa = factor(N), return factor(4*N) */
    1089             : static GEN
    1090          42 : fa_shift2(GEN fa)
    1091             : {
    1092          42 :   GEN P = gel(fa,1), E = gel(fa,2);
    1093          42 :   if (absequaliu(gcoeff(fa,1,1), 2))
    1094             :   {
    1095          21 :     E = shallowcopy(E);
    1096          21 :     gel(E,1) = addiu(gel(E,1), 2);
    1097             :   }
    1098             :   else
    1099             :   {
    1100          21 :     P = shallowconcat(gen_2, P);
    1101          21 :     E = shallowconcat(gen_2, E);
    1102             :   }
    1103          42 :   return mkmat2(P, E);
    1104             : }
    1105             : 
    1106             : /* P = prime divisors of N(E). Return the product of primes p in P, a_p != -1
    1107             :  * HACK: restrict to small primes since large ones won't divide our C-long
    1108             :  * discriminants */
    1109             : static GEN
    1110          42 : get_bad(GEN E, GEN P)
    1111             : {
    1112          42 :   long k, l = lg(P), ibad = 1;
    1113          42 :   GEN B = cgetg(l, t_VECSMALL);
    1114         147 :   for (k = 1; k < l; k++)
    1115             :   {
    1116         105 :     GEN p = gel(P,k);
    1117         105 :     long pp = itos_or_0(p);
    1118         105 :     if (!pp) break;
    1119         105 :     if (! equalim1(ellap(E,p))) B[ibad++] = pp;
    1120             :   }
    1121          42 :   setlg(B, ibad); return ibad == 1? NULL: zv_prod_Z(B);
    1122             : }
    1123             : 
    1124             : /* list of pairs [Q,N/Q], where Q | N and gcd(Q,N/Q) = 1 */
    1125             : static GEN
    1126          42 : find_div(GEN N, GEN faN)
    1127             : {
    1128          42 :   GEN listQ = divisors(faN);
    1129          42 :   long j, k, l = lg(listQ);
    1130             : 
    1131          42 :   gel(listQ, 1) = mkvec2(gen_1, N);
    1132        1624 :   for (j = k = 2; k < l; ++k)
    1133             :   {
    1134        1582 :     GEN Q = gel(listQ, k), NQ = diviiexact(N, Q);
    1135        1582 :     if (is_pm1(gcdii(Q,NQ))) gel(listQ, j++) = mkvec2(Q,NQ);
    1136             :   }
    1137          42 :   setlg(listQ, j); return listQ;
    1138             : }
    1139             : 
    1140             : static long
    1141        8652 : testDisc(GEN bad, long d) { return !bad || ugcdiu(bad, -d) == 1; }
    1142             : /* bad = product of bad primes. Return the NDISC largest fundamental
    1143             :  * discriminants D < d such that (D,bad) = 1 and d is a square mod 4N */
    1144             : static GEN
    1145          42 : listDisc(GEN fa4N, GEN bad, long d, long ndisc)
    1146             : {
    1147          42 :   GEN v = cgetg(ndisc+1, t_VECSMALL);
    1148          42 :   pari_sp av = avma;
    1149          42 :   long j = 1;
    1150             :   for(;;)
    1151             :   {
    1152        8652 :     d -= odd(d)? 1: 3;
    1153        8652 :     if (testDisc(bad,d) && unegisfundamental(-d) && Zn_issquare(stoi(d), fa4N))
    1154             :     {
    1155         420 :       v[j++] = d;
    1156         420 :       if (j > ndisc) break;
    1157             :     }
    1158        8610 :     set_avma(av);
    1159             :   }
    1160          42 :   set_avma(av); return v;
    1161             : }
    1162             : /* L = vector of [q1,q2] or [q1,q2,q2']
    1163             :  * cd = (b^2 - D)/(4N) */
    1164             : static void
    1165      166880 : listfill(GEN N, GEN b, GEN c, GEN d, GEN D, GEN L, long *s)
    1166             : {
    1167      166880 :   long k, l = lg(L);
    1168      166880 :   GEN add, frm2, a = mulii(d, N), V = mkqfb(a,b,c,D), frm = qfbred_i(V);
    1169      600089 :   for (k = 1; k < l; ++k)
    1170             :   { /* Lk = [v,frm] or [v,frm,frm2] */
    1171      597611 :     GEN Lk = gel(L,k);
    1172             :     long i;
    1173     1509515 :     for (i = 2; i < lg(Lk); i++) /* 1 or 2 elements */
    1174     1076306 :       if (gequal(frm, gel(Lk,i)))
    1175             :       {
    1176      164402 :         GEN v = gel(Lk, 1);
    1177      164402 :         if (cmpii(a, gel(v,1)) < 0) gel(Lk,1) = V;
    1178      164402 :         return;
    1179             :       }
    1180             :   }
    1181        2478 :   frm2 = qfbred_i(mkqfb(d, negi(b), mulii(c,N), D));
    1182        2478 :   add = gequal(frm, frm2)? mkvec2(V,frm): mkvec3(V,frm,frm2);
    1183        2478 :   vectrunc_append(L, add);
    1184        2478 :   *s += lg(add) - 2;
    1185             : }
    1186             : /* faN4 = factor(4*N) */
    1187             : static GEN
    1188         420 : listheegner(GEN N, GEN faN4, GEN listQ, GEN D)
    1189             : {
    1190         420 :   pari_sp av = avma;
    1191         420 :   const long kmin = 30;
    1192         420 :   long h = itos(quadclassno(D));
    1193         420 :   GEN ymin, b = Zn_sqrt(D, faN4), L = vectrunc_init(h+1);
    1194         420 :   long l, k, s = 0;
    1195       13020 :   for (k = 0; k < kmin || s < h; k++)
    1196             :   {
    1197       12600 :     GEN bk = addii(b, mulsi(2*k, N));
    1198       12600 :     GEN C = diviiexact(shifti(subii(sqri(bk), D), -2), N);
    1199       12600 :     GEN div = divisors(C);
    1200       12600 :     long i, l = lg(div);
    1201      179480 :     for (i = 1; i < l; i++)
    1202             :     {
    1203      166880 :       GEN d = gel(div, i), c = gel(div, l-i); /* cd = C */
    1204      166880 :       listfill(N, bk, c, d, D, L, &s);
    1205             :     }
    1206             :   }
    1207         420 :   l = lg(L); ymin = NULL;
    1208        2898 :   for (k = 1; k < l; k++)
    1209             :   {
    1210        2478 :     GEN t, Q, Lk = gel(L,k), f = gel(Lk,1);
    1211        2478 :     GEN y = lift_points(listQ, f, &t, &Q);
    1212        2478 :     gel(L, k) = mkvec3(t, stoi(lg(Lk) - 2), Q);
    1213        2478 :     if (!ymin || gcmp(y, ymin) < 0) ymin = y;
    1214             :   }
    1215         420 :   if (DEBUGLEVEL > 1)
    1216           0 :     err_printf("Disc %Ps : N*ymin = %Pg\n", D,
    1217             :                            gmul(gsqrt(ymin, DEFAULTPREC),N));
    1218         420 :   return gerepilecopy(av, mkvec3(ymin, L, D));
    1219             : }
    1220             : 
    1221             : /* Q | N, P = prime divisors of N, R[i] = local epsilon-factor at P[i].
    1222             :  * Return \prod_{p | Q} R[i] */
    1223             : static long
    1224         147 : rootno(GEN Q, GEN P, GEN R)
    1225             : {
    1226         147 :   long s = 1, i, l = lg(P);
    1227         581 :   for (i = 1; i < l; i++)
    1228         434 :     if (dvdii(Q, gel(P,i))) s *= R[i];
    1229         147 :   return s;
    1230             : }
    1231             : 
    1232             : static void
    1233          42 : heegner_find_disc(GEN *points, GEN *coefs, long *pind, GEN E,
    1234             :                   GEN indmult, long ndisc, long prec)
    1235             : {
    1236          42 :   long d = 0;
    1237             :   GEN faN4, bad, N, faN, listQ, listR;
    1238             : 
    1239          42 :   ellQ_get_Nfa(E, &N, &faN);
    1240          42 :   faN4 = fa_shift2(faN);
    1241          42 :   listQ = find_div(N, faN);
    1242          42 :   bad = get_bad(E, gel(faN, 1));
    1243          42 :   listR = gel(obj_check(E, Q_ROOTNO), 2);
    1244             :   for(;;)
    1245           0 :   {
    1246          42 :     pari_sp av = avma;
    1247          42 :     GEN list, listD = listDisc(faN4, bad, d, ndisc);
    1248          42 :     long k, l = lg(listD);
    1249          42 :     list = cgetg(l, t_VEC);
    1250         462 :     for (k = 1; k < l; ++k)
    1251         420 :       gel(list, k) = listheegner(N, faN4, listQ, stoi(listD[k]));
    1252          42 :     list = vecsort0(list, gen_1, 0);
    1253          56 :     for (k = l-1; k > 0; --k)
    1254             :     {
    1255          56 :       long bprec = 8;
    1256          56 :       GEN Lk = gel(list,k), D = gel(Lk,3);
    1257          56 :       GEN sqrtD = sqrtr_abs(itor(D, prec)); /* sqrt(|D|) */
    1258          56 :       GEN indmultD = heegner_indexmultD(faN, indmult, itos(D), sqrtD);
    1259             :       do
    1260             :       {
    1261             :         GEN mulf, indr;
    1262             :         pari_timer ti;
    1263          56 :         if (DEBUGLEVEL) timer_start(&ti);
    1264          56 :         mulf = ltwist1(E, D, bprec+expo(indmultD));
    1265          56 :         if (DEBUGLEVEL) timer_printf(&ti,"ellL1twist");
    1266          56 :         indr = mulrr(indmultD, mulf);
    1267          56 :         if (DEBUGLEVEL) err_printf("Disc = %Ps, Index^2 = %Ps\n", D, indr);
    1268          56 :         if (signe(indr)>0 && expo(indr) >= -1) /* indr >=.5 */
    1269             :         {
    1270             :           long e, i, l;
    1271          42 :           GEN pts, cfs, L, indi = grndtoi(sqrtr_abs(indr), &e);
    1272          42 :           if (e > expi(indi)-7)
    1273             :           {
    1274           0 :             bprec++;
    1275           0 :             pari_warn(warnprec, "ellL1",bprec);
    1276           0 :             continue;
    1277             :           }
    1278          42 :           *pind = itos(indi);
    1279          42 :           L = gel(Lk, 2); l = lg(L);
    1280          42 :           pts = cgetg(l, t_VEC);
    1281          42 :           cfs = cgetg(l, t_VECSMALL);
    1282         189 :           for (i = 1; i < l; ++i)
    1283             :           {
    1284         147 :             GEN P = gel(L,i), z = gel(P,2), Q = gel(P,3); /* [1 or 2, Q] */
    1285             :             long c;
    1286         147 :             gel(pts, i) = qfb_root(gel(P,1), sqrtD);
    1287         147 :             c = rootno(Q, gel(faN,1), listR);
    1288         147 :             if (!equali1(z)) c *= 2;
    1289         147 :             cfs[i] = c;
    1290             :           }
    1291          42 :           if (DEBUGLEVEL)
    1292           0 :             err_printf("N = %Ps, ymin*N = %Ps\n",N,
    1293           0 :                        gmul(gsqrt(gel(Lk, 1), prec),N));
    1294          42 :           *coefs = cfs; *points = pts; return;
    1295             :         }
    1296             :       } while(0);
    1297             :     }
    1298           0 :     d = listD[l-1]; set_avma(av);
    1299             :   }
    1300             : }
    1301             : 
    1302             : GEN
    1303         154 : ellanal_globalred_all(GEN e, GEN *cb, GEN *N, GEN *tam)
    1304             : {
    1305         154 :   GEN E = ellanal_globalred(e, cb), red = obj_check(E, Q_GLOBALRED);
    1306         154 :   *N = gel(red, 1);
    1307         154 :   *tam = gel(red,2);
    1308         154 :   if (signe(ell_get_disc(E))>0) *tam = shifti(*tam,1);
    1309         154 :   return E;
    1310             : }
    1311             : 
    1312             : static GEN
    1313          42 : vecelnfembed(GEN x, GEN M, GEN et)
    1314          91 : { pari_APPLY_same(gmul(M, etnf_to_basis(et, gel(x,i)))) }
    1315             : 
    1316             : static GEN
    1317          42 : QXQV_inv(GEN x, GEN T)
    1318          91 : { pari_APPLY_same(QXQ_inv(gel(x,i), T)) }
    1319             : 
    1320             : static GEN
    1321          42 : etnfnewprec(GEN x, long prec)
    1322         126 : { pari_APPLY_same(nfnewprec(gel(x,i),prec)) }
    1323             : 
    1324             : static GEN
    1325          49 : vec_etnf_to_basis(GEN et, GEN x)
    1326         154 : { pari_APPLY_same(etnf_to_basis(et,gel(x,i))) }
    1327             : 
    1328             : static GEN
    1329          42 : makenfA(GEN sel, GEN A, GEN cb)
    1330             : {
    1331          42 :   GEN etal = gel(sel,1), T = gel(etal,3);
    1332          42 :   GEN et = gel(etal,1), M = etnf_get_M(et);
    1333          42 :   long v = etnf_get_varn(et);
    1334          42 :   GEN al = vecelnfembed(A, M, et);
    1335          42 :   GEN th = gmul(M, etnf_to_basis(et, pol_x(v)));
    1336          42 :   return mkvec5(etal,QXQV_inv(A, T),cb,al,th);
    1337             : }
    1338             : 
    1339             : GEN
    1340          56 : ellheegner(GEN E)
    1341             : {
    1342          56 :   pari_sp av = avma;
    1343             :   GEN z, P, ht, points, coefs, s, om, indmult;
    1344             :   GEN sel, etal, et, cbb, A, dAi, T, Ag, At;
    1345             :   long ind, indx, lint, k, l, wtor, etor, ndisc, ltors2, selrank;
    1346          56 :   long bitprec = 16, prec = nbits2prec(bitprec)+1;
    1347             :   pari_timer ti;
    1348             :   GEN N, cb, tam, torsion, nfA;
    1349          56 :   E = ellanal_globalred_all(E, &cb, &N, &tam);
    1350          56 :   if (ellrootno_global(E) == 1)
    1351           7 :     pari_err_DOMAIN("ellheegner", "(analytic rank)%2","=",gen_0,E);
    1352          49 :   torsion = elltors(E);
    1353          49 :   wtor = itos( gel(torsion,1) ); /* #E(Q)_tor */
    1354          49 :   etor = wtor > 1? itou(gmael(torsion, 2, 1)): 1; /* exponent of E(Q)_tor */
    1355          49 :   sel = ell2selmer_basis(E, &cbb, prec);
    1356          49 :   etal = gel(sel,1); A = gel(sel,2); et = gel(etal,1); T = gel(etal,3);
    1357          49 :   ltors2 = lg(et)-2; selrank = lg(A)-1;
    1358          49 :   Ag = selrank > ltors2+1 ? pol_1(etnf_get_varn(et)): gel(A,selrank);
    1359          49 :   At = vecslice(A,1,ltors2);
    1360          49 :   dAi = gsupnorm(vec_etnf_to_basis(et,A),prec);
    1361             :   while (1)
    1362          42 :   {
    1363             :     GEN hnaive, l1;
    1364             :     long bitneeded;
    1365          91 :     if (DEBUGLEVEL) timer_start(&ti);
    1366          91 :     l1 = ellL1_bitprec(E, 1, bitprec);
    1367          91 :     if (DEBUGLEVEL) timer_printf(&ti,"ellL1");
    1368          91 :     if (expo(l1) < 1 - bitprec/2)
    1369           7 :       pari_err_DOMAIN("ellheegner", "analytic rank",">",gen_1,E);
    1370          84 :     om = ellR_omega(E,prec);
    1371          84 :     ht = divrr(mulru(l1, wtor * wtor), mulri(gel(om,1), tam));
    1372          84 :     if (DEBUGLEVEL) err_printf("Expected height=%Ps\n", ht);
    1373          84 :     hnaive = hnaive_max(E, ht);
    1374          84 :     if (DEBUGLEVEL) err_printf("Naive height <= %Ps\n", hnaive);
    1375          84 :     hnaive = gadd(shiftr(hnaive,-1),glog(dAi,prec));
    1376          84 :     bitneeded = itos(gceil(divrr(hnaive, mplog2(prec)))) + 12;
    1377          84 :     if (DEBUGLEVEL) err_printf("precision = %ld\n", bitneeded);
    1378          84 :     if (bitprec>=bitneeded) break;
    1379          42 :     bitprec = bitneeded;
    1380          42 :     prec = nbits2prec(bitprec)+1;
    1381             :   }
    1382          42 :   indmult = heegner_indexmult(om, wtor, tam, prec);
    1383          42 :   ndisc = maxss(10, (long) rtodbl(ht)/10);
    1384          42 :   heegner_find_disc(&points, &coefs, &ind, E, indmult, ndisc, prec);
    1385          42 :   if (DEBUGLEVEL) timer_start(&ti);
    1386          42 :   s = heegner_psi(E, N, points, bitprec);
    1387          42 :   if (DEBUGLEVEL) timer_printf(&ti,"heegner_psi");
    1388          42 :   l = lg(points);
    1389          42 :   z = mulsr(coefs[1], gel(s, 1));
    1390         147 :   for (k = 2; k < l; ++k) z = addrr(z, mulsr(coefs[k], gel(s, k)));
    1391          42 :   z = gsub(z, gmul(gel(om,1), ground(gdiv(z, gel(om,1)))));
    1392          42 :   if (DEBUGLEVEL) err_printf("z=%.*Pg\n",nbits2ndec(bitprec), z);
    1393          42 :   lint = wtor > 1 ? ugcd(ind, etor): 1;
    1394          42 :   indx = lint*2*ind;
    1395          42 :   if (vals(indx) >= vals(etor))
    1396          35 :     A = mkvec(Ag);
    1397             :   else
    1398           7 :     A = mkvec2(Ag, QXQ_mul(Ag, gel(At,1), T));
    1399          42 :   gmael(sel,1,1) = etnfnewprec(et, prec);
    1400          42 :   nfA = makenfA(sel, A, cbb);
    1401          42 :   P = heegner_find_point(E, nfA, om, ht, gmulsg(2*lint, z), indx, prec);
    1402          42 :   if (DEBUGLEVEL) timer_printf(&ti,"heegner_find_point");
    1403          42 :   if (cb) P = ellchangepointinv(P, cb);
    1404          42 :   return gerepilecopy(av, P);
    1405             : }
    1406             : 
    1407             : /* Modular degree */
    1408             : 
    1409             : static GEN
    1410          70 : ellisobound(GEN e)
    1411             : {
    1412          70 :   GEN M = gel(ellisomat(e,0,1),2);
    1413          70 :   return vecmax(gel(M,1));
    1414             : }
    1415             : /* 4Pi^2 / E.area */
    1416             : static GEN
    1417         140 : getA(GEN E, long prec) { return mpdiv(sqrr(Pi2n(1,prec)), ellR_area(E, prec)); }
    1418             : 
    1419             : /* Modular degree of elliptic curve e over Q, assuming Manin constant = 1
    1420             :  * (otherwise multiply by square of Manin constant). */
    1421             : GEN
    1422          70 : ellmoddegree(GEN E)
    1423             : {
    1424          70 :   pari_sp av = avma;
    1425             :   GEN N, tam, mc2, d;
    1426             :   long b;
    1427          70 :   E = ellanal_globalred_all(E, NULL, &N, &tam);
    1428          70 :   mc2 = sqri(ellisobound(E));
    1429          70 :   b = expi(mulii(N, mc2)) + maxss(0, expo(getA(E, LOWDEFAULTPREC))) + 16;
    1430             :   for(;;)
    1431           0 :   {
    1432          70 :     long prec = nbits2prec(b), e, s;
    1433          70 :     GEN deg = mulri(mulrr(lfunellmfpeters(E, b), getA(E, prec)), mc2);
    1434          70 :     d = grndtoi(deg, &e);
    1435          70 :     if (DEBUGLEVEL) err_printf("ellmoddegree: %Ps, bit=%ld, err=%ld\n",deg,b,e);
    1436          70 :     s = expo(deg);
    1437          70 :     if (e <= -8 && s <= b-8) return gerepileupto(av, gdiv(d,mc2));
    1438           0 :     b = maxss(s, b+e) + 16;
    1439             :   }
    1440             : }

Generated by: LCOV version 1.13