"""
This scripts seaches for totally real number fields K=Q(d_1^(1/2), ...,d_t^(1/2)) which are p-rational. Final result is printed in file examples-<p>-<t>.txt where p and t are to be replaced by their value.
"""

import itertools
load Schirokauer.sage

p = 97     # replace 97 with any other prime > 3
t_max = 10 # remplace with any other positive integer   

ds=[]
for t in range(1,t_max+1):
    print "searchind d_",t,"...",
    gd=open("examples-"+str(p)+"-"+str(t)+".txt","w")
    d = max(ds + [1])  
    found = false
    while not found:
        d += 1
        if gcd(d,prod(ds+[1])*p) != 1 or not d.is_squarefree():
            continue  # we must guarantee that Gal(K/Q) = (Z/2Z)^t 
                      # and p is unramified
        prop_2_10 = true
        for sbset in itertools.product([0,1],repeat=t):
            D=prod([ds[i]^sbset[i] for i in range(t-1)])*d^sbset[t-1] 
            D=D.squarefree_part()
            if D==1:
               continue
            eps=QuadraticField(D).units()[0]  # fundamental unit
            if legendre_symbol(D,p) == -1:
                E=p^2-1
            else:
                E=p-1
            if  Schirokauer(eps,p,E) == [0,0] or QuadraticField(D).class_number() % p == 0:
                prop_2_10 = false
                break         
        if prop_2_10:
            ds.append(d)
            print d
            found = true
    # at this point prop_2_10 == true so the Proposition 2.10 can be applied
    # to prove that K is p-rational
    gd.write(str(ds))
    gd.close()