/* 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)));
}