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

Generated by: LCOV version 1.13