def sdpoly(n, prec=53):
    v = primes_first_n(n)
    K = RealField(prec)
    R.<x> = K[]
    v = [K(a).sqrt() for a in v]
    C = cartesian_product_iterator([[-1,1]]*n)
    return prod( x + sum(s[i]*a for i,a in enumerate(v)) for s in C)

def sdpoly2(n, prec=53):
    v = primes_first_n(n)
    K = RealIntervalField(prec)
    R.<x> = K[]
    v = [K(a).sqrt() for a in v]
    C = cartesian_product_iterator([[-1,1]]*n)
    return prod( x + sum(s[i]*a for i,a in enumerate(v)) for s in C)

def sdpoly3(n, prec=53):
    """
    Use interval arithmetic -- this is provably correct no matter what prec is.
    """
    v = primes_first_n(n)
    C = cartesian_product_iterator([[-1,1]]*n)
    while True:
        K = RealIntervalField(prec)
        print prec
        R.<x> = K[]
        w = [K(a).sqrt() for a in v]
        f = prod( x + sum(s[i]*a for i,a in enumerate(w)) for s in C)
        z = []
        retry = False
        for a in f.list():
            t, i = a.is_int()
            if not t:
                prec *= 2
                retry = True
                break
            z.append(i)
        if not retry:
            return ZZ['x'](z)


# Code by Jeroen Demeyer to compute Swinnerton-Dyer polynomials very quickly using p-adics:
# Lift a padic `x` to ZZ, but centered around zero.
def centerlift(x):
    modulus = x.parent().prime_pow(x.precision_absolute())
    z = ZZ(x);
    if (2*z > modulus):
        z -= modulus
    return z

# bound = bound on the absolute value of the coefficients of
# the resulting polynomial.  NOT provably correct.
def sdpoly4(n, bound):
    L = primes_first_n(n)
    for p in Primes():
        if all([gcd(p,s) == 1 and is_square(Mod(s,p)) for s in L]):
            break
    prec = ceil(log(bound)/log(p))
    print "Using p =", p
    print "Precision:", p, "^", prec
    K = Qp(p, prec, print_mode="terse", print_pos=False)
    
    sqrts = [sqrt(K(s),extend=False) for s in L]
    
    n = len(L)
    padic_roots = []
    for k in range(0, 2^n):
        binary = ZZ(k).digits(base=2,padto=n)
        root = sum([sqrts[i]*(binary[i]*2-1) for i in range(0,n)])
        padic_roots.append(root)
    t = polygen(K)
    pol_padic = prod([t - r for r in padic_roots])
    coeffs_ZZ = [centerlift(c) for c in pol_padic.list()]
    max_coeff = max([abs(c) for c in coeffs_ZZ])
    print "Largest coefficient:", p, "^", ceil(log(max_coeff)/log(p))
    return PolynomialRing(ZZ, names='t')(coeffs_ZZ)

# Example:
#time swinnerton_dyer([2,3,5,7,11,13,17,19,23,29], 2^4000)   # bound is heuristic