July 2008 - mathfun@mailman.xmission.com

[math-fun] Smallest magic squares of powers
by Christian Boyer 01 Sep '08

01 Sep '08

Oct 7th, 2007: the smallest known magic square of 4th powers is now a 44x44 square constructed by Jaroslaw Wroblewski (Poland) and Hugo Pfoertner (Germany). Jaroslaw built the semi-magic square in June, then Hugo needed a very impressive computing job of 3 months, on 57 processors of a SGI computer http://www.sgi.com/products/remarketed/origin3000/, in order to find two magic diagonals in this square! A big progress, because the previous smallest known was a 243x243 square, coming from the 243x243 tetramagic square of Pan Fengchu (China). The best known results, and their first authors, are now: Magic square Smallest known ------------ --------------- Of squares 4x4 (Euler, 1770) Of cubes 9x9 (Boyer, 2006) Of 4th powers 44x44 (Wroblewski-Pfoertner, 2007) Of 5th powers 729x729 (from the pentamagic square of Li Wen, 2003) ON ANY POWER, WE DO NOT KNOW WHAT THE SMALLEST POSSIBLE MAGIC SQUARE IS! For example, as you know, we still do not know if a 3x3 magic square of squares is possible... 4x4, 5x5, 6x6, 7x7, and 8x8 magic squares of cubes are unknown... And the problem of the diagonals is extremely difficult, as proved by the impressive computing search of Hugo. It remains a big gap between the smallest known magic and semi-magic squares (=with or without magic diagonals) of 4th powers: 44x44 versus 4x4... Here is the 4x4 semi-magic square of 4th powers (Morgenstern-Boyer, 2006): 14101^4 938^4 931^4 37762^4 35866^4 20881^4 21038^4 13393^4 24806^4 30191^4 30418^4 9263^4 413^4 32026^4 31787^4 1106^4 Christian.

2 2

[math-fun] new inverter, new format : 2.459 billion mathematical constants.
by Simon Plouffe 24 Aug '08

24 Aug '08

Hello, I am downloading it now, the new version of my inverter. I changed the format to 64 digits and normalized. Well, there are 2.459 billion mathematical constants in it. This time I used Maple and Mathematica and in some occasion Pari-Gp to generate the numbers. The data is presently <dumped> into a work directory here : http://www.lacim.uqam.ca/~plouffe/ip/ the files are gziped, to take a look at one just click on the file and normaly the firefox browser will uncompress it on the fly. If you do not use firefox you may save it on your computer and uncompress it with winrar of gunzip (unix). The pc (actualy it is a mac) where the database will reside is currently being organized in MontrÃ©al, there are 90000 files for a total of 297 gigabytes. The disk is a 300 gigabytes, so I think it is time I stop generating constants! In the database there are something like 25 million constants that were generated by the OEIS database. For the moment, I intend to at least organize the data so that anybody interested in it may download (yes download) the database. This summer I intend to re-work the lookup script to take into account that new format. The database it of course what I think is an interesting constant. The base table has about 186 known constants http://www.lacim.uqam.ca/~plouffe/b000.txt and the rest are entries similar to the first and second versions of this project (the inverse symbolic calculator 1995-now) the plouffe inverter (1998-now). If you have any suggestions, of course they are welcome. note : the current plouffe_inverter program to analyze math constants will be compatible with the new format of course. It is still in Maple syntax, I should change it to mathematica, mathematica is faster, much faster than maple for many computations... For the nostalgic, the current inverter with 200 million constants is all here : http://www.lacim.uqam.ca/~plouffe/plouffe_inverter_tables/ Simon Plouffe

2 1

Re: [math-fun] Theta function halfangles
by rwg＠sdf.lonestar.org 05 Aug '08

05 Aug '08

