def testrange(type=0):
    """
    This is the driver function for factoring a range of integers.
    Factors will be found for each odd z between zaa and zb.

    The constants pmap, plist, sqlist, sqresidue, sqroot, are calculated
    in testrange so that they can be calculated only once for the entire range
    of z's to be factored.

    pmap is a biymap showing which positive integers < 10000 are prime.

    plist is the list of positive primes < 10000.

    
    """

# when you find out how,
# write the following tables to a file.
# then write code to have testrange check if the file exists,
# and if it exists, read the tables from the file.
# if the file does not exist, then generate the tables,
# and write the tables to the file.



    pmap = setpmap()
    plist = setplist(pmap)
    import time
    
    

    pow10 = 0

    loop = 0
    while loop == 0:

        if type == 0:
            
            print " \n Last number of digits  used is ",pow10,".\n"
            pow10 = int(raw_input("Enter the number of digits for numbers relatively difficult to factor.   "))
            zb = int(raw_input("Enter the number of z's to factor with that many digits."))

            pow10A = pow10/2 -1
            pow10B = pow10 - pow10A -2



            tenpowerA = pow(10,pow10A)
            tenpowerB = pow(10,pow10B)

            from random import randint
            startA = randint(tenpowerA,10*tenpowerA)
            startB = randint(tenpowerB,10*tenpowerB)



            q1 = 6 * (startA/6) + 6
            q2 = 6 * (startB/6) + 6

            kount = 0



            stop = 0
            while stop == 0:
                kount = kount + 1
                if kount > zb:
                    break
                q = 0
                while q == 0:
                    p1 = q1 -1
                    q = ifprime(p1,plist)

                    if q == 1:
                        q1 = q1 + 6
                        break
                    p1 = q1 + 1
                    q = ifprime(p1,plist)
                
                    if q == 1:
                        q1 = q1 + 6
                        break
                    q1 = q1 + 6

                q = 0
                while q == 0:
                    p2 = q2 - 1
                    q = ifprime(p2,plist)
                
                    if q == 1:
                        q2 = q2 + 6
                        break
                    p2 = q2 + 1
                    q = ifprime(p2,plist)
      
                    if q == 1:
                        q2 = q2 + 6
                        break
                    q2 = q2 + 6


                z = p1 * p2
                twolog2z = 2 * klog2(z)

            
                fac = twofactor(z,plist,pmap,twolog2z)
                if fac[0] == 0:
                    print "\n Could not factor ",z," within a reasonable time. "
                    break
                else:
                    print " z = ",z," x = ",fac[0]," y = ",fac[1]," difficulty = ",fac[2]


        if type == 1:
            
            print " \n Last number of digits  used is ",pow10,".\n"
            pow10 = int(raw_input("Enter the number of digits for numbers relatively difficult to factor.   "))
            zb = int(raw_input("Enter the number of z's to factor with that many digits."))

            pow10A = pow10/2
            pow10B = pow10 - pow10A

            tenpowerA = 10**pow10A
            tenpowerB = 10**pow10B
            fivetenpowerB = 5 * tenpowerB


            q1 = 6 * (tenpowerA/6)
            q2 = 6 * (fivetenpowerB/6)

            kount = 0



            stop = 0
            while stop == 0:
                kount = kount + 1
                if kount > zb:
                    break
                q = 0
                while q == 0:
                    p1 = q1 -1
                    q = ifprime(p1,plist)

                    if q == 1:
                        q1 = q1 + 6
                        break
                    p1 = q1 + 1
                    q = ifprime(p1,plist)
                
                    if q == 1:
                        q1 = q1 + 6
                        break
                    q1 = q1 + 6

                q = 0
                while q == 0:
                    p2 = q2 - 1
                    q = ifprime(p2,plist)
                
                    if q == 1:
                        q2 = q2 + 6
                        break
                    p2 = q2 + 1
                    q = ifprime(p2,plist)
      
                    if q == 1:
                        q2 = q2 + 6
                        break
                    q2 = q2 + 6


                z = p1 * p2
                twolog2z = 2 * klog2(z)
            
                fac = twofactor(z,plist,pmap,twolog2z)
                if fac[0] == 0:
                    print " Could not factor ",z," within a reasonable time. "
                    break
                else:
                    print " z = ",z," x = ",fac[0]," y = ",fac[1]," difficulty = ",fac[2]






def klogp(z,p):
    
    e = 0
    m = z
    while m > 1:
        m = m/p
        e = e + 1
    return e
    
