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

[Valentin Albillo]

 
Hi, all !   Happy San Valentín 2021 to all of you !  
 
Almost a year has elapsed since my previous S&SMC, exactly 6 years 6 since I registered in this new-style forum, and today it's San Valentin's Day so here you are, a brand-new Short & Sweet Math Challenge #25  "San Valentin's Special: Weird Math", once again intended to give you a chance to dust off your chosen HP calculator and actually use it, while also flexing your HP-programming muscle *and* your HP-assisted sleuthing abilities (*NOT* your Google Search proficiency).

Try all 6 Concoctions 6 below and see what you can do about them !

Rules:

• Any HP calc of your choice may be used (physical or virtual) but I suggest a Minimum Recommended Model (MRM) for each Concoction, which is the simplest model I deem capable of solving it more or less comfortably. Lesser models might, but not comfortably.
   
• Using anything other than a physical or virtual HP calculator is strictly disallowed, thus no Wolfram Alpha/Mathematica/Matlab/Pari, VBA, Excel. Pascal, C/C#/C++, Java, Python, Lua, etc., please go elsewhere for that. The idea here is that you actually use your HP calc, not your Lua on a laptop, so write your code in a language originally supported on some HP calc (RPN, RPL, 71B BASIC, Hppl, etc).
   
• Googling for the solutions is über-lame and by doing so you're admitting to yourself than you just can't cope and have to cheat to try and save face  .

Concoction the First: Weird limit
[MRM: HP-11C and up]

Many HP calculators have a function that returns a (pseudo-)random number uniformly distributed between 0 (inclusive) and 1 (exclusive), such as RND in the HP-11C, HP-15C, HP-71B and many other models, and RAN in the HP-42S.

Now get your calc and conduct a test where you generate a couple of random numbers and add them up, then if the sum is less than 1 you continue generating and adding up more random numbers one at a time until the sum eventually exceeds 1, while keeping count of just how many random numbers did you generate.

For instance, suppose your first test generates 0.87 and 0.65 and their sum is 1.52, which is greater than 1 already so the test is over and 2 random numbers were generated in all. Now you conduct another test and the sequence of random numbers is 0.21, 0.07, 0.16, 0.35, 0.19 and finally 0.58, which makes the sum (1.56) exceed 1, so the test is over and we had to generate 6 random numbers in all.


The Challenge:

Write a program that simulates the process for a given number of tests and outputs the average count of random numbers generated per test, and then (the sleuthing part) use the program to help you answer these questions:

• 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 ?
• b.  What would the average count be for the sum to exceed 2 ? To exceed 2.021 ? To exceed Pi ?

Surprising, isn't it ? But there's more surprises: now get a suitably fast HP model (physical or virtual) and conduct a sizable number of tests (say ~100,000) to find the average count for sums exceeding 3, 4, 5 and 5+1/6 (no kidding, try it) on the one hand, and for sums exceeding 10, 15 and 20 on the other, and analyze the results you get, in order to answer these additional questions:

• c.  What do you think is the asymptotic expression for the average count needed to exceed large values of the sum ?
• d.  Can you explain the constant component of said expression ?

Use a seed of 1 for the random number generator at the very start of your program (for instance, RANDOMIZE 1 on the HP-71B and 1, SEED on the HP-42S) and give your answers accurate to at least 2-3 digits. Please do not use/post any theoretical formulas to get the results for now, do it empirically by just generating and using actual random numbers.

I'll give my original solutions for both the HP-42S and the HP-71B, as well as my comments on the results.


Concoction the Second: Weird Sum
[MRM: HP-11C and up]

Consider the following infinite Albillo sum:



where, other than 2021, the coefficients are the prime numbers p_k in order: p₁ = 2, p₂ = 3, p₃ = 5, p₄ = 7, ...

Note: Observe that for k=1 in the sum above, the product in the numerator is the empty product, thus equal to 1 by definition.


The Challenge:

Write a program to compute and output the sum as accurately as possible, and then (the sleuthing part) use your HP calc (perhaps conduct some experiments) to try and attempt to answer this question: What's so weird about this sum ? (BTW, forget about Googling for it because I concocted it myself and it's nowhere else to be found AFAIK.)

I'll give my original solution for the HP-71B, as well as my comments on the result.


Concoction the Third: Weird Integral
[MRM: HP-15C and up]

Consider the following definite Albillo integral:



