Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

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

LCOV - code coverage report
Current view: top level - basemath - lfunquad.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25406-bf255ab81b) Lines: 356 367 97.0 %
Date: 2020-06-04 05:59:24 Functions: 44 45 97.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2018  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             : /**       L-functions: values at integers of L-functions           **/
      16             : /**             of primitive quadratic characters                  **/
      17             : /********************************************************************/
      18             : #include "pari.h"
      19             : #include "paripriv.h"
      20             : 
      21             : static GEN
      22         784 : RCpol(long k, long t, GEN c)
      23             : {
      24         784 :   GEN P = cgetg(k+3, t_POL);
      25             :   long l;
      26             : 
      27         784 :   gel(P,k+2) = c;
      28        2863 :   for(l = 0; l < k; l++)
      29             :   {
      30        2079 :     c = diviiexact(mulii(c, muluu(2*k-1 - 2*l, k-l)), mulss(l+1, 2*l-t));
      31        2079 :     gel(P,k-l+1) = c;
      32             :   }
      33         784 :   P[1] = evalsigne(1) | evalvarn(0); return P;
      34             : }
      35             : static GEN
      36         301 : vecRCpol(long r, long d)
      37             : {
      38         301 :   long k, K = d - 1, t = 2*r - 3;
      39         301 :   GEN v = cgetg(d + 1, t_VEC), c = int2n(2*K);
      40         784 :   for (k = 0; k <= K; k++)
      41             :   { /* c = 2^(2K) binomial(n/2,k), an integer */
      42         784 :     gel(v,k+1) = RCpol(k, t, c);
      43         784 :     if (k == K) break;
      44         483 :     c = diviuexact(muliu(c, t - 2*k), 2*k + 2);
      45             :   }
      46         301 :   return v;
      47             : }
      48             : /* D a t_INT */
      49             : static GEN
      50        2289 : RgXV_rescale(GEN v, GEN D)
      51             : {
      52             :   long j, l;
      53             :   GEN w;
      54        2289 :   if (equali1(D)) return v;
      55        2289 :   w = cgetg_copy(v, &l);
      56       15862 :   for (j = 1; j < l; j++) gel(w,j) = RgX_rescale(gel(v,j), D);
      57        2289 :   return w;
      58             : }
      59             : static GEN
      60       95599 : euler_sumdiv(GEN q, long v)
      61             : {
      62       95599 :   GEN u = addui(1, q);
      63      123354 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
      64       95599 :   return u;
      65             : }
      66             : 
      67             : /* [p^{k-1},p^{k-3},...,p^{k-2(d-1)-1}] * (s/p), s = 1 or -1 */
      68             : static GEN
      69        6447 : vpowp(long k, long d, long p, long s)
      70             : {
      71        6447 :   GEN v = cgetg(d + 1, t_VEC), p2 = sqru(p);
      72             :   long j;
      73        6447 :   gel(v, d) = powuu(p, k - 2*d + 1);
      74        6447 :   if (s == -1 && (p & 3L) == 3) togglesign_safe(&gel(v,d));
      75       95599 :   for (j = d-1; j >= 1; j--) gel(v, j) = mulii(p2, gel(v, j+1));
      76        6447 :   return v;
      77             : }
      78             : static GEN
      79         231 : usumdivk_0_all(long k, long d)
      80             : {
      81         231 :   GEN v = cgetg(d + 1, t_COL);
      82             :   long j;
      83         231 :   constbern(k >> 1);
      84         861 :   for (j = 1; j <= d; j++)
      85             :   {
      86         630 :     long n = k + 2 - 2*j;
      87         630 :     gel(v,j) = gdivgs(bernfrac(n), - (n << 1));
      88             :   }
      89         231 :   return v;
      90             : }
      91             : static GEN
      92        2961 : usumdivk_fact_all(GEN fa, long k, long d)
      93             : {
      94             :   GEN res, P, E, pow;
      95             :   long i, j, l;
      96        2961 :   res = cgetg(d + 1, t_COL);
      97        2961 :   P = gel(fa, 1); l = lg(P);
      98        2961 :   E = gel(fa, 2); pow = cgetg(l, t_VEC);
      99        8106 :   for (i = 1; i < l; i++) gel(pow, i) = vpowp(k, d, P[i], 1);
     100       45192 :   for (j = 1; j <= d; j++)
     101             :   {
     102       42231 :     GEN v = cgetg(l, t_VEC);
     103      134960 :     for (i = 1; i < l; i++) gel(v,i) = euler_sumdiv(gmael(pow,i,j), E[i]);
     104       42231 :     gel(res, j) = ZV_prod(v);
     105             :   }
     106        2961 :   return res;
     107             : }
     108             : 
     109             : /* Hadamard product */
     110             : static GEN
     111        3647 : RgV_mul(GEN a, GEN b)
     112             : {
     113        3647 :   long j, l = lg(a);
     114        3647 :   GEN v = cgetg(l, t_COL);
     115       57463 :   for (j = 1; j < l; j++) gel(v,j) = gmul(gel(a,j), gel(b,j));
     116        3647 :   return v;
     117             : }
     118             : static GEN
     119        1512 : RgV_multwist(GEN a, GEN P, long k, long dim, long d, long v2, long N4)
     120             : {
     121        1512 :   GEN v = cgetg(dim+1, t_COL);
     122             :   long j;
     123        4879 :   for (j = 1; j <= d; j++)
     124             :   {
     125             :     GEN z;
     126        3367 :     gel(v,j) = z = gmul(gel(a,j), gel(P,j));
     127        3367 :     if (j + d <= dim)
     128             :     {
     129        2086 :       if (N4 == 3) z = negi(z);
     130        2086 :       if (v2) z = shifti(z, (k - 2*j + 1)*v2);
     131        2086 :       gel(v, j + d) = z;
     132             :     }
     133             :   }
     134        1512 :   return v;
     135             : }
     136             : 
     137             : /* r = k - 2*j, 0<=j<d, factor s=an+b, 0<=s<lim. Check if n starts at 0 or 1
     138             :  * P(D,(an+b)^2), (D-s^2)/N = (D-b^2)/N - 2abn/N - a^2n^2/N and guarantee
     139             :  *  N | D-b^2, N | 2ab, and N | a^2 (except N=8, D odd):
     140             :  * N=4: a=2, b=0,1\equiv D: D = 0,1 mod 4.
     141             :  * N=8: a=4, b=2 if D/4 odd, 0 if D/4 even: D = 0 mod 4 or 1 mod 8
     142             :  * N=12: a=6, b=3 if D odd, 0 if D even: D = 0,1 mod 4
     143             :  * N=-12: a=6, b=5,1 if D odd, 4,2 if D even: D = 0,1 mod 4
     144             :  * N=16: a=8, b=7,1 if D = 1 mod 16, 5,3 if D = 9 mod 16: D = 1 mod 8 */
     145             : /* Cost: O( sqrt(D)/a d^3 log(D) ) */
     146             : static GEN
     147        1274 : sigsum(long k, long d, long a, long b, long D, long N, GEN vs, GEN vP)
     148             : {
     149             :   pari_sp av;
     150        1274 :   GEN S, keep0 = NULL, vPD = RgXV_rescale(vP, stoi(D));
     151        1274 :   long D2, n, c1, c2, s, lim = usqrt(labs(D));
     152             : 
     153        1274 :   D2 = (D - b*b)/N; c1 = (2*a*b)/N; c2 = (a*a)/N;
     154        1274 :   av = avma; S = zerocol(d);
     155        4921 :   for (s = b, n = 0; s <= lim; s += a, n++)
     156             :   {
     157        3647 :     long Ds = c2 ? D2 - n*(c2*n + c1) : D2 - ((n*(n+1)) >> 1);
     158        3647 :     GEN v, P = gsubst(vPD, 0, utoi(s*s));
     159        3647 :     if (vs)
     160        1232 :       v = gel(vs, Ds+1);
     161             :     else
     162        2415 :       v = Ds? usumdivk_fact_all(factoru(Ds), k, d)
     163        2415 :             : usumdivk_0_all(k,d);
     164        3647 :     v = RgV_mul(v, P);
     165        3647 :     if (!s) keep0 = gclone(v); else S = gadd(S, v);
     166        3647 :     if (gc_needed(av, 1)) S = gerepileupto(av, S);
     167             :   }
     168        1274 :   S = gmul2n(S, 1);
     169        1274 :   if (keep0) { S = gadd(S, keep0); gunclone(keep0); }
     170        1274 :   return S;
     171             : }
     172             : 
     173             : static GEN
     174         119 : sigsum4(long k, long d, long D, GEN vs, GEN vP)
     175         119 : { return sigsum(k, d, 2, odd(D), D, 4, vs, vP); }
     176             : 
     177             : /* D != 5 (mod 8) */
     178             : static GEN
     179         210 : sigsum8(long k, long d, long D, GEN vs, GEN vP)
     180             : {
     181         210 :   if (D&1L) return gmul2n(sigsum(k, d, 2, 1, D, 8, vs, vP), -1);
     182         210 :   return sigsum(k, d, 4, 2*odd(D >> 2), D, 8, vs, vP);
     183             : }
     184             : 
     185             : /* D = 0 (mod 3) */
     186             : static GEN
     187         273 : sigsum12(long k, long d, long D, GEN vs, GEN vP)
     188         273 : { return sigsum(k, d, 6, 3*odd(D), D, 12, vs, vP); }
     189             : 
     190             : /* D = 1 (mod 3) */
     191             : static GEN
     192          35 : sigsumm12(long k, long d, long D, GEN vs, GEN vP)
     193             : {
     194          35 :   long fl = odd(D);
     195          35 :   GEN res = sigsum(k, d, 6, 4 + fl, D, 12, vs, vP);
     196          35 :   res = gadd(res, sigsum(k, d, 6, 2 - fl, D, 12, vs, vP));
     197          35 :   return gmul2n(res, -1);
     198             : }
     199             : 
     200             : /* D = 1 (mod 8) */
     201             : static GEN
     202         301 : sigsum16(long k, long d, long D, GEN vs, GEN vP)
     203             : {
     204         301 :   long fl = (D&15L) == 1;
     205         301 :   GEN res = sigsum(k, d, 8, 5 + 2*fl, D, 16, vs, vP);
     206         301 :   return gadd(res, sigsum(k, d, 8, 3 - 2*fl, D, 16, vs, vP));
     207             : }
     208             : 
     209             : /* N = 4 (as above), 8 (factor (1+(D/2))), 12 (factor (1+(D/3))),
     210             :    16 (only D=1 mod 8). */
     211             : static GEN
     212         231 : Dpos(long d, long N, long B)
     213             : {
     214         231 :   GEN vD = cgetg(maxss(B, d) + 1, t_VECSMALL);
     215             :   long D, step, c;
     216         231 :   switch(N)
     217             :   {
     218          49 :     case 4:  D = 5;  step = 1; break;
     219          63 :     case 8:  D = 8;  step = 4; break;
     220          56 :     case 12: D = 12; step = 3; break;
     221          49 :     case 16: D = 17; step = 8; break;
     222          14 :     default: D = 13; step = 3; break; /* -12 */
     223             :   }
     224        1491 :   for (c = 1; c <= d || D <= B; D += step)
     225        1260 :     if (sisfundamental(D)) vD[c++] = D;
     226         231 :   setlg(vD, c); return vD;
     227             : }
     228             : 
     229             : typedef GEN (*SIGMA_F)(long,long,long,GEN,GEN);
     230             : static SIGMA_F
     231         231 : get_S_even(long N)
     232             : {
     233         231 :   switch(N) {
     234          49 :     case 4: return sigsum4;
     235          63 :     case 8: return sigsum8;
     236          56 :     case 12:return sigsum12;
     237          49 :     case 16:return sigsum16;
     238          14 :     default:return sigsumm12; /* -12 */
     239             :   }
     240             : }
     241             : 
     242             : static GEN Lfeq(long D, long k);
     243             : /* Euler numbers: 1, 0, -1, 0, 5, 0, -61,... */
     244             : GEN
     245         203 : eulerfrac(long n)
     246             : {
     247         203 :   pari_sp av = avma;
     248         203 :   if (n < 0) pari_err_DOMAIN("eulerfrac", "index", "<", gen_0, stoi(n));
     249         196 :   if (odd(n)) return gen_0;
     250         189 :   switch(n)
     251             :   {
     252          14 :     case 0: return gen_1;
     253          28 :     case 2: return gen_m1;
     254          14 :     case 4: return utoipos(5);
     255          14 :     case 6: return utoineg(61);
     256          14 :     case 8: return utoipos(1385);
     257          14 :     case 10:return utoineg(50521);
     258          14 :     case 12:return utoipos(2702765);
     259          14 :     case 14:return utoineg(199360981);
     260             :   }
     261          63 :   return gerepileuptoint(av, gmul2n(Lfeq(-4, n+1), 1));
     262             : }
     263             : 
     264             : static GEN
     265         301 : mfDcoefs(GEN F, GEN vD, long d)
     266             : {
     267         301 :   long l = lg(vD), i;
     268         301 :   GEN v = mfcoefs(F, vD[l-1], d), w = cgetg(l, t_COL);
     269         301 :   if (d == 4)
     270         252 :     for (i = 1; i < l; i++) gel(w, i) = gel(v, (vD[i]>>2)+1);
     271             :   else
     272         742 :     for (i = 1; i < l; i++) gel(w, i) = gel(v, vD[i]+1);
     273         301 :   return w;
     274             : }
     275             : 
     276             : static GEN
     277         301 : myinverseimage(GEN M, GEN R, GEN *pden)
     278             : {
     279         301 :   GEN c = Q_remove_denom(QM_gauss_i(M, R, 1), pden);/* M*res / den = R */
     280         301 :   if (!c) pari_err_BUG("theta brackets");
     281         301 :   return c;
     282             : }
     283             : 
     284             : static GEN
     285         301 : Hcol(GEN k, long r, GEN vD, long d, long N2)
     286             : {
     287         301 :   long i, l = lg(vD);
     288             :   GEN v;
     289         301 :   if (r < 5)
     290             :   {
     291         175 :     v = mfDcoefs(mfEH(k),vD,d);
     292         707 :     for (i = 1; i < l; i++)
     293         532 :       if (N2 != 1 && vD[i] % N2) gel(v,i) = gmul2n(gel(v,i), 1);
     294         175 :     return v;
     295             :   }
     296         126 :   v = cgetg(l, t_COL);
     297         896 :   for (i = 1; i < l; i++)
     298             :   {
     299         770 :     pari_sp av = avma;
     300         770 :     GEN c = Lfeq(odd(r)? -vD[i]: vD[i], r); /* fundamental */
     301         770 :     if (N2 != 1 && vD[i] % N2) c = gmul2n(c, 1);
     302         770 :     gel(v, i) = gerepileupto(av, c);
     303             :   }
     304         126 :   return v;
     305             : }
     306             : 
     307             : /***********************************************************/
     308             : /*   Modular form method using Half-Integral Weight forms  */
     309             : /*                      Case D > 0                         */
     310             : /***********************************************************/
     311             : static long
     312         231 : dimeven(long r, long N)
     313             : {
     314         231 :   switch(N)
     315             :   {
     316          49 :     case 4:  return r / 6 + 1;
     317          70 :     case 12: return r / 3 + 1;
     318         112 :     default: return r / 4 + 1;
     319             :   }
     320             : }
     321             : static long
     322         231 : muleven(long N) { return (N == 4)? 1: 2; }
     323             : 
     324             : /* L(\chi_D, 1-r) for D > 0 and r > 0 even. */
     325             : static GEN
     326         231 : modulareven(long D, long r, long N0)
     327             : {
     328         231 :   long B, d, i, l, N = labs(N0);
     329         231 :   GEN V, vs, R, M, C, den, L, vP, vD, k = sstoQ(2*r+1, 2);
     330         231 :   SIGMA_F S = get_S_even(N0);
     331             : 
     332         231 :   d = dimeven(r, N);
     333         231 :   B = muleven(N) * mfsturmNgk(N, k);
     334         231 :   vD = Dpos(d, N0, B);
     335         231 :   vP = vecRCpol(r, d);
     336         231 :   l = lg(vD); B = vD[l-1] / N; V = vecfactoru(1, B);
     337         231 :   vs = cgetg(B+2, t_VEC); gel(vs,1) = usumdivk_0_all(r, d);
     338         777 :   for (i = 1; i <= B; i++) gel(vs, i+1) = usumdivk_fact_all(gel(V,i), r, d);
     339         231 :   M = cgetg(l, t_MAT);
     340         938 :   for (i = 1; i < l; i++)
     341             :   {
     342         707 :     pari_sp av = avma;
     343         707 :     gel(M,i) = gerepileupto(av, S(r, d, vD[i], vs, vP));
     344             :   }
     345         231 :   M = shallowtrans(M);
     346         231 :   if (r == 2*d)
     347             :   { /* r = 2 or (r = 4 and N = 4) */
     348         126 :     GEN v = mfDcoefs(mfderiv(mfTheta(NULL), d+1), vD, 1);
     349         126 :     gel(M, d) = gadd(gel(M, d), gdivgs(v, N*(2*d - 1)));
     350             :   }
     351         231 :   R = Hcol(k, r, vD, 1, (N == 8 || N0 == 12)? N >> 2: 1);
     352             :   /* Cost is O(d^2) * bitsize(result) ~ O(d^3.8) [heuristic] */
     353         231 :   C = myinverseimage(M, R, &den);
     354             : 
     355             :   /* Cost: O( sqrt(D)/c d^3 log(D) ), c from findNeven */
     356         231 :   L = RgV_dotproduct(C, S(r, lg(C)-1, D, NULL, vP));
     357         231 :   return den? gdiv(L, den): L;
     358             : }
     359             : 
     360             : /***********************************************************/
     361             : /*   Modular form method using Half-Integral Weight forms  */
     362             : /*                      Case D < 0                         */
     363             : /***********************************************************/
     364             : 
     365             : static long
     366          70 : dimodd(long r, long kro, long N)
     367             : {
     368          70 :   switch(N)
     369             :   {
     370           0 :     case 1: switch (kro)
     371             :     {
     372           0 :       case -1:return (r + 3) >> 2;
     373           0 :       case 0: return (r + 2)/3;
     374           0 :       case 1: return (r + 1) >> 2;
     375             :     }
     376           7 :     case 3: return kro? (r + 1) >> 1: ((r << 1) + 2)/3;
     377          28 :     case 5: switch (kro)
     378             :     {
     379           0 :       case -1:return (3*r + 2) >> 2;
     380          28 :       case 0: return r;
     381           0 :       case 1: return (3*r - 1) >> 2;
     382             :     }
     383           7 :     case 6: return kro == 1 ? (r + 1) >> 1 : r;
     384          28 :     default: return r;
     385             :   }
     386             : }
     387             : 
     388             : static GEN
     389          70 : Dneg(long n, long kro, long d, long N)
     390             : {
     391          70 :   GEN vD = cgetg(maxss(n, d) + 1, t_VECSMALL);
     392          70 :   long D, c, step, N2 = odd(N)? N: N>> 1;
     393          70 :   switch(kro)
     394             :   {
     395          21 :     case -1: D = -3; step = 8; break;
     396          14 :     case 1:  D = -7; step = 8; break;
     397          35 :     default: D = -8; step = 4; break;
     398             :   }
     399        1694 :   for (c = 1; D >= -n || c <= d; D -= step)
     400        1624 :     if (kross(-D, N2) != -1 && sisfundamental(D)) vD[c++] = -D;
     401          70 :   setlg(vD, c); return vD;
     402             : }
     403             : 
     404             : static GEN
     405          35 : div4(GEN V)
     406             : {
     407          35 :   long l = lg(V), i;
     408          35 :   GEN W = cgetg(l, t_VECSMALL);
     409         329 :   for (i = 1; i < l; i++) W[i] = V[i] >> 2;
     410          35 :   return W;
     411             : }
     412             : 
     413             : static GEN
     414        1498 : usumdivktwist_fact_all(GEN fa, long k, long d)
     415             : {
     416        1498 :   GEN V, P, E, pow, res = cgetg(d + 1, t_VEC);
     417             :   long i, j, l;
     418             : 
     419        1498 :   P = gel(fa, 1); l = lg(P);
     420        1498 :   E = gel(fa, 2);
     421        1498 :   if (l > 1 && P[1] == 2) { l--; P++; E++; } /* odd part */
     422        1498 :   pow = cgetg(l, t_VEC);
     423        2800 :   for (i = 1; i < l; i++) gel(pow, i) = vpowp(k, d, P[i], -1);
     424        1498 :   V = cgetg(l, t_VEC);
     425        4795 :   for (j = 1; j <= d; j++)
     426             :   {
     427        6167 :     for (i = 1; i < l; i++) gel(V,i) = euler_sumdiv(gmael(pow,i,j), E[i]);
     428        3297 :     gel(res, j) = ZV_prod(V);
     429             :   }
     430        1498 :   return res;
     431             : }
     432             : 
     433             : static long
     434          70 : mulodd(long N, long kro)
     435             : {
     436          70 :   if (N == 1 || N == 2) return 1;
     437          56 :   if (kro != 1) return kro? 5: 7;
     438           0 :   if (N == 3) return 4;
     439           0 :   if (N == 5) return 5;
     440           0 :   return 2;
     441             : }
     442             : 
     443             : /* Cost: O( sqrt(D)/a d^3 log(D) ) */
     444             : static GEN
     445        1015 : sigsumtwist(long k, long dim, long a, long b, long Da, long N, GEN vs, GEN vP)
     446             : {
     447        1015 :   GEN vPD, S = zerocol(dim), keep0 = NULL;
     448        1015 :   long D2, n, c1, c2, s, lim = usqrt(Da), d;
     449             :   pari_sp av;
     450             : 
     451        1015 :   if (N > 2 && kross(Da, N == 6 ? 3 : N) == -1) return S;
     452        1015 :   d = (dim + 1) >> 1;
     453        1015 :   vPD = RgXV_rescale(vP, stoi(Da));
     454        1015 :   D2 = (Da - b*b)/N; c1 = (2*a*b)/N; c2 = (a*a)/N;
     455        1015 :   av = avma;
     456        2527 :   for (s = b, n = 0; s <= lim; s += a, n++)
     457             :   {
     458        1512 :     long v2, D4, Ds2, Ds = D2 - n*(c2*n + c1); /* (Da - s^2) / N */
     459             :     GEN v, P;
     460        1512 :     if (!Ds) continue;
     461        1512 :     v2 = vals(Ds); Ds2 = Ds >> v2; D4 = Ds2 & 3L; /* (Ds/2^oo) mod 4 */
     462        1512 :     if (vs)
     463        1323 :       v = gel(vs, Ds+1);
     464             :     else
     465         189 :       v = usumdivktwist_fact_all(factoru(Ds2), k, d);
     466        1512 :     P = gsubst(vPD, 0, utoi(s*s));
     467        1512 :     v = RgV_multwist(v, P, k, dim, d, v2, D4);
     468        1512 :     if (!s) keep0 = gclone(v); else S = gadd(S, v);
     469        1512 :     if (gc_needed(av, 1)) S = gerepileupto(av, S);
     470             :   }
     471        1015 :   S = gmul2n(S, 1);
     472        1015 :   if (keep0) { S = gadd(S, keep0); gunclone(keep0); }
     473        1015 :   return gmul2n(S, -2*(d-1));
     474             : }
     475             : 
     476             : /* Da = |D|; [sum sigma_r^(1)(Da-s^2), sum sigma_r^(2)(Da-s^2)], N = 1 */
     477             : static GEN
     478           0 : sigsumtwist11(long k, long dim, long Da, long N, GEN vs, GEN vP)
     479           0 : { return sigsumtwist(k, dim, 1, 0, Da, N, vs, vP); }
     480             : 
     481             : /* Da = |D| or |D|/4 */
     482             : /* [sum sigma_r^(1)((Da-s^2)/N), sum sigma_r^(2)((Da-s^2)/N)] */
     483             : /* Case N|Da; N not necessarily prime. */
     484             : static GEN
     485         161 : sigsumtwist12p0(long k, long dim, long Da, long N, GEN vs, GEN vP)
     486         161 : { return sigsumtwist(k, dim, N, 0, Da, N, vs, vP); }
     487             : 
     488             : /* [sum sigma_r^(1)((Da-s^2)/p), sum sigma_r^(2)((Da-s^2)/p)] */
     489             : /* Case p\nmid Da */
     490             : /* p = 3: s = +-1 mod 3;
     491             :  * p = 5: s = +-1 mod 5 if Da = 1 mod 5, s = +-2 mod 5 if Da = 2 mod 5;
     492             :  * p = 7: s=+-1, +-2, +-3 if Da=1,4,2 mod 7;
     493             :  * p = 6: s=+-1, +-2, +-3 if Da=1,4,3 mod 6 */
     494             : static GEN
     495         504 : sigsumtwist12pt(long k, long dim, long Da, long N, GEN vs, GEN vP)
     496             : {
     497         504 :   long t = Da%N, e = 0;
     498             :   GEN res;
     499         504 :   if (t == 1) e = 1;
     500         210 :   else if (t == 4) e = 2;
     501          63 :   else if (t == 2 || t == 3) e = 3;
     502         504 :   res = sigsumtwist(k, dim, N, N-e, Da, N, vs, vP);
     503         504 :   if (N-e != e) res = gadd(res, sigsumtwist(k, dim, N, e, Da, N, vs, vP));
     504         504 :   return res;
     505             : }
     506             : 
     507             : static GEN
     508          63 : sigsumtwist12_6(long r, long dim, long Da, long N, GEN vs, GEN vP)
     509             : {
     510          63 :   if (Da%12 == 6) return sigsumtwist12p0(r, dim, Da, N, vs, vP);
     511          42 :   return sigsumtwist12pt(r, dim, Da, N, vs, vP);
     512             : }
     513             : static GEN
     514         602 : sigsumtwist12_N(long r, long dim, long Da, long N, GEN vs, GEN vP)
     515             : {
     516         602 :   if (Da%N == 0) return sigsumtwist12p0(r, dim, Da, N, vs, vP);
     517         462 :   return sigsumtwist12pt(r, dim, Da, N, vs, vP);
     518             : }
     519             : 
     520             : typedef GEN (*SIGMA_Fodd)(long,long,long,long,GEN,GEN);
     521             : static SIGMA_Fodd
     522          70 : get_S_odd(long N)
     523             : {
     524          70 :   if (N == 1) return sigsumtwist11;
     525          70 :   if (N == 6) return sigsumtwist12_6;
     526          63 :   return sigsumtwist12_N;
     527             : }
     528             : 
     529             : /* L(\chi_D, 1-r) for D < 0 and r > 0 odd. */
     530             : static GEN
     531          70 : modularodd(long D, long r, long N0)
     532             : {
     533          70 :   long B, d, i, l, dim, kro = kross(D, 2), Da = labs(D), N = labs(N0);
     534          70 :   GEN V, vs, R, M, C, den, L, vP, vD, vD4, k = sstoQ(2*r+1, 2);
     535          70 :   SIGMA_Fodd S = get_S_odd(N);
     536             : 
     537          70 :   dim = dimodd(r, kro, N); d = (dim + 1) >> 1;
     538          70 :   vP = vecRCpol(r, d);
     539          70 :   B = mulodd(N, kro) * mfsturmNgk(4*N, k);
     540          70 :   vD = Dneg(B, kro, dim + 5, N);
     541          70 :   vD4 = kro ? vD : div4(vD);
     542          70 :   l = lg(vD); B = vD4[l-1] / N; V = vecfactoru(1, B);
     543          70 :   vs = cgetg(B+2, t_VEC); gel(vs,1) = NULL; /* unused */
     544        1379 :   for (i = 1; i <= B; i++) gel(vs,i+1) = usumdivktwist_fact_all(gel(V,i), r, d);
     545          70 :   M = cgetg(l, t_MAT);
     546         665 :   for (i = 1; i < l; i++)
     547             :   {
     548         595 :     pari_sp av = avma;
     549         595 :     gel(M,i) = gerepileupto(av, S(r, dim, vD4[i], N, vs, vP));
     550             :   }
     551          70 :   M = shallowtrans(M);
     552          70 :   R = Hcol(k, r, vD, kro? 1: 4, odd(N)? N: N >>1);
     553             :   /* Cost O(d^2) * bitsize(result) ~ O(d^3.7) [heuristic] */
     554          70 :   C = myinverseimage(M, R, &den);
     555             : 
     556          70 :   if (!kro) Da >>= 2;
     557             :   /* Cost: O( sqrt(D)/c d^3 log(D) ), c from findNodd */
     558          70 :   L = RgV_dotproduct(C, S(r, lg(C)-1, Da, N, NULL, vP));
     559          70 :   if (N0 < 0 && (N0 != -6 || Da%3)) den = den? shifti(den,1): gen_2;
     560          70 :   return den? gdiv(L, den): L;
     561             : }
     562             : 
     563             : /********************************************************/
     564             : /*        Using the Full Functional Equation            */
     565             : /********************************************************/
     566             : /* prod_p (1 - (D/p)p^(-k))
     567             :  * Cost O( D/log(D) (k log(kD))^mu ), mu = multiplication exponent */
     568             : static GEN
     569        1911 : Linv(long D, long k, ulong den)
     570             : {
     571             :   pari_sp av;
     572        1911 :   long s, bit, lim, Da = labs(D), prec;
     573        1911 :   double km = k - 1, B = (k-0.5) * log(km*Da/17.079) + 12; /* 17.079 ~ 2Pi e */
     574             :   forprime_t iter;
     575             :   ulong p;
     576             :   GEN P, Q;
     577        1911 :   if (den) B += log((double)den);
     578        1911 :   bit = maxss((long)(B * k)/(M_LN2 * km), 32) + 32;
     579        1911 :   prec = nbits2prec(bit);
     580        1911 :   lim = (long)exp( (B-log(km)) / km ); /* ~ D / (2Pi e) */
     581        1911 :   u_forprime_init(&iter, 3, lim); av = avma;
     582        1911 :   s = kross(D, 2);
     583        1911 :   if (!s) P = real_1(prec);
     584             :   else
     585             :   {
     586        1113 :     Q = real2n(-k, nbits2prec(bit - k));
     587        1113 :     P = (s == 1)? subir(gen_1, Q): addir(gen_1, Q);
     588             :   }
     589      106386 :   while ((p = u_forprime_next(&iter)))
     590             :   {
     591             :     long bitnew;
     592             :     GEN Q;
     593      104475 :     s = kross(D, p); if (!s) continue;
     594      102305 :     bitnew = (long)(bit - k * log2(p));
     595      102305 :     Q = divrr(P, rpowuu(p, k, nbits2prec(maxss(64, bitnew))));
     596      102305 :     P = s == 1? subrr(P, Q): addrr(P, Q);
     597      102305 :     if (gc_needed(av,1)) P = gerepileuptoleaf(av, P);
     598             :   }
     599        1911 :   return P;
     600             : }
     601             : 
     602             : static GEN
     603        1911 : myround(GEN z, ulong d)
     604             : {
     605             :   long e;
     606        1911 :   if (d) z = mulru(z, d);
     607        1911 :   z = grndtoi(z, &e); if (e >= -4) pari_err_BUG("lfunquad");
     608        1911 :   return d? Qdiviu(z, d): z;
     609             : }
     610             : 
     611             : /* D != 1, k > 2; L(\chi_D, 1-k) using func. eq. */
     612             : static GEN
     613        1911 : Lfeq(long D, long k)
     614             : {
     615             :   GEN z, res;
     616        1911 :   long Da, prec, den = 0;
     617             : 
     618        1911 :   if ((D > 0 && odd(k)) || (D < 0 && !odd(k))) return gen_0;
     619        1911 :   Da = labs(D);
     620        1911 :   if (Da & 3)
     621             :   {
     622        1113 :     long d = (Da - 1) >> 1, kd = k / d;
     623        1113 :     if (odd(kd) && !(k % d) && uisprime(Da)) den = kd * Da;
     624             :   }
     625         798 :   else if (Da == 4) den = 2;
     626        1911 :   z = Linv(D, k, den); prec = lg(z);
     627        1911 :   z = mulrr(z, powrs(divru(Pi2n(1, prec), Da), k));
     628        1911 :   if (Da != 4) { z = mulrr(z, sqrtr_abs(utor(Da,prec))); shiftr_inplace(z,-1); }
     629        1911 :   res = divrr(mpfactr(k-1, prec), z);
     630        1911 :   if (odd(k/2)) togglesign(res);
     631        1911 :   return myround(res, den);
     632             : }
     633             : 
     634             : /* heuristic */
     635             : static long
     636        1379 : usefeq(long D, long k, double c)
     637             : {
     638        1379 :   if (k == 2) return 0;
     639        1281 :   if (D < 0) { k = 2*k; D = -D; }
     640        1281 :   return sqrt(D*c) <= k;
     641             : }
     642             : 
     643             : static long
     644         616 : findNeven(long D, double *c)
     645             : {
     646         616 :   long r = D%3;
     647         616 :   if (!r) { *c = 3; return 12; }
     648         525 :   if ((D&7L) == 1) { *c = 2; return 16; }
     649         476 :   if (!odd(D)) { *c = 2; return 8; }
     650         280 :   if (r == 1) { *c = 1.5; return -12; }
     651         189 :   *c = 1; return 4;
     652             : }
     653             : static long
     654         763 : findNodd(long D, long k, double *c)
     655             : {
     656         763 :   long Dmod8 = D&7L, r;
     657         763 :   if (log(k) > 0.7 * log((double)-D)) { *c = 1; return odd(D)? 2: 1; }
     658         343 :   if (D%7 == 0 && Dmod8 == 5) { *c = 3.5; return 7; }
     659         343 :   if (D%6 == 0) { *c = 3; return 6; }
     660         315 :   if (D%5 == 0) { *c = 2.5; return 5; }
     661         294 :   if (D%3 == 0) { *c = 1.5; return 3; }
     662         245 :   if (Dmod8 == 5)
     663             :   {
     664          63 :     r = smodss(D, 7);
     665          63 :     if (r!=1 && r!=2 && r!=4) { *c = 7./6; return -7; }
     666             :   }
     667         182 :   if (smodss(D, 3) != 1 && !odd(D)) { *c = 1.5; return -6; }
     668         182 :   r = smodss(D, 5); if (r != 2 && r != 3) { *c = 5./4; return -5; }
     669          70 :   *c = 1; return 2;
     670             : }
     671             : 
     672             : /* k <= 0 */
     673             : static GEN
     674        2674 : lfunquadneg_i(long D, long k)
     675             : {
     676             :   double c;
     677             :   long N;
     678             : 
     679        2674 :   if (D == 1) return k == 0 ? gneg(ghalf) : gdivgs(bernfrac(1-k), k-1);
     680        2597 :   if (!sisfundamental(D)) pari_err_TYPE("lfunquad [D not fundamental]",stoi(D));
     681        2597 :   if (k == 0) return D < 0? hclassno(stoi(-D)): gen_0;
     682        2548 :   if ((D > 0 && !odd(k)) || (D < 0 && odd(k))) return gen_0;
     683        1470 :   if (D == -4) return gmul2n(eulerfrac(-k), -1);
     684        1379 :   k = 1 - k;
     685        1379 :   N = D < 0? findNodd(D, k, &c): findNeven(D, &c);
     686        1379 :   if (usefeq(D, k, c)) return Lfeq(D, k);
     687         301 :   return D < 0? modularodd(D,k,N): modulareven(D,k,N);
     688             : }
     689             : /* need k <= 0 and D fundamental */
     690             : GEN
     691        2674 : lfunquadneg(long D, long k)
     692        2674 : { pari_sp av = avma; return gerepileupto(av, lfunquadneg_i(D, k)); }

Generated by: LCOV version 1.13