def abcdquad(z,test=0):
    steps = 0
    pmap = setpmap()
    plist = setplist(pmap)
    log2z = klog2(z)
    twolog2z = 2 * log2z
    if test == 0:
        fac = checktriv(z,plist,twolog2z,pmap)
        if fac[0] != 0:
            return fac
    print "\n Begin quadratic factoring."

    root = ksqrt(z)
    roothalf = ksqrt(z/2)
    steps = steps + 1
    x = gcd(z,root)
    if x > 1:
        print " Found factor by integer square root of z has a factor in common with z."
        y = z/x
        print " z = ",z," x = ",x," y = ",y
        return [x,y,steps]



    oldw2 = 0
    w2 = 0
    j = 0
    jcount = 0
    fourlog2z = 2 * twolog2z
    while jcount < fourlog2z:
        steps = steps + 1
        jcount = jcount + 1
        j = j + 1
        t = root + j
        tsq = t * t
        diff = tsq - z
        diffsqrt = ksqrt(diff)

        if test > 0:
            print "\n Testing for difference of squares."
            print " z = ",z," root = ",root," t = ",t," tsq = ",tsq," diff = ",diff," diffsqrt = ",diffsqrt
            
        if diffsqrt * diffsqrt == diff:
            x = t - diffsqrt
            y = t + diffsqrt
            print "\n Found factors by difference of squares."
            print " z = ",z," x = ",x," y = ",y
            return [x,y,steps]
            

        steps = steps + 1

        j = j - 1
        while w2 == oldw2:
            j = j + 1
            v = roothalf + j
            v2 = v * v
            t1 = z - v2
            B = t1/v
            C = t1 - B * v
            u1 = B * v
            usum = v2 + u1 + C
            w2 = B*B - 4 * C

        
        oldw2 = int(w2)
        if w2 < 0:
            if test > 0:
                print "\n z = ",z," j = ",j," w2 = ",w2
            continue
        w = ksqrt(w2)
        wtest = w * w
        if test > 0:
            print "\n"
            print " z = ",z," j = ",j," v = ",v," v2 = ",v2
            print " t2 = ",v2," t1 = ",t1," t0 = ",C
            print " u2 = ",v2," u1 = ",u1," u0 = ",C," usum = ",usum
            print " B = ",B," C = ",C
            print " w2 = ",w2," w = ",w," wtest = ",wtest

            raw_input(" Press Enter when ready to continue. ")

        if wtest == w2:
            root1 = (-B -w)/2
            root2 = (-B+w)/2
            x = gcd(z,v -root1)
            y = z/x
            if test > 0:
                print "\n Test1: root1 = ",root1," z = ",z," x = ",x," y = ",y, " x * y = ",x*y
            if x*y != z:
                x = gcd(z,v -root2)
                y = z/x
                if test > 0:
                    print "\n Test2: "," root2 = ",root2," z = ",z," x = ",x," y = ",y, " x * y = ",x*y
            if 1 < x < z:
                print " z = ",z," x = ",x," y = ",y
                return [x,y,steps]

        for kk in range(1,fourlog2z):
            k = root - kk
            steps = steps + 1
            g = k*k + B * k + C
            x = gcd(z,g)
            if 1 < x < z:
                y = z/x
                print "\n Found factors by using z is in ring a + b sqrt(",w2,")."
                print " z = ",z," x = ",x," y = ",y
                return [x,y,steps]


    print " Could not factor by quadratic factoring."
    return [0,0,0]
    




    
