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

Generated by: LCOV version 1.13