import math
try:
    reload (WindowsBMP)
except NameError:
    import WindowsBMP

try:
    reload (ColourMap)
except NameError:
    import ColourMap

#==============================================================
# Pythagoras Triples
#==============================================================

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def hcf (n1, n2):
    "Highest common factor of two integers"

    n1 = abs (int (n1))
    n2 = abs (int (n2))

    while 1:
        if n1 > n2:
            t = n1
            n1 = n2
            n2 = t

        if n1 == 0:
            return n2

        if n1 == 1:
            return 1

        n2 = n2 % n1        
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def series (n, k):
    """
    Return kth triple (a, b, b+n) that satisfies
            a^2 + b^2 = (b+n)^2

    => a^2 = 2nb + n^2         (sub a = 2nk + n)
    =>(2nk + n)^2 = 2nb + n^2
    => 4.n^2.k^2 + 4.n^2.k + n^2 = 2nb + n^2
    => 2n.k^2 + 2nk = b

    if n is an odd square (m**2)
    
       a^2 = 2nb + n^2         (sub a = 2mk + m, m = sqrt(n))
    => 4.m^2.k^2 + 4.m^2.k + m^2 = 2nb + n^2
    => 4.k^2 + 4.k + 1 = 2b + n
    => 2.k^2 + 2.k + (1-n)/2 = b    (n is odd)

    if 2n is a square number (m**2)

       a^2 = 2nb + n^2         (sub a = mk, m = sqrt(2n))
    => m^2.k^2 = 2nb + n^2
    => k^2 = b + n^2/m^2
    => k^2 - n/2 = b           (n must be even if 2n is square)

    Note: only the series where n = (2a+1)^2 or 2n == m^2 are interesting.
    All the others generate triples which are multiples of smaller ones.
    """

    ret = None

    m = int(math.sqrt (2*n))

    if m * m == 2 * n:
        a = m * k
        b = k * k - n / 2
        return (a,b,b+n)

    m = int(math.sqrt (n))

    if m * m == n and (n%2) == 1:
        a = 2 * m * k + m
        b = 2 * k * (k + 1) + (1-n)/2 
        return (a,b,b+n)
    
    a = n * (2 * k + 1)
    b = 2 * n * k * (k + 1)
    return (a,b,b+n)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def primes (n, limit):
    """
    Returns a list of the prime triples generated by the seed 'n',
    (see 'series')

    Prime triples are where there are no common factors and A < B < C
    """

    ret = []
    
    k = 1
    
    while (1):
        tpl = series (n,k)

        a = tpl [0]
        b = tpl [1]
        c = tpl [2]
        
        if a > limit or b > limit:
            return ret
        
        if a < b and b < c:
            h = hcf (a,b)

            if h == 1:
                ret.append (tpl)

        k += 1        

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def seeds (limit):
    """
    Returns a list seeds for generating prime triples
    seeds are either odd squares or half even squares
    """

    seeds = []

    lim = int (math.sqrt (limit) - 1) / 2
    seeds += [(2*i+1)**2 for i in range (lim+1)]
    lim = int (math.sqrt (limit / 2))
    seeds += [2*i*i for i in range (1,lim+1)]
    seeds.sort()

    return seeds

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def allPrimes (limit):
    """
    Returns a list of the prime triples where a and b are less
    than limit
    """

    ret = []
    sds = seeds(limit)

    for i in sds:
        ret += primes (i, limit)
        
    return ret

#==============================================================
# Pythagorean Triples
#==============================================================

class PythagoreanTripleDrawer:
    """
    Draws a bitmap showing the location of Pythagorean Triples
    """
#----------------------------------------------------------------------------
    def __init__(self):
        self.bm_width = 640
        self.bm_height = 640
        self.x_origin = -319    # Centre image at (0,0)
        self.y_origin = -319    
        self.bin = 1            # one square per pixel
        
#----------------------------------------------------------------------------
    def getLimits(self):
        """
        Work out the range of values required to draw the picture
        """

        self.x_max = self.x_origin + self.bin * self.bm_width -1
        self.y_max = self.y_origin + self.bin * self.bm_height -1
        xmax       = max (-self.x_origin, self.x_max)
        ymax       = max (-self.y_origin, self.y_max)
        self.limit = max (xmax, ymax)
#----------------------------------------------------------------------------
    def draw (self, fname, inc_derived=0):
        """
        Draws a bitmap showing the location of all the triples
        that fit in a limit x limit square.
        """

        self.getLimits ()
        points = allPrimes (self.limit)
        self.bitmap = WindowsBMP.WindowsBMP ()
        self.bitmap.create24Bit (self.bm_width, self.bm_height, ColourMap.vdk_blue)

    # Draw derived first so that the primes will overwrite them
    
        if inc_derived:
            self.drawDerived(points)

    # Draw the points in the list

        for pt in points:
            self.drawPoint(pt[0],pt[1],ColourMap.cyan)
    # Save

        self.bitmap.writeFile(fname)
    
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def drawDerived(self, points):
        """
        Derived points are those obtained by multiplying the prime
        triples by an integer.
        """
        for pt in points:
            n = 2
            while 1:
                x = pt [0] * n
                y = pt [1] * n

                if x > self.limit or y > self.limit:
                    break

                self.drawPoint(x,y,ColourMap.blue)

                n += 1
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def drawPoint(self,x,y,colour):
        """
        Draw a point (plus it's 7 reflections)
        """
        points = [[x,y],[y,x],[-x,y],[y,-x],[x,-y],[-y,x],[-x,-y],[-y,-x]]

        for pt in points:
            if self.mapPoint (pt):
                self.bitmap.setPixel (pt[0], pt[1], colour)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def mapPoint(self, pt):
        """
        Convert a point to bitmap co-ordinates.
        Returns false if the point can't be mapped.
        """
        if pt [0] < self.x_origin or pt [0] > self.x_max:
            return 0

        if pt [1] < self.y_origin or pt [1] > self.y_max:
            return 0

        pt [0] = (pt [0] - self.x_origin) / self.bin
        pt [1] = (pt [1] - self.y_origin) / self.bin
        return 1