where Γ is the Gamma function, ln is the natural logarithm (i.e., base e) and φ is the Golden Ratio = (1+ √5)/2.


The Challenge:

Use your HP calc to compute (either manually or writing a program to do it) and output the value of the definite integral as accurately as possible, and then (the sleuthing part) use your HP calc (perhaps conduct some experiments) to try and attempt to answer this nagging question: What's so weird about this integral ?. (Again, forget about Googling for it because I concocted it myself and it's nowhere else to be found either.)

I'll give my original solution for the HP-71B, as well as my comments on the result.


Concoction the Fourth: Weird Graph
[MRM: HP-PRIME and other graphing models]

Consider the following polynomial in two real variables x, y:

      P(x, y) = 9 x⁸ + 9 y⁸ + 36 x² y⁶ + 54 x⁴ y⁴ + 36 x⁶ y² - 100 x⁶ - 4 y⁶ - 108 x² y⁴ - 204 x⁴ y² + 182 x⁴ - 10 y⁴ - 84 x² y² - 100 x² - 4 y² + 9


The Challenge:

Use your graphing calc (remember: no Wolfram Alpha, etc), either by writing and running some program or manually (but then succintly specify the operations performed) to accurately plot the resulting 2D graph for P(x, y) = 0, and somewhat describe what you see in the graph you get, giving also relevant parameters (say ranges for x and y, or maybe things like centers or focii or radii or asymptotes or intersections or zeros/poles, whatever.

If you think it might help, you may also attempt to factorize the polynomial, but in any case the main question is: What's so weird about this graph ?

I wont post code or manual operations as I don't own any graphing calculators but I'll give the resultant graphic you should get, as well as extensive commentary.


Concoction the Fifth: Weird Primes
[MRM: any fast physical or virtual HP calc]

In Milos Forman's 1984 "Amadeus" film (winner of 8 Academy Awards {aka Oscars}, including Best Picture) Salieri comments on the perfection of Mozart's music:

      "Displace one note and there would be diminishment, displace one phrase and the structure would fall."

Now let's bring that observation to the realm of prime numbers and consider a prime number so 'Perfectly Prime' (a PP for short, pronounced "Pepe") that changing any single digit would diminish its primeness by turning it into a composite number. Note: We're talking about base-10 digits here.


The Challenge:

Write a program (the faster & shorter, the better) for your HP calc to compute: (a) the 5 smallest PP, (b) the first PP greater than 500 million, (c) the first PP greater than 777,777,777 and only for very fast programs/devices, the second PP greater than 666,666,666.

I'll give my original solution for the HP-71B, as well as my comments on the results.


Concoction the Sixth: Weird Year
[MRM: HP-11C and up]

Note: All that follows is stated utterly tongue-in-cheek, not to be taken seriously at all. No offence whatsoever is meant to anyone who's been negatively affected during 2020.

Unless you've been hiding under a rock last year, surely you're sorely aware that 2020 was a catastrophic year and many of you might wonder why did it came that way. I know I did, and being fully convinced that this Universe of ours is a mathematical entity subject to mathematical rules, I have been analyzing the matter exhaustively using my trusty HP calculators and have finally succeeded in unraveling the mystery !! At long last, I now know the reason why 2020 was a catastrophic year and of course the reason is of a mathematical nature, as expected.

To wit, the reason is that 2020 shares a very striking numeric property with many other catastrophic years such as, e.g.: the year 662 (the Damghan earthquake killed 40,000 people), the year 458 (the Antioch earthquake killed 80,000), the year 1348 (the Black Death plague, which killed up to 200 million, was at its apogee), the year 1556 (the Shaanxi earthquake killed 830,000) or the year 1849 (the Great Irish Famine killed ~1,500,000), to name but a few.

That can't be a mere coincidence !  Moreover and just in case this wasn't evidence enough, the number 666 (the infamous Number of the Beast of apocalyptic fame) also shares that very property as well.


The Challenge:

Use your trusty HP calculator to assist you in your sleuthing to try and discover what simple but highly remarkable (striking, in fact) numeric property all the aforementioned numbers have in common, and then write a program to find out and output a listing of all years between AD 4 and AD 5000 (both included) which have this property (hint: less than 100). Of course the listing should include all mentioned past years as well as future years thus predicted to be potentially catastrophic, up to that limit.

The questions are: (a) How many years will be listed in the output ? (b) What will be the next predicted potentially catastrophic year after 2020 ?, and (c) Should we be concerned ?

As an additional hint to help finding the remarkable shared property, remember Occam's Razor: the property can be unambiguously stated by saying that the year's number is a "(five words)".

I'll give my solution, two short programs (6 and 7 lines resp.) for the HP-71B which produce the listing and also accept a given year in range and demonstrate^* whether it has the required numeric property (thus, if it indeed was/might be catastrophic) or not.

• * E.g.: For property "The year's number is a factorial" and year 720 you would output "720 = 1x2x3x4x5x6" demonstrating the property, while for year 721 you would ouput "721 = not a factorial".

And last but certainly not least, a couple' important caveats:

• Please do NOT include CODE panels in your replies to this thread, as it makes it difficult for me to generate the online PDF document which will include the whole thread. I expect you'll kindly comply with this requirement but otherwise you'll risk your carefully crafted code appearing truncated or not at all in the final PDF and thus being irretrievably lost from the online document and making your posting it moot. With no CODE panels, to prevent solutions posted before yours spoiling your fun, I advise to simply avoid reading any of them before posting your own. Thank you.
   
• Designing, testing and formatting these Challenges and their solutions takes and awful lot of time and effort. Hence, if you do enjoy them and would like to see more posted in the future, consider participating or commenting on them so that I get feedback of your appreciation.

I'll post my original solutions in a week or so, for you to have enough time (and no excuses) to try them all. That is, if you've got what it takes ...

That's all. Enjoy ! ... and that's an order !

V. 

[Vincent Weber]

Hi Valentin,

Is the Python syntax available on the HP Prime CAS allowed ?

I kind of fell in love with Python, now that it is the mandatory language on educational calculators in France (starting in 2017 with the awesome Numworks, then Casio, then TI, with HP supporting its syntax along with HPPL on the Prime). I am teaching my kids programming with both Scratch and Python, and I have to admit that I like it a lot - Scratch gives the youngest good ideas of structural programming with friendly drag and drop graphical blocks, and Python is quite easy to catch up for older kids.

I find Python more elegant and expressive than any Basic, even the massively enhanced 71B one (but granted, the 71B can interact seamlessly with advanced mathematical features such as matrices. To compete that, Python calculators should provide Numpy, which is not the case yet).

On the Numworks (did I say how much I love this tiny, stylish, fast, well designed machine with awesome software ? Valentin you should love it, as you love beautiful objects ) the Python menus allow you to generate skeletons of code with utmost ease.

Anyway... I already got the (surprising indeed !) answer for #1. If you ban Python altogether, I shall try with RPN on Free42

Cheers,

Vincent

[J-F Garnier]

I don't know the reason but I'm particularly attracted by the integral problems (like this previous challenge).
So I tested the "weird" integral:

(02-14-2021 08:58 PM)Valentin Albillo Wrote:  

where Γ is the Gamma function, ln is the natural logarithm (i.e., base e) and φ is the Golden Ratio = (1+ √5)/2.

The expression to integrate looks complicate, although without any trap, it is not defined at the integral boundaries but that's not a problem with the Romberg algorithm.
I took my old trusted HP-32S and entered the program:

LBL I
RCL X
LN
1
-
x!
RCL F
x²
RCL X
-
LN
1
-
x!
+
LASTx
x<>y
/
RTN

Then:
FIX 11
put the golden ratio in F: 5 SQRT 1 + 2 / STO F
FN= I
1
RCL F
∫FN dX
and quickly got the answer up to at least 11 places.

So what's special? It was fast and easy, without any problem.
Oh wait, it was too easy, too fast. That's weird.

[...Investigating a bit...]

Now that I think I found what is "weird", I can even do the calculation by hand, and get the symbolic result without having to identify the numeric result:
[... hidden for now ...]

J-F

[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

[J-F Garnier]

(02-15-2021 08:51 PM)robve Wrote:  Challenge three
Evaluated on the HP PRIME displays 5.61670944148E-2.

Can you double check?
I got .30901699438 on my 32S.

And what do you think is weird?

J-F

[edited: ok I see your mistake, you missed the LN in the numerator]

[ramon_ea1gth]

I tried the weird limit in my physical HP 50g in User-RPL, reaching 100,000 iterations. When I saw the outcome going to e, some kind of law of large numbers came to my mind.

Input stack:
2: number of tests, CNT
1: limit to exceed, XCD (e.g., 1, pi)

My code:
« 0 → CNT XCD TOT
« 1 RDZ
1 CNT FOR I
RAND 1 → TMP NR
« DO RAND 'TMP' STO+ 1 'NR' STO+
UNTIL TMP XCD >
END
NR 'TOT' STO+
»
NEXT
TOT CNT /
»
»

[Valentin Albillo]

.
Thanks for your interest and for participating, some quick answers to your questions:

Vincent Weber Wrote:Is the Python syntax available on the HP Prime CAS allowed ? Smile [...] If you ban Python altogether, I shall try with RPN on Free42 Smile

If it's available for the HP Prime then it's technically allowed but I'd rather have a program in Hppl for the Prime or in RPN/RPL/BASIC/FORTH/Assembler/etc. for any other HP calc. Anyway, your choice.

robve Wrote: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?

Yep. If changing any single digit to any other digit results in a prime number, the original one is not a PP.

Regards.
V.

[Vincent Weber]

Thanks Valentin.
Anyway it seems that similar solutions to mine have already been posted. This one is actually fairly easy, from a programming perspective, if not from a theorical perspective - but you didn't allow theory here

Cheers,

Vincent

[Werner]

Free42 code for challenges 1,2 and 3:

VA#1: average nr of random numbers in X runs to sum over Y

e.g. 2 ENTER 100000 XEQ"VA1" -> 4.67827

00 { 39-Byte Prgm }
01▸LBL "VA1"
02 LSTO "X"
03 LSTO "C"
04 SIGN
05 SEED
06 CLX
07 STO ST Z
08▸LBL 02
09 ISG ST Z
10 X<>Y
11 RAN
12 +
13 X≤Y?
14 GTO 02
15 CLX
16 DSE "C"
17 GTO 02
18 R^
19 RCL÷ "X"
20 END

results identical to robve's post. I would say that for 'larger' n, the limit is 2(n+1/3)..

00 { 73-Byte Prgm }
01▸LBL "VA2"
02 1
03 LSTO "P" @ term
04 CLX
05 LSTO "S" @ sum
06 2
07 XEQ 14
08 3
09▸LBL 03
10 XEQ "P?"
11 X=0?
12 GTO 00
13 R↓
14 XEQ 14
15 X=Y?
16 RTN
17 R↓
18▸LBL 00
19 R↓
20 2
21 +
22 GTO 03
23▸LBL 14
24 RCL "P"
25 2021
26 RCL ST Z
27 +
28 ÷
29 ×
30 STO "P"
31 X<> ST L
32 X<> "S"
33 STO+ "S"
34 RCL "S"
35 END

The P? primality testing program returns X and 0 when not prime, 1 when prime. You may use eg

00 { 114-Byte Prgm }
01>LBL "P?"
02 ENTER
03 STO ST Z
04 2
05 X=Y?
06 GTO 01
07 X>Y?
08 SIGN
09 MOD
10 X=0?
11 RTN
12 CLX
13 3
14 X=Y?
15 GTO 01
16 MOD
17 X=0?
18 RTN
19 CLX
20 5
21 X=Y?
22 GTO 01
23 MOD
24 X=0?
25 RTN
26 CLX
27 7
28>LBL 03
29 X^2
30 X>Y?
31 GTO 01
32 SQRT
33 MOD
34 X=0?
35 RTN
36 CLX
37 4
38 LASTX
39 +
40 MOD
41 X=0?
42 RTN
43 CLX
44 2
45 LASTX
46 +
47 MOD
48 X=0?
49 RTN
50 CLX
51 4
52 LASTX
53 +
54 MOD
55 X=0?
56 RTN
57 CLX
58 2
59 LASTX
60 +
61 MOD
62 X=0?
63 RTN
64 CLX
65 4
66 LASTX
67 +
68 MOD
69 X=0?
70 RTN
71 CLX
72 6
73 LASTX
74 +
75 MOD
76 X=0?
77 RTN
78 CLX
79 2
80 LASTX
81 +
82 MOD
83 X=0?
84 RTN
85 CLX
86 6
87 LASTX
88 +
89 GTO 03
90>LBL 01
91 SIGN
92 END

@robve: you made the same mistake I did at first.. read the formula carefully, the nominator goes to pk-1 and the denominator to pk.

The sum converges quickly to 4.948045522...E-04, which is exactly 1/2021. Why? I'm sure Valentin will enlighten us in a few days ;-)

VA#3: integral
Program to be used with INTEG

00 { 37-Byte Prgm }
01▸LBL "VA3"
02 MVAR "X"
03 1.25
04 SQRT
05 1.5
06 +
07 RCL- "X"
08 LN
09 GAMMA
10 RCL "X"
11 LN
12 GAMMA
13 RCL+ ST Y
14 ÷
15 END

Making use of the fact that φ^2 = φ+1
Same observations as J-F.. though the fact that it's easy to see what the value is, does not explain (to me) why it converges so quickly.

Cheers, Werner

[rprosperi]

(02-16-2021 09:56 AM)Werner Wrote:  The sum converges quickly to 4.948045522...E-04, which is exactly 1/2021. Why? I'm sure Valentin will enlighten us in a few days ;-)

The first S&S Math Challenge in 2021?

[Namir]

Regarding the weird sum that uses primes, I used Matlab to calculate each term and the updated sum. The sum seems to reach the asymptotic value of 9.90099737e-04.

The value of 2021, which I will call the Albillo Parameter helps the summation to reach the assumptotic value at the prime of 7. The smaller values of the Albillo Parameter will delay the summation from reach an assymptotic value (which depends on the Albillo Parameter). For example if the Albillo Parameter is, say, 101, the summation reaches the assymptotic value of 2.00059191e-02 at prime of 23.

My main focus is on the sumamtion in general and the role that the Albillo Parameter plays.

My own question is what is the relationship between the summation's assymptotic value and the value of the Albillo Parameter?

Cool summation!!!!

:-)

Namir

[robve]

For challenge five, which I have finished, it does not seem possible to implement efficient sets on the HP PRIME because lists are limited to only 20K elements. There is no bitvector or efficient set type that we can use.

Because of this limitation I rewrote the program in C (sorry for Cheating...). With this I found:

perfect 294001
perfect 505447
perfect 584141
perfect 604171
perfect 971767

Note that the program checks for PP by changing one digit of a prime at a time to check if that number is prime or not. 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.

I updated my post with EDITs to add my answer as well as make a correction to the weird sum (I misread k-1).

- Rob

[J-F Garnier]

(02-16-2021 09:56 AM)Werner Wrote:  VA#3: integral
[..]
Making use of the fact that φ^2 = φ+1
Same observations as J-F.. though the fact that it's easy to see what the value is, does not explain (to me) why it converges so quickly.

I thought the other way around: I understood that it's so fast because each sample set used by the algorithm gives the exact value of the integral, so the algorithm terminates after just the first 2 estimations (7 samples used).
Then I thought that even one sample is enough to give the exact value, using the midpoint that contributes to 1/2, hence the result.

And yes, φ^2 = φ+1 is the key :-)

