triplecorrelation.com

Did you know that

and that

These are neat facts that I picked up from reading Borwein & Borwein’s paper in the recent American Math Monthly. What other functions have that property?

That is the beautiful title of a paper by Vinay Deolalikar at HP Labs Palo Alto.  His proof is being reviewed by experts, so it is probably too early to get excited about the result. However, I found it fascinating that he used statistical arguments from Markov random field theory that I remembered from Geman & Geman’s famous 1984 paper.  It’s inspiring that he worked out the proof on his own, and on his own free time apparently.  HP Labs is not like AT&T, where pure research is encouraged.

Next week I am going to the International Congress of Mathematicians, where the Fields Medal will be awarded.  They still haven’t announced the winners, perhaps saving it for the conference.  It seems amazing to say this, but Deolalikar’s quiet announcement may have upstaged the Fields medal. The Fields medals have been upstaged before: by Perlman, who refused to show up despite proving the Poincare conjecture, and Wiles, who was stiffed by the arcane rules that the winner be under 40.

By the way if you write “P”, then “&”, “#”, followed by “8800 NP” without the quotes, you get the title of Deolalikar’s paper in HTML.

Update (14 August): Deolalikar has posted a summary of his proof on his HPL homepage.   He shows the essence of the proof is that k-SAT for k>9 grows exponentially, and thus there is no polynomial time algorithm.  K-SAT is the problem of checking whether a boolean expression ( X1X2(~X3)+X2(~X1)X3+…) with at most K literals per clause.

When I looked up P=NP on Google’s image search, I found this

On Homer’s right is P=NP.  On his left is Euler’s result, which when written

is the most beautiful formula ever.

I read in Malcom Line’s book (thanks, Google Book Search!)  a nice discussion of repeating digit primes.  It goes like this.

We know that 11 is a prime number. Are there any other primes with only a single repeating digit, i.e., of the form 111…11? Since 11*101=1111, 1111*101=111111 (six digits), we can see that any even length repeating digit number is not a prime. Likewise, since 111 is not prime, neither is 111111111 (nine digits) or any number with 3n digits, since we can use a construction like 111111*1001001=111111111. Using this line of reasoning we can see that any repeating digit number whose length is not prime cannot be prime either. So that leaves only prime lengths. We know by using MATLAB’s “factor” command that (using monetary notation to group digits into thousands)

11,111 = 41*271 % 5 digits is not prime

1,111,111 = 239*4,629 % neither is 7 digits

The “factor” command is limited to 2^32, or 10 digits, which is not enough for the 11 digit repeater. However, we can write a simple MATLAB script to find that

11,111,111,111 = 21,649 * 513,239

Similarly, we can find that the 13 digit repeater factors into three primes

1,111,111,111,111 = 53 * 79 * 26,5371,653

The 17 digit repeater requires more work, since it exceeds the 4 byte integer math capability of MATLAB. We have to devise a long division routine. We can take advantage of the 15 digit precision of a IEEE 754 “double”, which MATLAB uses for floating point calculations, and divide the 17 digit string starting from the left most 15 digits. This gives

11,111,111,111,111,111 = 2,071,723 * 10,726,419,557

One way to check this is to use the “multiplystrings” script below.

How about the 19 digit repeater? It is said to be prime in Lines’s book  (see above).  I don’t have a script to check it yet, though I know it has no prime factors with less than 9 digits. If you know of a good primality script in matlab or octave, please let me know.

We all know that convolution and multiplication are Fourier transforms of each other as operations.  It may be less well known that multiplication of numbers is essentially convolution of their digits.  For example, if we multiply “123″ by “321″, the result “39483″ is also obtained by convolving vectors [1,2,3] with [3,2,1].  In MATLAB or Octave, we can see this by

result =  conv([1,2,3],[3,2,1])

We get that “result” is [3,8,14,8,3].   If we scan the vector from right to left and carry any digits in the 10-s or higher places leftward, we get the answer [3,8+1,4,8,3] = [3,9,4,8,3].   

There is really no better way than convolution to multiply large integers, particularly those that exceed the bit depth available for standard arithmetic.  Consider mutiplying ‘123456789′ with ‘987654321′.  In MATLAB, we get

123456789*987654321 = 1.2193e+017

If we instead multiply by convolving the digits with carry, we get the actual result:

 121932631112635269.

 Here is a MATLAB script that implements multiplication of arbitrary length strings. 

function result = multiplystrings(a,b)
% uses convolution to multiply the two strings a and b
% for example
% result = multiplystrings(’123′,’321′)
% gives
% result = ‘39843′;
% This is a great way to multiply large integers
% r kakarala
na = length(a); nb = length(b);

for k = 1:na
anum(k) = a(k)-’0′;
end; % make “a” a vector

for k = 1:nb
bnum(k) = b(k)-’0′;
end; % and “b” as well