def abcdpower(z):
    trace = 0
    print "\n Factoring ",z," by abcdpower."
    
    if z < 0:
        print "\n z is negative."
        return [-1,-z,0]
    if z < 4:
        print "\n z is < 4"
        return [1,z,0]
    if z%2 == 0:
        print "\n z is even."
        return [2,z,1]
    pmap = setpmap()
    plist = setplist(pmap)
    q = ifprime(z,plist)
    if q == 1:
        print "\n z is probably prime."
        return [1,z,2]

    
    root = ksqrt(z)
    if root * root == z:
        print "\n Found factors by z is square."
        print " z = ",z," x = ",root," y = ",root
        return [root,root,3]

    steps = 3
    for g in range(1,2):
        steps = steps + 1

        high = root + g
        highsq = high * high
        diff = highsq - z
        diffsqrt = ksqrt(diff)
        if diffsqrt * diffsqrt == diff:
            print "\n Found factors by difference of squares."
            x = high - diffsqrt
            y = high + diffsqrt
            if x * y != z:
                print "\n Error in theory of routine abcdpower."
                print " z = ",z," k = ",k," j = ",j," B = ",B," C = ",C," x = ",x," y = ",y
                return [ x,y,steps]
            return [x,y,steps]
            
        p = plist[g]
        if z%p == 0:
            print "\n Found factors by direct division by prime."
            x = p
            y = z/p
            print " z = ",z," x = ",x," y = ",y
            return [x,y,steps]
        logpz = klogp(z,p) + 1
        ppowk = 1

        for k in range(1,logpz):
            
            ppowk = ppowk * 2

            t1 = z - ppowk*ppowk
            BpC = t1/ppowk
            t3 = t1 - ppowk * BpC

            if trace > 1:
                print "\n Trace:  z = ",z," logpz = ",logpz," ppowk = ",ppowk," t1 = ",t1," t3 = ",t3," BpC = ",BpC

            discrim = BpC * BpC - 4 * t3
            if discrim == 0:
                steps = steps + 1
                x = ppowk + BpC/2
                y = z/x
                if x * y != z:
                    print "\n Error in theory of routine abcdpower."
                    print " z = ",z," k = ",k," j = ",j," B = ",B," C = ",C," x = ",x," y = ",y
                    return [ x,y,steps]
                if 1 < x < z:
                    return [x,y,steps]

                    
            if discrim > 0:
                steps = steps + 1
                w = ksqrt(discrim)
                B = (w + BpC)/2
                C = BpC - B
                if B * C == t3:
                    x = ppowk + C
                    y = ppowk + B
                    if x * y != z:
                        print "\n Error in theory of routine abcdpower."
                        print " z = ",z," k = ",k," j = ",j," B = ",B," C = ",C," x = ",x," y = ",y
                        return [ x,y,steps]
                    if 1 < x < z:
                        return [x,y,steps]

            expon = logpz - ppowk - 2
            if expon > -1:
                ppowj = pow(p,logpz - ppowk - 2)
                steps = steps + 1
                t1 = z - ppowk * ppowj
                C = t1/ppowj
                t2 = t1 - ppowj * C
                B = t2/ppowk
                t3 = t2 - ppowk * B



                if trace > 1:
                    print "\n Trace: ppowk = ",ppowk," ppowj = ",ppowj," t1 = ",t1," C = ",C," t2 = ",t2," B = ",B," t3 = ",t3
                
               
                if t3 == B * C:
                    x = ppowk + C
                    y = ppowj + B
                    if x * y != z:
                        print "\n Error in theory of routine abcdpower."
                        print " z = ",z," k = ",k," j = ",j," B = ",B," C = ",C," x = ",x," y = ",y
                        return [ x,y,steps]
                    if 1 < x < z:
                        return [x,y,steps]



    print " abcdpower failed to find factors."
    return [0,0,0]




def checktriv(z,plist,twolog2z,pmap):
    
    
    if z < 0:
        print "\n z is negative."
        return [-1,-z,0]
    if z < 4:
        print "\n z is < 4"
        return [1,z,0]

    if z < 10000:
        q = pmap[z]
        if q == z:
            print "\n z is prime < 10000."
            return [1,z,1]
    else:
        q = ifprime(z,plist)
        if q == 1:
            print "\n z is prime."
            return [1,z,1]


    limit = int(twolog2z)
    if limit > 1230:
        limit = 1230
    for j in range(1,limit):
        p = plist[j]
        if z%p == 0:
            print "\n Found factor by direct division by prime < ",plist[limit],"."
            return [p,z/p,j]
    return [0,0,0]