J-F

[Werner]

The key to the value of the integral in #3 is to realize that the integral is symmetric in the integral boundaries. Basically it is

integral(a,b,f(a+b-x)/(f(x) + f(a+b-x)),dx)

if you substitute y = a+b-x, dy=-dx, you get

integral(b,a,-f(y)/((y) + f(a+b-y)),dy) = integral(a,b,f(y)/((y) + f(a+b-y)),dy)

Meaning the value of the integral is the sum of both divided by two, or (b-a)/2, whatever f(x) is. So here it is (φ-1)/2.

Cheers, Werner

[Gerson W. Barbosa]

This should provide all the answers to #6:

20

« { } 1 DUP2 ROT 5 ROLL
START NEXTPRIME + LASTARG NIP
NEXT DROP SQ DUP SIZE 2 SWAP
FOR i 1 OVER SIZE 1 + i -
FOR j DUP j DUP i + 1 - SUB ∑LIST ROT + SWAP
NEXT
NEXT + SORT
»

EVAL

« REVLIST 1 119
START TAIL
NEXT SORT
»

EVAL

Perhaps 19 or less would suffice as the argument for the first program, but I have not enough time to check this out right now.

Best regards,

Gerson.

[Nihotte(lma)]

(02-14-2021 08:58 PM)Valentin Albillo Wrote:   
Hi, all !   Happy San Valentín 2021 to all of you !  
 
