I ran some more experiments; below are my conclusions. All perfect squares of Gaussian integers are Pythagorean triples; are all Pythagorean triples squares of Gaussian integers? So, if we are interested in 'converting' a complex floating point number z=x+iy into a Pythagorean triple (px+ipy)/ph, here's one way to do it: 1. Normalize the floating point complex number z: z <- z/|z| 2. Multiply by a large radius R>0: z <- R*z (# bits of R ~= final size) 3. Take the _complex square root_ of z: z <- sqrt(z) 3. Round the real & imaginary parts: z <- round(realpart(z))+i*round(imagpart(z)) 4. Square this Gaussian integer: z <- z*z (|z| is now an integer) 5. Result is the complex rational z/|z|. The result of step 4 is a perfect square of a Gaussian integer, hence Pythagorean. The problem is that the minimum # bits required to simply _represent_ N distinguishable rational points on the unit circle leads to substantial _angular errors_ in these points. For example, we need 30 bits in the denominator (|z|) to distinguish 100,000 rational points around the unit circle. However, the positional error can be as large as 69% of the angle 2*pi/100000. If we add more bits, we can reduce this angular error. Going to 36 bits gets the angular error under 10%, going to 44 bits gets the angular error under 1%, and going to 58 bits gets the angular error under 1%. Clearly, minimizing the angular error is very expensive in the number of bits. Below is my test program: (defvar *fudge* 0) ; number of additional bits to add (setq *fudge* 0) (defun ratcis (phi &optional (n (ash 1 15))) "Return a Pythagorean rational approximation to cis(phi)." (declare (type double-float phi) (type unsigned-byte n)) (let* ((phid2 (/ phi 2.0d0)) ; need phi/2. (r (round (* (/ n pi) (expt 2.0d0 *fudge*)))) (c (cos phid2)) (s (sin phid2)) (x (round (* r c))) (y (round (* r s))) (z (complex x y)) (z2 (* z z)) (nphi (phase z2)) (z2norm (* z2 (conjugate z2))) (iabsz2 (isqrt z2norm))) (declare (type double-float phid2 c s) (type single-float nphi) (type unsigned-byte r z2norm iabsz2) (type integer x y)) ; (print `(ratcis phi ,phi phid2 ,phid2 r ,r x ,x y ,y z ,z z2 ,z2)) (values z2 (= z2norm (expt iabsz2 2)) (integer-length x) (integer-length y) (integer-length iabsz2) (/ (abs (- phi nphi)) (/ (* 2 pi) 100000))))) (defun tratcis (n &aux (maxerr 0.0d0) (maxerri 0) (maxlen 0)) (declare (type unsigned-byte n maxlen maxerri) (type double-float maxerr)) (dotimes (i (ash n -2)) (unless (zerop i) (multiple-value-bind (res pyth-p xl yl hl err) (ratcis (/ (* 2 pi i) n) n) (declare (type unsigned-byte hl) (ignore res xl yl) (type double-float err)) (unless pyth-p (print `(tratcis bad pyth-p i ,i))) (when (> hl maxlen) (setq maxlen hl)) (when (> err maxerr) (setq maxerr err maxerri i))))) (print `(tratcis maxlen ,maxlen maxerr ,maxerr "1/maxerr" ,(floor (/ maxerr)) maxerri ,maxerri))) At 04:48 PM 8/14/2013, Henry Baker wrote:
If I understand you correctly, you're recommending the following sequence of operations:
To find a Pythagorean approximation for floating point (x+iy)/|x+iy|.
Denote atan2(y,x) by phi.
So tan(phi)=y/x, cos(phi)=x/sqrt(x^2+y^2), sin(phi)=y/sqrt(x^2+y^2).
What we really want is a rational approximation to tan(phi/2).
tan(phi/2) = sin(phi)/(1+cos(phi)) = y/(sqrt(x^2+y^2)+x), or = (1-cos(phi))/sin(phi) = (sqrt(x^2+y^2)-x)/y,
depending upon the sign of x, so we get the most accurate answer.
We then use continued fractions to produce a '1/2-length' rational approximation to t=tan(phi/2).
The Pythagorean result is ((t^2-1)+2it)/(t^2+1), which is 'full length'. --- I ran this procedure on the test case:
cos(phi)+i*sin(phi), where phi=2*%pi/100000.
The best my previous procedure was able to achieve was 506,605,917 + 31,831*i, which is _not_ Pythagorean.
However, the above procedure produces (using one of the rational continued fraction approximations for tan(phi/2)):
506,606,280 + 31,831*i, which _is_ Pythagorean, with hypotenuse 506,606,281 !!
So maybe there is hope for a 'Pythagorean number system' after all !
I now need to check all 100,000 points -- phi=2pik/100000, k=0..99999 -- to make sure that they can all be represented with reasonable size denominators.
At 01:46 PM 8/14/2013, Warren D Smith wrote:
Henry Baker wanted to approximate reals x by "pythagorean rationals" x =approx= a/b with squareroot(a^2+b^2)=integer.
Obviously, yes, any such pythagorean (a,b) does have a continued fraction (namely just use the CF for a/b); and as Asimov remarked the set of rational points on the unit circle is dense, hence any x can be approximated arbitrarily closely by a pythagorean rational.
The pythagorean rationals form a group in the sense that any two rational points on the unit circle can be multiplied (viewed as complex numbers). That is, if v+iw has v,w rational with v^2+w^1, then view that (v,w) as the rational (v*LCM(vdenominator,wdenominator)) / (w*LCM(vdenominator,wdenominator)).
Euclid found all pythagorean triples (a,b,c) namely a = (m-n)*(m+n), b=2*m*n, c=m^2+n^2 which enables easily finding arbitrarily many arbitrarily close pythagorean rational approximations to x by this method: 1. use CFs to find large a,b so x =approx= a/b 2. solve for m,n satisfying above two simultaneous quadratic equations (as real numbers) 3. round off m,n to nearest integers 4. find resulting pythagorean a,b.