Threaded Mode | Linear Mode

robve · (This post was last modified: 03-01-2021 06:15 PM by robve.)

Challenge one

Note: code indent replaced by U+00A0 (NO-BREAK SPACE) by hand and may not be copyable/usable.

HP PRIME
--------
EXPORT MC()
BEGIN
  LOCAL SUM, N, TRIALS;
  RANDSEED(1);
  N := 0;
  FOR TRIALS FROM 1 TO 100000 DO
    SUM := 0;
    WHILE SUM <= 1 DO
      SUM := SUM + RANDOM();
      N := N + 1;
    END;
  END;
  MSGBOX(N/100000);
END;

BASIC
-----

randomize 1
n=0
for t=1 to 100000
sum=0
while sum<=1
  sum=sum+rnd(0)
  n=n+1
wend
next t
print n/100000

Python
------

import random

def MC():
  random.seed(1)
  n = 0
  for t in range(100000):
    sum = 0
    while sum <= 1:
      sum = sum + random.random()
      n += 1
  print(n/100000)

if __name__ == "__main__":
  MC()

a. What do you think is, in the limit, the average count of generated random numbers for their sum to exceed 1 ? Can you recognize what the exact value would be?

2.71959 when trials approaches infinity the exact value is e

b. What would the average count be for the sum to exceed 2 ? To exceed 2.021 ? To exceed Pi?

2 -> 4.67827
2.021 -> 4.71806
3 -> 6.66808
pi -> 6.95027
4 -> 8.66601
5 -> 10.66641
5+1/6 -> 10.99947
10 -> 20.65914
15 -> 30.66700
20 -> 40.66927

c. What do you think is the asymptotic expression for the average count needed to exceed large values of the sum?

I have an idea what kind of inequality in probability theory applies here, which uses EXP() but cannot recall from the top of my head without looking it up. The value appears to be converging to a value between 2k+1/6 and 2(k-1)+e.

d. Can you explain the constant component of said expression?

Working on this...

Challenge two

HP PRIME
--------

EXPORT WS2021()
BEGIN
  LOCAL I,J,F,P,SUM,PROD;
  L0 := [2,3,5,7,11,13,17,19];
  FOR P FROM 23 TO 99999 STEP 2 DO
    F := 0;
    FOR I FROM 1 TO SIZE(L0) DO
      IF P MOD L0[ I] = 0 THEN
        F := 1;
        BREAK;
      END;
    END;
    IF F = 0 THEN
      L0 := append(L0, P);
    END;
  END;
  SUM := 0;
  PROD := 1;
  FOR I FROM 1 TO SIZE(L0) DO
    PROD := PROD*L0[ I]/(2021+L0[ I]);
    SUM := SUM + PROD;
  END;
  MSGBOX(SUM);
END;

This program displays the value 9.90099737136E-4 for primes up to 99999 (and also for smaller bounds than 99999, such as 97 which is weird).

EDIT
Oops, double check: I read the sum wrong, for the term's numerator has k-1 while the denominator has k. Here is the simple correction. This gives 4.94804552201E-4 and does not change for primes >11 (this sum should converge quickly):

EXPORT WS2021()
BEGIN
  LOCAL I,J,F,P,SUM,PROD;
  L0 := [2,3];
  FOR P FROM 5 TO 17 STEP 2 DO
    F := 0;
    FOR I FROM 1 TO SIZE(L0) DO
      IF P MOD L0[ I] = 0 THEN
        F := 1;
        BREAK;
      END;
    END;
    IF F = 0 THEN
      L0 := append(L0, P);
    END;
  END;
  SUM := 0;
  PROD := 1;
  FOR I FROM 1 TO SIZE(L0) DO
    SUM := SUM + PROD/(2021+L0[ I]);
    PROD := PROD*L0[ I]/(2021+L0[ I]);
  END;
  MSGBOX(SUM);
END;

Interesting to observe is that the reciprocal of the sum 1 / 4.94804552201E-4 = 2021.

/EDIT

Challenge three

Evaluated on the HP PRIME displays 0.309016994375 (updated to correct a typo in the input).

I don't have the HP-71B (love it), but could use my other BASIC pocket computers since I wrote my own Romberg integration program and Gamma Lanczos approximation in BASIC some time ago. But alas, non-HP machines are not allowed.

EDIT: anyway, let me add this SHARP PC-1350 BASIC program (Valentin, your favorite sharp calc) that produces the answer 0.3090169944 which is perfect to the full 10 digits (updated to correct a typo):

' ROMBERG QUADRATURE WITH HIGH PRECISION - Numerical recipes qromb + polint p.134
' see also https://en.wikipedia.org/wiki/Romberg's_...ementation

' Functions to integrate are defined with label "F1", "F2", etc.

' VARIABLES
' A,B range
' F$,F function label (or number) to integrate
' Y result
' E relative error: integral = Y with precision E (attempts E = 1E-10)
' H step size
' N max number of Romberg steps (=10)
' I iteration step
' U current row
' O previous row
' J,S,X scratch
' A(27..26+2*N) scratch auto-array

100 INPUT "a=";A
110 INPUT "b=";B
120 INPUT "F";F
130 E=1E-9,N=10,F$="F"+STR$ F,X=A: GOSUB F$: S=Y,X=B: GOSUB F$
140 H=B-A,O=27,U=O+N,A(O)=H*(S+Y)/2,I=1
150 H=H/2,S=0
160 FOR J=2 TO 2^I STEP 2: X=A+J*H-H: GOSUB F$: S=S+Y: NEXT J
170 A(U)=H*S+A(O)/2,S=4
180 FOR J=1 TO I: A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1),S=4*S: NEXT J
190 IF I>1 IF ABS(A(U+I)-A(O+I-1))<E*A(O+I-1) LET Y=A(U+I): PRINT Y: END
200 S=O,O=U,U=S,I=I+1: IF I<N GOTO 150
210 Y=A(O+N-1),S=U+N-2,E=ABS((Y-A(S))/A(S)): PRINT Y,E: END