def twofactor(z,plist,pmap,twolog2z):
    """

    Apply each of several factoring methods to factor z.
    plist is the list of the positive primes < 10000.
    pmap is the bitmap of primes < 10000.
    Actually, pmap is not quite a bitmap.
    pmap(p) = p if and only if p is prime < 10000.
    pmap(p) = 0 if and only if p is nonprime < 10000.

    twolog2z is twice the integer part of log base 2 of z.


    """
    print "\n\n\n"
    import time
    tbegin = time.time()
    if z < 0:
        print "\n z is negative."
        tend = time.time()
        return [-1,-z,0,tend - tbegin ]
    if z < 4:
        print "\n z is < 4"
        tend = time.time()
        return [1,z,0,tend - tbegin ]

    if z < 10000:
        q = pmap[z]
        if q == z:
            print "\n z is prime < 10000."
            tend = time.time()
            return [1,z,1,tend - tbegin ]
    else:
        q = ifprime(z,plist)
        if q == 1:
            print "\n z is prime."
            tend = time.time()
            return [1,z,1,tend - tbegin]


    tbegin = time.time()
    limit = int(twolog2z)
    if limit > 1230:
        limit = 1230
    for j in range(1,limit):
        p = plist[j]
        if z%p == 0:
            tend = time.time()
            print "\n Found factor in ",tend - tbegin, " seconds, by direct division by prime < ",plist[limit],"."
            return [p,z/p,j,tend - tbegin ]

    tend = time.time()
    print "\n Direct division by prime gave up after ",tend - tbegin," seconds."



    tbegin = time.time()
    fac = abcdmn(z,twolog2z,plist)
    if fac[0] != 0:
        tend = time.time()
        print "\n abcdmn found factors in ", tend - tbegin, " seconds."
        return fac
    tend = time.time()
    print "\n abcdmn gave up after ",tend - tbegin," seconds."


    
    tbegin = time.time()
    fac = abcdquad(z)
    if fac[0] != 0:
        tend = time.time()
        print "\n abcdquad found factors in ",tend - tbegin," seconds."
        return fac
    tend = time.time()
    print "\n abcdquad gave up after ",tend - tbegin," seconds."

    tbegin = time.time()
    fac = abcdpower(z)
    if fac[0] != 0:
        tend = time.time()
        print "\n abcdpower found factors in ",tend - tbegin," seconds."
        return fac
    tend = time.time()
    print "\n abcdpower gave up after ",tend - tbegin," seconds."




    tbegin = time.time()
    fac = algebrafermat(z,twolog2z)
    if fac[0] != 0:
        tend = time.time()
        print "\n algebrafermat found factors in ", tend - tbegin, " seconds."
        return fac
    tend = time.time()
    print "\n algebrafermat gave up after ",tend - tbegin," seconds."


    tbegin = time.time()
    fac = abcd(z,plist,pmap,twolog2z)
    if fac[0] != 0:
        tend = time.time()
        " Found factor by abcd method in ", tend - tbegin," seconds.",  fac ,"= ",fac
        return fac
    print " Could not find factor by abcd method."
    tend = time.time()
    print "\n abcd gave up after ",tend - tbegin," seconds."
    


    tbegin = time.time()
    fac = pminus1(z,twolog2z,plist)
    if fac[0] != 0:
        tend = time.time()
        print "\n Found factor by pminus1 in ",
        return fac
    tend = time.time()
    print "\n p minus 1 gave up after ",tend - tbegin," seconds."

    tbegin = time.time()
    fac = pminus1v2(z,twolog2z)
    if fac[0] != 0:
        tend = time.time()
        print "\n pminus1 version 2 found factor in ", tend - tbegin," seconds."
        return fac
    tend = time.time()
    print "\n p minus 1 version 2 gave up after ",tend - tbegin," seconds."

    
    tbegin = time.time()
    fac = brent(z,twolog2z)
    if fac[0] != 0:
        tend = time.time()
        print "\n Brent's method found factor in ",tend - tbegin," seconds."
        return fac
    tend = time.time()
    print "\n Brent's method gave up after ",tend - tbegin," seconds."

    tbegin = time.time()
    fac = brentpminus1(z,twolog2z)
    if fac[0] != 0:
        tend = time.time()
        print "\n Brent in colaboration with p minus 1 found factor in ",tend - tbegin," seconds."
        return fac
    tend = time.time()
    print "\n p minus 1 using Brent's method gave up after ",tend - tbegin," seconds."


    return [0,0,0]