resultnum = conv(anum,bnum); % multiply!
nr = length(resultnum);
carry = 0;

for k = nr:-1:1 % convert it back to a string
v = resultnum(k) + carry;
d = rem(v,10); % separate digit and higher place digits
carry = (v-d)/10;  % this handles carries greater than a single digit as well
result(k) = char(d+’0′); %convert back to a character
end;
if (carry > 0)
result = [num2str(carry) result];
end;

Here are some non-obvious things I learned about Octave.

1) It is easy to generate all N! permutations of the N elements of a vector X. Use

>P = perms(X)

The permutations don’t come out in lexicographic order unfortunately.

2) Typing “whos” gives a list of all variables and built-in functions. By default, Octave uses the UNIX “less” pager to show output. This has the unfortunate effect of clearing the screen after the last item has been displayed (and you press “f”). You can avoid using “less” by typing

> more off

This can go into your startup script “octaverc”, which resides in a Windows installation in the folder

C:\Program Files\Octave\share\octave\site\m\startup

Installing Octave for Linux is easy if there is a RPM or similar package. However, I couldn’t find one for my Linux distrubtion, OpenSuSe (10.2).

I decided to roll up my sleeves and compile it. It took me a while to figure out how. Here’s what I ended up doing.

1) Install “compat-g77″ to allow Fortran 77 compilation. Download the rpm from this ftp site.

Then, as root, type “rpm -i rpmfilename”

This is the only step that was not obvious, at least not to me. I tried various other ways to install a Fortran compiler, namely “gfortran”, “g95″, and also rebuilding “gcc” from scratch. None of that is necessary. Just install the RPM, and have Fortran 77 on your system.

2) Install Gnuplot 4.2 from

http://www.gnuplot.info/download.html

The procedure for this installation, and for the others described below, is the same in each case.

• tar zxvf tarfilename
• cd directorynameofpackage
• configure
• make
• make check
• make install # as root

3) Install “FFTW”, “GMP”, “MPFR”, and “GLPK”. Look each of these up on Google and download the tarballs. Install these as in step 2.

4) Download the Octave source tarball, and install as in step 2. This will take about a hour, depending on your machine.

Octave is a public domain version of MATLAB. It is probably 70% compatible with MATLAB, which is close enough for educational purposes. Version 3.0 is a major update to Octave, with many new graphics and image processing capabilities.

The windows installer is simple and effective. I selected the JHandles options, meaning graphics is handled by the JAVA and Open-GL based “JHandles” package.

Here is the sombero plot that serves as Octave’s logo. To generate it, use

>sombrero(32),title(’sombrero(32)’);

The screenshot above (click to enlarge) shows how close to Infinity you can get in different programs.  In Matlab, the largest floating point number appears to be

1.79769313486231580e308

In Google Calculator, the largest number I can get displayed is

e^(709)

On my new HP 35S, the largest number displayed is 1.0 e 500. The closest I can come to this in a computation is ask for 253! 1.9329 x (using RPN), which returns 9.9999e499.

Windows Calculator, that surprisingly powerful tool, can handle much larger numbers. I don’t know what the limit is. The largest number I can generate without much effort is

7.6150458132720644624236074295991e+1042494

which is an exponent of more than one million. I’m sure there is room to go higher. If anyone knows what the maximum representable number is in Windows Calculator, please let me know. If you have another tool that allows you to get closer to Infinity, I’d love to find out.

I recently bought a HP 35S calculator. Though I like the calculator, I found something surprising: when asked for sin(pi), it does not report zero, but rather -2.067e-13. Similarly, for cos(pi/2), it reports -5.1034e-12. However, some results are still exact: cos(pi) returns -1.0 and sin(pi/2) = 1.0. Finite numerical precision is no doubt behind these results, but not all calculators use the same algorithm. For example, my Casio FX115, which cost less than half as much as the HP, gets a perfect score on the four values: sin(pi) = 0, cos(pi/2) = 1.0, sin(pi/2) = 1.0 and cos(pi/2) = 0. Furthermore, the Calculator built into Windows is also four for four.

Now consider MATLAB on a standard PC: despite being a high-end numerical analysis program, it reports sin(pi)=1.22e-16 and cos(pi/2)=6.12e-17; though it reports cos(pi) = -1.0 and sin(pi/2)=1.0. In the product documentation , under “doc sin”, it explains

“The expression sin(pi) is not exactly zero, but rather value the size of the floating-point accuracy eps, because ‘pi’ is only a floating-point approximation to the exact value of pi.

Under “algorithm”, MATLAB’s documentation says:

“sin uses FDLIBM, which was developed at SunSoft,a Sun Microsystems, Inc. business, by Kwok C. Ng, and others. For informationabout FDLIBM, see http://www.netlib.org.”

I also tried Google calculator. It too gets a perfect score

sin(pi)=0

cos(pi/2)=0

cos(pi) = - 1

sin(pi/2)=1