July 2013 - mathfun@mailman.xmission.com

Re: [math-fun] computing non-rigorously defined quantities
by Steve Witham 24 Jul '13

24 Jul '13

>From: Dan Asimov <dasimov(a)earthlink.net> > >What method or formula are you using to get the Taylor coefficients >of expm1^t at x = 1/2, where expm1(x) := e^x - 1 ??? (I don't >know what sympy is or what query you have submitted to it.) I've sent a more detailed answer to Dan. The short story is, I started with N (around 18) terms of the Taylor series of f(x) = e^((x+e)/e) - e P = x + x^2/(2! e) + x^3/(3! e^2) ... I generated M polynomials, for 0 through M-1 iterations of P, each up to N terms. (I'm forgetting whether M needs to be N+1 or N-1 right now.) Visualize an array of the coefficients, with the iterations as rows. Say t is the row number, starting at zero. The column for the constant term is all zero. the coefficient of x is all 1, constant down the column x^2, a linear progression x^3, a quadratic x^4, a cubic ... I inferred polynomials of t from the columns, giving a polynomial in both x and t. I substituted in x = .5, reducing the whole thing to a polynomial in t. The general way to calculate the function is to use the polynomial for t mod 1, then apply either f or its inverse repeatedly as needed. >I ask particularly because I have read that e^x - 1 (and hence >e^(x/e)) does not have any non-integer analytic iterates. See I. N. >Baker, Zusammensetzungen ganzer Functionen, Math. Zeitschr. Bd. 69, >pp. 121-163 (1958). I haven't read the article, haven't even tried complex values of t. For reals less than about 13, it seems well-behaved. --Steve

2 1

[math-fun] Dolphins have and use audible individual names
by Warren D Smith 24 Jul '13

24 Jul '13

http://www.guardian.co.uk/environment/2013/jul/22/bottlenose-dolphins-names…

4 4

[math-fun] series convergence problem
by Warren D Smith 24 Jul '13

24 Jul '13

I guess the answer to this is a known "folk theorem" but have not seen it explicitly stated before PUZZLE: Given a finite-dimensional vector V=(v0,v1,v2,...,vk) define the function F(n) via F(n) = n^v0 * ln(n)^v1 * lnln(n)^v2 * lnlnln(n)^v3 * ... Then: for which vectors V is it true that some N>1 exists such that the series sum(n>N) 1/F(n) converges? -- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)

2 1

Re: [math-fun] 'Rationalize' a vector of real numbers ?
by Adam P. Goucher 23 Jul '13

23 Jul '13

> The question now is how to 'rationalize' a _vector_ of given numbers, but do so with a single denominator. I think that you can use the Lenstra^2-Lovasz lattice basis reduction algorithm. Here's one I made earlier: http://cp4space.wordpress.com/2012/09/12/lunisolar-calendars/ Sincerely, Adam P. Goucher

1 0

[math-fun] Atrocious function plots
by Fred lunnon 23 Jul '13

23 Jul '13

A short while back, I described a purported plot of the classical Beta function at http://mathworld.wolfram.com/BetaFunction.html --- or maybe it was the the very similar offering at https://en.wikipedia.org/wiki/Beta_function --- as "passable identikit". Given that identikit reconstructions are notorious for bearing at best only a remote resemblance to their subjects, the verdict may well have been more appropriate than I could known at the time. The sole noteworthy feature of these hilarious travesties proves to be merely a lattice of vertical planar artifacts, resulting from a misguided attempt to plot across lines of discontinuity, and completely obliterating any structural detail of the actual surface. In the course of drafting an account of the recent discussion of binomial coefficients in this list, I took up Gosper's challenge to produce an accurate plot of the associated function of two real variables; and while I was at it knocked out a plot of the (closely related) Beta function, which I have now posted at https://www.dropbox.com/s/rzcojf8wibd0mpw/beta5.png Compare the competing versions; and weep. Fred Lunnon

1 0

[math-fun] 'degrees of separation' hopscotch
by Henry Baker 22 Jul '13

22 Jul '13

According to the NYTimes, http://well.blogs.nytimes.com/2011/11/22/four-degrees-of-separation/ "Using data on the links among 721 million Facebook users, a team of scientists discovered that the average number of acquaintances separating any two people in the _United States_ was *** 4.37, *** and that the number separating any two people *** in the _world_ was 4.74. *** As John Markoff and Somini Sengupta report in todayÂs New York Times, the findings highlight the growing power of the emerging science of social networks: "The original Âsix degreesÂ finding, published in 1967 by the psychologist Stanley Milgram, was drawn from 296 volunteers who were asked to send a message by postcard, through friends and then friends of friends, to a specific person in a Boston suburb. "The new research used a slightly bigger cohort: 721 million Facebook users, more than *** one-tenth *** of the worldÂs population. The findings were posted on FacebookÂs site Monday night.Â" ... -------- Here's the new twist from yesterday: http://www.theatlanticwire.com/politics/2013/07/nsa-admits-it-analyzes-more… "Chris Inglis, the [NSA] agency's deputy director, was one of several government representativesÂincluding from the FBI and the office of the Director of National IntelligenceÂtestifying before the House Judiciary Committee this morning.... But Inglis' statement was new. Analysts look *** "two or three hops" *** from terror suspects when evaluating terror activity, Inglis revealed. Previously, the limit of how surveillance was extended had been described as *** two *** hops. This meant that if the NSA were following a phone metadata or web trail from a terror suspect, it could also look at the calls from the people that suspect has spoken withÂone hop. And then, the calls that second person had also spoken withÂtwo hops. Terror suspect to person two to person three. Two hops. And now: *** A third hop. *** " ... "For a sense of scale, researchers at the University of Milan found in 2011 that everyone on the Internet was, on average, 4.74 steps away from anyone else. The NSA explores relationships up to *** three *** of those steps." "Inglis' admission didn't register among the members of Congress present, ..." --- This reminds me of a possibly apocryphal story about a debate in Congress over whether Reagan's 'Star Wars' program would work. When one of the scientists testified that some number would have to be something like 10^20 for the system to work, and the best they'd been able to achieve by that date was 10^10, one of the Congressmen remarked "then we're halfway there!"

3 5

[math-fun] Exact coordinate arithmetic for graphics
by Henry Baker 20 Jul '13

20 Jul '13

Although modern graphical systems -- e.g., GPU's -- now have full HW floating point arithmetic, even double precision isn't enough for certain processes -- e.g., 'constructive solid geometry' (CSG), where it's possible to get inconsistent topological answers based on different arithmetic calculations of questions like 'inside/outside', 'intersect/not intersect', etc. A clever 1989 paper by Sugihara & Iri suggested a very interesting solution for CSG by noticing that _all of the planes_ involved in the final answer to a CSG operation were _already present_ in the initial objects; some might be missing, but _no new planes need to be created_. This means that if you start by representing the initial data objects using only planes, _then you never need to compute any new coordinates_, and all CSG operations involve only the shuffling of existing planes. Sugihara, K., and Iri, M.. "A Solid Modelling System Free from Topological Inconsistency". J. Info. Proc., v12#4, 1989, 380-393. Of course, you have to compute one heck of a lot of determinants, all of which are computed only for their sign bit (!), and many of these may require multiple precision to get the answer topologically correct, but you are guaranteed not to lose any precision in the _representation_ of the final answer. Sugihara & Iri represent planes by their classical algebra Ax+By+Cz+D=0 representation, except that they require D to have twice the precision of A,B,C. The overall sign of this equation indicates which side of the plane is being represented. Sugihara & Iri point out that you might as well clear fractions and use integers in this representation, because all of that floating point machinery isn't helping you, since you're going to have to do multiple precision determinant operations anyway. 'Lines' are now just a pair of planes (you don't have to do any arithmetic at all), and 'points' are a triple of planes (still no arithmetic). [There is a technical issue with vertices whose degree isn't exactly 3. In these cases, the 'planes' representation doesn't worry about vertex degree, because the topology only worries about planes and intersections. If/when you convert from this 'planes' representation into 'vertex' representation, there may be some vertices which are incredibly close together. But this situation is completely dual to the issue of faces with more than 3 edges and which are slightly non-planar as a result of imperfect arithmetic.] Sugihara & Iri also point out that their representation is the _dual_ of the more traditional _vertex_ representation, where points are the primary data, lines are pairs of points, and planes are triples of points. However, S&I don't carry through on the complete duality, because they still treat vertex points as triples, whereas most graphics systems now represent points as [x y z 1] in 'homogeneous' coordinates, which handle both rotational and translational transformations more elegantly. But if you have homogeneous coordinates anyway, then you might as well clear fractions, and represent points as [D*x D*y D*z D], where D is the denominator. But now that you have integers, you can convert to/from the 'planes' representation without loss of precision. Some operations, like CSG ops, seem better suited for the 'planes' representation, while other graphics ops are better suited for the 'vertices' representation; it is nice to be able to move back & forth without loss of precision. One question: what is the appropriate meaning for a negative denominator in [x y z D] ? In the 'planes' representation, the sign tells us which side of the plane we're on, but I'm stumped about what the overall sign of the point [x y z D] should mean. (Points and anti-points ??) An answer to this question would solve a problem I've had for 50 years: what is the meaning of a negative denominator in a rational number ? Second question: GCD's. Typical rational arithmetic, as found in Common Lisp/Maxima/etc., involves huge numbers of GCD operations, which are the limiting factor in the use of these numbers. With _four_ numbers instead of just two, GCD's are even less likely to be unequal to 1, yet GCD=1 is the worst case scenario for standard GCD implementations. So you probably want to defer GCD's until the numbers get to be perhaps twice as big as you would expect from a precision analysis.

