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 16375-9f41ae0) Lines: 675 698 96.7 %
Date: 2014-04-19 Functions: 55 55 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 334 393 85.0 %

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

Generated by: LCOV version 1.9