def algebrafermat(z,twolog2z):
    """

    Extend the concept of factoring by difference of squares.

    Classical Fermat factoring factors z by findint
    t and s such that t * t - s * s = z.
    then the factors of z are ( t - s) and (t + s).

    This routine extends the concept.

    Choose t > sqrt(z).

    if  d = t*t - z is not a perfect square,
    step through some range of integers of the form  k = h*h - d * g * g
    and test if k has any factor in common with z.

    """
    
    
    newlimit = 100 * twolog2z
    steps = 0
    print "\n Trying to factor using Algebraic Fermat."
    steps = steps + 1
    r = ksqrt(z)
    if z == r*r:
        print "\n Found factors by z is square."
        fac = [r,r,1]
        return fac
    steps = steps + 1
    x = gcd(z,r)
    if 1 < x < z:
        print "\n Found factor by integer sqrt of z has a factor in common with z."
        fac = [x,z/x,[steps]]
        return fac

    for j in range(1,twolog2z):
        t = r + j
        d = t*t-z
        steps = steps + 1
        rd = ksqrt(d)
        if rd * rd == d:
            print "\n Found factor by difference of squares."
            x = t - rd
            y = t + rd
            fac = [x,y,[steps,j]]
        steps = steps + 1
        x = gcd(z,d)
        if 1 < x < z:
            print "\n Found factor by gcd with difference of z from square."
            return [x,z/x,[steps,j]]


        BpC = d/t
        BC = d - t * BpC

        w2 = BpC*BpC + 4 * BC
        w = ksqrt(w2)
        if w * w == w2:
            B = (BpC - w)/2
            C = (BpC + w)/2
            x = t - C
            y = t - B
            if x*y != z:
                print "\n Error in algebra fermat theory."
                print " z = ",z," t = ",t," d = ",d," BpC = ",BpC," BC = ",BC," B = ",B," C = ",C," x = ",x," y = ",y
                return [0,0,0]
            print "\n Found factor by Quadratic factoring of z."
            print " z = ",z," x = ",x," y = ",y
            return [x,y,steps]


        for jd in range(1,newlimit):
            td = rd + jd
            td2 = td*td
            for m in range(twolog2z):
                xd = td2 - d * m*m
                if xd < 2:
                    break
                steps = steps + 1
                x = gcd(z,xd)
                if 1 < x < z:
                    print "\n Found factor in ring  a + b sqrt(",d,")."
                    fac = [x,z/x,[steps,j,jd,m]]
                    return fac
    print "\n Could not find factors by Algebraic Fermat in reasonable time."
    return [0,0,0]


def abcdmn(z,twolog2z,plist):
    method = 2
    trace = 0
    print "\n Trying to find factors by abcdmn"


    limitB = 1230
    if limitB > twolog2z:
        limitB = twolog2z

    limitN = 2 * twolog2z

    steps = 0
    for m in range(1,twolog2z):
        if trace > 1:
            print " Outer loop.  m = ",m

        
        for n in range(m,limitN):
            if trace > 1:
                print " Second  nested loop. n = ",n," m = ",m


            mn = m * n
            #print " m*n = ",mn

            steps = steps + 1
            r = ksqrt(z/mn)

            #print " Square root of z/mn = ",r

            
            limit = twolog2z

            if limit > r:
                limit = r

            #print " Range of A is ",r-limit," to ",r,"."

            for j in range(limit):
                A = r - j
                t1 = z - mn * A * A
                t2 = t1/A
                t3 = t1 - A * t2

                #print " A = ",A,"\n t3 = BC = ",t3,"\n t2 = B m + C n = ",t2," t1 = z - mn * A * A = ",t1

                if t3 == 0:
                    x = A
                    y = z/A
                    print "\n Found factors by A is factor of Z"
                    return [x,y,[steps,j]]

                steps = steps + 1


                if method == 2:
                    t2sq = t2 * t2
                    discrim1 = t2sq - 4 * mn * t3
                    if discrim1 > 0:
                        
                        k1 = ksqrt(discrim1)
                        if k1 * k1 == discrim1:
                            sum = (t2 + k1)/2
                            diff = (t2 - k1)/2
                            Csum = sum/n
                            Cdiff = diff/n
                            if n * Csum == sum:
                                x = A * m + Csum
                                y = z/x
                                if x * y == z:
                                    print "\n Found factors by alternative calculation in abcdmn."
                                    print "\n z = ",z," x = ",x," y = ",y
                                    return [x,y,steps]
                            if n * Cdiff == diff:
                                x = A * m + Cdiff
                                y = z/x
                                if x * y == z:
                                    print "\n Found factors by alternative calculation in abcdmn."
                                    print "\n z = ",z," x = ",x," y = ",y
                                    return [x,y,steps]

                    

                    discrim2 = t2sq + 4 * mn * t3
                    k2 = ksqrt(discrim2)
                    if k2 * k2 == discrim2:
                        sum = (t2 + k2)/2
                        diff = (t2 - k2)/2
                        Csum = sum/n
                        Cdiff = diff/n
                        if n * Csum == sum:
                            x = A * m + Csum
                            y = z/x
                            if x * y == z:
                                print "\n Found factors by alternative calculation in abcdmn."
                                print "\n z = ",z," x = ",x," y = ",y
                                return [x,y,steps]
                            
                        if n * Cdiff == diff:
                            x = A * m + Cdiff
                            y = z/x
                            if x * y == z:
                                print "\n Found factors by alternative calculation in abcdmn."
                                print "\n z = ",z," x = ",x," y = ",y
                                return [x,y,steps]
                    
                    

                if method == 1:
                    
                    for jB in range(1,limitB):
                        B = plist[jB]
                        C = t3/B

                        steps = steps + 1
                    
                        if B * C == t3:
                            if t2 == B * m + C * n:
                                x = A * m + C
                                y = A * n + B
                                if x * y != z:
                                    print " Error.  Found factors incorrectly. "
                                    return [x,y,[steps,j,B] ]
                                else:
                                    print "\n Found factors by function abcdmn."
                                    return [x,y,[ steps,j,B]]


    return [0,0,0]


    
    