300 "F1" V=X,X=1.618033989,X=LN(X*X-V): GOSUB "GAMMA": W=Y
310 X=LN V: GOSUB "GAMMA": Z=Y,X=1.618033989,X=LN(X*X-V): GOSUB "GAMMA"
320 Y=W/(Z+Y): RETURN

' gamma(X) -> Y using Lanczos approximation
400 "GAMMA" IF X<=0 LET Y=9E99: RETURN
410 Y=1+76.18009173/(X+1)-86.50532033/(X+2)+24.01409824/(X+3)
420 Y=Y-1.231739572/(X+4)+1.208650974E-3/(X+5)-5.395239385E-6/(X+6)
430 Y=EXP(LN(Y*2.506628275/X)+(X+.5)*LN(X+5.5)-X-5.5): RETURN

/EDIT

Challenge four

Plot with HP PRIME shows two (cross) eyes. Very funny!

I used the advanced graphing app to enter the X-Y equation. The "eyes" radii are close to the root of 17/6 (updated) and centered at (±4/3,0).

The two "pupils" centered at (±1,0) take some time to converge. They may never computationally converge exactly to a single point.

Now on to the last two challenges. But some information is missing about the PP definition: do you mean that changing any digit to any other digit ALWAYS produces a composite number? So for example 17 is not a PP because 7->9 gives 19?

EDIT

Challenge five

PP #1: 294001
PP #2: 505447
PP #3: 584141
PP #4: 604171
PP #5: 971767

I wrote a Sieve of Eratosthenes program with an addition to filter perfect primes.

The big problem is that HPPL does not permit lists lengths exceeding 20K. So unfortunately I had to rewrite the program in C and run it. I wish HPPL has bitvectors or efficient sets!

EXPORT PP()
BEGIN
  LOCAL I,J,K,N,M;
  LOCAL F,D,P,Q,R;
  N := 9999;
  // INIT LIST OF 0 (COMPOSITE) OR 1 (PRIME)
  L0 := MAKELIST(odd(I),I,1,N,1)
  // SIEVE PRIMES
  I := 3;
  WHILE 1 DO
    WHILE I < N AND L0[ I] = 0 DO
      I := I+1;
    END;
    IF I = N THEN
      BREAK;
    END;
    FOR J FROM 2*I TO N STEP I DO
      L0[J] := 0;
    END;
    I := I+2;
  END;
  // FILTER PERFECT PRIMES
  P := 3;
  WHILE 1 DO
    WHILE P < N AND L0[P] = 0 DO
      P := P+1;
    END;
    IF P = N THEN
      BREAK;
    END;
    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 L0[Q+K*I] THEN
          F := 0;
        END;
      END;
    END;
    IF F THEN
      PRINT("PERFECT "+P);
    END;
    P := P+2;
  END;
END;

I spent a lot of time typing this program in on my HP PRIME (and the other ones). This challenge is meant to enter programs on the HP machine!

Note that the program checks for PP by changing digits. This includes changing the leading digit to any digit 0 to 9, including 0. So for example, 199 -> 099 is not PP. It seems odd to me to change the leading digit to 0 to check for PP.

The second HPPL program is less efficient, but permits a lower bound B and upper bound E for the search:

EXPORT WPP()
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;

The first method (sieving) is more efficient. The asymptotic running time is linear in the max perfect prime p we're hunting for, but requires memory of size p. The asymptotic running time of the second method is roughly square in p but requires constant memory. It takes a few minutes to find the first five perfect primes.

My Cheat program gives an answer in a fraction of a second:

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

int main()
{
  long i, j, k, n = 1000000;
  // sieve (0 = composite, 1 = prime)
  char *s = (char*)malloc(n);

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

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

  // sieve for perfect primes
  for (i = 3; i < n; i += 2)
  {
    while (i < n && s[i] == 0)
      ++i;
    if (i < n)
    {
      long m = floor(log10(i));
      long p = 1; // p = 10^j
      char h = 1; // h = 1 if i is perfect
      // for each digit in prime i, from least to most significant
      for (j = 0; j <= m; ++j, p *= 10)
      {
        // q = prime i with jth digit zero
        long q = i%p + i/p/10*p*10;
        // d = jth digit of prime i
        long d = i/p%10;
        // twiddle jth digit and check for prime
        for (k = 0; k <= 9; ++k)
          if (k != d && s[q+k*p])
            h = 0;
      }
      if (h)
        printf("perfect %ld\n", i);
    }
  }
  free(s);
}

b. 500004469 (first PP > 500000000)

c. 777781429 (first PP > 777777777)

d. 666999929 (second PP > 666666666)

For small primes with k digits such as k=1, k=2, …, k=5 digits there are no perfect primes in base 10, because primes with k digits are "too close" (as in a Hamming distance kind of way). The prime number theorem tells us that primes are distributed roughly as N/log(N) so that a randomly picked integer not greater than N has a probability of 1/log(N) to be prime. A perfect prime of k digits can be perturbed by its definition to generate 9k distinct composite integers. Roughly, the chance that an integer of k digits is prime is 1/log(10^k)=0.43/k. The chance that the 9k integers are all composite is (1-0.4343/k)^(9k), assuming k is sufficiently large. It turns out that the chance approaches 2% for large k and is half that for small k (though the log constant is somewhat arbitrary). More importantly, there are also far more integers to pick as potential perfect primes for large k. Based on this, it seems reasonable to see perfect primes for large k and there are infinitely many of them.

/EDIT