> From: rwg(a)sdf.lonestar.org [mailto:rwg@sdf.lonestar.org] > Sent: Wed 7/9/2008 9:22 AM > To: math-fun(a)mailman.xmission.com > Cc: "rwg(a)sdf.lonestar.orgrwg"@sdf.lonestar.org > Subject: Theta function halfangles > > (I promised Cheny Xu, proprietor of Jade Palace Restaurant, www.jadepalace.us > that I'd solve a problem with the help of some fish he fed me last night.) > > It's not obvious how to extract a tractable halfangle formula from the many > double angle formulae in W&W. After you eliminate the three undesired > species, the resultant polynomial is too hairy. Actually, it would > probably simplify if Macsyma or Mma only knew how, but I can't even get Mma > 6.0 to EllipticTheta[1, Pi, q] -> 0 ! (Hey guys, need a programmer?) > > Fortunately, I (re?)conjectured a family of theta determinant identities, > http://www.tweedledum.com/rwg/thet.PNG (too tall for IE 7), > arxiv.org/pdf/math/0703470, the first few (including the one used here) proved by > Rich (and Hilarie?). These are nice when actually typeset: > > x > theta (-, q) = > 1 2 > > 3 3 3 > - theta (0, q) theta (x, q) + theta (0, q) theta (x, q) - theta (0, q) theta (x, q) > 4 4 3 3 2 2 1/4 > (-----------------------------------------------------------------------------------) , > 2 > > x > theta (-, q) = > 2 2 > > 3 3 3 > - theta (0, q) theta (x, q) + theta (0, q) theta (x, q) + theta (0, q) theta (x, q) > 4 4 3 3 2 2 1/4 > (-----------------------------------------------------------------------------------) , > 2 > > x > theta (-, q) = > 3 2 > > 3 3 3 > theta (0, q) theta (x, q) + theta (0, q) theta (x, q) + theta (0, q) theta (x, q) > 4 4 3 3 2 2 1/4 > (---------------------------------------------------------------------------------) , > 2 > > x > theta (-, q) = > 4 2 > > 3 3 3 > theta (0, q) theta (x, q) + theta (0, q) theta (x, q) - theta (0, q) theta (x, q) > 4 4 3 3 2 2 1/4 > (---------------------------------------------------------------------------------) . > 2 > > E.g., > pi pi > theta (--, q) = theta (--, q) = > 1 4 2 4 > > 2 2 > theta (0, q) theta (0, q) (theta (0, q) - theta (0, q)) > 3 4 3 4 1/4 > (-------------------------------------------------------) , > 2 > > pi pi > theta (--, q) = theta (--, q) = > 3 4 4 4 > > 2 2 > theta (0, q) theta (0, q) (theta (0, q) + theta (0, q)) > 3 4 4 3 1/4 > (-------------------------------------------------------) . > 2 But, duh, %pi %pi 4 theta (---, q) = theta (---, q) = theta (0, q ) 4 4 3 4 4 > Finally, this must be in Tannery & Molk, but I don't recall seeing a Lambert type > series come out as a pure theta instead of a log derivative theta: > > inf 2 > ==== j 2 j - 1 theta (0, q) - 1 > \ (- 1) q 4 > > --------------- = ----------------. > / 2 j - 1 4 > ==== q + 1 > j = 1 There seems to be a lot of these. theta (0, q) 2 log(-------------------) inf 1/4 ==== n n 2 theta (0, q) q \ (- 1) q 3 > ---------- = ------------------------. / n 2 ==== n (q + 1) n = 1 > --rwg >

5 19

[math-fun] Superfluous Multinomials
by wouter meeussen 04 Aug '08

04 Aug '08

dear funsters, when posting stuff like this, I presume it is generally known, or at least could be. Motivation: just fun and beauty. 1/. Intro consider the expansion of (a+b+c+d)^3 it is the sum over M(u,v,w,x) a^u b^v c^w d^x with positive or zero integers u+v+w+x=3 and M(u,v,w,x) is the multinomial u! v! w! x! /(u+v+w+x)! When we drop the multinomial coefficients, then we get the 'compositions' of 4 elements in groups of three. in Mma: Expand[(a+b+c+d)^3 /3! ]/. q_^n_ ->q^n n! gives a^3 + a^2*b + a*b^2 + b^3 + a^2*c + a*b*c + b^2*c + a*c^2 + b*c^2 + c^3 + a^2*d + a*b*d + b^2*d + a*c*d + b*c*d + c^2*d + a*d^2 + b*d^2 + c*d^2 + d^3 this is equivalent to Plus @@ Apply[Times,{a,b,c,d}^#& /@Compositions[3,4],{1}] 2/. Generalise replace the (a+b+c+d) above by Sum(j=1..n, t[ j ] ) , t stands for 'temporary', now replace t[j] with E^(2 Pi I j/n) and it's obvious that this sum equals zero. And that any power of zero is also zero. But now, dropping multinomials, we find: Table[ Plus@@ Apply[Times,(Exp[2 Pi I/n]^Range[n])^#& /@ Compositions[k,n],{1}],{n,5},{k,n}] and we find that it's also zero, except it's 1 when k equals n, that is, we take all n roots of (-1) and put it to the n-th (multinomial-less) power. 3/. an instructive typo : Exp[2 Pi I/k] instead of Exp[2 Pi I/n]: Table[ Plus@@ Apply[Times,(Exp[2 Pi I/k]^Range[n])^#& /@ Compositions[k,n],{1}],{n,5},{k,n}] equals Table[Ceiling[n/m],{n,9},{m,n}] alias A120885 That leads to this silly remark: Table[ Plus@@ (Exp[2 Pi I/k]^Range[n])^k ,{n,16},{k,n}] is integer, except for the following couples of {n,k} : {7, 5}, {8, 5}, {9, 7}, {10, 7}, {10, 8}, {11, 7}, {11, 8}, {11, 9}, {12, 5}, {12, 7}, {12, 8}, {12, 9}, {12, 10}, {13, 5}, {13, 8}, {13, 9}, {13, 10}, {13, 11}, {14, 8}, {14, 9}, {14, 10}, {14, 11}, {14, 12}, {15, 9}, {15, 10}, {15, 11}, {15, 12}, {15, 13}, {16, 7}, {16, 9}, {16, 10}, {16, 11}, {16, 12}, {16, 13}, {16, 14}, ... and, where's the system in *this* madness? Wouter.

4 6

Re: [math-fun] Presentation on floating point implementation
by Henry Baker 29 Jul '08

29 Jul '08

Here's some additional info on quadratics, including trig solutions: http://home.pipeline.com/~hbaker1/sigplannotices/sigcol05.pdf At 06:06 AM 7/29/2008, James Buddenhagen wrote: >I saw this back around 1966 in a paper by George Forsythe, some form >of which may have been in the Monthly (or maybe not?). The original >technical report seems to be here > >ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/66/40/CS-TR-66-40.pdf > >Jim Buddenhagen > >On Mon, Jul 28, 2008 at 11:16 PM, Mike Speciner <ms(a)alum.mit.edu> wrote: >> Those were nice. One of the slides reminded me that for solving quadratics >> (-b +- sqrt(b^2-4ac))/2a = 2c/(-b -+ sqrt(b^2-4ac)), >> so better precision for the two roots can be obtained by picking one root >> from each side so that both can have the sign of the sqrt same as -b. (And >> then it's easy to see the RHS root approaching -c/b as a->0.) >> >> I don't remember where I first saw this trick, but I've certainly made use >> of it many times. >> >> --ms >> >> -----Original Message----- >> From: math-fun-bounces(a)mailman.xmission.com >> [mailto:math-fun-bounces@mailman.xmission.com]On Behalf Of Jason >> Sent: Sunday, July 27, 2008 14:21 >> To: math-fun(a)mailman.xmission.com >> Subject: [math-fun] Presentation on floating point implementation >> >> I thought these slides were well-paced and interesting: >> >> http://www.research.scea.com/gdc2003/fast-math-functions_p1.pdf >> http://www.research.scea.com/gdc2003/fast-math-functions_p2.pdf

1 0

[math-fun] Presentation on floating point implementation
by Jason 29 Jul '08

29 Jul '08

I thought these slides were well-paced and interesting: http://www.research.scea.com/gdc2003/fast-math-functions_p1.pdf http://www.research.scea.com/gdc2003/fast-math-functions_p2.pdf

4 3

Re: [math-fun] dilog puzzle
by Dan Asimov 27 Jul '08

27 Jul '08

<< > (Yes, Henry, ln phi = asinh 1/2 .-) I'm not even managing to clamber in on the ground floor here. _Why_ is it true? >> Lessee: sinh(x) := (e^x - e^-x)/2, so one branch (the usual one) of the inverse function is asinh(y) = ln(x + sqrt(x^2 + 1). So asinh(1/2) = ln( 1/2 + sqrt(5/4) ) = ln(phi). --Dan _____________________________________________________________________ "It don't mean a thing if it ain't got that certain je ne sais quoi." --Peter Schickele

4 3

[math-fun] Zipf, codes & cumulative probability
by Steve Witham 25 Jul '08

25 Jul '08

On 2006.06.15, Henry Baker started a thread on binary encodings for Zipfian distributions. That is, distributions like { 1/1, 1/2, 1/3, 1/4, 1/5 ... 1/m } He asked (at various points in the thread): 1) This must be in the literature, can anyone give a ref? 2) What efficient encoding(s) exist(s)? 3) Although a Zipf distribution must be finite, are there universal encodings (i.e. for infinite ranges) that approach 100% efficiency for { 1/1 .. 1/m } as m gets big? 4) Universal codes give about log2(n) bits to the number n; are there inefficient codes that give, say, sqrt(n) bits? This letter promotes cumulative probability functions as a handle on all this. 1) References As far as I recall, no one found any for Zipf distributions in the first round of this thread. I notice this sentence in http://en.wikipedia.org/wiki/Universal_code_(data_compression) : "With universal codes, the implicit distribution is approximately a power law such as 1 / n^2 (more precisely, a Zipf distribution). " 2) Encoding Zipf for 1 through finite m. Say you have a relative probability distribution f'( n ) for n = 1..m. This will go easier for both of us if f'(n) is sorted in decreasing order. Now, if there is a simple function f''(x) that approximates f'(n) (& a couple more conditions), we're in business. e.g., for the Zipf distribution, f'(n) = 1/n; f''(x) = 1/x One condition is that you need to be able to integrate f''(n), and also find the inverse of the resulting function pretty easily. e.g. F''(x) = log(x) + c; G''(y) = exp(y - c) This is the key point: f'(n) is a relative probability distribution; f''(x) is similar and almost a probability density function; F''(x) is a corresponding, almost cumulative probability function. Why settle for almost? Let F(x) be F''(x) adjusted so that F(.5) = 0 and F( m + .5 ) = 1. Let f(x) be d F(x)/dx and G(y) is the inverse of F(x). e.g., F(x) = k log( x ) + c; G(x) = exp( (y-c)/k ); f(x) = k/x (The slight fudge here is that instead of giving n a probability proportional to f'(n) in our code, we're really going to give it the definite integral from n-.5 to n+.5 of f(x) dx. This difference approaches zero for larger numbers and can be re-fudged for smaller numbers, so we'll ignore the issue from now on.) Now you can pretty much encode the value n as the binary representation of F(n), and decode codeword w with G(w). The length of the codeword for n is approximately -log2( f(n) ) (unrelated to the fact that F() is log-like in the example here). How exactly though? Here's an existence proof of a code. Say you've got F, G and f now. To encode n, generate F(n-1) (if n>1), F( n ), and F(n+1) (if n < m) to at least a couple more bits of precision than -log2(f(n)) would call for. Use just enough bits of F(n) to distinguish from both neighbors (or the one neighbor if n = 1 or n = m). I believe the code you get that way uses the whole code space. It'll be fairly efficient, but not as good as a Huffman code. I've been thinking about tricks to make it more efficient, but they would obscure the point that it's the binary representation of the cumulative probability function. Decoding is a little trickier since what you're given is a message without knowing where the first codeword ends. Nevertheless you treat the whole message as a binary number w, n = G(w) +/- 1, and you consume about -log2(f(n)) bits after looking at the neighbors. 3) Universal codes and their efficiency for Zipf distributions A "universal" code is one that can encode numbers over an infinite range, e.g. any n >= 1. The whole method above works for infinite ranges if you have a distribution (a positive function whose integral is finite) over the range. The only problem is you need unlimited-precision arithmetic. In fact, for decoding you need to increase your precision until you know you've gone far enough. Typical "mixed" methods for universal codes break the number into an exponent and a mantissa, and use binary for the mantissa and some other method for the exponent. But with this method you just look at the values of a single function. One advantage to this is that it can be smoother and therefore more efficient at the low end. What kind of function, though? One trick is to pick a function F() that, over the range in question (e.g. .5 to infinity) goes from zero to one, and let f() be F's derivative. One F() I was looking at way back when was: F(x) = 1 - ( 1 / log2( x + 1.5 ) ) Using my soggy brain plus a table of derivatives I get f(x) = d/dx F(x) = -1/( (x+1.5) (ln 2) (log2(x+1.5))^2 ) This means that the length of the codeword for n, -log2(f(n)) is about log2(ln 2) + log2(n+1.5) + 2 log2( log2( n+1.5 ) ) In other words, a small negative constant, plus enough bits for a mantissa, plus twice enough bits for an exponent...only creamy and smooth. As m, the largest n in a Zipf distribution, gets large the constant and the log(log) terms become less significant, so the code based on this function approaches 100% efficiency (slowly though). It's slightly strange that having the length of the codeword be the log of the number is the natural limit for a universal code. The same is true of the way we write decimal integers. But this corresponds to the probabilities of a Zipf distribution, which is the distribution none of us knew or found a reference to in the coding literature. Another interesting thing to think about is that no finite-length code can grow more slowly than log2(n) in the limit. I guess you could optimize for a distribution of distributions, e.g., Zipf distributions for a distribution of m. By the way, there are known ways to generate infinite Huffman codes, although they mostly work for variations on geometric distributions. 4) Non-optimal and funky universal codes If you want a universal code whose wordlength grows at a more-than- log rate, say like sqrt(n), then you take -log2(f''(x)) = sqrt(x) f''(x) = 2 ^ (-sqrt(x)) integrate that to get F''(x), then adjust to get F(), f() and G(). An interesting puzzle might be to take a given binary code (e.g, unary, Zeckendorf/Fibonacci, ternary, septary, various mixed universal codes) and see whether the implied distribution resembles a simple function. I think some universal codes relate to the "log*" function, but I don't know whether a natural continuous version of log* has been found. --Steve p.s. none of this needed to mention my *previous* pet universal tool, histograms of codeword-lengths.

