[math-fun] New (?) Euclidean algorithm
I've been hacking a new version of the Common Lisp & Maxima "rationalize" function which produces the worst rational approximation to a complex number that falls within epsilon of that complex number. In the case of Maxima, the epsilon is set by the global variable "ratepsilon". The traditional Euclidean algorithm for this process focuses on the number itself (which I call c), and only later worries about the epsilon and the radius of convergence. My algorithm starts with the _circle_ whose center is the given complex number c and whose real radius is r; we can denote this circle by (c,r). Since translation and inversion in the complex plane both take circles to circles (more accurately, "generalized circles" to "generalized circles"), we can perform a Euclidean procedure on the _circle_ (c,r) itself. At each stage we keep track of the transformations (translations & inversions), so that we can recover the original circle, if necessary. We first look at the given circle (c,r). If there is a (Gaussian) integer point within this circle, then we are done, because such a point is rational, so this is our answer. If we compute cr=round(c) (= round(realpart(c))+i*round(imagpart(c))), then this number cr is (one of) the closest integer points to c. If cr is "inside" the circle (c,r), then we are done. Now if r>sqrt(2)/2, then cr _has_ to be inside the circle (c,r). So we now consider the case where cr is not inside (c,r). After translating the origin to cr, we know that this new center is within sqrt(2)/2 from the origin -- i.e., inside the unit circle. We then invert the given circle in the unit circle, which will put the center of that new circle outside the unit circle. Furthermore, the radius of that new circle will grow bigger. We can now compute the closest integral point to the center of that new circle, and so forth. At each of these Euclidean steps, the circle radius grows larger, and must eventually enclose an integral point, at which time we invert the sequence of transformations to compute the final rational answer. This answer will be rational, because the integral point was transformed by inversions and translations by (Gaussian) integers. This discussion is not a _proof_, but shows how the algorithm was conceived. Here's the algorithm in Maxima: euclid(zc,m):= block([q,iqmat], /* zc is input circle in Schwerdtfeger's representation */ /* m accumulates transforms; m starts out as the 2x2 identity matrix */ /* Original circle is always transpose(m^^-1).zc.conjugate(m^^-1) */ q:cround(center(zc)), if is(inside(q,zc)) then m.matrix([q],[1]) else (iqmat:matrix([1,q],[0,1]), euclid(transpose(iqmat.imat).zc.conjugate(iqmat.imat),m.iqmat.imat))); This produces is a column vector matrix([num],[denom]); the result is num/denom. euclid(gcircle(%pi,1.0e-7),ident(2)) produces [ - 122453 ] [ ] [ - 38978 ] = 122453/38978, which is within 10^-7 of pi. (%i13) abs((-122453)/(-38978)-%pi),numer; (%o13) 3.9724384226502707E-8 We use Hans Schwerdtfeger's representation for generalized circles as 2x2 Hermitian matrices. See this Wikipedia article for more details. http://en.wikipedia.org/wiki/Generalised_circle /* We test for "inside" by substituting the point into the equation. */ inside(w,c):=is(equation(c,w)<=0); Although very elegant, there is a problem with Schwerdtfeger's representation for circles: the squared radius of the circle is computed by the difference of two large quantities, so when we approach the limit of machine precision, we may inadvertently end up with a negative number, and hence an _imaginary_ radius! Perhaps this problem can be ameliorated with a slightly different representation.
That is absolutely brilliant! It seems to resemble the process of generating a continued fraction approximation of a complex number: 1. Subtract [z] from z, where [] indicates the closest Gaussian (or Eisenstein, or whatever) integer. 2. Raise z to the power of -1 (i.e. invert in the unit circle and conjugate). 3. Repeat the first two steps ad infinitum. This differs from the ordinary CF algorithm in that the closest integer, rather than the floor, is chosen at each stage (as floor is meaningless in an unordered set such as the complex numbers). Let's see how this version produces a different CF for pi: Ordinary: [3;7,15,1,292,1,1,1,2,1,...] Modified: [3;7,16,-294,3,-4,5,-15,-3,2,...] I decided to search the OEIS for this sequence: http://oeis.org/A133593 -------- What about exp(1)? According to Mathematica, it is: [3; -4, 2, 5, -2, 7, -2, 9, -2, -11, ...] This is only slightly more complicated than the regular continued fraction for e. And exp(i) using Gaussian integers? [1+i; -2+i, 1+3i, -2, 1-5i, -2, 1+7i, -2, 1-9i, -2, 1+11i, ...] Very similar. What about exp(2)? [7; 3, -2, -3, -18, -6, 2, 6, 30, 9, -2, -9, -42, -12, 2, 12, 54, 15, -2, -15, -66, -18, 2, 18, 78, 21, -2, -21, -90, -24, ...] Excluding the obvious sign alteration, the columns are all basic arithmetic progressions. exp(3)? [20; 12, -3, -4, -4, 7, -3, -17, 2, 16, 2, 13, 14, 4, 6, 3, -2, -2, -2, -2, -3, -6, 5, -2, -68, -7, -6, 5, 3, -3, ...] Hmm, that appears completely irregular! Maybe it is still a simple interleaving of arithmetic progressions, but with a very large period. Let's try exp(1/2): [2; -3, 7, -2, -10, 2, 14, -2, -18, 2, 22, -2, -26, 2, 30, ...] Again, very regular. -------- Note that these continued fractions (of irrational numbers) do not permit the appearance of ones, in the same way that an ordinary CF cannot contain zeroes. Phi has the simplest CF; I wonder what has the 'simplest' modified CF of [2; 2, 2, 2, 2, 2, 2, ...]. Well, there's the obvious recurrence x = 2 + 1/x, which leads to the quadratic equation x^2 - 2x - 1 = 0 with positive real root 2 + sqrt(2). -------- Sincerely, Adam P. Goucher
Using the closest integer ("round" = least absolute remainder) isn't new. I don't have my Knuth V.I handy, but I think that using "round" instead of "floor" is several hundred years old. Rounding is also necessary for convergence in the complex plane -- e.g., for Gaussian integer gcd's. See this discussion: http://home.pipeline.com/~hbaker1/Gaussian.html My new (?) idea is focusing on the circle itself, rather than the point, BTW, the inverse of a circle (C,r) is the circle (C',r)/(CC'-r^2), and as you can imagine, things get very hairy when r^2=CC' -- i.e., when the circle goes directly through the origin. If the circle to be inverted encloses the origin, then its inverse will be "inside out" -- i.e., it will enclose infinity (oo), but exclude the origin. For small r, one could approximate 1/(CC'-r^2) with a Taylor expansion: 1/(CC'-r^2) ~ 1/(CC')(1+r^2/(CC')+...) So if we choose a smaller radius to take account of the approximation for the center point, then we might be able to use this Taylor expansion until r^2>1/2, which is all we will ever need. ---- Note that the sequence of partial quotients is slightly different from the classical algorithm. In particular, my algorithm might be "too" good, in the sense that it might not return the _worst_ approximation that fits within the original circle, but one that is "one click" better. If you are really looking for the worst approximation -- i.e., the one with the smallest denominator, my algorithm may not produce it. ---- It should be possible to get an estimate on the roundoff error involved in these procedures by carrying out 2 parallel computations: the traditional one using floating point numbers, and mine using rational numbers; we just have to agree on the partial quotients. Since the Mobius transform is rational and hence exact, if we ever arrive at a situation where the point to be approximated no longer lies "inside" its own circle, then we have clear evidence of roundoff error, and we must stop -- this is the best approximation we're going to be able to get within the limits of machine precision. At 01:04 PM 7/8/2012, Adam P. Goucher wrote:
That is absolutely brilliant! It seems to resemble the process of generating a continued fraction approximation of a complex number:
1. Subtract [z] from z, where [] indicates the closest Gaussian (or Eisenstein, or whatever) integer.
2. Raise z to the power of -1 (i.e. invert in the unit circle and conjugate).
3. Repeat the first two steps ad infinitum.
This differs from the ordinary CF algorithm in that the closest integer, rather than the floor, is chosen at each stage (as floor is meaningless in an unordered set such as the complex numbers). Let's see how this version produces a different CF for pi:
Ordinary: [3;7,15,1,292,1,1,1,2,1,...] Modified: [3;7,16,-294,3,-4,5,-15,-3,2,...]
I decided to search the OEIS for this sequence:
--------
What about exp(1)? According to Mathematica, it is:
[3; -4, 2, 5, -2, 7, -2, 9, -2, -11, ...]
This is only slightly more complicated than the regular continued fraction for e. And exp(i) using Gaussian integers?
[1+i; -2+i, 1+3i, -2, 1-5i, -2, 1+7i, -2, 1-9i, -2, 1+11i, ...]
Very similar. What about exp(2)?
[7; 3, -2, -3, -18, -6, 2, 6, 30, 9, -2, -9, -42, -12, 2, 12, 54, 15, -2, -15, -66, -18, 2, 18, 78, 21, -2, -21, -90, -24, ...]
Excluding the obvious sign alteration, the columns are all basic arithmetic progressions. exp(3)?
[20; 12, -3, -4, -4, 7, -3, -17, 2, 16, 2, 13, 14, 4, 6, 3, -2, -2, -2, -2, -3, -6, 5, -2, -68, -7, -6, 5, 3, -3, ...]
Hmm, that appears completely irregular! Maybe it is still a simple interleaving of arithmetic progressions, but with a very large period. Let's try exp(1/2):
[2; -3, 7, -2, -10, 2, 14, -2, -18, 2, 22, -2, -26, 2, 30, ...]
Again, very regular. --------
Note that these continued fractions (of irrational numbers) do not permit the appearance of ones, in the same way that an ordinary CF cannot contain zeroes. Phi has the simplest CF; I wonder what has the 'simplest' modified CF of [2; 2, 2, 2, 2, 2, 2, ...].
Well, there's the obvious recurrence x = 2 + 1/x, which leads to the quadratic equation x^2 - 2x - 1 = 0 with positive real root 2 + sqrt(2). --------
Sincerely,
Adam P. Goucher
It seems like this same algorithm ought to work on a regular hyperbolic tiling, too: every point lies in some regular polygon labeled by an element of the Coxeter group. We apply the inverse of that element to the point and get a point within the polygon centered at the origin, then do inversion through the circumcircle (using the hyperbolic metric) and repeat. On Sun, Jul 8, 2012 at 10:23 AM, Henry Baker <hbaker1@pipeline.com> wrote:
I've been hacking a new version of the Common Lisp & Maxima "rationalize" function which produces the worst rational approximation to a complex number that falls within epsilon of that complex number. In the case of Maxima, the epsilon is set by the global variable "ratepsilon".
The traditional Euclidean algorithm for this process focuses on the number itself (which I call c), and only later worries about the epsilon and the radius of convergence.
My algorithm starts with the _circle_ whose center is the given complex number c and whose real radius is r; we can denote this circle by (c,r).
Since translation and inversion in the complex plane both take circles to circles (more accurately, "generalized circles" to "generalized circles"), we can perform a Euclidean procedure on the _circle_ (c,r) itself. At each stage we keep track of the transformations (translations & inversions), so that we can recover the original circle, if necessary.
We first look at the given circle (c,r). If there is a (Gaussian) integer point within this circle, then we are done, because such a point is rational, so this is our answer. If we compute cr=round(c) (= round(realpart(c))+i*round(imagpart(c))), then this number cr is (one of) the closest integer points to c. If cr is "inside" the circle (c,r), then we are done. Now if r>sqrt(2)/2, then cr _has_ to be inside the circle (c,r).
So we now consider the case where cr is not inside (c,r). After translating the origin to cr, we know that this new center is within sqrt(2)/2 from the origin -- i.e., inside the unit circle. We then invert the given circle in the unit circle, which will put the center of that new circle outside the unit circle. Furthermore, the radius of that new circle will grow bigger. We can now compute the closest integral point to the center of that new circle, and so forth.
At each of these Euclidean steps, the circle radius grows larger, and must eventually enclose an integral point, at which time we invert the sequence of transformations to compute the final rational answer. This answer will be rational, because the integral point was transformed by inversions and translations by (Gaussian) integers.
This discussion is not a _proof_, but shows how the algorithm was conceived.
Here's the algorithm in Maxima:
euclid(zc,m):= block([q,iqmat], /* zc is input circle in Schwerdtfeger's representation */ /* m accumulates transforms; m starts out as the 2x2 identity matrix */ /* Original circle is always transpose(m^^-1).zc.conjugate(m^^-1) */ q:cround(center(zc)), if is(inside(q,zc)) then m.matrix([q],[1]) else (iqmat:matrix([1,q],[0,1]), euclid(transpose(iqmat.imat).zc.conjugate(iqmat.imat),m.iqmat.imat)));
This produces is a column vector matrix([num],[denom]); the result is num/denom.
euclid(gcircle(%pi,1.0e-7),ident(2)) produces
[ - 122453 ] [ ] [ - 38978 ] =
122453/38978, which is within 10^-7 of pi.
(%i13) abs((-122453)/(-38978)-%pi),numer; (%o13) 3.9724384226502707E-8
We use Hans Schwerdtfeger's representation for generalized circles as 2x2 Hermitian matrices. See this Wikipedia article for more details.
http://en.wikipedia.org/wiki/Generalised_circle
/* We test for "inside" by substituting the point into the equation. */
inside(w,c):=is(equation(c,w)<=0);
Although very elegant, there is a problem with Schwerdtfeger's representation for circles: the squared radius of the circle is computed by the difference of two large quantities, so when we approach the limit of machine precision, we may inadvertently end up with a negative number, and hence an _imaginary_ radius!
Perhaps this problem can be ameliorated with a slightly different representation.
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
-- Mike Stay - metaweta@gmail.com http://www.cs.auckland.ac.nz/~mike http://reperiendi.wordpress.com
On Sun, Jul 8, 2012 at 4:23 PM, Mike Stay <metaweta@gmail.com> wrote:
It seems like this same algorithm ought to work on a regular hyperbolic tiling, too: every point lies in some regular polygon labeled by an element of the Coxeter group.
Sorry, the 2n triangles in each regular n-gon are labeled by elements of the Coxeter group.
We apply the inverse of that element to the point and get a point within the polygon centered at the origin, then do inversion through the circumcircle (using the hyperbolic metric) and repeat.
On Sun, Jul 8, 2012 at 10:23 AM, Henry Baker <hbaker1@pipeline.com> wrote:
I've been hacking a new version of the Common Lisp & Maxima "rationalize" function which produces the worst rational approximation to a complex number that falls within epsilon of that complex number. In the case of Maxima, the epsilon is set by the global variable "ratepsilon".
The traditional Euclidean algorithm for this process focuses on the number itself (which I call c), and only later worries about the epsilon and the radius of convergence.
My algorithm starts with the _circle_ whose center is the given complex number c and whose real radius is r; we can denote this circle by (c,r).
Since translation and inversion in the complex plane both take circles to circles (more accurately, "generalized circles" to "generalized circles"), we can perform a Euclidean procedure on the _circle_ (c,r) itself. At each stage we keep track of the transformations (translations & inversions), so that we can recover the original circle, if necessary.
We first look at the given circle (c,r). If there is a (Gaussian) integer point within this circle, then we are done, because such a point is rational, so this is our answer. If we compute cr=round(c) (= round(realpart(c))+i*round(imagpart(c))), then this number cr is (one of) the closest integer points to c. If cr is "inside" the circle (c,r), then we are done. Now if r>sqrt(2)/2, then cr _has_ to be inside the circle (c,r).
So we now consider the case where cr is not inside (c,r). After translating the origin to cr, we know that this new center is within sqrt(2)/2 from the origin -- i.e., inside the unit circle. We then invert the given circle in the unit circle, which will put the center of that new circle outside the unit circle. Furthermore, the radius of that new circle will grow bigger. We can now compute the closest integral point to the center of that new circle, and so forth.
At each of these Euclidean steps, the circle radius grows larger, and must eventually enclose an integral point, at which time we invert the sequence of transformations to compute the final rational answer. This answer will be rational, because the integral point was transformed by inversions and translations by (Gaussian) integers.
This discussion is not a _proof_, but shows how the algorithm was conceived.
Here's the algorithm in Maxima:
euclid(zc,m):= block([q,iqmat], /* zc is input circle in Schwerdtfeger's representation */ /* m accumulates transforms; m starts out as the 2x2 identity matrix */ /* Original circle is always transpose(m^^-1).zc.conjugate(m^^-1) */ q:cround(center(zc)), if is(inside(q,zc)) then m.matrix([q],[1]) else (iqmat:matrix([1,q],[0,1]), euclid(transpose(iqmat.imat).zc.conjugate(iqmat.imat),m.iqmat.imat)));
This produces is a column vector matrix([num],[denom]); the result is num/denom.
euclid(gcircle(%pi,1.0e-7),ident(2)) produces
[ - 122453 ] [ ] [ - 38978 ] =
122453/38978, which is within 10^-7 of pi.
(%i13) abs((-122453)/(-38978)-%pi),numer; (%o13) 3.9724384226502707E-8
We use Hans Schwerdtfeger's representation for generalized circles as 2x2 Hermitian matrices. See this Wikipedia article for more details.
http://en.wikipedia.org/wiki/Generalised_circle
/* We test for "inside" by substituting the point into the equation. */
inside(w,c):=is(equation(c,w)<=0);
Although very elegant, there is a problem with Schwerdtfeger's representation for circles: the squared radius of the circle is computed by the difference of two large quantities, so when we approach the limit of machine precision, we may inadvertently end up with a negative number, and hence an _imaginary_ radius!
Perhaps this problem can be ameliorated with a slightly different representation.
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
-- Mike Stay - metaweta@gmail.com http://www.cs.auckland.ac.nz/~mike http://reperiendi.wordpress.com
-- Mike Stay - metaweta@gmail.com http://www.cs.auckland.ac.nz/~mike http://reperiendi.wordpress.com
participants (3)
-
Adam P. Goucher -
Henry Baker -
Mike Stay