/* Auxiliary functions. */ ZetaDN(P, s) = zeta(s) * prod(j = 1, #P, 1 - P[j]^(-s)); LchiN(L, Ebad, s) = { my([P, E] = Ebad); lfun(L, s) * prod(j = 1, #P, subst(E[j], 'x, P[j]^(-s))); } LchiNinit(D, P) = { my(Ebad = [], Pbad = []); for (j = 1, #P, my(p = P[j], s = kronecker(D, p)); if (s, Ebad = concat(Ebad, 1 - s*'x); Pbad = concat(Pbad, p))); return ([Pbad, Ebad]); } Oddpart(n) = n >> valuation(n,2); /* The real work; D is fundamental. */ HLW2(D, N) = { my(B = getlocalbitprec(), lim, S1, S2, L, P, v, Ebad); localbitprec(32); lim = ceil(B*log(2)/log(N/2)); localbitprec(B + lim + exponent(lim)); L = lfuninit(D, [1/2, lim, 0]); v = vector(lim); forfactored(X = 1, lim, my([n] = X, S = 0); \\ FIXME: loop over odd divisors? fordivfactored(X, Y, my([d] = Y); if (d % 2, S += moebius(Y) << (n/d))); v[n] = S / (2*n); ); P = setunion(factor(abs(D))[,1]~, primes([2, N])); Ebad = LchiNinit(D, P); S1 = sum(n = 1, lim, v[n] * log(LchiN(L, Ebad, n))); S2 = sum(n = 2, lim, (v[n] - if (n%2 == 0, v[n/2])) * log(ZetaDN(P, n))); return (S1 + S2); } /* Compute the Hardy-Littlewood constant of aX^2+bX+c. */ HardyLittlewood2(A, N = 50) = { my(D = poldisc(A), S, P); if (poldegree(A) != 2, error("polynomial of degree != 2")); my([a, b, c] = Vec(A)); if (issquare(D) || gcd([2 * a, a + b, c]) > 1, return (0)); N = max(N, 3); /* Take care of the prime p = 2. */ S = if ((a + b) % 2, 1., 2.); /* Take care of odd primes dividing a. */ P = factor(Oddpart(a))[,1]; for (j = 1, #P, my(p = P[j]); S *= if (b % p, (p - 1) / (p - 2), p / (p - 1)) ); /* Take care of odd primes dividing the index f. */ my([D0, f] = coredisc(D, 1)); P = factor(Oddpart(f))[,1]; S /= prod(j = 1, #P, my(p = P[j]); 1 - kronecker(D0, p) / (p - 1); ); /* Take care of the primes p <= N. */ S *= prodeuler(p = 3, N, 1 - kronecker(D0, p) / (p - 1)); /* Do the real work */ return (S * exp(-HLW2(D0, N))); }