3 3

[math-fun] Zipf, codes & cumulative probability
by Steve Witham 21 Jul '08

21 Jul '08

>f'(n) is a relative probability distribution; >f''(x) is similar and almost a probability density function; >F''(x) is a corresponding, almost cumulative probability function. (F'' is the integral of f''.) >Let F(x) be F''(x) adjusted so that F(.5) = 0 and F( m + .5 ) = 1. >Let f(x) be d F(x)/dx and G(y) is the inverse of F(x). There's an argument for making f(x) = F(n+.5) - F(n-.5) since that's the approximate chunk of code space this method gives to n. One version of my method needs to invert f(), though, and I'm not sure that's always so easy. (If I thought this would all be as easy (& applicable to as many kinds of distributions) with differences and sums instead of derivatives and integrals, I'd do that since it's more correct.) >Here's an existence proof of a code. [...] To encode n, generate > F(n-1) (if n>1), > F( n ), and > F(n+1) (if n < m) >to at least a couple more bits of precision than -log2(f(n)) would call >for. Use just enough bits of F(n) to distinguish from both neighbors >(or the one neighbor if n = 1 or n = m). > >I believe the code you get that way uses the whole code space. No, it does get skips (although never overlaps). Example: 1 00 2 01 3 1010 << 100* is skipped 4 1011 5 110 6 11100 7 11101 8 1111 >It'll >be fairly efficient, but not as good as a Huffman code. Yes, particularly when there isn't a skip! Here are some numbers for small m: 1 .. 2 : 2.15066010309 uniq_bin 0.901524769097 huffman 0.901524769097 1 .. 3 : 2.15066010309 uniq_bin 0.982262472148 huffman 0.982262472148 1 .. 4 : 2.15066010309 uniq_bin 0.996698112389 huffman 0.996698112389 1 .. 5 : 2.15066010309 uniq_bin 0.973892050742 huffman 0.985497700151 1 .. 6 : 2.15066010309 uniq_bin 0.981575348361 huffman 0.981575348361 1 .. 7 : 2.15066010309 uniq_bin 0.981487394329 huffman 0.981487394329 1 .. 8 : 2.15066010309 uniq_bin 0.888319512382 huffman 0.97949565162 ^^^ this is the code above with a skip ...and some large m's.: 1 .. 16 : uniq_bin 0.978314000279 huffman 0.988963192851 1 .. 256 : uniq_bin 0.99063060087 huffman 0.993030526135 1 .. 4096 : uniq_bin 0.994895289226 huffman 0.9972730518 1 .. 65536 : uniq_bin 0.994403026895 huffman 0.996899550792 1 .. 1048576 : uniq_bin 0.991023380916 huffman 0.997858670277 These are "comparing apples to apples" by using F(n+.5) - F(n-.5) as the probabilities for both methods. >Another interesting thing to think about is that no finite-length code >can grow more slowly than log2(n) in the limit. It has to be ever-so-slightly bigger than log2(n). Something like size(n) = log2(n) + size( log2(n) ) Notice that that goes berserk if you recur enough times, so it can't be right. It looks like Elias might have been the one to have thought this through, but maybe not. The whole idea of "optimum" codes is a little strange: Huffman's proof that his method is optimum (among "prefix-free codes") is a sort of touchstone for entropy coders, but the differences in efficiency are like the .991 vs. .998 above. It's just that when you're analyzing an algorithm, once you know it's optimal, you don't have to think about that issue anymore. The difference between one universal code that encodes a number like 2^(2^(2^(2^(2^(2^2))))) with 2^(2^(2^(2^(2^2)))) + 2^(2^(2^(2^2))) + 2^(2^(2^2)) + 100 bits compared to one that uses 2^(2^(2^(2^(2^2)))) + 2^(2^(2^(2^2))) + 2^(2^(2^2)) + 50 bits is hypernegligible. But fascinating. But, granted that it's a perverse question when applied to coding efficiency, what's the mathematical question equivalent to asking for the absolutely most bit-squeezingly efficient code for super-super-super...astronomical numbers, after stripping away the veils of bits and efficiency? Something like, the slowest way to approach an asymptote, or the most nearly non-converging sequence? --Steve

1 0

Re: [math-fun] becoming a 21st century mathematician
by Steve Witham 19 Jul '08

19 Jul '08

Jim Propp asked for essential computer tools for mathematicians. I like the X Screen Saver collection: http://www.jwz.org/xscreensaver/ Probably all Linux systems have this, but it's also available for non-X Mac OS machines. I like it just for the "Ant" screen saver, which does random combinations of number of ants, number of states, number of sides for the cells. It keeps doing things I didn't expect and making me wonder what's next. Some tools are important, but this is an essential timewaster. --Steve

1 0