2550.sagews

def intlist(n):
    l=[0]
    i=1
    for i in range(1,n):
        l.append(i)
        l.append(-i)
    return l

def frep(a,b,c,n):
    l=[0]
    for x in intlist(n):
        for y in intlist(n):
            m=a*x^2 +b*x*y +c*y^2
            i=0
            while i<len(l) and l[i]!=m:
                if i==len(l)-1:
                    l.append(m)
                i=i+1
    return l

def period(d):
    l=continued_fraction_list(sqrt(d))
    p=[0,l[0]]
    q=[1,d-(l[0])^2]
    i=1
    while p[i]!=p[i-1] and q[i]!=q[i-1]:
        p.append(l[i]*q[i] -p[i])
        q.append((d-(p[i+1])^2)/q[i])
        i=i+1
    if p[i-1]==p[i]:
        return 2*(i-1)
    if q[i-1]==q[i]:
        return 2*i+1

def seventest(d):
    l=continued_fraction_list(sqrt(d),bits=5*period(d))
    c=continued_fraction(sqrt(d),bits=5*period(d))
    p=[0,l[0]]
    q=[1,d-(l[0])^2]
    i=1
    while i<period(d):
        p.append(l[i]*q[i] -p[i])
        q.append((d-(p[i+1])^2)/q[i])
        i=i+1
    for k in range(period(d)+1):
        if((c.pn(k))^2 -d*(c.qn(k))^2)==-7:
            print "pk/qk=",convergent(l,k),"for k=",k

def pqlist(d):
    l=continued_fraction_list(sqrt(d),bits=5*period(d))
    p=[0,l[0]]
    q=[1,d-(l[0])^2]
    i=1
    while i< period(d):
        p.append(l[i]*q[i] -p[i])
        q.append((d-(p[i+1])^2)/q[i])
        i=i+1
    for x in range(period(d)+1):
        print "P_",x,"is",p[x]," ","Q_",x,"is",q[x]