Sums of all Squares

[CRGreathouse]

I'm looking at your posts on the other forum. You write:

#include <math.h>
#include <stdio.h>
#include <limits.h>
#pragma precision=log10l(ULLONG_MAX)/2

typedef enum { FALSE=0, TRUE=1 } BOOL;

BOOL is_prime( int p ){
if( p == 2 ) return TRUE;
else if( p <= 1 || p % 2 == 0 ) return FALSE;
else {
BOOL prime = TRUE;
const int to = sqrt(p);
int i;
for(i = 3; i <= to; i+=2) 
if (!(prime = p % i))break;
return prime;
}
}

BOOL is_mersenne_prime( int p ){
if( p == 2 ) return TRUE;
else {
const long long unsigned m_p = ( 1LLU << p ) - 1;
long long unsigned s = 4;
int i;
for (i = 3; i <= p; i++){
s = (s * s - 2) % m_p;
}
return s == 0;
}
}

int main(int argc, char **argv){

const int upb = log2l(ULLONG_MAX)/2; 
int p;
printf(" Mersenne primes:\n");
for( p = 2; p <= upb; p += 1 ){
if( is_prime(p) && is_mersenne_prime(p) ){
printf (" M%u",p);
}
}
printf("\n");
}

This code looks fine for p up to 31, but it can't work any higher. It uses multiplication mod 2^p - 1, so the numbers used need to be able to fit (2^p - 2)^2 at the very least. For that section they use 64-bit unsigned integers, which are barely enough for 31. For the next Mersenne prime 2^61 - 1, you need at least 122 bits (practically, 128); for the following, at least 178 bits (practically, 192).

Also, I would be wary of the call to sqrt; it may not round properly for sufficiently large numbers.

Just for kicks, here's how I'd code the basic trial division prime tester:

int isprime(int p ) {
  if(!(p&1))
    return p == 2;
  if(p%3 == 0)
    return p == 3;
  if(p < 25)
    return p > 1;
  int i, to = (int)sqrt(p + 0.5);
  for(i = 5; i <= to; i += 6) {
    if (p%i == 0)
      return 0;
    if (p%(i+2) == 0)
      return 0;
  }
  return 1;
}

[davar55]

Another possibly interesting variation:

2^2 + 3^3 + 5^5 + ... + p^p = 10^mK

What is the smallest prime p such that
the sum of squares of all primes up to p
is a multiple of 10 (or 100 or 1000 and so on).

(I think somewhere in these series we'll get
numeric pointers to empirical evidence
connecting some number theoretical conjectures.)

[CRGreathouse]

Originally Posted by davar55

Another possibly interesting variation:

2^2 + 3^3 + 5^5 + ... + p^p = 10^mK

What is the smallest prime p such that
the sum of squares of all primes up to p
is a multiple of 10 (or 100 or 1000 and so on).

11, 751, 1129, then something big.

Originally Posted by davar55

(I think somewhere in these series we'll get
numeric pointers to empirical evidence
connecting some number theoretical conjectures.)

Why?

[davar55]

Something big? Now wouldn't that be cool to see ...

[CRGreathouse]

Originally Posted by davar55

Something big? Now wouldn't that be cool to see ...

Nah, it turns out it's only 361,649. My code had been storing the full numbers (which are large) and was running slowly as a result.

[davar55]

So far, this (sum p^p) sequence is:

11
751
1129
361649

Couldn't have done that without a computer.

[NBtarheel_33]

Originally Posted by davar55

So far, this (sum p^p) sequence is:

11
751
1129
361649

Couldn't have done that without a computer.

2^2 = 4
2^2 + 3^3 = 4 + 27 = 31
2^2 + 3^3 + 5^5 = 4 + 27 + 3125 = 3156
2^2 + 3^3 + 5^5 + 7^7 = 826699

The next term brings the sum to 285 billion and change (i.e. 285e9), and the term after that brings the sum to over 303 trillion (i.e. 303e12).

Looking at the sum mod 10, we see that 2^2 = 4, 3^3 = 7, 5^5 = 5, 7^7 = 3, 11^11 = 1. This adds up to 0 (mod 10). Hence,

2^2 + 3^3 + 5^5 + 7^7 + 11^11 = 285,312,497,210 would be the first time the sum is divisible by 10.

Likely, the way to attack the problem for mod 10^n, would be to ask questions about the residues of p^p (mod 10^n) in general.

[davar55]

I suppose my batteries would have been sufficient to get
that first term, but for the second I would have needed my
home-built multi-precision (computer-resident) calculator,
written in C but currently in mothballs.

Does the sum for 361649 really end in four zeros?
Maybe those sums should be in the same table.

This could get interesting.

[CRGreathouse]

Originally Posted by davar55

Does the sum for 361649 really end in four zeros?

Five, actually -- it will be fourth and fifth on your list:

11, 751, 1129, 361649, 361649, 12462809, 12462809, 1273183931, 1273183931.

Would someone check me on these? The repeating entries are freaking me out.

[davar55]

Well, so far the sequence of sum(p^p)mod(10^n)
(I think it needs a better name) is:

11
751
1129
361649
361649
12462809
12462809
1273183931
1273183931

Any bets that the tenth and eleventh terms are equal,
and whether the twelfth breaks THAT pattern?

[CRGreathouse]

Nothing to 10^10.