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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ellanal.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17236-ec2ebc9) Lines: 672 731 91.9 %
Date: 2014-12-19 Functions: 56 58 96.6 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 303 371 81.7 %

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

Generated by: LCOV version 1.9