Re: [math-fun] pythagorean rationals & approximations by them
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.
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.
It's interesting, though very easy to prove, that the set G of rational points QxQ \int S^1 on the unit circle forms a group under ordinary complex multiplication. (This has been discussed here many years ago, but:) Whereas QxQ seems, like ZxZ, to have a squarish structure, the subset G of QxQ is rotationally homogeneous: Given any two elements of G, there is (of course) some group element that carries one onto the other. Also interesting is its group structure. I carelessly used to think G was isomorphic to Q/Z, which is easily seen to be the direct sum of the groups H_p := {k/p^n : (k,n) in Z^2}/Z for each prime p. But instead, G is unexpectedly isomorphic to the direct sum of a countable number of infinite cyclic groups and one copy of Z/4Z. --Dan
An interesting question is how best to represent G's elements on a computer. The usual representations I can think of blow up extremely quickly. Perhaps you just bite the bullet & factor the number into its primes, and then keep a list of non-zero prime exponents, while treating 2 specially. At 12:24 PM 8/15/2013, Dan Asimov wrote:
It's interesting, though very easy to prove, that the set G of rational points QxQ \int S^1 on the unit circle forms a group under ordinary complex multiplication.
(This has been discussed here many years ago, but:) Whereas QxQ seems, like ZxZ, to have a squarish structure, the subset G of QxQ is rotationally homogeneous: Given any two elements of G, there is (of course) some group element that carries one onto the other.
Also interesting is its group structure. I carelessly used to think G was isomorphic to Q/Z, which is easily seen to be the direct sum of the groups H_p := {k/p^n : (k,n) in Z^2}/Z for each prime p.
But instead, G is unexpectedly isomorphic to the direct sum of a countable number of infinite cyclic groups and one copy of Z/4Z.
--Dan
On 15/08/2013 16:48, Henry Baker wrote:
All perfect squares of Gaussian integers are Pythagorean triples; are all Pythagorean triples squares of Gaussian integers?
If I'm understanding the question right (maybe not) the answer is "yes, of course, up to multiplying everything by an integer". The general Pythagorean triple is k.(m^2-n^2), k.2mn, k.(m^2+n^2) where k can be any integer. If we ignore k, then m^2-n^2 + 2mn.i = (m+in)^2 so we have the square of a Gaussian integer. If we don't ignore k, we can fail to have the square of a Gaussian integer; e.g., take any primitive triple, e.g. m=2,n=1 giving 3,4,5, and then let k=7 giving 21,28,35. 21+28i isn't the square of a Gaussian integer because if it were we'd have 7|w|^2=|z|^2, hence 7=|z|^2/|w|^2, but both |w|^2 and |z|^2 have to have an even number of factors of 7 because 7 = -1 (mod 4). But this sort of thing -- having a common factor with an odd power of a prime that's -1 mod 4 -- is the only way it can happen. -- g
OK, so I guess the most general statement is that if all common prime factors of a,b,c which are =3 (mod 4) have even exponents, then a+bi is the square of a Gaussian integer. In particular, if (a,b)=1, then a+bi must be the square of a Gaussian integer. Correct? Since I'm primarily trying to come up with _small_ representations, this constraint isn't a problem, because additional factors of such primes would only help with the representation of numbers on the axes, and we already have excellent representations of numbers on the axes! At 03:14 PM 8/15/2013, Gareth McCaughan wrote:
On 15/08/2013 16:48, Henry Baker wrote:
All perfect squares of Gaussian integers are Pythagorean triples; are all Pythagorean triples squares of Gaussian integers?
If I'm understanding the question right (maybe not) the answer is "yes, of course, up to multiplying everything by an integer".
The general Pythagorean triple is k.(m^2-n^2), k.2mn, k.(m^2+n^2) where k can be any integer. If we ignore k, then
m^2-n^2 + 2mn.i = (m+in)^2
so we have the square of a Gaussian integer.
If we don't ignore k, we can fail to have the square of a Gaussian integer; e.g., take any primitive triple, e.g. m=2,n=1 giving 3,4,5, and then let k=7 giving 21,28,35. 21+28i isn't the square of a Gaussian integer because if it were we'd have 7|w|^2=|z|^2, hence 7=|z|^2/|w|^2, but both |w|^2 and |z|^2 have to have an even number of factors of 7 because 7 = -1 (mod 4). But this sort of thing -- having a common factor with an odd power of a prime that's -1 mod 4 -- is the only way it can happen.
participants (3)
-
Dan Asimov -
Gareth McCaughan -
Henry Baker