#!/usr/bin/env python

import sys
import getopt
from spigot import Spigot

def gcd(a,b):
    while b != 0:
        a, b = b, a % b
    return a

def spig_decimal(spig, length=20):
    s = ""
    spig.sp.base()
    while len(s) < length:
        s = s + spig.sp.readfmt()
    return s

tan_expr = "tan(theta/2)"
angle_expr = "atan(height/base)"
limit = None
print_angle = print_gaussian = False

opts, args = getopt.gnu_getopt(sys.argv[1:], "dl:as")
for opt, val in opts:
    if opt == "-d":
        tan_expr = "tand(theta/2)"
        angle_expr = "atand(height/base)"
    elif opt == "-l":
        limit = int(val)
    elif opt == "-a":
        print_angle = True
    elif opt == "-s":
        print_gaussian = True

if len(args) != 1:
    usage = """\
usage:   pythangle.py [options] <desired-angle>
options: -d               interpret angle in degrees rather than radians
         -l LIMIT         output only the first LIMIT approximations
         -a               also print the actual angle of each approximation
         -s               also print the Gaussian integer we squared to get it
"""
    if len(args) == 0:
        sys.stdout.write(usage)
        sys.exit(0)
    else:
        sys.stderr.write(usage)
        sys.exit(1)

# The basic idea of this program is as follows.
#
# A Pythagorean triple is obtained by squaring a Gaussian integer.
# Proof: any Gaussian integer u+vi has a norm sqrt(u^2+v^2) which is
# the square root of an integer. So its square, (u+vi)^2, must have a
# norm which _is_ an integer, i.e. that norm together with (u+vi)^2's
# real and imaginary parts form a Pythagorean triple.
#
# Squaring a Gaussian integer doubles its argument. So if we want a
# Pythagorean triple (square Gaussian integer) with an argument near
# theta, we should start by finding a Gaussian integer with an
# argument near _half_ of theta, and then square it. And of course
# finding Gaussian integers near to a given argument is the same
# problem as finding pairs of integers near to a given ratio - it's
# just that the argument / ratio are in a different representation -
# so what we really do is look for rational approximations v/u to
# tan(theta/2), each of which corresponds to a Gaussian integer u+vi
# with argument close to theta/2, whose square is a Pythagorean triple
# with argument close to theta itself.

angle = Spigot(args[0])
tan_halfangle = Spigot(tan_expr, {"theta":angle})

# However, this isn't quite the whole story. A Gaussian integer u+vi
# only gives a _primitive_ Pythagorean triple when squared if it is
# not only in lowest terms itself (that is, gcd(u,v)=1) but also
# exactly one of u,v is even. (Clearly 'lowest terms' rules out u,v
# _both_ even, but u,v both odd causes trouble as well because the
# real and imaginary parts both work out even in that situation.)
#
# But that seems odd. Suppose one of our continued-fraction
# convergents does have u,v both odd (which is of course completely
# plausible); then squaring it gives a Pythagorean triple with all
# three values even, which means that there's a way to scale it down
# by a factor of two. What's _that_ the square of? Surely that, also
# being a Pythagorean triple, must be the square of _some_ Gaussian
# integer - and yet it surely has to be sqrt(2) times (u+vi), which
# isn't integral at all!
#
# The answer to this apparent paradox is that if squaring (u+vi) gives
# some triple A+Bi with A,B both even, then it's not (A/2)+(B/2)i
# which is the square of a smaller Gaussian integer: it's the number
# you get by _swapping_ the real and imaginary parts, i.e.
# (B/2)+(A/2)i. Specifically, this will turn out to be the square of
# (u-v)/2 + (u+v)/2 i. (Which is a Gaussian integer, if u,v are both
# odd, because u-v and u+v are both even.) In argument terms, swapping
# the parts of the Pythagorean triple subtracts its argument from
# pi/2, and hence the corresponding transformation to its square root
# u+vi is to subtract _its_ argument from pi/4, which we achieve by
# dividing (1+i) by (u+vi). (with an abs() to arrange that large
# angles don't end up negative)
#
# So we also seek rational approximations to tan(pi/4 - theta/2) as
# well as to tan(theta/2), and between both of those, we hope to find
# a larger set of usefully small primitive Pythagorean triples.

tan_other_halfangle = Spigot("abs((1-x)/(1+x))", {"x":tan_halfangle})

# FIXME. Finish up this other_halfangle business in some way. Known
# problems: sign of tan_other_halfangle can go wrong (in which case,
# surely, also sign of tan_halfangle too?); need to comment the
# reasoning behind what we're even doing here.

def iterator(spig, swap=False):
    for v, u in spig.convergents():
        base, height, hypot = (u*u-v*v, 2*u*v, u*u+v*v)

        if base == 0 or height == 0:
            continue

        d = reduce(gcd, [base, height, hypot])
        base /= d
        height /= d
        hypot /= d

        if swap:
            base, height = height, base

        yield base, height, hypot, u, v

def merge(iter1, iter2):
    val1 = next(iter1)
    val2 = next(iter2)
    while True:
        minhypot = min(val1[2], val2[2])
        if val1[2] == minhypot:
            yield val1
        else:
            yield val2
        while val1[2] <= minhypot:
            val1 = next(iter1)
        while val2[2] <= minhypot:
            val2 = next(iter2)

for base, height, hypot, v, u in merge(
        iterator(tan_halfangle), iterator(tan_other_halfangle, swap=True)):
    output = "%d %d %d" % (base, height, hypot)
    if print_angle:
        angle = Spigot(angle_expr, {"height":Spigot(height),
                                    "base":Spigot(base)})
        output += " angle=" + spig_decimal(angle)
    if print_gaussian:
        output += " source=%d+%di" % (u, v)
    print output

    if limit is not None:
        limit -= 1
        if limit == 0:
            break