From: rusin@washington.math.niu.edu (Dave Rusin) Newsgroups: sci.math Subject: Re: Something like continued fractions? Date: 24 Oct 1995 21:58:03 GMT In article <46jb37$h6p@news.ox.ac.uk>, Thomas Womack wrote: >Is there a way, like the continued fraction algorithm for finding numbers >such that ax+b is close to zero for irrational x, of finding numbers such >that ax^2+bx+c is close to zero (for a,b,c integers)? > Given how good the continued fraction approximations are, I'd expect it >to be possible to find very good quadratic approximations with small a,b,c; >OK, what I'm really looking for is a short expression which is approximately >equal to pi, and closer to pi than 355/113 is. In an incredible display of chance events, this post arrived just moments after I was trying to remember the same thing myself (hoping to make test questions for calculus students who would calculate the root to a quadratic or cubic with Newton's method and jump to the assumption that the root was really pi). One way to describe what you want is to have smallish numbers a, b, c making M*(ax^2+bx+c) also smallish, when M is really big (and x is fixed). The trick is to realize that this is the same as finding integers a, b, and c so that the vector a(1,0,0,Mx^2)+b(0,1,0,Mx)+c(0,0,1,M) has all its coordinates as small as possible (with x given and M fixed but large and otherwise arbitrary). So what you need is an algorithm which finds elements in lattices (subgroups of R^n) which are small in, say, the euclidean norm. There is, after all, a minimum to the non-zero distances between lattice points and the origin if, as in this example, the lattice is discrete. Well, you're in luck. There is an algorithm which can rapidly determine, if not the smallest element in a lattice, at least a fairly small one -- indeed it can find a basis of smallish lattice elements. It's usually called the LLL routine, (Lenstra Lenstra and Lovasz, Math Ann 261 (1982) pp. 515-534) and it's encoded into maple and other CA packages. In maple you just do readlib(lattice); x:=evalf(Pi); M:=1000000000; lattice([[1,0,0,M*x^2], [0,1,0,M*x], [0,0,1,M]]); and voil`a, you have three small vectors in this lattice. One of them, for example, is [780, -2255, -614, 350.], which must of course correspond to 780X^2-2255*X-614, a polynomial with a root which is deucedly close to Pi. I have also been shown this polynomial, which I assume was found by similar means: 242 x^4 + 113 x - 23928. Its root is quite close to Pi. I know there's a lot of material relating the size of the polynomial and its degree to the separation of its root from Pi (or any other fixed transcendental) just like the 1/q^2 result for continued fractions, but I don't have references for it to hand. Perhaps these examples will suffice for your purposes. dave ============================================================================== From: bobs@mathworks.com (Bob Silverman) Newsgroups: sci.math Subject: Re: More of the pseudo-continued fractions Date: 25 Oct 1995 08:36:13 -0400 In article <46l45b$a8e@news.ox.ac.uk>, Thomas Womack wrote: >I'm afraid my last post wasn't very clear. What I'm looking for is a way, >short of exhaustive search, of finding near-identities like the following : > >5*e^2+14*e-75 = 0.001 > You want to look up integer relation finding algorithms. (such as one by Forcade). David Bailey wrote some articles on these in Math. Comp. in the late 80's early 90's. -- Bob Silverman The MathWorks Inc. 24 Prime Park Way Natick, MA ============================================================================== From: nsinger@eos.hitc.com (Nick Singer) Newsgroups: sci.math Subject: Re: More of the pseudo-continued fractions Date: 25 Oct 1995 19:37:19 GMT In article <46l45b$a8e@news.ox.ac.uk>, mert0236@sable.ox.ac.uk (Thomas Womack) says: >What I'm looking for is a way, short of exhaustive search, >of finding near-identities like 5*e^2+14*e-75 = 0.001. > >I can think of one way - producing a large array of the fractional parts of >n*e, sorting it, and then looking through it with binary search for the >fractional part of m*e^2. Obviously (since e is transcendental), you don't get >an exact hit, but if you find the closest number located you can use this to >show that m*e^2-n*e has a small fractional part. Alternatively, look for >1-the fractional part of m*e^2 to find solutions to m*e^2+n*e is about 0. > Recall e's continued fraction expansion: e=[2,1,2,1,1,4,1,1,6,1,1,8,...]. Now recall (e.g. from reading G. Chrystal, Textbook of Algebra, chap. XXXIII; or Hardy and Wright, Intro. to Theory of Numbers) that every recurring continued fraction expansion is a "simple quadratic surd" [and vice versa], i.e. of the form (P + sqrt(Q))/R, where P, Q, and R are integers. But these simple quadratic surds are just the sort of things that are roots of quadratics with integer coefficients. So try to make a fit. Try a=[2,1,2,1,2,...] = 1+sqrt(3), which satisfies a^2 -2a -2=0. If you substitute e for a, you get -0.0475 on the right hand side. You could try b=[2,1,2,1,1,4,1,1,2,1,1,4,...] =(7+sqrt(15))/4, which satisfies 8b^2-28b+17=0. When you stick in e for b, the right hand side is 0.00056. Or you might try c=[2,1,2,1,1,4,1,1,4,1,1,4,1,1,...]=(35-sqrt(26))/11, which satisfies 11c^2-70c+109=0. Substituting e for c gives a right hand side of -0.00011. The example you gave corresponds to d=(-7+sqrt(424))/5 =[2,1,2,1,1,4,1,1,2,1,7,1,1,13,...], which is not such a great fit. Hope this is what you were looking for. ============================================================================== From: Thomas Womack Newsgroups: sci.math Subject: Summary of pseudo-continued-fractions Date: Mon, 30 Oct 1995 13:44:13 +0000 Approximating numbers quadratically and higher ---------------------------------------------- I had quite a lot of replies to my search for information on this, most of which proposed one of the four following methods. 1) LLL's lattice-reduction routine This routine looks for points close to the origin in a lattice. I used it very much as a 'black box', although it is described in Matt Ann 261 pp 515-534 and in Cohen's 'Computational Algebraic Number Theory'. I don't think, however, that it was designed for the task I'm setting it. This is because, when I asked it to find the close points in the lattice [1,0,0,0,10^9*Pi^3],[0,1,0,0,10^9*Pi^2],[0,0,1,0,10^9*Pi], [0,0,0,1,10^9], it returned three vectors with a zero in the last component and one with a one there; the polynomials derived from these had roots fairly close to Pi, but not staggeringly close given the size of their coefficients. 2. PSOS algorithm This is described in a paper by D Bailey and H Ferguson in Math Comp 53, pp 649-656. What it is used for there is to produce empirical tests for the transcendence of certain numbers, by trying to find a polynomial which the number satisfies. The paper gives an algorithm for it, which I have not yet been able to implement, and no discussion as to why it works. There are no indications in the paper as to whether it produces good intermediate results (OK, there may be, but I haven't found them yet). 3. Continuing continued fractions periodically A periodic continued fraction (eg [1,7,9,11,4,3,6,3,6,3,6,...]) represents the solution to a quadratic equation. Several people have suggested that I could get good quadratic approximations to a number (the method does not extend further) by taking a good rational approximation's continued fraction and extending it periodically. This works very well for e, where the continued fraction is virtually periodic anyway; for pi, where there is no such obvious pattern, it is less successful. 4. Variants on Euclid's GCD algorithm This seems the most promising approach. You start off with an array of the powers of the number you're trying to approximate, then perform the following algorithm : 1. find the smallest entry 2. subtract integer multiples of it from the other entries to make the absolutes value of the other entries as small as possible 3. record what you've done algebraically 4. repeat until bored (or until your computer runs out of precision, which is usually earlier). 5. Brute-force [I haven't implemented 1..4 yet] These go up to (coefficient of x^2) + (coefficient of x) = 20000 Each set took about 2.5 minutes to calculate on an Alpha mini Pi^2 Pi Result 0 1 3.1415926535897931000000000 1 1 13.0111970546791510000000000 0 7 21.9911485751285520000000000 1 8 35.0023456298077060000000000 12 4 131.0016234274314700000000000 4 39 162.0005310943593600000000000 15 35 257.9998088919830900000000000 49 31 580.9999879146621400000000000 34 674 2452.9999981565588000000000000 1146 1013 14493.0000017348650000000000000 200 1971 8166.0000004433541000000000000 445 1900 10361.0000003053720000000000000 690 1829 12556.0000001673890000000000000 935 1758 14751.0000000294060000000000000 8417 251 83861.0000000201540000000000000 6201 3606 72529.9999999998980000000000000 e^2 e Result 0 1 2.7182818284590451000000000 1 1 10.1073379273896950000000000 2 3 22.9329576832384350000000000 1 5 20.9804652412258750000000000 4 2 34.9927880526406910000000000 11 1 83.9978989166961870000000000 5 14 75.0012260930798790000000000 29 1 217.0009086974479000000000000 16 15 158.9991250097760700000000000 21 29 234.0003511028559700000000000 45 16 376.0000337072239600000000000 92 159 1111.9999718266079000000000000 137 175 1488.0000055338319000000000000 283 153 2506.9999957516075000000000000 74 553 2050.0000024587198000000000000 420 328 3995.0000012854398000000000000 766 103 5940.0000001121589000000000000 442 2190 9219.0000000526561000000000000 118 4277 12497.9999999931530000000000000 4283 3636 41530.9999999970610000000000000 8448 2995 70564.0000000009600000000000000 w = 1.234567891011121314151617181920 ... (known transcendental) w^2 w Result 0 1 1.2345678910111213000000000 1 2 3.9932936595378905000000000 11 1 18.0003045436832490000000000 0 81 99.9999991719008250000000000 6040 43 9258.9999995079925000000000000 5959 224 9358.9999997022369000000000000 5878 405 9458.9999998964813000000000000 5797 586 9559.0000000907276000000000000 5473 1391 10059.0000000396120000000000000 5149 2196 10558.9999999884930000000000000 1342 11594 16359.0000000089400000000000000 7868 10389 24818.0000000076540000000000000 I've also got some cubic and quartic ones (again, brute force searched) for e Cubic approximations ( ae^3+be^2+ce=d) a,b,c<1000 a b c d 0 0 1 20.0855369231876680000000000 0 0 11 220.9409061550643400000000000 0 0 12 241.0264430782520200000000000 0 0 35 702.9937923115684300000000000 0 0 152 3053.0016123245255000000000000 0 0 491 9861.9986292851445000000000000 0 0 643 12915.0002416096700000000000000 0 1 506 10165.9999649614180000000000000 0 21 151 3089.9999937989778000000000000 0 104 249 5284.0000040334698000000000000 0 125 400 8373.9999978324468000000000000 0 229 649 13658.0000018659180000000000000 0 546 220 5903.0000014399257000000000000 0 671 620 14276.9999992723730000000000000 1 312 41 1679.0000004288472000000000000 1 983 661 15955.9999997012200000000000000 3 73 332 6889.0000002726083000000000000 6 463 235 6023.0000001192238000000000000 8 224 526 11232.9999999629840000000000000 9 853 138 5156.9999999658394000000000000 35 932 958 22034.0000000001890000000000000 Quartic approximations (ae^4+be^3+ce^2+de = f) a b c d f 0 0 0 1 54.5981500331442360000000000 0 0 0 2 109.1963000662884700000000000 0 0 0 5 272.9907501657211800000000000 0 0 1 45 2477.0022884146783000000000000 0 0 4 78 4338.9978502780014000000000000 0 0 5 16 973.9980851462461300000000000 0 0 6 61 3451.0003735609243000000000000 0 0 25 75 4596.9996755655093000000000000 0 0 26 13 1231.9999104337544000000000000 0 1 27 15 1364.0000292516895000000000000 0 2 1 66 3629.0000027676256000000000000 0 11 23 37 2512.0000005727024000000000000 0 42 91 82 6418.9999995231847000000000000 0 59 29 32 2489.9999997121413000000000000 0 70 52 69 5002.0000002848446000000000000 3 86 35 72 4890.0000002422230000000000000 10 45 99 44 4587.0000001238886000000000000 13 61 82 47 4475.0000000812688000000000000 16 77 65 50 4363.0000000386472000000000000 19 93 48 53 4250.9999999960264000000000000 82 70 29 34 3235.0000000037930000000000000 I'll post some results got with more elaborate method when I've learnt enough Maple programming to write the code for them. Any comments? Tom Womack ============================================================================== Date: Tue, 31 Oct 1995 15:34:33 +0500 From: nsinger@eos.hitc.com (Nick Singer) Subject: Re: More of the pseudo-continued fractions To: rusin@math.niu.edu (Dave Rusin) Thanks for sending me your earlier posting. There is no question that the approach you outline is more general and more powerful. The continued-fraction- expansion-fitting made sense for quadratic approximations to e, which has a predictable C.F. expansion, but is useless for higher-order polynomials, and is of questionable utility for other real numbers. I don't "have" the LLL algorithm, but I created something 15 yrs ago that does the same thing (i.e. lattice reduction) when I implemented the spectral test from Knuth vol. 2. You say that my method has the advantage of being hand-calculable, but so is the lattice reduction technique, especially if you are looking only at quadratics. I blew a few hours cheerfully (mindlessly) carrying out the algorithm by "hand"--just me and MS Excel, sans macros. My method (and LLL's ?) does (basically) integer-multiple Gram-Schmidt orthogonalization, reducing the length of (at least) one of the basis vectors at each iteration, until it can't anymore. Since we include the error as the last component of the basis vectors, the (final) results depend strongly on M. The "best" answer for one value of M may not even show up at all in the reduction process for a larger M. In any case, we seem to get |x-x0| on the order of 1/M. Here are some of my hand-cranked solutions (one of which seems "better" than the example of 780*x^2-2255*x-614 for pi.) pi: M x^2 x 1 |x-pi| 2,648,617: 49 -82 -226 8*10^-8 10^10 245 -71 -2195 9*10^-11 10^14 6201 3606 -72530 2*10^-15 e: M x^2 x 1 |x-e| 10^14 11345 -29775 -2892 5*10^-15 10^11 324 -2087 3279 2*10^-10 10^11 442 2190 -9219 1*10^-11 10^11 3929 -9195 -4037 1*10^-12 10^7 54 -181 93 2*10^-7 I also tried a quartic fit for pi; I came up with 9x^4+65x^3-336x^2+262x-399. The error |x-pi| seems to be about 1*10^-14. >Nice idea. I already wrote the poster with the next-to-last quadratic >you posted, and made the same comment you did about his original one. >I found it using the lattice reduction technique I posted in a previous >article, which admits a generalization to cubics etc., but your method >has the advantage of being hand-calculable. > >dave ============================================================================== From: rusin@math.niu.edu (Dave Rusin) Newsgroups: sci.math.num-analysis,sci.math.research,sci.math Subject: Re: anybody see this formula? [algebraic approximants to Pi] Date: 11 Jun 1998 06:43:49 GMT In article <6l4n4r$lr8@mcmail.CIS.McMaster.CA>, Zdislav V. Kovarik wrote: >Here is something more spectacular (and I'd like to see how this example >was discovered): > >Find the positive root of > > p(x) = 242 * x^4 + 113 * x - 23928 I cannot say how this example was first discovered, but such examples are easy to find. Is the above more remarkable than p(x) = 5 * x^4 + 5 * x^3 + 69 * x^2 - 389 * x - 101 of the same degree, with much smaller coefficients, whose root is even closer to Pi? One can find increasingly good rational approximants to Pi whose minimal polynomials have fairly small coefficients; see e.g. http://www.math.niu.edu/~rusin/known-math/index/11YXX.html dave [ moderator's note: Maple even has a command for this, called minpoly. Some other nice ones for Pi: 2+29 X-22 X^2+4 X^3 has root 3.141592 9956 -72530+3606 X+6201 X^2 has root 3.14159265358979 54182 46-104 X+2 X^2+2 X^3+79 X^4+68 X^5-64 X^6+11 X^7 has root 3.14159265358979323 355 -2223+285 X+983 X^2-4373 X^3+1306 X^4 has root 3.1415926535897932384 958 ]