2 1

Re: [math-fun] computing non-rigorously defined quantities
by Steve Witham 20 Jul '13

20 Jul '13

> > At 08:41 AM 7/2/2013, James Propp wrote: >>> For some thoughts about the enterprise of computing quantities one >>> hasn't rigorously defined, see >>> >>>http://mathoverflow.net/questions/135536/procedure-based-as-opposed-to-defi… > >> (which I hope will elicit interesting responses over the next few days). >Date: Tue, 02 Jul 2013 14:16:11 -0600 >From: rcs(a)xmission.com > >Some of the funsters' flow work fits here: Hah! "Fits," but comfortably? >Another, even more surprising example, is f(x) = sqrt2^x. This has fixed >points at x = 2 & 4, and is well-behaved (nearly linear, smooth, etc.) >near the fixed points and in the interval between. sqrt(2) is close to B = e^(1/e) = ~1.444667 whose exponential has one ("double?") fixed point, B^e = e and, you know, B^x is an honest exponential. If you scale and translate B^x a bit, (B^(ex + e) - e)/e = e^x - 1 = expm1(x) = x + x^2/2! + x^3/3! + ... expm1(0) = 0 If you set the initial x to (say) 1/2, then the flow, expm1^t(1/2) has a Taylor series (of t) with nice rational coefficients and an interval around zero where it converges. Except I'm getting the coefficients from sympy and don't see a pattern beyond knowing what I asked it. expm1 is not a true exponential, but the further up the flow you *start*, the closer it fits the "towers" of e. (log1p(x) = ln(1 + x), inverse of expm1.) + -inf ^ + 0.61301989 ^ | exp ln | | expm1 log1p | v -0.00000000 ^ v 0.84599770 ^ | exp ln | | expm1 log1p | v 1.00000000 ^ v 1.33030161 ^ | exp ln | | expm1 log1p | v 2.71828183 ^ v 2.78218395 ^ | exp ln | | expm1 log1p | v 15.15426224 ^ v 15.15426250 ^ | exp ln | | expm1 log1p | v 3814279.10476017 +----v 3814279.10476017 + | | +------------------------------------------+ The outer and inner U's, including turning around short of infinity, are procedures to map between, not only the two sequences, but between any point on the expm1^t flow and the implied exp^t flow. The difference between the two columns at x is like e^-x, so you get like e^^n bits of precision in n steps, which you can carry down the other side (er up, in the diagram). There isn't a simple mapping between flows of exponentials of different bases, but this technique makes flows for e and e^(1/e), and I think with some tweaks, for any base (including yet another answer for sqrt(2)!). So, I start with a function with a single fixed point at zero and extract something from that, but then go almost out to infinity and back a different way to cook it. And the mapping must make a complete hash of complex numbers. And even if there's a way to fix that... Personally I'm not happy with just these procedures or numbers or curves, although they're something I've been able to do and they're pretty. I wish they were somehow close or related to something with a real reason to live. I wonder about Beauty in math, the smell that's supposed to tell you you're on to something. --Steve + -inf ^ + 0.61301989 ^ exp | | ln expm1 | | log1p v -0.00000000 ^ v 0.84599770 ^ exp | | ln expm1 | | log1p v 1.00000000 ^ v 1.33030161 ^ exp | | ln expm1 | | log1p v 2.71828183 ^ v 2.78218395 ^ exp | | ln expm1 | | log1p v 15.15426224 ^ v 15.15426250 ^ exp | | ln expm1 | | log1p v 3814279.10476017 +------------------v 3814279.10476017 + | | +--------------------------------------------------------+

