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

Generated by: LCOV version 1.14