def abcd(z,plist,pmap,twolog2z):
    import time
    print "\n Attempting to find factors by A * D + B * C = z "
    tbegin = time.time()
    limit = twolog2z * 100
    if limit > 10000:
        limit = 10000
        sqrtlim = ksqrt(limit)

    if z < 0:
        print "\n z is negative."
        tend = time.time()
        print "\n abcd time is ",tend - tbegin," seconds."
        return [-1,-z,0]
    if z < 4:
        print "\n z is < 4."
        tend = time.time()
        print "\n abcd time is ",tend - tbegin," seconds."        
        return [1,z,0]

    if z < limit:
        q = pmap[z]
        if q == z:
            print "\n z is prime < ",limit,"."
            tend = time.time()
            print "\n abcd time is ",tend - tbegin," seconds."

            return [1,z,1]
    else:
        q = ifprime(z,plist)
        if q == 1:
            print "\n z is probably a prime."
            tend = time.time()
            print "\n abcd time is ",tend - tbegin," seconds."
            return [1,z,1]
    
    BC = range(limit)
    ADabove = range(limit)
    ADbelow = range(limit)
    steps = 1
    for j in range(limit):
        BC[j] = [j,1]
        ADabove[j] = [z+j,1]
        ADbelow[j] = [z-j,1]
    for k in range(2,limit):
        zmodk = z%k
        if zmodk == 0:
            print "\n Found factor by division by integer < ",limit,"."
            tend = time.time()
            print "\n abcd time is ",tend - tbegin," seconds."
            return [k,z/k,k]
        for j in range(k-zmodk,limit,k):
            ADabove[j].append(k)
        for j in range(zmodk,limit,k):
            ADbelow[j].append(k)

    
    for j in range(limit):
        lenbelow = len(ADbelow[j])
        for ja in range(  1,lenbelow  ):
            A = ADbelow[j][ja]
            D = ADbelow[j][0]/A
            if A > D:
                break
            for jb in range(1,len(BC[j])):
                B = BC[j][jb]
                C = BC[j][0]/B
                if B > C:
                    break
                steps = steps + 1
                fac = testabcd(z,A,B,C,D,twolog2z,steps,j,ja,jb)
                if fac[0] != 0:
                    print "\n Found factor by ad + bc = z"
                    tend = time.time()
                    print "\n abcd time is ",tend - tbegin," seconds."
                    return fac
        lenabove = len(ADabove[j])
        for ja in range(1,lenabove):
            A = ADabove[j][ja]
            D = ADabove[j][0]/A
            if A > D:
                break
            for jb in range(1,len(BC[j])):
                B = BC[j][jb]
                C = BC[j][0]/B
                if B > C:
                    break
                steps = steps + 1
                fac = testabcd(z,A,B,C,D,twolog2z,steps,j,ja,jb)
                if fac[0] != 0:
                    print "\n Found factor by ad - bc = z"
                    tend = time.time()
                    print "\n abcd time is ",tend - tbegin," seconds."
                    return fac
    print "\n Could not factor z = ",z," by abcd method within ",steps," steps."
    tend = time.time()
    print "\n abcd time is ",tend - tbegin," seconds."
    return  [0,z,steps]




