September 2010 - mathfun@mailman.xmission.com

[math-fun] A real, live, practical use for Lambert's W-function
by Bill Gosper 20 Sep '10

20 Sep '10

Kepler's equation (c1) 2*\a = alpha-epsilon*sin(alpha) (d1) 2 A = alpha - sin(alpha) epsilon was solved by Lagrange with an infinite series of Bessel functions (c6) alpha = 2*sum(bessel_j[n](epsilon*n)*sin(2*n*\a)/n,n,1,inf)+2*\a inf ==== bessel_j (epsilon n) sin(2 n A) \ n (d6) alpha = 2 > ------------------------------- + 2 A / n ==== n = 1 which is equivalent to the double sum alpha = 2*'sum((-1)^k*('sum(epsilon^(n+k)*(n-k)^(n+k-1)*sin(2*(n-k)*\a)/(2^(n+k)*n!),n,0,k))/k!,k,1,inf)+2*\a k-1 ==== n + k n + k - 1 k \ epsilon (n - k) sin(2 (n - k) A) (- 1) > ---------------------------------------------- inf / n + k ==== ==== 2 n! \ n = 0 alpha = 2 > ----------------------------------------------------------- / k! ==== k = 1 + 2 A This converges very poorly as |epsilon| nears 1. Although Lambert_w seems inadequate to solve Kepler's directly, we have alpha = 2*\a+im((w+1)*'sum(\p[k](w)*z^k/k!,k,0,inf)) inf k ==== P (w) z \ k alpha = 2 A + im((w + 1) > --------) / k! ==== k = 0 where w:=lambert_w(-(epsilon*%e^(-2*%i*\a)/2)) - 2 %i A epsilon %e w := lambert_w(- ------------------) 2 and z:=-epsilon^2/(w*(2*w+2)^2) 2 epsilon z = - ------------, 2 w (2 w + 2) w P (w) = -----, 0 w + 1 and P[>0](w) are (not obviously) polynomials \p[k]:=funmake('lambda,['([w]),factor(subst(lambert_w=lambda([[y]],w),exp(-2*k*w)*diff(exp(2*k*lambert_w(-epsilon*exp(2*%i*y)/2)),epsilon,k)/2/k/w^k*epsilon^k*(w+1)^(2*k-1)))]) P (w) = 1 1 P (w) = 3 w + 2 2 2 P (w) = 20 w + 26 w + 9 3 3 2 P (w) = 210 w + 404 w + 273 w + 64 4 4 3 2 P (w) = 3024 w + 7692 w + 7672 w + 3524 w + 625 5 5 4 3 2 P (w) = 55440 w + 175272 w + 230422 w + 156374 w + 54505 w + 7776 6 6 5 4 3 2 P (w) = 1235520 w + 4667952 w + 7605084 w + 6801904 w + 3507828 w 7 + 985830 w + 117649 ... Unfortunately, the radius of convergence is even worse, due to the k^(k-1) growth of the P[k] coefficients. But if we take exactly the same number of terms of twice the Lambert sum minus the Bessel sum, alpha = -'sum((-1)^k*('sum(epsilon^(n+k)*(n-k)^(n+k-1)*sin(2*(n-k)*\a)/(2^(n+k-1)*n!),n,0,k))/k!,k,1,kmax)+2*\a+2*im((w+1)*'sum(\p[k](w)*z^k/k!,k,0,kmax)) k-1 ==== n + k n + k - 1 k \ epsilon (n - k) sin(2 (n - k) A) (- 1) > ---------------------------------------------- kmax / n + k - 1 ==== ==== 2 n! \ n = 0 alpha = - > ----------------------------------------------------------- / k! ==== k = 1 kmax k ==== P (w) z \ k + 2 A + 2 im((w + 1) > --------) / k! ==== k = 0 the errors magically cancel, and we need only kmax ~ 13 terms for double precision, even when epsilon = 1, which, with alpha=y-x+pi/2, A = pi/4-x, lets us solve the descending tilted cosine wave y + x = cos(y - x) for y(x), an interesting involution. --rwg ----------- hbaker>FYI -- For those who 1) have a mortgage; that 2) isn't underwater; and 3) still have a job & can qualify for refinancing: "Optimal Mortgage Refinancing: A Closed Form Solution" http://www.economics.harvard.edu/faculty/laibson/files/Mortgage%20refinanci… http://www.nber.org/papers/w13487 >An actual mortgage refinancing calculator based on the formulae from this paper: http://zwicke.nber.org/refinance/index.py

1 0

[math-fun] Root-factorial representation?
by Marc LeBrun 18 Sep '10

18 Sep '10

Can every integer be represented by a composition of integer square root ($) and factorial (!) applied to 3? Some example representations: 3 = 3 6 = 3! 1 = 3$ 2 = 3!$ 720 = 3!! 26 = 3!!$ 5 = 3!!$$ 120 = 3!!$$! 10 = 3!!$$! ... If not, what is the smallest counterexample?

5 5

Re: [math-fun] domino networks... correction
by Michael Beeler 16 Sep '10

16 Sep '10

Yesterday I wrote: > If we allow holes, I like the 5x6 with two holes. There's only > one distinct such pattern, [...] Nonsense. There are three. My program says: With the two holes at the corners on a long side, 0 solutions. With the two holes at diagonally opposite corners, 1 distinct solution, in yesterday's email. With the two holes at the corners on a short side, 9 distinct solutions: 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 3 3 4 4 5 3 3 4 4 5 3 3 4 5 1 3 3 4 5 1 3 3 4 5 1 6 7 8 6 1 6 7 8 6 1 6 2 4 6 6 6 2 4 6 6 6 7 4 6 6 6 2 5 5 7 6 5 2 2 7 7 2 8 8 9 7 2 8 8 9 8 7 2 2 9 9 2 3 9 7 9 5 3 9 7 1 9 3 7 9 7 9 3 7 1 1 9 3 8 9 1 8 8 9 4 1 8 8 9 4 8 5 5 7 4 4 9 5 5 8 7 5 5 8 4 1 1 2 1 2 3 1 2 3 1 2 3 3 3 4 5 1 2 4 5 3 4 2 4 5 3 4 3 4 2 5 5 6 7 4 6 6 6 4 7 8 8 6 4 7 8 8 6 4 7 7 8 8 7 2 2 9 6 9 1 1 5 8 9 1 1 5 6 5 1 6 2 8 9 3 8 1 8 9 3 6 5 2 9 3 6 5 9 9 1 8 9 4 9 5 5 7 2 2 7 7 9 2 7 7 6 9 7 3 3 8 4 -- Mike

2 2

Re: [math-fun] Cubic eqn with 3 real zeros, again
by Henry Baker 16 Sep '10

16 Sep '10

About 20 years ago I went through much of Henrici's work on "circular arithmetic", and while it was quite elegant, I didn't find in my experiments that it gained in either speed or stability. I may have missed something, because I wasn't able to fully grok all of his analysis. BTW, "circular arithmetic" is a version of interval arithmetic in which you use a center (in the complex plane) & radius instead of an upper & lower bound. The cool thing is that 1/z converts circles to circles; unfortunately, it also turns them inside out! Since circular arithmetic is a version of interval arithmetic, it shares many of the same problems: i.e., x^2 is (or should be) a different function from x*x, etc. At 10:29 AM 9/16/2010, Fred lunnon wrote: >On 9/16/10, Joerg Arndt <arndt(a)jjj.de> wrote: >> ... >> This may not be what you have in mind but there are methods to >> find _all_ roots simultaneously, see e.g. >> http://en.wikipedia.org/wiki/Durand-Kerner_method >> and the references at the bottom of the page > >A very elegant method which computes all roots simultaneously, >yielding circles in which they are guaranteed to lie, is buried >(I don't remember exactly where --- anybody else know?) >in the strangely neglected > >Henrici P., Applied and Computational Complex Analysis (Wiley). >[Three volumes: 1974, 1977, 1986.] > >It's interesting that these methods actually gain in speed and >stability by simultaneity, as opposed to attempting to exclude the >other roots. > >WFL

1 0

Re: [math-fun] Two more damn "Why are these integral?" puzzles
by Bill Gosper 16 Sep '10

16 Sep '10

Joerg> Interesting. * Bill Gosper <billgosper(a)gmail.com <http://gosper.org/webmail/src/compose.php?send_to=billgosper%40gmail.com>> [Sep 15. 2010 09:04]: > [...] The factors... > 2 2 2 2 > z = 128 (y - 1) y (y + 1) (2 y - 1) (2 y + 1) (2 y - 1) (4 y - 3) > 4 > > 4 2 2 4 2 > (8 y - 8 y + 1) (16 y - 16 y + 1) > > 8 6 4 2 > (256 y - 512 y + 320 y - 64 y + 1) ...look somewhat Chebyshev-oid(?) Bingo. z[n](y)=T[2^(n+1)](x)-T[2^n](x), and w[n](x)=2*(T[2^n](x)+1)*(2*T[2^n](x)-1) . I was looking for Pollard rho recurrences that directly compute slow minus fast so you need only run one sequence. --rwg

1 0

[math-fun] Isodiametric problem
by victor miller 16 Sep '10

16 Sep '10

A few days ago, at the AMS meeting, I saw a neat talk by Michael Mossinghoff about "isodiametic" problems -- e.g. Fix the diameter of a conv ex n-gon and ask for the ones with maximum perimeter. He reduced the problem to finding univariate polynomials whose non-zero coefficients alternate between +/- 1 and are divisible by a particular cyclotomic polynomial.. He was able to calculate all such up to about 500. Each such polynomial is compactly encoded by th length of the gaps between the non-zero terms. For small n almost all were periodic. He showed how a certain family of periodic constrctions give solutions, but as n got bigger the "sporadic" ones swamped them. He conjectued that these are eventually all of them. He had lots of pictures. The paper is not available yet. Victor

7 9

Re: [math-fun] Cubic eqn with 3 real zeros, again
by Henry Baker 15 Sep '10

15 Sep '10

Yes, I've played with a version of this idea over the years, but have never gotten a satisfactory algorithm out of it, yet. The idea is an old one, but I don't recall off the top of my head what the name usually attached to this is. One of the problems is that the factor (x-r)/(x-r'), for constants r~r', but r/=r', doesn't act as much like "1" as one might like. So this "shielding" of a factor that you're sufficiently close to doesn't work as well as one might hope. From an algorithmic efficiency point of view, you're better off finding the root you're already close to, because convergence is quadratic. Then, you divide out that root so you have a smaller polynomial to work with for the next step. At 03:43 PM 9/15/2010, Allan Wechsler wrote: >Henry's post got me thinking in a different direction, maybe not the one he >had intended. I was thinking about Newton's method and the problem that it >is vulnerable to instabilities; even when it works perfectly and converges >fast, it gets preoccupied with one root and ignores the others unless you >start it over. > >So I started wondering whether anybody had played with a sort of adaptive >version of Newton's algorithm. While one "thread" is converging on a root, >a "kibitzer" thread divides the original function by x-r, as if r were the >actual root, and gets a modified function that can be the host for another >round of Newton's algorithm that would be scared away from the first root >and be more likely to converge on the others. Perhaps I'm barking up a >squirrel-less tree, but I was imagining different threads converging on >different roots, and dividing out by each other's approximations in order to >"isolate" their own target root. > >Has anyone heard of an approach like this? > >On Wed, Sep 15, 2010 at 6:14 PM, <rcs(a)xmission.com> wrote: > >> This suggests a meta-problem: If you find your numerical algorithm is >> converging slowly, can you introduce an extra dimension, perhaps >> imaginary, to speed things up? >> >> Rich >> >> ------------- >> >> Quoting Henry Baker <hbaker1(a)pipeline.com>: >> >> Has anyone seen the following treatment of the real cubic? >>> >>> A cubic equation with 3 real zeros can be reduced to either (your choice): >>> >>> 4*c^3-3*c=c3 or 4*s^3-3*s=-s3 >>> >>> where c,c3 reminds us of cos(a),cos(3*a), and s,s3 reminds us of >>> sin(b),sin(3*b), respectively. >>> >>> Now Newton's method _can_ be used to compute one of the zeros, but an >>> unfortunate choice of starting point will attempt to converge on a double >>> (or almost double) root, thereby seriously slowing the convergence. >>> >>> For example, Newton on the cosine version gives >>> >>> c' = (c3+8*c^3)/[3*(2*c-1)*(2*c+1)] >>> >>> which leads to very unfortunate behavior for c near +- 1/2. >>> >>> But if we recognize that we are attempting to solve the 1-dimensional >>> projection of a 2-dimensional problem, then we can avoid this convergence >>> problem. >>> >>> Specifically, if we are given 4*c^3-3*c=c3, then we first compute >>> s3=sqrt(1-c3^2), and then instead solve the problem >>> >>> z^3 = z3 = >>> = c3+i*s3 = >>> (c+i*s)^3 = c3+i*s3 >>> = c^3 - 3*c*s^2 - i*(s^3 - 3*s*c^2) >>> = c^3 - 3*c*(1-c^2) - i*(s^3 - 3*s*(1-s^2)) >>> = c^3 - 3*c +3*c^3 - i*(s^3 - 3*s + 3*s^3) >>> = 4*c^3 - 3*c - i*(4*s^3 - 3*s) >>> = c3 - i*(-s3) >>> = c3 + i*s3 >>> >>> But Newton's method on the complex problem z^3=z3 is much more robust, so >>> long as we stay away from z=0: >>> >>> z' = (z3+2*z^3)/(3*z^2) >>> >>> (Throw away s=imagpart(z) & keep c=realpart(z) for the answer. But >>> before doing so, we can trivially produce the other two roots with complex >>> multiplications: w*z and w^2*z, where w is cube root of unity: w^3=1.) >>> >>> But we know that |z|^2 must be 1, so we can easily avoid the region near >>> z=0. >>> >>> On modern hardware with multiple floating point units, complex arithmetic >>> may have no more latency than real arithmetic, so this method may be no >>> slower, and (due to faster convergence) may actually be faster.

1 0

[math-fun] Two more damn "Why are these integral?" puzzles
by Bill Gosper 15 Sep '10

15 Sep '10

This time, only 2nd order: (c44) (DONE,Z[N] := FACTOR(2*((Z[N-1]-2)*Z[N-1]^2*(Z[N-1]+1)-8*(Z[N-2]-2)*Z[N-2]^2*(Z[N-2]+1)^3*(2*Z[N-2]^2-2*Z[N-2]-1))/(Z[N-1]*(2*Z[N-1]+1)+(-2*Z[N-2]*(Z[N-2]+1)*(2*Z[N-2]+1)*(2*Z[N-2]+5))))); 2 (d44) z := factor(2 ((z - 2) z (z + 1) n n - 1 n - 1 n - 1 2 3 2 - 8 (z - 2) z (z + 1) (2 z - 2 z - 1)) n - 2 n - 2 n - 2 n - 2 n - 2 /(z (2 z + 1) + (- 2 z (z + 1) (2 z + 1) n - 1 n - 1 n - 2 n - 2 n - 2 (2 z + 5)))) n - 2 (c45) [Z[0] : (Y-1)*(2*Y+1),Z[1] : 2*(Y^2-1)*(4*Y^2-1)]; 2 2 (d45) [(y - 1) (2 y + 1), 2 (y - 1) (4 y - 1)] (c46) BLOCK([FANCY_DISPLAY : FALSE],FOR K FROM 0 THRU 5 DO PRINT(ARRAYMAKE(Z,[K]) = Z[K]))$ z = (y - 1) (2 y + 1) 0 2 2 z = 2 (y - 1) (4 y - 1) 1 2 2 z = 8 (y - 1) y (y + 1) (2 y - 1) (2 y + 1) (4 y - 3) 2 2 2 2 2 z = 32 (y - 1) y (y + 1) (2 y - 1) (2 y + 1) (2 y - 1) (4 y - 3) 3 4 2 (16 y - 16 y + 1) 2 2 2 2 z = 128 (y - 1) y (y + 1) (2 y - 1) (2 y + 1) (2 y - 1) (4 y - 3) 4 4 2 2 4 2 (8 y - 8 y + 1) (16 y - 16 y + 1) 8 6 4 2 (256 y - 512 y + 320 y - 64 y + 1) 2 2 2 2 z = 512 (y - 1) y (y + 1) (2 y - 1) (2 y + 1) (2 y - 1) (4 y - 3) 5 4 2 2 4 2 (8 y - 8 y + 1) (16 y - 16 y + 1) 8 6 4 2 2 (128 y - 256 y + 160 y - 32 y + 1) 8 6 4 2 (256 y - 512 y + 320 y - 64 y + 1) 16 14 12 10 8 (65536 y - 262144 y + 425984 y - 360448 y + 168960 y 6 4 2 - 43008 y + 5376 y - 256 y + 1) (c38) W[N] := FACTOR((W[N-1]*(W[N-1]^3+2*W[N-1]^2-6*W[N-1]-16)-W[N-2]*(W[N-2]+2)*(W[N-2]^6-2*W[N-2]^5-10*W[N-2]^4+8*W[N-2]^3+26*W[N-2]^2-4*W[N-2]-24))/(2*(W[N-1]^2+3*W[N-1]-W[N-2]*(W[N-2]+2)*(W[N-2]^2+6*W[N-2]+7)))); 3 2 (d38) w := factor((w (w + 2 w - 6 w - 16) n n - 1 n - 1 n - 1 n - 1 6 5 4 3 - w (w + 2) (w - 2 w - 10 w + 8 w n - 2 n - 2 n - 2 n - 2 n - 2 n - 2 2 2 + 26 w - 4 w - 24))/(2 (w + 3 w n - 2 n - 2 n - 1 n - 1 2 - w (w + 2) (w + 6 w + 7)))) n - 2 n - 2 n - 2 n - 2 (c39) [W[0] : 2*(2*X^2+X-1),W[1] : 4*X^2*(4*X^2-3)]; 2 2 2 (d39) [2 (2 x + x - 1), 4 x (4 x - 3)] (c40) FOR K FROM 0 THRU 5 DO PRINT(ARRAYMAKE(W,[K]) = W[K])$ 2 w = 2 (2 x + x - 1) 0 2 2 w = 4 x (4 x - 3) 1 2 2 4 2 w = 4 (2 x - 1) (16 x - 16 x + 1) 2 4 2 2 8 6 4 2 w = 4 (8 x - 8 x + 1) (256 x - 512 x + 320 x - 64 x + 1) 3 8 6 4 2 2 w = 4 (128 x - 256 x + 160 x - 32 x + 1) 4 16 14 12 10 8 (65536 x - 262144 x + 425984 x - 360448 x + 168960 x 6 4 2 - 43008 x + 5376 x - 256 x + 1) 16 14 12 10 w = 4 expt(32768 x - 131072 x + 212992 x - 180224 x 5 8 6 4 2 + 84480 x - 21504 x + 2688 x - 128 x + 1, 2) 32 30 28 (4294967296 x - 34359738368 x + 124554051584 x 26 24 22 - 270582939648 x + 392586854400 x - 401311006720 x 20 18 16 + 297124495360 x - 161296154624 x + 64266436608 x 14 12 10 8 - 18627952640 x + 3852599296 x - 550371328 x + 51597312 x 6 4 2 - 2924544 x + 87040 x - 1024 x + 1) --rwg

2 1

Re: [math-fun] mechanical realization of a 3-torus
by George W. Hart 14 Sep '10

14 Sep '10

Jim, I made a simple mockup of the idea by printing linkages onto a laser-printable transparency, cutting them out, punching some holes, and pinning it together with small nuts and bolts on a scrap of cardboard. Here are images with two different choices for which pair of sums is inside the unit circle: http://www.georgehart.com/MoMath/3-torus-A.jpg http://www.georgehart.com/MoMath/3-torus-B.jpg The links with solid lines are your u,v,w directions and the ones with dashed lines are just to maintain parallelism. In this form, there are some interlocking issues that make it difficult to quickly move to any desired configuration. Links will hang up on bolts instead of passing by them. I'm not sure how to get around that. As to your puzzle, my feeling from playing with it is that there is just one connected component to the space when one of the sums (red dots) is fixed with magnitude almost 3. But this prototype is not at all rigid, so it is hard to be confident. As to developing this idea into a possible exhibit for the Museum of Mathematics, I'd be interested to continue discussing that. George http://momath.org/ http://georgehart.com/ On 9/10/2010 5:05 PM, James Propp wrote: > I have an idea for a physical mathematical demo, and I'm wondering > if anyone has any thoughts about it (or any interest in actually > building it). > > Spitzer's "random flight theorem" says that a sum of n independent > random unit vectors in the plane has a probability of exactly > 1/(n+1) of having magnitude less than 1. In particular, the > chance that the sum of 3 independent random unit vectors has > magnitude less than 1 is exactly 1/4. It so happens that this > particular fact can be deduced from the following geometrical > fact, which may or may not be new: If u, v, and w are unit > vectors, then exactly two of the eight vectors > +u+v+w, +u+v-w, +u-v+w, +u-v-w, -u+v+w, -u+v-w, -u-v+w, -u-v-w > have magnitude less than 1 (except in the case when some of > these vectors have magnitude exactly 1). > > There's a nice way to physically model this "two out of eight" > fact with a linkage consisting of bars of length 2, so that > eight points move jointly so as to trace out the configuration > space {(+u+v+w, +u+v-w, +u-v+w, +u-v-w, -u+v+w, -u+v-w, -u-v+w, > -u-v-w): |u|=|v|=|w|=1}. Now imagine that the bars are made of > a transparent material like plexiglass, and that the model has > a backing that shows a disk of radius 1 centered at 0. As you > move around in configuration space, there are always exactly > two points that lie inside the disk, except for the borderline > case where four (or more) points lie on the boundary of the disk. > > I'm wondering if this might make a good object for the Museum > of Mathematics, or maybe just a fun object to own and to display > in my office. I think I know what happens to the configuration > space if you pin down one of the eight points, but it'd be cool > to really feel what happens. > > Note that one could describe this linkage simply as a mechanical > realization of a 3-torus (hence the title of this post). > > Jim > > _______________________________________________ > math-fun mailing list > math-fun(a)mailman.xmission.com > http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun

2 2

Re: [math-fun] mechanical realization of a 3-torus
by Dan Asimov 13 Sep '10

13 Sep '10

Seb, That's a very nice Java application! Just in case you're feeling ambitious, would you be willing to also build a Java program that displays three rotatable rods, U, V, W, of equal length L in the plane, each fixed at one end -- at points P, Q, R, respectively, that form the vertices of an equilateral triangle of side S. The lines connecting the variable ends of U, V, W -- i.e., P + U(x), Q + V(y), R + W(z) -- could be thick and colored red, blue, and green. (Here 0 <= x,y,z < 2pi are the angles of U, V, W.) This would display each point of the configuration space as a distinct colored triangle -- among triangles whose vertices lie on the three circles of radius L about P, Q, R. Since triangles are familiar objects, this might be a palatable way to digest the idea of a configuration space. Finally, it would be interesting to experiment with several values for the ratio S/L -- I'd say especially S/L = 0, 1, 2 -- to see which mechanism looks more interesting or informative. (Of course S/L = 0 means P = Q = R.) If you're feeling *really* ambitious, there could also be an automatic mode in which the user can push U, V, W off to keep moving by itself along some geodesic. Since almost all ideal such geodesics are dense, this will give a "grand tour" of the configuration space T^3. --Dan _____________________________________________________________________ "It don't mean a thing if it ain't got that certain je ne sais quoi." --Peter Schickele

1 0