#!/usr/bin/python

# Copyright Mike Doty 2007 <kingtaco@gentoo.org>
# This code is provided as-is!

# This code is some of my thoughts on factoring composites after
# coffee with my father.  It's not good coding, don't take it as
# such.  It's merely me playing around with numbers.

import pprint
from math import sqrt
import array

def isprime(c):
    x = 3
    while x < (sqrt(c) + 1):
        if (c % x) == 0:
#            print str(c)+" is not prime, factor is "+str(x)
            return False
        x += 2
#    print str(c)+" is prime"
    return True

def build_prime_list(n,m):
    x = n
    l = array.array('L')
    c = 0
    f = False
    while x < m:
        if isprime(x):
            f = True
            c += 1
            l.append(x)
        if f:
            f = False
            if c%500 == 0:
                print c
        x += 2
    return l

def build_prime_list2(filename):
    f = file(filename,'r')
    l = array.array('L')
    try:
        l.fromfile(f,2000)
    except:
        pass
    print "Number of primes is " + str(len(l))
    return l

def build_comp_list(l):
    l_a = array.array('L')
    l_b = array.array('L')
    l_c = array.array('L')
    cnt = 0
    for a in l:
        for b in l:
            cnt+=1
            c = a * b
            l_a.append(a)
            l_b.append(b)
            l_c.append(c)
            if cnt % 500 == 0:
                print cnt
    return (l_a,l_b,l_c);

def slow_sort(l): #bubble sort, bleh
    do_loop = True
    while do_loop:
        do_loop = False
        for n in xrange(len(l)-1):
            (a1,b1,c1) = l[n]
            (a2,b2,c2) = l[n+1]
            if c1 > c2: #swap
                do_loop = True
                l[n] = (a2,b2,c2)
                l[n+1] = (a1,b1,c1)
    return l
            


            
def do_magic(l):
    for (a,b,c) in l:
        print "%08d * %08d = %012d sqrt = %8.12f psqr = %08d pdiff=%d" % (a,b,c,sqrt(c),int(sqrt(c)) ** 2, c - (int(sqrt(c)) ** 2) )

def bruteprime(c):
    x = 3
    while x < (sqrt(c) + 1):
        if (c % x) == 0:
            return x
        x += 2
    return 1

def elegantxxxx(n):
    x = sqrt(n)
    l = array.array('L')
    r = array.array('L')
    while x >= 1.0:
        l.append(int(x))
        x = float(int(x)) * (x - float(int(x)))
    q = array.array('L')
    for z in xrange(len(l)):
        tmp = array.array('L')
        all = array.array('L')
        if z == 0:
            q.append(l[0])
            all = q
        else:
            for c in q:
                #tmp.append(abs(c-l[z]))
                #tmp.append(abs(c+l[z]))
                tmp.append((c-l[z]))
                tmp.append((c+l[z]))
            q = tmp
            #all += tmp
            #print q
    for z in q:
        try:
            if abs(z) > 1:
                if n % abs(z) == 0:
                    print "Factor found(1st pass): " + str(abs(z))
                    return (z,q,1)
        except:
            pass
#    for z in all:
#        try:
#            if abs(z) > 1:
#                if n % abs(z) == 0:
#                    print "Factor found(2nd pass): " + str(abs(z))
#                    return (z,all,2)
#        except:
#            pass
    print "Factor not found for :" + str(n)
    return (0,q,0)


def elegant(n):
    x = sqrt(n)
    l = []#array.array('L')
    r = []#array.array('L')
    co = sqrt(n) - float(int(sqrt(n)))
    while x > 1.0:
        l.append(int(x))
        x = float(int(x)) * (x - float(int(x)))
    q = []#array.array('L')
    for z in xrange(len(l)):
        tmp = []
#        all = []
        if z == 0:
            q.append(l[0])
            #all = q
        else:
            for c in q:
                #tmp.append(abs(c-l[z]))
                #tmp.append(abs(c+l[z]))
                tmp.append((c-l[z]))
                tmp.append((c+l[z]))
            q = tmp
            #all += tmp
            #print q
    for z in q:
        try:
            if abs(z) > 1:
                if n % abs(z) == 0:
                    print "Factor found(1st pass): " + str(abs(z))
                    return (z,q,1)
        except:
            pass
#    for z in all:
#        try:
#            if abs(z) > 1:
#                if n % abs(z) == 0:
#                    print "Factor found(2nd pass): " + str(abs(z))
#                    return (z,all,2)
#        except:
#            pass
    print "Factor not found for :" + str(n)
    return (0,q,0)

def psquare(n):
    return (int(sqrt(n))+1) ** 2

def do_magic2(l):
    good = [0,0]
    close = [0,0,0,0]
    psq = [0,0]
    bad = 0
    ctotal = 0
    cdiff_max = 0
    (l_a,l_b,l_c) = l
    for mark in xrange(len(l_c)):
        a=l_a[mark]
        b=l_b[mark]
        c=l_c[mark]
        print "Testing %d * %d = %d" % (a,b,c)
        (x,y,p) = elegant(c)
        if x != 0:
            good[p-1] += 1
        else:
            bad += 1
            for n in y: # see how close we got
                for t in (a,b):
                    ctotal += 1
                    if abs(t-n) == (psquare(c)-c):
                        psq[0] += 1
                    if abs(t-n) < (psquare(c)-c):
                        psq[1] += 1
                    if cdiff_max < (c - psquare(c)):
                        cdiff_max = c - psquare(c)
                    z = float(abs(t - n))/float(t)
                    if z <= 0.25:
                        close[0] += 1
                    elif z > 0.25 and z <= 0.50:
                        close[1] += 1
                    elif z > 0.50 and z <= 0.75:
                        close[2] += 1
                    else: # if z > 0.75:
                        close[3] += 1
    print "Found %d of %d keys using elegant method. %d(%2.2f%%) 1st pass, %d(%2.2f%%) 2nd pass." % (good[0] + good[1],
                                                                                                     len(l_c),
                                                                                                     good[0],
                                                                                                     100 * (float(good[0]) / float(len(l_a))),
                                                                                                     good[1],
                                                                                                     100 * (float(good[1]) / float(len(l_a))))
    print "Of %d mismatches, %2.2f%% tests were within 1-25%%, %2.2f%% within 25-50%%, %2.2f%% within 50-75%%, and %2.2f%% others." %(ctotal,
                                                                                                                                      100*(float(close[0])/float(ctotal)),
                                                                                                                                      100*(float(close[1])/float(ctotal)),
                                                                                                                                      100*(float(close[2])/float(ctotal)),
                                                                                                                                      100*(float(close[3])/float(ctotal)))
    print "Of %d mismatches, %2.2f%% tests were exactly offset of c - psquare(c) and %2.2f%% were less than c - psquare(c).  Max diff was %d." % (ctotal,
                                                                                                                                                  100*(float(psq[0])/float(ctotal)),
                                                                                                                                                  100*(float(psq[1])/float(ctotal)),
                                                                                                                                                  cdiff_max)
        
    
a1 = build_prime_list(3,50)
a2 = build_prime_list(101,150)
a = a1 + a2
#a = build_prime_list2('prime_array')#.tolist()
b = build_comp_list(a)
#print len(b)
#c = slow_sort(b)
#do_magic(c)
do_magic2(b)