...
That's all. Enjoy ! ... and that's an order !

V. 

UPDATED 2021.02.19 - 04.10PM to comply with the mention : Please do NOT include CODE panels in your replies to this thread

Hi,

Just to say, I've taken part of the challenge too !
However, I've progressed more slowly than you all...

For a. and b. points

I began with my HP35s :

** LBL V
INPUT S // 1 for stop
INPUT L // 100000 for loops
CLΣ
1
SEED
RCL L
STO I
009 CLx
STO J
STO C
012 1
STO + J
RANDOM
STO + C
RCL C
RCL S
x > y
GTO V012
0
RCL J
Σ+
DSE I
GTO V009
RCL L
RCL S
x-mean
RTN


XEQ V001
with S=1 and L=1E5 gave 2.71959 about 6 hours later

Then, I've brought out my 48G

<<
DEPTH DROPN
0 'SUM' STO
SEED RDZ
1 LOOPS
FOR i
0 DUP
DO
SWAP
1 +
SWAP
RAND +
DUP
TOP
UNTIL >= END
DROP 'SUM' STO+
NEXT
TOP LOOPS SUM OVER /
>>

with :
VAW is the program above
1 'TOP' STO // 1, 2 .. 20
1 'SEED' STO
100000 'LOOPS' STO
And the run of VAW gives
1
100000
2.71959
in the stack
and SUM is 271959

