July 2017 - mathfun@mailman.xmission.com

[math-fun] MROB now lets you build ries to include LambertW
by Bill Gosper 15 Jul '17

15 Jul '17

In[486]:= FixedPoint[2 + E^-# &, .0] Out[486]= 2.12002823898764 billgosper$ ries/ries -l6 -EW 2.12002823898764 Your target value: T = 2.12002823898764 mrob.com/ries x^x = 5 for x = T + 0.00934424 {53} x-2 = 1/9 for x = T - 0.00891713 {59} ln(x)^2 = W(1) for x = T + 0.00352165 {62} x^4 = e^3 for x = T - 0.00302822 {65} x/3 = 1/sqrt(2) for x = T + 0.0012921 {64} ln(x)^2 = 1/sqrt(pi) for x = T - 0.000644103 {67} e"/x = 1/pi+1 for x = T - 0.000489232 {73} x-1/x = sqrt(e) for x = T + 0.000314826 {67} sqrt(x-2) = ln(sqrt(2)) for x = T + 8.50145e-05 {77} sqrt(x/pi) = sinpi(ln(2)) for x = T + 2.54297e-05 {82} x/sqrt(phi) = 5/3 for x = T + 4.5102e-06 {84} sinpi(cospi(1/x)) = e/pi^2 for x = T + 1.73152e-06 {92} x"/(x-1) = 9"/phi for x = T - 3.21311e-07 {96} x-2 = W(1/e^2) for x = T + 1.33227e-15 {80} (Stopping now because best match is within 1.88e-15 of target value.) W(x) = LambertW(x) = inverse(x=w*e^w) phi = the golden ratio, (1+sqrt(5))/2 sqrt(x) = square root A"/B = Ath root of B pi = 3.14159... (I.e., 2.12002823898764 = 2 + W(e^-2).) But, thanks to El Capitan (which also broke my Macsyma)-:, billgosper$ gcc -o ries ries.c -DRIES_USE_SA_M64 -lm xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun until: billgosper$ xcode-select --install --rwg Now if ries just had γ and the Γ function . . . it would probably be too slow!-(

1 0

[math-fun] Can someone please explain how Somos4,
by Bill Gosper 14 Jul '17

14 Jul '17

initialized to 1,1,-1,x, Nest[Append[Rest[#], Factor[(#[[4]] #[[2]] + #[[3]]^2)/prnt@#[[1]]]] &, {1, 1, -1, x}, 17]; 1 1 -1 x 1+x -1-x+x^2 -1-x-x^3 -x (2+3 x) 1+3 x+3 x^2-x^4+x^5 -(1+x) (-1-2 x+2 x^2+3 x^3-x^4+x^5) -1-3 x-3 x^2-5 x^3-9 x^4-3 x^5+2 x^6-x^7 -x (-1-x+x^2) (3+9 x+9 x^2+5 x^3+2 x^4+x^6) 1+6 x+15 x^2+16 x^3-6 x^5+12 x^6+13 x^7+x^9+x^10 -(1+x+x^3) (1+5 x+4 x^2-15 x^3-30 x^4-13 x^5+6 x^6+x^7-3 x^8+x^9) -(1+x) (1+5 x+10 x^2+20 x^3+50 x^4+63 x^5+12 x^6-19 x^7+8 x^8+5 x^9-2 x^10+5 x^11-2 x^12+x^13) x (2+3 x) (-2-12 x-27 x^2-27 x^3-6 x^4+18 x^5+23 x^6+3 x^7-2 x^8+14 x^9+6 x^10-3 x^11+2 x^12) 1+10 x+45 x^2+110 x^3+135 x^4+45 x^5-2 x^6+211 x^7+429 x^8+311 x^9+105 x^10+69 x^11+56 x^12+5 x^13-9 x^14+6 x^15+4 x^16-3 x^17+x^18 . . . never manages to produce a denominator? (Very) similarly for 1,i,1,x: Nest[Append[Rest[#], Factor[(#[[4]] #[[2]] + #[[3]]^2)/prnt@#[[1]], Extension -> I]] &, {1, I, 1, x}, 17] 1 i 1 x 1+i x -i+x-i x^2 1+i x-i x^3 -i x (-2 i+3 x) 1+3 i x-3 x^2-x^4+i x^5 (1+i x) (i-2 x+2 i x^2-3 x^3+i x^4+x^5) 1+3 i x-3 x^2-5 i x^3+9 x^4+3 i x^5+2 x^6-i x^7 -x (1+i x+x^2) (-3-9 i x+9 x^2+5 i x^3-2 x^4+x^6) 1+6 i x-15 x^2-16 i x^3-6 i x^5-12 x^6-13 i x^7+i x^9-x^10 (1+i x-i x^3) (-i+5 x+4 i x^2+15 x^3+30 i x^4-13 x^5+6 i x^6-x^7+3 i x^8+x^9) -(-i+x) (-i+5 x+10 i x^2-20 x^3-50 i x^4+63 x^5+12 i x^6+19 x^7-8 i x^8+5 x^9-2 i x^10-5 x^11+2 i x^12+x^13) x (-2 i+3 x) (-2 i+12 x+27 i x^2-27 x^3-6 i x^4-18 x^5-23 i x^6+3 x^7-2 i x^8-14 x^9-6 i x^10-3 x^11+2 i x^12) 1+10 i x-45 x^2-110 i x^3+135 x^4+45 i x^5+2 x^6-211 i x^7+429 x^8+311 i x^9-105 x^10-69 i x^11+56 x^12+5 i x^13+9 x^14-6 i x^15+4 x^16-3 i x^17-x^18 --rwg

3 2

[math-fun] Any cool hacks for Gaussian elimination over GF(2) ?
by Bill Gosper 13 Jul '17

13 Jul '17

Date: 2014-03-31 02:46 The Symbolics Lisp Machine (and many of its contemporaries) had a BITBLT (Bit block transfer) instruction capable of Boolean functions on rectangular bit arrays. Our usual test subject was a 1 bit/pixel frame of a Zippy the Pinhead strip. I wrote a little demo called Thrash-Zippy which iteratively subjected this to various BITBLTcherous flips and rotations, the most interesting of which was a bitmatrix transpose (swap-xy), which took k1 lg n for n×n, and used lg n n×n masks, each of which I'm pretty sure took k2 lg n to compute. The kth mask designated all those bits which needed to switch with their counterparts 2^k to their northwest. I recall it startling Leo Guibas, who was writing a computational geometry book at the time, so maybe he wrote it up. Funster Stephen Jones has copies of my old files. I think the functions were in a file named bltfun. --rwg [ On 2017-07-02 17:42, Henry Baker wrote: > The 3 shears procedure is more exactly: > > 1. Shear the bitmap horizontally 45 degrees to the right; the X-axis > stays the same, the Y-axis tilts 45 degrees to the right. > 2. Shear the bitmap vertically 45 degrees down; the Y-axis drops down > 45 degrees to horizontal, the X-axis drops down 45 degrees below > horizontal. > 3. Shear the bitmap horizontally 45 degrees to the left; the Y-axis > remains horizontal, the X-axis shifts 45 degrees to become vertical. > > At 01:51 PM 7/2/2017, Henry Baker wrote: >>If you have high speed raster graphics hardware, you may already have the capability for 2D transpose. >> >>There's also the famous (at least in computer graphics circles) "three shears" algorithm for image rotation by 90 degrees without any loss of information: >> >>1. Shear the image 45 degrees to the right. >>2. Shear the image 90 degrees down. >>3. Shear the image 45 degrees to the left. >> >>The problem is that step #2 is usually just as hard in most memory systems as transposition, unless the hardware is set up to pick up word slices -- e.g., give me the 13th bit each of 32 words arranged as a 32-bit word. >> >>At 01:10 PM 7/2/2017, rcs(a)xmission.com wrote: >>>Transposition can be done in log2(N) passes over the matrix. >>>Imagine cutting your matrix into quarters, and exchanging the upper right >>>quarter with the lower left quarter. Next, do the same thing within each >>>quarter. Etc. If N isn't a power of 2, you may have to round it up; unsure. >>> >>>Rich >>> >>>----- >>>Quoting Henry Baker <hbaker1(a)pipeline.com>: >>>>Thanks, Adam. >>>> >>>>I'm already doing some of what you suggest. >>>> >>>>Is there any cleverness for *transposing* such bit matrices ? >>>> >>>>(As an aside, do *least squares* have any relevance over GF(2)?) >>>> >>>>Note that Gaussian elimination on extremely long vectors (much >>>>longer than the rank of the system) is faster if you invert the >>>>small square (rank x rank) matrix by itself, and then propagate that >>>> info to the rest of each row. >>>> >>>>Is there any cleverness that might be achieved if one had access to >>>>a *polynomial* multiply (and/or divide) operation over GF(2)[x] ? I >>>> seem to recall that this sort of HW is finally starting to show up >>>>as a "carryless multiply" operation. >>>> >>>>At 11:06 AM 7/2/2017, Adam P. Goucher wrote: >>>>>Yes, if the matrix has <= 64 rows, then you can store the columns as >>>>>64-bit integers and perform column operations using bitwise XOR. Also, >>>>>you can determine the uppermost and lowermost non-zero elements in a >>>>>column using the __builtin_clzll() instruction and variants thereof. >>>>> >>>>>This generalises, of course; with N rows, you can store each column >>>>>as a vector of ceil(N/64) 64-bit integers and perform column operations >>>>>with as many instructions. With SIMD, you can actually get away with as >>>>>few as ceil(N/256) instructions per column operation. >>>>> >>>>>The fact that there's only one non-zero element of GF(2) is also >>>>>helpful since there's no need to multiply the columns -- add/subtract >>>>>(i.e. XOR) is the only operation you ever need to perform! >>>>> >>>>>I'm using exactly this procedure in my PhD research on topological data >>>>>analysis. In particular, computing the persistent homology barcode from >>>>>the boundary map is (a slightly modified version of) Gaussian elimination >>>>>over GF(2). >>>>> >>>>>Best wishes, >>>>> >>>>>Adam P. Goucher >>>>> >>>>>>Sent: Sunday, July 02, 2017 at 6:29 PM >>>>>>From: "Henry Baker" <hbaker1(a)pipeline.com> >>>>>>To: math-fun(a)mailman.xmission.com >>>>>>Subject: [math-fun] Any cool hacks for Gaussian elimination over GF(2) ? >>>>>> >>>>>>Other than Strassen? >>>>>> >>>>>>I.e., are there any bit-twiddling hacks that are particular to >>>>>>GF(2) for Gaussian elimination? >

1 0

[math-fun] What is the probability that
by Bill Gosper 13 Jul '17

13 Jul '17

[a 32x32 array of random bits over GF(2) is singular?] Responding to Adam Goucher's praise of a sequence of Baez posts expositing q-Hypergeometric applications, I said One goodie that Baez skipped: If an m×n bitmatrix has 1s with probability p and 0s with probability q :=1-p, the matrix will have (mod 2) rank k with probability P[m,n,k] := p^k q^((k - m) (k - n)) QBinomial[m, k, q] QBinomial[n, k, q] QFactorial[k, q] . How could Baez resist presenting yet this other pun on q? I *must* have said this to math-fun because I later(?) said that (eavesdropper) George Andrews reminded me that the identity Sum[P[m,n,k],{k,Min[m,n]}] = 1 is a special case of q Chu-Vandermonde. But I can't find where I said it. m - k + 1 k + 1 (m - k) (n - k) (q ; q) (q ; q) q k n - k P(m, n, k) = --------------------------------------------------. (q; q) n - k -m -n m n + k (q , q ; q) q k = ----------------------- (q; q) k = P(n, m, q). On 2017-07-12 09:31, Eugene Salamin via math-fun wrote: > Consider the n×n matrices over GF(q). There are q^(n^2) of them. > Lets count how many are nonsingular. The first row must not be all > zero, so (q^n - 1) possibilities. Suppose we've selected the first k > rows. Being linearly independent, they span a k-dimensional space. > So there are (q^n - q^k) choices for the (k+1)-th row. Thus the > number of nonsingular matrices is > > (q^n - 1) (q^n - q) (q^n - q^2) ... (q^n - q^(n-1)). > > There are n factors in this product. We divide each factor by q^n to > get the probability that the matrix is nonsingular. > > (1 - 1/q) (1 - 1/q^2) (1 - 1/q^3) ... (1 - 1/q^n) > > The product converges nicely as n goes to infinity. For Henry's > question, q = 2, n = 32, but n hardly matters, and the probability of > being nonsingular is > > (1/2) (3/4) (7/8) ... = 0.288788, and 0.711212 for being singular. > > q P(nonsing) P(sing)2^2 0.688538 0.311462 > 2^3 0.859406 0.140594 > 2^4 0.933595 0.0664053 0.560126 0.439874 > 3^2 0.876560 0.123440 > 3^3 0.961591 0.0384095 0.760333 0.239667 > 5^2 0.958400 0.041600 > 7 0.836795 0.163205 > > -- Gene > > On Wednesday, July 12, 2017, 7:44:36 AM PDT, Henry Baker > <hbaker1(a)pipeline.com> wrote: > > a 32x32 array of random bits over GF(2) is singular? > > It appears to be non-negligible -- i.e., perhaps 10% ?? > > (This is mere eyeballing; I haven't done a serious Monte Carlo estimate.) > > My initial intuition was that it should have been negligible. >

1 0

[math-fun] normal bases for GF(2^32)
by Henry Baker 12 Jul '17

12 Jul '17

I'm using the industry standard crc32 polynomial as my reduction polynomial. There appear to be only 6 monomials that generate normal bases: x^9, x^15, x^22, x^23, x^30, x^31. Is this number 6 approx. O(log 32), or is this number 6/32=19% "approximately" 100%-71%=29%. We can't worry about "optimal normal bases" here, because GF(2^32) can't have an optimal basis: * 32 can't have an optimal normal basis of type 1, since 2 is not a primitive element of GF(33). * 32 can't have an optimal normal basis of type 2, since 2*32+1=65 isn't prime.

1 0

[math-fun] What is special about 323 and squares?
by Ed Pegg Jr 12 Jul '17

12 Jul '17

The minimal number of squares for rectangles up to longest side 380 is known. http://int-e.eu/~bf3/squares/view.html#319,323 needs a minimum of 18 squares. f(323,319)=18, f(323,293)=17, f(323,317)=16, f(323,283)=16, f(323,281)=16, Those are all of the cases up to larger side 380 that need more than 15 squares if the aspect ratio is less than 2. For 15 squares, add the value 352 as hard. f(323,X)=15, with X in {256, 271, 277, 307, 313} f(352,X)=15, with X in {283, 289, 293, 299, 307, 311, 317, 325, 329, 331, 333, 343, 347, 349, 351} For rectangles needing 14 squares, more than half the values have larger side 323 or 352. What is special about 323 (and 352) and squares? --Ed Pegg Jr

2 7

[math-fun] What is the probability that
by Henry Baker 12 Jul '17

12 Jul '17

a 32x32 array of random bits over GF(2) is singular? It appears to be non-negligible -- i.e., perhaps 10% ?? (This is mere eyeballing; I haven't done a serious Monte Carlo estimate.) My initial intuition was that it should have been negligible.

3 2

Re: [math-fun] What is the probability that
by Henry Baker 12 Jul '17

12 Jul '17

Indeed. Out of 10 million trials, I just got 7,112,737 singular. My intuition is badly out of whack wrt matrices over GF(2). Random matrices over larger fields -- e.g., reals -- are almost never singular. So now what is the probability that a random non-zero element over GF(2^32) does NOT generate a normal basis? I guess what I'm asking is whether the problem in generating a normal basis is getting a basis at all, or whether it has to do with something else. In other words, what's the difference between foo^(2^i), i=0..31 forming a basis and 32 randomly selected elements of GF(2^32) forming a basis ? At 08:14 AM 7/12/2017, rcs(a)xmission.com wrote: >Assuming you mean "singular over GF[2]", about 71%. > >The chance of being non-singular is 29%, or > >1/2 * 3/4 * 7/8 * 15/16 * ... > >cutting off the infinite product after 1 - 1/2^32. > >The idea: >Assume the space spanned by the first 31 rows is full-rank, >as large as possible. This is 2^31 vectors. What's the chance >that the bottom row is in this set, making the matrix singular? >50%. >What's the chance that row 31 is in the space spanned by the >first 30 rows (assuming the first 30 rows are full-rank)? >25% >Etc. > >This is a standard infinite product, equal to > >1 - 1/2 - 1/4 + 1/32 + 1/128 - 1/4096 - 1/32768 + ... > >The exponents of 2 are pentagonal numbers, n(3n+-1)/2, and the >signs are alternating pairs (after the initial 1). >The sum is nice in binary. >There's a slightly messier expression for the partial products. >The same idea works in GF[p]. Or mod P, for us lowbrow folks. > >Rich > >--- >Quoting Henry Baker <hbaker1(a)pipeline.com>: >>a 32x32 array of random bits over GF(2) is singular? >> >>It appears to be non-negligible -- i.e., perhaps 10% ?? >> >>(This is mere eyeballing; I haven't done a serious Monte Carlo estimate.) >> >>My initial intuition was that it should have been negligible.

1 0

[math-fun] square roots in GF(2^k) ?
by Henry Baker 10 Jul '17

10 Jul '17

I did a quick Google search on discrete square roots, but most refs were for GF(p). Is the complexity of taking a square root in GF(2^k) any less than taking an arbitrary discrete log? How about the simple question of whether an element *has* a square root? The reason I ask this question is due to the fact that: (a+b)^2 (mod pp) = a^2 + 2ab + b^2 = a^2+b^2 (mod pp), where pp is the irreducible polynomial for the field GF(2^k), and a,b are representative polynomials whose degree is < deg(pp)=2^k. So, (a*x^j+b)^2 = (a*x^j)^2 + b^2 = a^2*x^2j + b^2 So Karatsuba is even simpler/faster in these fields, hence, it might seem that finding such a square root might be easier than in some other fields.

2 1

[math-fun] Algebraic number question
by David Wilson 10 Jul '17

10 Jul '17

By happenstance, I noticed that v8 + v5 = 5.064495+ is just over 5, so I computed 1 / (v8 + v5 - 5) = 1/4 * (5v10 + 7v5 + 11v2 + 15) = 15.505053+ I noticed that the terms 5v10 = v250 = 15.811388+ 7v5 = v245 = 15.652475+ 11v2 = v242 = 15.556349+ 15 = v225 = 15.000000 are all in the same ball park. Is there some reason this should be the case? Can other similar examples be formulated?

2 4