1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
#!/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
|