The result appears more rapidly than with the HP35s, of course
6x faster, perhaps

Unsurprisingly, my results match those of robve's post for the successive limit values ​​1, 2 .. 20 retained


------


The result for a limit of 1 seems to be close to e.
But, I've tested something on my HP10BII+
I decided to enter all couple of result by Σ+ in the calculator

C STAT
1 INPUT 2.71959 and Σ+
2 INPUT 4.67827 and Σ+
2.021 INPUT 4.71806 and Σ+
3 INPUT 6.66808 and Σ+
pi INPUT 6.95027 and Σ+
4 INPUT 8.66601 and Σ+
5 INPUT 10.66641 and Σ+
5+1/6 INPUT 10.99947 and Σ+
10 INPUT 20.65914 and Σ+
15 INPUT 30.66700 and Σ+
20 INPUT 40.66927 and Σ+
then,
[BLUE] REGR and [-] to select 0 - bESt Fit [INPUT]
and,
2 [ORANGE] ŷ,m displaying bESt Fit with the choice of 1 - LinEAr and the result of 4.6773856...
followed by [ORANGE] ^x,r resulting in 2 and [SWAP] giving 0,999999... as a correlation coefficient to describe the goodness of the fit.

Now here is where I am going with my initial idea :
1000 [ORANGE] ŷ,m 1999.68225...
100000 [ORANGE] ŷ,m 199900.966545...
200000 [ORANGE] ŷ,m 399801.25371...
1E9 [ORANGE] ŷ,m 1999002872.33...
9999999999 [ORANGE] ŷ,m 19990028715,2...
I don't know what to think but, it's giving a result near of 2*x decreased by something that is proportional to x (based on m, in fact) !!!