def brent(z,twolog2z):
    steps = 0
    for start in range(100,100+twolog2z):
        a = start
        b = start
        for j in range(twolog2z):
            steps = steps + 1
            a = (a * a + 1)%z
            b = (b * b + 1)%z
            b = (b * b + 1)%z
            x = gcd(z,kabs(b-a))
            if x > 1:
                fac = [x,z/x,[steps,start,j]]
                print "\n Found factors by Brent's method."
                return fac
    return [0,0,0]
    

def brentpminus1(z,twolog2z):
    steps = 0
    kpow = 2
    for start in range(100,100+twolog2z):
        a = start
        b = start
        for j in range(twolog2z):
            steps = steps + 1
            a = (a * a + 1)%z
            b = (b * b + 1)%z
            b = (b * b + 1)%z
            d = kabs(b-a)
            kpow2 = pow(kpow,d,z)
            x = gcd(z,kpow2-1)
            kpow = int(kpow2)
            if x > 1:
                fac = [x,z/x,[steps,start,j]]
                print "\n Found factors by Brent's method applied to p-1."
                return fac
    return [0,0,0]


def gcd(a,b):
    """
    Returns the greatest common divisor of a and b.

    
    """
    if a * b == 0:
        return 0
    r1 = a
    r2 = b
    flag = 0
    while flag == 0:
        r3 = r1%r2
        if r3== 0:
            return r2
        r1 = r2
        r2 = r3





def ifprime(z,plist):
    """
    
    Returns 1 if z is probably prime.
    
    """

    lenp = len(plist)
    for j in range(1,lenp):
        w = plist[j]
        if w == z:
            q = 1
            return q
        q = strong(z,w)
        if q == 0:
            return q
    return 1

def kabs(z):
    if z < 0:
        return -z
    return z


def klog2(z):

    """

    Returns  integer part of log base 2 of z.

    """
    
    
    e = 0
    m = z
    while m > 1:
        m = m/2
        e = e + 1
    return e



def ksqrt(z):

    """

    Returns integer part of square root of z.

    """
    
    if z < 1:
        return 0
    if z == 1:
        return 1
    j = 1
    k = z
    flag = 0
    while flag == 0:
        j = k
        k = (j + z/j)/2
        if k == j + 1:
            return j
        if k == j:
            return j
        

def pminus1(z,twolog2z,plist):
    kpow = 2
    for j in range(2,twolog2z):
        p = plist[j]
        kpow = pow(kpow,p-1,z)
        x = gcd(z,kpow-1)
        if x > 1:
            print "\n Found factors by p minus 1."
            return [x,z/x,j]
    return [0,0,0]


def pminus1v2(z,twolog2z):
    kpow = 2
    for j in range(twolog2z):
        kpow = pow(kpow,z-j,z)
        x = gcd(z,kpow-1)
        if x > 1:
            print "\n Found factors by p minus 1 version 2."
            return [x,z/x,j]
    return [0,0,0]




def setplist(pmap):

    """

    Given the bitmap pmap for primes < 10000,
    generate the prime list of primes < 10000.

    """
    
    plist = range(1230)
    n = 0
    if pmap[3] != 3:
        print "\n pmap has not been calculated! "
        return
    for j in range(2,10000):
        if pmap[j] != -1:
            n = n + 1
            plist[n] = j
    plist[0] = n
    return(plist)




def setpmap():

    """

    Generate the bitmap of primes < 10000.

    """

    
    pmap=range(10001)
    pmap[0] = 0
    pmap[1] = 0
    pmap[2] = 2
    for j in range(3,10000,2):
        pmap[j] = j
        pmap[j+1] = -1
    for k in range(3,100,2):
        if pmap[k] != -1:
            jstart = k*k
            jend = 10000
            jinc = 2 * k
            for j in range(jstart,jend,jinc):
                pmap[j] = -1
    return(pmap)


def strong(z,w):

    """

    Tests whether or not z is strongly  probably prime according to witness w.


    """
    
    
    trace = 0
    if z < 2:
        return 0
    t = z - 1
    a = 0
    while t%2 == 0:
        t = t/2
        a = a + 1
    test = pow(w,t,z)
    if test ==1:
        q = 1
        return q
    if test + 1 == z:
        q = 1
        return q
    for j in range (1,a):
        test = (test * test)%z
        if test + 1 == z:
            q = 1
            return q
    q = 0
    if trace > 0:
        print "\n found ",z," is composite, with witness = ",w
    return q



