/* In the cyclic cubic case, we need z^3/z_K=z^2/(z_K/z) */ /* In the noncyclic cubic case, we need z^2/z_k=z/L_D and zz_k/z_K=z_k/(z_K/z)=zL_D/(z_K/z). */ /* Auxiliary functions. */ /* ZetaDN and LchiN are defined in Hardylittlewood2.gp. */ LchiKNinit(nf, P) = { my(Ebad = [], Pbad = [], D = nf.disc); for (j = 1, #P, my(E, p = P[j], v = idealprimedec(nf, p)); E = if (#v == 1, if (D % p == 0, next); 1 + 'x + 'x^2, #v == 2, if (D % p == 0, 1 - 'x, 1 - 'x^2), #v == 3, (1 - 'x)^2); Ebad = concat(Ebad, E); Pbad = concat(Pbad, p); ); return ([Pbad, Ebad]); } HLW3(A, N, P, D) = { my(flnoncyc = (D != 1), nf, lim); my(LK, va, vb, vc, S1, S2, S3); my(B = getlocalbitprec(), EKbad); localbitprec(32); lim = ceil(B*log(2)/log(N/3)); nf = nfinit(A); localbitprec(B + ceil(1.585*lim) + exponent(lim)); LK = lfundiv(lfuncreate(A), lfuncreate(1)); LK = lfuninit(LK, [1/2, lim, 0]); va = vector(lim); vb = vector(lim); vc = vector(lim); forfactored(X = 1, lim, my([n] = X); S2 = S3 = 0; fordivfactored(X, Y, my([d] = Y, n3 = 3^(n/d), mob = moebius(Y)); if (mob, if (d % 2 && flnoncyc, S2 += mob * (n3 - 1)); if (d % 3, S3 += mob * n3) ) ); va[n] = -S3 / (3*n); if (flnoncyc, vb[n] = -(va[n] + S2 / (2*n))); vc[n] = va[n] - vb[n] - if (n%3 == 0, va[n/3], 0); ); EKbad = LchiKNinit(nf, P); S1 = sum(n = 1, lim, va[n] * log(LchiN(LK, EKbad, n))); S2 = 0; if (flnoncyc, my(LD, Echibad = LchiNinit(D, P)); LD = lfuninit(D, [1/2, lim, 0]); S2 = sum(n = 1, lim, vb[n] * log(LchiN(LD, Echibad, n))) ); S3 = sum(n = 2, lim, vc[n] * log(ZetaDN(P, n))); return (S1 + S2 + S3); } /* Compute the Hardy-Littlewood constant of aX^3+bX^2+cX+d. */ HardyLittlewood3(A, N = 50) = { my(DA = poldisc(A), S = 1., Da6, P, v = variable(A)); if (poldegree(A) != 3, error("polynomial of degree != 3")); my([a, b, c, d] = Vec(A)); if (!polisirreducible(A) || gcd([6 * a, 2 * b, a + b + c, d]) > 1, return (0)); N = max(N, 5); /* Take care of bad primes and p <= N. */ Da6 = abs(6 * a * DA); P = setunion(factor(Da6)[,1]~, primes([5, N])); for (j = 1, #P, my(p = P[j]); S *= (p - #polrootsmod(A, p)) / (p - 1) ); /* Do the real work. */ if (a != 1, A = a^2 * subst(A, v, v / a)); return (S * exp(HLW3(A, N, P, coredisc(DA)))); } HardyLittlewood(A, N = 50) = { my(d = poldegree(A)); if (d <= 0, return (0), d == 1, my(a = abs(polcoeff(A, 1))); return (a / eulerphi(a)), d == 2, return (HardyLittlewood2(A, N)), d == 3, return (HardyLittlewood3(A, N)), error("HardyLittlewood not implemented for d >= 4") ); }