[Gerson W. Barbosa]

(02-18-2021 01:41 PM)Gerson W. Barbosa Wrote:  Perhaps 19 or less would suffice as the argument for the first program...

Indeed 19 is enough:

19

« { } 1 DUP2 ROT 5 ROLL
START NEXTPRIME + LASTARG NIP
NEXT DROP SQ DUP SIZE 2 SWAP
FOR i 1 OVER SIZE 1 + i -
FOR j DUP j DUP i + 1 - SUB ∑LIST ROT + SWAP
NEXT
NEXT + SORT
»

EVAL

« REVLIST 1 99
START TAIL
NEXT SORT
»

EVAL

→

{ 4 9 13 25 34 38 49 74 83 87 121 169 170 195 204 208 289 290 339 361 364 373 377 458 529 579 628 650 653 662 666 819 841 890 940 961 989 1014 1023 1027 1179 1348 1369 1370 1469 1518 1543 1552 1556 1681 1731 1802 1849 2020 2189 2209 2310 2330 2331 2359 2384 2393 2397 2692 2809 2981 3050 3150 3171 3271 3320 3345 3354 3358 3481 3530 3700 3721 4011 4058 4061 4350 4489 4519 4640 4689 4714 4723 4727 4852 4899 }



P. S.:

2020 is in the list because 2020 = 17² + 19² + 23² + 29². In fact every number in the list is either a squared prime or the sum of an ordered sequence thereof.

The number of consecutive prime squares starting with 4 (2²) required to safely producing a list of all years that share this property, up to a certain year, can be mathematically determined, but that’s way beyond my math skills. The following program uses a borrowed formula.

5000