def testabcd(z,A,B,C,D,twolog2z,steps,j,ja,jb):
    x = gcd(z,A)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,1]]
    x = gcd(z,B)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,2]]
    x = gcd(z,C)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,3]]
    x = gcd(z,D)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,4]]
    x = gcd(z,A+B)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,5]]
    x = gcd(z,A+C)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,6]]
    x = gcd(z,D+B)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,7]]
    x = gcd(z,D+C)
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,8]]
    x = gcd(z,kabs(A-B))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,9]]
    x = gcd(z,kabs(A-C))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,10]]
    x = gcd(z,kabs(D-B))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,11]]
    x = gcd(z,kabs(D-C))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,12]]
    x = gcd(z,kabs(A+B+C-D))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,13]]
    x = gcd(z,kabs(A+B-C+D))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,14]]
    x = gcd(z,kabs(A-B+C+D))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,15]]
    x = gcd(z,kabs(-A+B+C+D))
    if 1 < x < z :
        return [x,z/x,[steps,j,ja,jb,16]]
    AD = A * D
    BpC = B + C
    mmax = (-B-C + ksqrt( BpC * BpC + 4 * A *D) )/(2 * A)
    mid = mmax/2
    
    if AD == z:
        return [A,z/A,[steps,j,ja,jb,17]]
    if AD < z:
        for m in range(twolog2z):
            top = D - B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,18]]
        for n in range(twolog2z):
            top = D - C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,19]]

        for joy in range(twolog2z):
            m = mmax - joy
            top = D - B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,20]]
        for joy in range(twolog2z):
            n = mmax - joy
            top = D - C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,21]]
        for joy in range(twolog2z):
            m = mid + joy
            top = D - B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,22]]
        for joy in range(twolog2z):
            n = mid + joy
            top = D - C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,23]]
    

    else:
        for m in range(twolog2z):
            top = D + B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,24]]
        for n in range(twolog2z):
            top = D + C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,25]]

        for joy in range(twolog2z):
            m = mmax - joy
            top = D + B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,26]]
        for joy in range(twolog2z):
            n = mmax - joy
            top = D + C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,27]]
        for joy in range(twolog2z):
            m = mid + joy
            top = D + B * m
            x = A * m + C
            if x == 0:
                continue
            n = top/x
            if n * x == top:
                y = z/x
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. M = ",m,"."
                    return [x,z/x,[steps,j,ja,jb,28]]
        for joy in range(twolog2z):
            n = mid + joy
            top = D + C * n
            y = A * n + B
            if y == 0:
                continue
            m = top/y
            if y * m == top:
                x = z/y
                if x* y != z:
                    print "\n Error in abcd. z = ",z," A = ",A," B = ",B," C = ",C," D = ",D," m = ",m," x = ",x," y = ",y
                    return [0,0,0]
                if 1 < x < z:
                    print "\n Found factors by testing trial values of M and N in A * M * N + B * M + C * N = D. N = ",n,"."
                    return [y,z/y,[steps,j,ja,jb,29]]


    fib1 = -B
    fib2 = D - B
    for j in range(twolog2z):
        fib3 = fib1 + fib2
        while fib3 > z:
            fib3 = fib3 - z
        x = gcd(z,fib3)
        if 1 < x < z:
            print "\n Found factors by Fibonacci of A * D + B * C = z. sequence number = [",j,",1]."
            return [x,z/x,[steps,j,ja,jb,30]]
        fib1 = fib3 + fib2
        while fib1 > z:
            fib1 = fib1 - z
        x = gcd(z,fib1)
        if 1 < x < z:
            print "\n Found factors by Fibonacci of A * D + B * C = z. sequence number = [",j,",2]."
            return [x,z/x,[steps,j,ja,jb,31]]
        fib2 = fib1 + fib3
        while fib2 > z:
            fib2 = fib2 - z
        x = gcd(z,fib2)
        if 1 < x < z:
            print "\n Found factors by Fibonacci of A * D + B * C = z. sequence number = [",j,",3]."
            return [x,z/x,[steps,j,ja,jb,32]]


    return [0,0,0]



def vmod(b,p):
    w1 = 1
    w2 = 0
    w3 = b
    v1 = 0
    v2 = 1
    v3 = p
    while v3 != 1:
        q = w3/v3
        v1,w1 = w1 - q * v1,v1
        v2,w2 = w2 - q * v2,v2
        v3,w3 = w3 - q * v3,v3
#        print " vmod:  b = ",b," p = ",p," v1,v2,v3 = ",v1,v2,v3," w1,w2,w3 = ",w1,w2,w3," q = ",q

    return (p + v1)%p