2 1

Re: [math-fun] exact arithmetic for geometry
by Henry Baker 20 Jul '13

20 Jul '13

Warren, you are correct, but this is far easier said than done. See a good summary here: Yap, Chee, and Dube, Thomas. "The Exact Computation Paradigm". 1994? (Chapter in "Computing in Euclidean Geometry", by Du and Hwang, 1994.) Available in pdf on the Internet. At 03:04 PM 7/19/2013, Warren D Smith wrote: >a good solution to the geometry consistency problem is to use exact arithmetic. >More precisely, you use ordinary lo-precision arithmetic except when it gets >dicey, then you go for more precision... the net effect is, hi-precision calculations are >rarely needed, so it is efficient. > >This is a pain, but it solves everything.

1 0

Re: [math-fun] primes - {2} + {1}? Not.
by Bill Gosper 20 Jul '13

20 Jul '13

How the heck did RWG come across it? WFL On 7/19/13, Bill Gosper < billgosper(a)gmail.com> wrote: What the heck use is A006005, "The odd prime numbers together with 1"? Claim: the presence of 1 and the absence of 2 suggest these are actually "Numbers 2n+1 such that 0=Mod[Binomial[2n,n]-(-1)^n,2n+1]" See Concrete Math (G,K&P) Exercise 5.110 I.e., A065091 ∪ A163209 ∪ {1}, suggesting that A163209 should really be *1,5907*, *1194649*, 12327121, "Catalan pseudoprimes: (odd) nonprime integers n=2*m+1 satisfying A000108< http://oeis.org/A000108>(m) == (-1)^m * 2 (mod n). " A176997, " Odd numbers n such that (2^(n-1) mod n) = (4^(n-1) mod n) = (8^(n-1) mod n)= .. = (k^(n-1) mod n) for k=2, 4, 8, .., smallest power of 2>n" should be extended through 341 to avoid matching A006005. (Is this not equivalent to "Numbers n such that 1 = (2^(n-1) mod n) = (4^(n-1) mod n) = (8^(n-1) mod n)= .. = (k^(n-1) mod n) for k=2, 4, 8, .., smallest power of 2>n "?) Likewise, 0, 1, 2, 3, 5, 6, 8, 9, 11, 14, 15, 18, 20, 21, 23, 26, 29, 30, 33,... A102781 "Number of even numbers less than the n - th prime." is probably less useful than "Numbers n such that 0=Mod[Binomial[2n,n]-(-1)^n,2n+1]" which also include 2953, (1093^2 - 1)/2, and (3511^2 - 1)/2. --rwg See math-fun-archive-931123-950622.txt excerpt below. Anybody know how far the Catalan pseudoprime search has been pushed? E.g., past the Wieferich search? http://gradelle.educanet2.ch/christian.aebi/.ws_gen/9/catalan.pdf gives two stronger Catalan pseudoprime tests: 0 == Mod[(-1)^((p - 1)/2)* CatalanNumber[(p - 1)/2] - 2 (-1 + 2^p) (1 - p), p^2] 0 == Mod[(-1)^((p - 1)/2)*CatalanNumber[(p - 1)/2] - (1 - p + p^2)*2^(2*p - 1), p^3] AFAIK, *no* false positives are known! Can the sexy trick for computing binomial(n,k) mod p by looking at the carries when adding k to n-k in base p be turned into an efficient probable prime test? Can the trick be extended to mod p^2 and mod p^3 to facilitate the Catalans in the above identities? Failing that, how about just 0 == Mod[(-1)^((p - 1)/2)* CatalanNumber[(p - 1)/2] - 2^p , p] (after checking for 5907, 1093^2, and 3511^2)? --rwg Why doesn't Mma Binomial take a Modulus option? Macsyma has binomial_modulus. I forget if it can Chinese remainder for composite modulus. >From rwg(a)NEWTON.Macsyma.COM Wed May 31 01:11:00 1995 Received: from optima.CS.Arizona.EDU by cheltenham.CS.Arizona.EDU; Wed, 31 May 1995 01:20:23 MST Received: from newton.macsyma.com by optima.cs.arizona.edu (5.65c/15) via SMTP id AA24896; Wed, 31 May 1995 01:20:21 MST Received: from SWEATHOUSE.macsyma.com ([192.233.166.105]) by NEWTON.Macsyma.COM via INTERNET with SMTP id 320962; 31 May 1995 04:11:10-0400 Date: Wed, 31 May 1995 01:11-0700 From: Bill Gosper <rwg(a)NEWTON.macsyma.com> Subject: q-binom mod 1-q^p To: math-fun(a)cs.arizona.edu Cc: mrahman(a)math.carleton.ca, george(a)math.nwu.edu, askey(a)math.wisc.edu, math-fun(a)cs.arizona.edu Message-Id: <19950531081123.8.RWG(a)SWEATHOUSE.macsyma.com> n Recall that ( ) mod p can [be] very quickly found from the carries when adding k k to n-k in base p, when p is a prime. (In fact, doing this for several p and then Chinese remaindering is a good way to compute binomial coeffs.) Recall also the q-binomial coeff: k /===\ n - k + j n | | 1 - q [ ] := | | --------------, k | | j j = 1 1 - q which becomes the vanilla one when q -> 1. 2 p n What is known about the remainder mod 1-q^p ? E.g., for [ ] p n it's 2 n p 2 n ( ) - ( ) p n p n 1 - q 2 n --------------- ------ + ( ), p 1 - q n which gives the sizes of the residue classes when p n distinct positive integers <= 2 p n are summed mod p. E.g., the probability that the sum of 1000 distinct integers from [1,2000] will be divisible by 5 is 2 n ( ) 1 n 1 (1 - -) --------- + - p 2 n p p ( ) n p with p = 5 and n = 200, roughly 1/5 + 4.021284314257178521246*10^-482. The other four residues are equally probable, very roughly 1/5 - 10^-482. >From 75240.125(a)compuserve.com Wed May 31 05:44:10 1995 Received: from optima.CS.Arizona.EDU by cheltenham.CS.Arizona.EDU; Wed, 31 May 1995 02:45:25 MST Received: from dub-img-2.compuserve.com by optima.cs.arizona.edu (5.65c/15) via SMTP id AA27641; Wed, 31 May 1995 02:45:21 MST Received: by dub-img-2.compuserve.com (8.6.10/5.950515) id FAA11025; Wed, 31 May 1995 05:45:07 -0400 Date: 31 May 95 05:44:10 EDT From: "Fred W. Helenius" <75240.125(a)compuserve.com> To: math-fun <math-fun(a)cs.arizona.edu> Subject: Re: doubling mod P Message-Id: <950531094409_75240.125_HHB42-1(a)CompuServe.COM> > I think that the search for numbers like 1093 & 3511 has reached > at least 1G, but don't know details. > > I think I've heard these numbers named after Wieferich. Even more > exclusive is 3*11*179, the only known pseudoprime satisfying > > binom(p-1,(p-1)/2) == (-1)^((p-1)/2) mod p. I've found two other members of this exclusive club. They turn out to be the squares of the two Wieferich primes: 1093^2 = 1194649 and 3511^2 = 12327121. An outline proof: (1) If q is a Wieferich prime, then binom(q-1,(q-1)/2) == (-1)^((q-1)/2) mod q^2. This is an immediate consequence of Hardy & Wright's Theorem 133. (2) If q is an odd prime, then binom(q^2-1,(q^2-1)/2) == (-1)^((q-1)/2) binom(q-1,(q-1)/2) mod q^2. This follows from writing the first binomial as the product of terms (q^2 - k)/k, 1 <= k <= (q^2-1)/2. For q(q-1)/2 of the terms, the numerator and denominator are free of q, so each of these terms reduces to -1 mod q^2; the remaining (q-1)/2 terms form the second binomial. Combining (1) and (2), for a Wieferich prime q, binom(q^2-1,(q^2-1)/2) == 1 mod q^2, which is the original condition with p replaced by q^2. >From rwg(a)NEWTON.Macsyma.COM Thu Jun 1 00:39:00 1995 Received: from optima.CS.Arizona.EDU by cheltenham.CS.Arizona.EDU; Thu, 1 Jun 1995 01:06:08 MST Received: from newton.macsyma.com by optima.cs.arizona.edu (5.65c/15) via SMTP id AA17460; Thu, 1 Jun 1995 01:06:05 MST Received: from SWEATHOUSE.macsyma.com ([192.233.166.105]) by NEWTON.Macsyma.COM via INTERNET with SMTP id 321147; 1 Jun 1995 03:39:30-0400 Date: Thu, 1 Jun 1995 00:39-0700 From: Bill Gosper <rwg(a)NEWTON.macsyma.com> Subject: Re: doubling mod P To: 75240.125(a)compuserve.com Cc: math-fun(a)cs.arizona.edu In-Reply-To: <950531094409_75240.125_HHB42-1(a)CompuServe.COM> Message-Id: <19950601073915.9.RWG(a)SWEATHOUSE.macsyma.com> > I think I've heard these numbers named after Wieferich. Even more > exclusive is 3*11*179, the only known pseudoprime satisfying > > binom(p-1,(p-1)/2) == (-1)^((p-1)/2) mod p. I've found two other members of this exclusive club. They turn out to be the squares of the two Wieferich primes: 1093^2 = 1194649 and 3511^2 = 12327121. Gack, retroactive deja vu. I think this may have been in Ilan Vardi's Mathematica book.

1 0