« DUPDUP 4 * SWAP LN SQ / 3 XROOT 3 * CEIL { } 1 DUP2 ROT 5 ROLL
START NEXTPRIME + LASTARG NIP
NEXT DROP SQ DUP SIZE 2 SWAP
FOR i 1 OVER SIZE 1 + i -
FOR j DUP j DUP i + 1 - SUB ∑LIST ROT + SWAP
NEXT
NEXT + SORT REVLIST
WHILE DUP HEAD PICK3 >
REPEAT TAIL
END REVLIST NIP
»

EVAL

Checksum: # 180Ch
Length: 213.5 bytes

This will return a list with 91 elements, same as above.

[Nihotte(lma)]

(02-14-2021 08:58 PM)Valentin Albillo Wrote:  Hi, all ! Happy San Valentín 2021 to all of you !


Concoction the Third: Weird Integral
[MRM: HP-15C and up]

Consider the following definite Albillo integral:



where Γ is the Gamma function, ln is the natural logarithm (i.e., base e) and φ is the Golden Ratio = (1+ √5)/2.


The Challenge:

Use your HP calc to compute (either manually or writing a program to do it) and output the value of the definite integral as accurately as possible, and then (the sleuthing part) use your HP calc (perhaps conduct some experiments) to try and attempt to answer this nagging question: What's so weird about this integral ?. (Again, forget about Googling for it because I concocted it myself and it's nowhere else to be found either.)

That's all. Enjoy ! ... and that's an order !

V.

Here is the program I used on my HP35s for the Concoction the Third and the Weird Integral:

LBL W
RCL P
x²
RCL X
-
LN
!
ENTER
ENTER // 2 times because ENTER has disabled moves in the stack...
RCL X
LN
!
+
/
RTN

With FN= W
1 and RCL P (where P is 1.61803398875)
∫FN d X
INTEGRATING gives 3.09016994375E-1 and it seems to be 1/(2*φ)

[Valentin Albillo]

      
Hi, all:

First of all, thanks a lot for your interest and your participation in this challenge, really much appreciated. As stated in my OP, I'll post my original solutions and extensive comments in a few days but as there's still a number of questions left unanswered so far, I'll give you all a

       Last Chance

to address and complete them before I do. This includes the following:

Concoction the First: Weird limit

      Point d is still unaddressed:  d. Can you explain the constant component of said [asymptotic] expression ?

Concoction the Second: Weird Sum

      The main question is wholly unaddressed: What's so weird about this sum ?. Perhaps a little sleuthing (i.e.: conducting some experiments) would be of help.

Concoction the Fourth: Weird Graph

      This is mostly unaddressed. In particular no program or description of the operations required to produce the graph have been posted so far and worse, no one has produced and posted the graph itself, i.e.: an actual image of it. This is the kind of functionality for which graphics calculators are intended and I included this part specifically to give HP Prime's or RPL-models' users some opportunity to show off their calculator's worthiness for this challenge.

      Also, no attempts to factorize the polynomial have been reported, either successful or unsuccessful. This is the kind of functionality for which CAD is intended. Why don't you give it a try ? Hint: It would help to answer the main question.

      Finally, as for the main question proper, What's so weird about this graph?, it's still left unanswered. Apart from its so-described "funny" appearance, there's more to the graph than it seems at first sight. Some sleuthing would surely help.

Concoction the Sixth: Weird Year

      Both RPL code and the resulting list of the years have been produced (without comments or explanation), but nearly all the questions have been left unanswered, i.e.:

• What is the "simple but highly remarkable (striking, in fact) numeric property" ?
• (a) How many years will be listed in the output ? ,
• (b) What will be the next predicted potentially catastrophic year after 2020 ?,
• (c) Should we be concerned ?

      Also, although not explicitly asked, no program has been produced to accept a given year in range and demonstrate* whether it has the required numeric property (thus, if it indeed was/might be catastrophic) or not, which would be nice as I state that my original solution does exactly this.

      * E.g.: For property "The year's number is a factorial" and year 720 you would output "720 = 1x2x3x4x5x6" demonstrating the property, while for year 721 you would ouput "721 = not a factorial"

      Finally, programs written in other than RPL would be welcome for variety and to let readers better understand what the code does and how their RPN calcs (say) would deal with the problem.


As stated, Last Chance. Thanks and best regards.  
V.

[robve]

For challenge two (weird sum), 1 / 4.94804552201E-4 = 2021.

Pretty neat!

- Rob