[VA] Short & Sweet Math Challenge #25 "San Valentin's Special: Weird Math"

[robve]

Some additional comments on concoction five.

First, the HP PRIME HPPL program with the Sieve of Eratostenes is the smart way to approach this problem efficiently. Unfortunately, sieving for primes above 20,000 or so becomes tricky, due to HPPL list length restrictions to 20,000. One could populate the list with 64 bit integers and use bit banging to increase the prime sieving space 64 fold, which is just enough to find the first 5 perfect primes.

Without sieving is straightforward to implement with a few loops, but slow and takes a few minutes rather than seconds (e.g. compared to C). The program hunts for perfect primes within designated bounds B and E:


EXPORT SLOWPP()
BEGIN
  LOCAL B,D,E,F,I;
  LOCAL J,K,P,Q;
  // BEGIN SEARCH AT B
  B := 11;
  // END SEARCH AT E
  E := 9999999;
  FOR P FROM B TO E STEP 2 DO
    IF isprime(P) THEN
      M := FLOOR(LOG(P));
      F := 1;
      FOR J FROM 0 TO M DO
        I := 10^J;
        // Q = P WITH JTH DIGIT SET TO 0
        Q := FLOOR(P/I/10)*I*10+ROUND(FP(P/I)*I);
        // D = JTH DIGIT OF P (0 to 9)
        D := FLOOR(P/I) MOD 10;
        FOR K FROM 0 TO 9 DO
          IF K <> D AND isprime(Q+K*I) THEN
            F := 0;
            BREAK;
          END;
        END;
        IF F=0 THEN
          BREAK;
        END;
      END;
      IF F THEN
        PRINT("PERFECT "+P);
      END;
    END;
  END;
END;

Hence, hereby my alternative submission in HPPL.

Second, I stated that theoretically more perfect primes should exist for bases smaller than 10 and also gave a list of the first base-perfect primes for base 2 to 20. Running additional experiments corroborates this claim empirically verified for primes up to 1,000,000,000 (with a small modification made to the C program):


base | PP count |
2    | 7179981  |
3    | 10070244 |
4    | 1521240  |
5    | 7627392  |
6    | 12456    |
7    | 3762873  |
8    | 98926    |
9    | 400411   |
10   | 3167     |
11   | 1018816  |
12   | 37       |
13   | 558553   |
14   | 403      |
15   | 2821     |
16   | 743      |
17   | 153894   |
18   | 0        |


An interesting pattern emerges, which shows something else, besides the frequency of base-perfect primes is higher for smaller bases as expected. Note that prime base perfect primes are more prevalent.

EDIT
Try it yourself: my C program to find PPs for any base for the given max:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char **argv)
{
  if (argc > 2)
  {
    long i, j, k;
    long base = atoi(argv[1]); // perfect prime base
    double lnbase = log(base);
    long max = atoi(argv[2]);  // primes up to max (exclusive)
    long maxd = 0;             // primes of maxd digits
    long num = 0;              // number of perfect primes found
    long len = 0;              // length of the array
    char *sieve;               // the sieve array (0 = composite, 1 = prime)

    // make max odd, e.g. 1000000000 -> 99999999 has 9 digits to perturb, not 10
    max -= !(max & 1L);

    // len = base^maxd for maxd>0 such that len >= max
    do
      len = pow(base, ++maxd);
    while (len < max);

    sieve = (char*)malloc(len);

    // init sieve, we should keep only odd values
    for (i = 0; i < len; ++i)
      sieve[i] = (i & 1L);

    // sieve for primes
    for (i = 3; i < len; i += 2)
    {
      while (i < len && sieve[i] == 0)
        ++i;
      for (j = 2*i; j < len; j += i)
        sieve[j] = 0;
    }

    // sieve for perfect primes
    for (i = 3; i < len; i += 2)
    {
      while (i < len && sieve[i] == 0)
        ++i;
      if (i > max)
        break;
      if (i < len)
      {
        long m = floor(log(i)/lnbase); // m = number of digits of i
        long p = 1;                    // p = base^j
        char h = 1;                    // h = 1 if i is perfect
        // for each jth digit in prime i, from least to most significant
        for (j = 0; j <= m; ++j, p *= base)
        {
          // q = prime i with jth digit zero
          long q = i%p + i/p/base*p*base;
          // d = jth digit of prime i
          long d = i/p%base;
          // twiddle jth digit and check for prime
          for (k = 0; k < base; ++k)
            if (k != d && sieve[q+k*p])
              h = 0;
        }
        if (h)
          ++num;
      }
    }

    printf("num=%ld\n", num);

    free(sieve);
  }
  else
  {
    printf("Usage: pps <BASE> <MAX>\n");
  }
}

/EDIT

- Rob