Threaded Mode | Linear Mode

C.Ret · (This post was last modified: 02-28-2021 09:37 AM by C.Ret.)

Ah! Ah!

Many Thanks to Valentin for this good puzzle and the good explanation and pointing out the fact that the graph is not exactly a circle.

I didn't notice it and was wandering why I was unable to factorize Po(x,y) even with the correct equation of the expected, but outfitted circle equation.

By the fact, checking whenever the graph is or isn’t a true circle is so easy by entering its equation in the Advanced Graphic Application of the HP Prime and zoom enough to get a good view of the disparity:

Fig.: Po(x,y) black Circle green.

How Can I have miss verifying this weird fact ?

J-F Garnier · (This post was last modified: 02-28-2021 11:09 AM by J-F Garnier.)

Now that the challlenge is closed, let me add some comments:

(02-28-2021 02:18 AM)Valentin Albillo Wrote: My original solution for "Concoction the Third: Weird Integral"

What's so weird about this integral?

Quote:The Sleuthing

As before, I'll describe a typical sleuthing procedure to try and answer the above question. First of all, let's compute the integral's value to full 12-digit accuracy using the HP-71B, right from the command line:

>DESTROY ALL @ P=(1+SQR(5))/2 [ENDLINE]
>INTEGRAL(1,P,0,GAMMA(LN(P^2-IVAR))/(GAMMA(LN(IVAR))+GAMMA(LN(P^2-IVAR)))) [ENDLINE]

First, would you have used the Enhanced Math ROM, you could have just typed:

>INTEG(1,P,0,GAMMA(LN(P^2-IX))/(GAMMA(LN(IX))+GAMMA(LN(P^2-IX)))) [ENDLINE]
After all, you are the one who suggested the use of such shortcuts Smile

As I mentioned in my solution, the first aspect I noticed was the speed; I wasn't expecting to get a fast answer on a physical calculator of the 32S class for such complicate-looking equation combining the gamma and log functions.
I also noticed that the speed and the answer were independent of the target accuracy (set by the display mode on the 32S - from FIX 1 to FIX 11). That was really weird.
BTW, using a slow machine sometimes highlights some aspects that can be missed on faster systems such as Free42.

Moving to the HP-71/Emu71, I investigated more:

30 Y=0 @ A=0 @ F=(1+SQR(5))/2
60 DEF FNF(X)
70 A=GAMMA(LN(F*F-X)) @ Y=A/(GAMMA(LN(X))+A)
90 PRINT X;Y
100 FNF=Y
110 END DEF
130 E=.0001
140 DISP INTEGRAL(1,F,E,FNF(IVAR))
150 END

>RUN
1.30901699438 .500000000007
1.09656781074 .169413871978
1.52146617801 .830586128022
1.19554981675 .328227901097
1.422484172 .671772098903
1.02655614795 4.82413106951E-2
1.5914778408 .951758689306

It turned out that only 7 samples were used, whatever was the target accuracy set by the E value.
Having some knowledge of the Romberg algorithm, I recognized the first sample to be the midpoint, then the next two samples used for the 1st estimation and the last four samples used for the 2nd estimation.
It was quickly obvious that the sample pairs (2-3, 4-5, 6-7) gave values that sum to exactly 1.
That gave me the key for the explanation of the speed: the summation for any sample set was the same (and the exact one), so the algorithm reached its termination condition after the first two estimations.
Extrapolating to a single sample at the midpoint, I got the symbolic answer: (φ-1) x 1/2

This is of course a consequence of the symmetry of the function to integrate as Werner and Valentin explained.
It is interesting to note that my observations had also to do with the HP implementation of the Romberg algorithm: although the samples are non-uniform, they are symmetric relative to the midpoint, and that gave the effect I reported.

Thanks for the challenge that was at the intersect of my interests for math puzzles and numerical algorithms (such as used in HP machines).

J-F

Albert Chan · 02-28-2021, 01:21 PM

(02-18-2021 09:28 AM)Werner Wrote: 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)

Without noticing the symmetry, if we "fold" the integral, we get the same result.

[Image: TEST6-DISREGARD.jpg]

Let f(x) be integrand of above integral.

I = ∫(f(x), x = 1 .. φ)
= ∫(f(x) + f(1+φ-x), x = 1 .. φ) / 2
= ∫(f(x) + f(φ²-x), x = 1 .. φ) / 2
= ∫(1, x = 1 .. φ) / 2
= (φ-1)/2

robve · (This post was last modified: 03-02-2021 04:03 PM by 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

robve · 03-01-2021, 03:03 AM

(02-28-2021 02:18 AM)Valentin Albillo Wrote: [*]robve computed the integral using the HP Prime but alas, the result he posted is wrong. He also posted a BASIC program for the SHARP PC-1350 which produced the same wrong result. Anyway, thanks for trying ...

Ah, that explains why I was getting nowhere with this one! Thanks for debugging my answer now that this challenge is over. Garbage in, garbage out: by putting the wrong expression in the integral the integration area did not offer any insights in what is weird about this integration. Somehow I missed the LN in the numerator of the integrand. The rest of the integrand expression was correct. Go figure...

With the LN correction applied to the numerator, I noticed that the Romberg integration ran much faster this time on the SHARP PC-1350 compared to the incorrect integrand. When viewing the value of I (the iterator) I noticed that I=2 when the integration converges. Two successive trapezoidal approximations at alternating points are indistinguishable within MachEps. This means that the curve should be sufficiently close to linear between 1<=x<=phi. Plotting the integrand on the HP PRIME shows that this hunch is indeed the case:

Observing 0<=f(x)=(x-1)/(phi-1)<=1 for 1<=x<=phi. Integrating this f gives the same result 0.309016994375. Analytically, the value of the integral of f between 1 and phi is:

$$\frac{\phi^2-2\phi+1}{2\phi-2}=\frac{(1-\phi)^2}{2\phi-2}=-\frac{1}{2}\frac{(1-\phi)^2}{1-\phi}=\frac{\phi-1}{2}$$

Noticing this behavior of the quadrature algorithm is now (a bit late) helpful to answer what is weird about this integral. I am glad to use my trusty (albeit crusty) Romberg on the PC-1350.

Also, what is the fun of doing math and calc exercises if we don't implement numerical integration ourselves?

The Romberg quadrature program on the PC-1350 produces 0.3090169944 (10 digits exact) with the change X=LN(X*X-V) to line 300 to correct the integrand:

300 "F1" V=X,X=1.618033989,X=LN(X*X-V): GOSUB "GAMMA": W=Y

- Rob

robve · 03-01-2021, 07:05 PM

One more observation that may be of interest. Based on the visual clue that the integrant behaves as (x-1)/(phi-1), but is not the same, it is prudent to numerically verify that the function (x-1)/(phi-1) equals or closely approximates the integrand

$$ \frac{\Gamma\ln(\phi^2-x)}{\Gamma\ln x+\Gamma\ln(\phi^2-x)} $$

This can be done by listing the results of the two functions side-by-side for x between 1 and phi in increments of say 0.1 to compute the absolute or relative error of the difference between the two. We can also take the absolute difference of the two functions as the integrand and numerically evaluate the integral:

$$ \int_1^\phi \left|\frac{\Gamma\ln(\phi^2-x)}{\Gamma\ln x+\Gamma\ln(\phi^2-x)}-\frac{x-1}{\phi-1}\right| dx \approx 5.66308603939\times10^{-3}$$

The value is the result on the HP PRIME. My PC-1350 Romberg integrator gives the same 5.66308E-3. Therefore, the total area of the difference between the two functions is small, but not zero. Based on the numerical results, it does not appear convincing to state that

$$ \frac{\Gamma\ln(\phi^2-x)}{\Gamma\ln x+\Gamma\ln(\phi^2-x)} \stackrel{?}{=} \frac{x-1}{\phi-1} $$

They are approximate, but not equal. Function (x-1)/(phi-1) under estimates the integrand on the interval 1<=x<=phi/2 and overestimates the integrand on interval phi/2<=x<=phi by about 0.014 symmetrically:

Because of this symmetry, one could argue that

$$ \int_1^\phi \frac{\Gamma\ln(\phi^2-x)}{\Gamma\ln x+\Gamma\ln(\phi^2-x)} dx = \frac{\phi-1}{2} \approx 0.309016994375 $$

Furthermore, the integral of the difference is near zero

$$ \int_1^\phi \frac{\Gamma\ln(\phi^2-x)}{\Gamma\ln x+\Gamma\ln(\phi^2-x)}-\frac{x-1}{\phi-1} dx \approx -1.13498070108\times10^{-13} $$

Suggesting that the over and under errors of the (x-1)/(phi-1) curve nicely cancel out over 1<=x<=phi.

- Rob

Valentin Albillo · 03-02-2021, 02:55 AM

Hi all,

These are my original solutions for "Concoction the Fifth: Weird Primes" and "Concoction the Sixth: Weird Year", which concludes all my original solutions for the six Concoctions in this S&SMC #25. You can follow these links to see my previously posted original solutions to:

      "Concoction the First: Weird Limit"
      "Concoction the Second: Weird Sum".
      "Concoction the Third: Weird Integral" and "Concoction the Fourth: Weird Graph"

Note: My HP-71B code might use keywords from the JPC ROM, MATH ROM, HP-IL ROM and STRINGLX LEX file, executed on go71b, while RPN code is for the HP-42S, executed on a DM42.

My original solution for "Concoction the Fifth: Weird Primes"

      "Consider a prime number so 'Perfectly Prime' (a PP for short) that changing any single [base 10] digit would turn it into a composite number. Write a program 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 the second PP greater than 666,666,666."

What's so weird about these primes?

The Sleuthing

The first thing to do is to write a program which will find these PP starting from any given integer, and this straightforward 4-line HP-71B program will nicely do:

      1  DESTROY ALL @ INPUT K,N @ FOR P=1 TO K
      2  N=FPRIM(N+2) @ N$=STR$(N) @ FOR I=1 TO LEN(N$) @ M$=N$ @ FOR D=0 TO 9
      3  M$[I,I]=STR$(D) @ IF M$=N$ THEN 4 ELSE IF NOT PRIM(VAL(M$)) THEN 2
      4  NEXT D @ NEXT I @ DISP P;N @ NEXT P

It asks how many PP to output and the lower limit to begin the search from. Let's run it to output the 5 smallest PP, as asked. We begin at 11 as obviously no single-digit prime can be a PP:

      [RUN]  ? 5,11 [ENDLINE]

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

As for the remaining questions, we proceed likewise:

      [RUN]  ? 1,5E8       [ENDLINE]  ->   1  500004469    { PP #1318 }
      [RUN]  ? 1,777777777 [ENDLINE]  ->   1  777781429    { PP #2259 }
      [RUN]  ? 2,666666666 [ENDLINE]  ->   1  666850699    { PP #1845 }
                                           2  666999929    { PP #1846 }

The Results

Once the sleuthing's over, the results are as follows:

The 5 smallest PP are: 294001, 505447, 584141, 604171 and 971767.
The first PP > 500 million is 500004469.
The first PP > 777,777,777 is 777781429.
The second PP > 666,666,666 is 666999929.

The Comments

Though for a prime being a PP seems a rare occurrence (the very first one is 294001; that any do actually exist is weird), it's been proved that in fact there's an infinity of PP and what's more, a positive proportion (in terms of asymptotic density) of the primes are PP. That's weirder !

We can try to (very roughly) "guesstimate" said proportion: there are 5,761,455 primes and 334 PP less than 100 million, so the rate comes out as 0.00580 %, i.e.: one PP per ~ 17,200 primes. Going for the next order of magnitude, there are 50,847,534 primes and 3,167 PP less than 1 billion, so the rate now comes out as 0.00623 %, i.e.: one PP per ~ 16,000 primes, hence the proportion, though very small, seems to be increasing (and in any case, it's asymptotically > 0 .)

Furthermore, if that wasn't weird enough, there's also a positive proportion of primes which are PPP (Perfect Prime Plus), which have the property that if any single digit is changed (including any zero among the prime's infinitely many leading zeros), then the resulting number is always composite. That's weirdest !!

For instance, notice that the first PP, 294001, is not a PPP since changing, say, the second leading zero of ...000294001 results in ...010294001, which is prime. Matter of fact, even though a (much smaller, but still) positive proportion of the primes are PPP, none are yet known.

Let's end these Comments by listing some peculiar PP you may find useful to check or optimize your code:

The first five PP found above are all the PP less than 1 million. The next five are 1062599, 1282529, 1524181, 2017963 and 2474431.
PP #24 and PP #25 are consecutive and extremely close: 7469789 and 7469797, respectively.
PP #51 and PP #52 are consecutive and both look like nearby date ranges: 1985_1991 and 2021_2327, resp.
PP #1404 to PP #1414 are consecutive and end in 1,1,1,1,7,7,1,1,7,7,1, resp.
PP #1846 and PP #2539 have only 3 different digits each: 666999929 and 844448333, resp.
PP #3048 = 969094909, has the odd digit 9 in all odd positions and even digits in all even positions.

So much for Perfect Primes but there's more: you might want to try your hand at these four lingering questions (for which I won't provide answers, you'll be on your own):

Are there any Twin PP, i.e: two PP whose difference is 2 ?
Are there any PP which remain so (i.e.: the result is always composite) when changing any two digits ?
Find a PPP. There's an infinity of them and you'll get worldwide fame in the math field if you do,
What about Composites ? What's the first composite (not divisible by 2 or 5) so "Perfectly Composite" [PC] that it remains composite when you change any single [base 10] digit ? (hint: it begins with 2 and ends with 9)

The Hall of Fame

This time the one expert which dealt with this Concoction the Fifth: Weird Primes is:

robve, who wrote HPPL code for the HP Prime and found the correct answers to all four questions asked. He also wrote a "Cheat program" (his own words, mine too) in C which ran much faster (D'oh!), but more importantly, he gave some theoretical considerations to estimate the chance that a prime number is a PP and correctly concluded that "it seems reasonable to see perfect primes for large k [digits] and there are infinitely many of them".

Afterwards he posted an alternative submission in HPPL and augmented his sleuthing above and beyond the call of duty, producing a list of the first PP for every base from 2 to 19, then another list with the PP counts for the same bases, commenting on the results and the interesting pattern obtained. That's what sleuthing is all about !

Finally he also said that "It seems odd to me to change the leading digit to 0 to check for PP" but the reason for that is made clearer in the light of the existence of PPP, as stated in the Comments above.

My original solution for "Concoction the Sixth: Weird Year"

      "2020 shares a very striking numeric property with many other catastrophic years [...] try and discover what simple numeric property {which can be unambiguously stated by saying that the year's number "is a (five words)"} is shared by all the aforementioned numbers, and then write a program to output a listing of all years between AD 4 and AD 5000 (both included) which have this very property.

      Additional 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 ?

The Sleuthing

Here, the very first thing to do is to discover the "striking numeric property" in question, which isn't as difficult as it seems because there are only so many such properties which can be unambiguously stated in just five words, e.g.: the number is a "sum of some prime numbers" does fit, but this is far too generic because every integer > 1 has that property. Many variations are possible (say "5" instead of "some" or replacing "primes" by "squares", "cubes", "factorials", etc.)

One very useful heuristic is to analyze possible properties of the smallest numbers given: we first discover some property they have in common and then check if it applies to the bigger numbers as well. The numbers stated as sharing the sought-for property are 458, 662, 666, 1348, 1556, 1849 and 2020, so the smallest number is 458 and the sleuthing process begins with it.

After discarding many sums of squares and cubes (e.g.: "sum of four square numbers", as every integer is a sum of four square numbers, including 0²) and "sum of some prime numbers" and variations, we eventually discover that 458 = 13² + 17², a sum of exactly two non-zero squares, but regrettably 662 is not. We also notice that the numbers 13 and 17 are primes so we refine the property to "sum of some squared primes" but again that's too generic and gets us nowhere.

Cutting to the chase, eventually we also notice that 13 and 17 are consecutive primes so the property becomes "sum of consecutive squared primes" (5 words!) and this time we hit the jackpot: 458 = 13² + 17² and 662 = 3² + 5² + 7² + 11² + 13² + 17², and checking the bigger numbers we readily get:

      666   = 2² + 3² + 5² + 7² + 11² + 13² + 17²
      1348 = 13² + 17² + 19² + 23²
      1556 = 2² + 3² + 5² + 7² + 11² + 13² + 17² + 19² + 23²
      1849 = 43² { a sum with just one summand is also a sum; the empty sum is 0 }
      2020 = 17² + 19² + 23² + 29²

so bingo !. Now it's just a matter of writing a program to find all the years in the given interval having that property, and this 7-line program for the HP-71B will do:

      1  DESTROY ALL @ OPTION BASE 1 @ DIM P(19) @ K=2 @ L=0 @ C=0
      2  REPEAT @ L=L+1 @ P(L)=K*K @ K=FPRIM(K+1) @ UNTIL K>70
      3  FOR N=4 TO 5000 @ J=L @ WHILE P(J)>N @ J=J-1 @ END WHILE
      4  M=N @ I=J
      5  M=M-P(I) @ IF NOT M THEN C=C+1 @ DISP N; @ GOTO 7
      6  IF M<0 THEN J=J-1 @ GOTO 4 ELSE I=I-1 @ IF I THEN 5
      7  NEXT N @ DISP @ DISP "Total:";C

      [RUN]
                 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

              Total: 91

where the years given as examples appear in bold and the next one after 2020 appears in bold red, which is 2189 = 13² + 17² + 19² + 23² + 29².

The Results

Based on the data obtained by the sleuthing process above and the results from running the program, the answers are:

The simple numeric property is: The years' numbers are a "sum of consecutive squared primes" (five words)
The listing of all years between AD 4 and AD 5000 (both included) which have this property is the output above, 91 years in all.
What will be the next predicted potentially catastrophic year after 2020 ?: It will be 2189, as seen in the output.
Should we be concerned ?: Well, with ever-advancing technology one can never say for sure, but I very much doubt any people reading this in 2021 will be kicking and alive by 2189, so no worries !

The Comments

The joke explanation of why 2020 was such a catastrophic year is an original idea of mine, and implementing it was just a matter of finding some nice but simple numeric property of 2020. After some sleuthing I found it to be that 2020 = 17² + 19² + 23² + 29², which are the squares of consecutive primes (pretty nice indeed), and then I wrote the program to compute all other years between AD 4 and AD 5000 sharing this property. So far, so good.

Then I searched Wikipedia for a list of big catastrophes, and cross-checking the years found there with the years listed by my program I finally selected the ones most remarkable while ignoring numbers too low, to avoid reducing the difficulty too much (e.g.: if chosing AD 13 it would be instantly recognizable as 2² + 3², all too easy!) and le voilà, Concoction 6 ready !

I also wrote the following 6-line program for the HP-71B which accepts a given year in range and demonstrates whether it has the required numeric property (thus, if it indeed was/might be catastrophic) or not.).

      1  DESTROY ALL @ OPTION BASE 1 @ DIM P(19),S$[80] @ K=2 @ L=0
      2  REPEAT @ L=L+1 @ P(L)=K*K @ K=FPRIM(K+1) @ UNTIL K>70 @ INPUT N
      3  J=L @ WHILE P(J)>N @ J=J-1 @ END WHILE
      4  S$="" @ M=N @ I=J
      5  M=M-P(I) @ S$=STR$(SQR(P(I)))&"^2+"&S$ @ IF NOT M THEN S$[LEN(S$)]="" @ DISP S$ @ END
      6  IF M<0 THEN J=J-1 @ GOTO 4 ELSE I=I-1 @ IF I THEN 5 ELSE DISP "Not a sum"

      [RUN]
            ?  2020   ->   17^2+19^2+23^2+29^2                 { VAL(S$) = 2020 }
            ?  666    ->   2^2+3^2+5^2+7^2+11^2+13^2+17^2      { VAL(S$) = 666 }
            ?  1849   ->   43^2                                { VAL(S$) = 1849 }
            ?  555    ->   Not a sum

The Hall of Fame

Again, the one and only expert who dared to deal with this Concoction the Sixth: Weird Year is no other but

Gerson W. Barbosa, who posted RPL code for such models as the HP50G and BASIC code for the HP-75C, both of which correctly produced the list of "weird" years, and he also gave their total number (91) and stated explicitly the remarkable property all these numbers share. Alas, he didn't post any details on his sleuthing, particularly how he managed to find the correct "remarkable and striking numerical property", which would be fascinating for sure but never mind, Well Done !

That's all, this concludes my original solutions for all 6 Concoctions 6 of this S&SMC #25 of mine. Thank you very much to those who contributed their solutions or at least posted some comments, much appreciated.

I really hope you enjoyed it as well as all the readers in general. Over and out.

Note: For any comments not directly related to the math or code here but to ancillary matters such as this or that opinion on the rules or "Halls of Fame" or whatever, please PM me instead of posting them here. Let's keep this thread strictly mathematical and algorithmical in nature. Thanks.

Best regards.
V.

EdS2 · (This post was last modified: 03-03-2021 09:14 AM by EdS2.)

Very nice, thanks again for the motivations and explanations.

Just one query: I'd interested to know the runtimes of some of these 71B programs. (Even if approximate, from memory.) For example, the time to find PP#1 and then PP#2. Another example, the most expensive calculation relating to the PP challenge - how long did that run for?

(If I had a 71B I could find this out directly!)

(Edit: answered in PM. VA (re)noting the use of go71b, about 128x faster than 71b, and taking 186s and 337s for PP's 1 and 2, and 4117s for PP# 1846. Noting also that a real 71b is about 10x or more the speed of a 41C. HP85 about 5x faster than 71. And noting that none of the code is written for speed.)

robve · 03-02-2021, 09:32 PM

Very nice last challenge and results!

The combinatorial space of possible methods to apply seems quite large, perhaps too large for some to consider attacking this challenge when given limited time to work on this, despite helpful hints such as “five words”.

(03-02-2021 02:55 AM)Valentin Albillo Wrote: After discarding many sums of squares and cubes (e.g.: "sum of four square numbers", as every integer is a sum of four square numbers, including 0²)

I for sure could have guessed this challenge had something to do with a sum of squares for a specific year^*.

(03-02-2021 02:55 AM)Valentin Albillo Wrote: Cutting to the chase, eventually we also notice that 13 and 17 are consecutive primes so the property becomes "sum of consecutive squared primes" (5 words!)

OK (duh), it happens that with many such methods other “special years” will coincidently fall in place too, with the often sought after number of the beast (which according to Numberphile is debatable whether it is 666 or 616 because both spell out Nero, the most evil ruler of the Roman era). If the solutions to the equation are very limited, then it won’t be likely. But if there are a couple of hundred solutions to an equation in a limited integer range and the definition of "special year" is less specific, then the odds are in favor of matching some special years.

The result is nice because of the series of primes is consecutive. Not only does this reduce the number of solutions reasonably, and thus the years that match, it also saves time not having to check all subsets of a certain set of integers (such as squares of primes). For example, checking all subsets of integers between 1 and sqrt(Y) for a given year Y is an O(2^Y/2) time process to compute. Reducing it to linear O(Y) makes it not only fast, but much simpler to implement (on a calculator). The method is linear in Y to generate primes up to sqrt(Y) and to check if the sum of squares of consecutive primes 2 to sqrt(2) is equal to Y (i.e. in time O(sqrt(Y²))=O(Y)).

Now that the results are in, I will share an effective approach to tackle the kinds of questions that require an exhaustive search among subsets. In general, we can approach these type of problems using the classic generate-and-test paradigm. Essentially brute force, but a bit smarter to generate combinations with a classic combination generating algorithm.

Listed below is a generic generate-and-test HPPL program to generate and check all subsets of a set of integers that sum up to a given integer, with a modification to sum up squares of primes to match a given year. First, the values are generated with GENVALUES and stored in list L1. Second, all subsets of the list are then tested to verify if their sum of squares match the given number NUM:

EXPORT GENVALUES(N)
BEGIN
LOCAL I;
  L1:=[2];
  FOR I FROM 3 TO N DO
    IF isprime(I) THEN
      L1:=append(L1,I);
    END;
  END;
END;

EXPORT SUMSQ(NUM)
BEGIN
  LOCAL N=FLOOR(SQRT(NUM)),K,I,F=0,X,Y;
  GENVALUES(N);
  FOR K FROM 1 TO SIZE(L1) DO
    L0:=MAKELIST(I,I,1,K);
    REPEAT
      L2:=MAKELIST(L1[L0[ I]],I,1,K);
      IF ΣLIST(L2^2)=NUM THEN
        PRINT(L2);
        F:=1;
      END;
      PNXCB(SIZE(L1),K,X,Y);
    UNTIL L0(K)=K;
    // IF F THEN
    //  RETURN;
    // END;
  END;
END;

Note that L0:=MAKELIST(I,I,1,K) is an index vector with values 1 to K. This vector is used to generate all subsets of values though indirect indexing, with elements of L0 pointing to elements of L1. For example, L0=[1,2,3] is the subset of the first three primes L1=[2,3,5,...] whereas L0=[1,4,7] is the subset of three primes L1=[2, , ,7, , ,17,…]. The list L2:=MAKELIST(L1[L0[ I]],I,1,K) produces this list of subset of primes, which we then check if their sum of squares equals NUM with ΣLIST(L2^2)=NUM. If so, the list L2 satisfies the conditions and is displayed.

Also note that we start with the smallest subsets first, a set of size 1, then a set of size 2, and so on by populating L0:=MAKELIST(I,I,1,K) for K from 1 to the size of L1.

The program displays the smallest subsets only. If all sets should be displayed, then uncomment the three lines in SUMSQ. To check with another function rather than squaring, just replace L2^2 with another function and adjust sqrt(NUM) to assign to N accordingly (we limit the iteration space to sqrt(NUM) to avoid iterating over values who’s square exceeds NUM and therefore are never part of the solution).

The workhorse for this program is the clever PNXCB combination generator which I’ve rewritten below in HPPL:

EXPORT PNXCB(N,K,X,Y)
BEGIN
  LOCAL J=1;
  WHILE 1 DO
    X:=L0(J);
    Y:=X+1;
    IF J=K THEN
      IF Y>N THEN
        Y:=K;
        L0(J):=Y;
        RETURN;
      END;
      L0(J):=Y;
      IF J>1 THEN
        L0(J-1):=X;
        X:=J-1;
      END;
      RETURN;
    END;
    IF Y<L0(J+1) THEN
      L0(J):=Y;
      IF J>1 THEN
        L0(J-1):=X;
        X:=J-1;
      END;
      RETURN;
    END;
    J:=J+1;
    X:=L0(J);
    Y:=X-1;
    IF Y>=J THEN
      L0(J):=Y;
      IF J>1 THEN
        Y:=J-1;
        L0(J-1):=Y;
      END;
      RETURN;
    END;
    IF J=K THEN
      Y:=N;
      L0(J):=Y;
      RETURN;
    END;
    J:=J+1;
  END;
END;

We start with an index list L0 containing the values 1 to K. PNXCB takes L0 as input. The next combination of the binomial C(N,K) (N choose K) is generated by PNXCB and is stored in L0.

PNXCB generates exactly C(N,K) selections K of N in L0. But we do not need to count them to stop the REPEAT loop! Here is the clever part: if we start with a list of values 1 to K then the PNXCB cycles through all combinations to arrive back to the initial list, exactly after C(N,K) calls to PNXCB. Therefore, the C(N,K) cycle is complete when the last value in the list is K, a useful property of PNXCB.

The PNXCB routine rewritten for classic BASIC calculators:

' Algorithm PNXCB Combination Generators
' PNXCB(P,N,K) cycles through "pointers" P(1:K) to N items
' Calling PNXCB NCR(N,K) times produces all combinations
' If initiallly P(I)=I, then NCR(N,K) cycle is complete when P(K)=K

' TEST PROGRAM
10 CLEAR: N=5,K=2: INPUT "N=";N,"K=";K
20 DIM P(K): FOR I=1 TO K: P(I)=I: NEXT I
30 FOR I=1 TO K: PRINT P(I);: NEXT I: PRINT
40 GOSUB 100
' all pointers were left-most initialized, if P(K)=K then end
50 IF P(K)=K END
60 GOTO 30

' PNXCB(P,N,K) -> P, X (old), Y (new), changes J
100 J=1: REM GOTO 180 FOR DOWN-UP SEQUENCE

110 X=P(J),Y=X+1
120 IF J<>K GOTO 160
130 IF Y>N LET Y=K,P(J)=Y: RETURN

140 P(J)=Y: IF J>1 LET P(J-1)=X,X=J-1
150 RETURN

160 IF Y<P(J+1) GOTO 140
170 J=J+1
180 X=P(J),Y=X-1
190 IF Y>=J GOTO 230
200 IF J=K LET Y=N,P(J)=Y: RETURN
210 J=J+1
220 GOTO 110

230 P(J)=Y: IF J>1 LET Y=J-1,P(J-1)=Y
240 RETURN

It's all basic, isn't it?

- Rob

^*) the year 2021 has 17 minimal sets of 3 squares that sum up to 2021. Run SUMSQ with GENVALUES changed to L1:=MAKELIST(I,I,1,N) and uncomment the IF F THEN clause to display these minimal solutions.

Gerson W. Barbosa · (This post was last modified: 03-04-2021 02:59 PM by Gerson W. Barbosa.)

(03-02-2021 02:55 AM)Valentin Albillo Wrote: Gerson W. Barbosa, who posted RPL code for such models as the HP50G and BASIC code for the HP-75C, ...
...
Alas, he didn't post any details on his sleuthing, particularly how he managed to find the correct "remarkable and striking numerical property", which would be fascinating for sure but never mind, Well Done !

Hello, Valentín,

Thank you for the compliment and sorry for the belated reply (I’ve been busier than usual lately).

Unlike you, didactics is not one of my strengths, but I will try to demonstrate how I arrived at the solution to your sixth concoction. At first I was not intending to participate in any of these as I was 400 kilometers away from home and I had brought no calculator along. When I visited this thread again, all the other concoctions had been solved. It was late in the evening already, I had neither pencil nor paper but I found a copybook and a pen somewhere and decided to give it a go. As I am not smart enough to look at a small set of numbers and spot its common property right away, I submitted all your examples but one to Wolfram Alpha, which told me that
2020 = 16² + 42² = 24² + 38²;
1348 = 18² + 32²;
1556 = 20² + 34²;
458 = 13² + 17²;
1849 = 43² + 0² and
666 = 15² + 21².
“Quite easy!”, I thought. These are just numbers that can be represented as a sum of two squares (although that didn’t look like a ‘striking’ property to me). Then I wrote a simple RPN program on the emulated 49G+ on my smartphone running m48+:

« { } 0 71
FOR x x 71
FOR y x SQ y SQ + +
NEXT
NEXT SORT DUP SIZE 1 - 1 SWAP
FOR k k GETI UNROT GETI NIP ROT == { k { 0 } REPL } IFT
NEXT SORT
WHILE DUP HEAD NOT
REPEAT TAIL
END
»

After a rather long time, it returned a list with way more elements than the one specified by you (less than 100 elements), even when discarding the ones greater than 5000. Also, one of your examples were missing: 662. It turns out that was the only one I had left out.

I was disappointed of course, but not everything was lost. I noticed that two of your examples involved either the some of the squares of consecutive primes or the square of a prime: 458 and 1849. The Wikipedia article on 666 only confirmed that 666 = 2² + 3² + 5² + 7² + 11² + 13² + 17². Likewise, it wasn’t difficult to verify that all the remaining examples could be expressed as the sum of the squares of consecutive prime numbers. Now, I only needed to figure out a suitable algorithm. But at daft o’clock my reasoning was starting to fail, so I went to sleep and early next morning after I woke up I soon came out with a working program. The argument to the program in post #15 is the number of elements in the list of the squares of the first consecutive prime numbers, { 4, 9, 25, 49, 121, ... }. I started with a 15-element list, but I noticed the resulting list of ‘weird years’ would not include some of your examples. Then I increased it to 16, then to 20, when the list seemed to be complete (the list up to 5000 would not be changed by further increases). I thought it would be nice if the program calculated the size of the basic list automatically, so I researched at OEIS, which led me to an interesting recent article on the subject. That’s where I ‘borrowed’ the formula I used in my next programs, which has to do with the size of the list of the first squares of primes necessary to produce the list of all years up to certain limit (hopefully I have simplified the correct formula – the article looks complicated to me). Next day I was tired after driving all the way back home, but I took the time to write an HP-75C version per your suggestion of using programming languages other than RPL for more clarity. But it payed off, as I noticed the BASIC program was more optimized than the RPL program (the output list never would grow longer than it needed to be). As a result, I was able to write a faster RPL program, but really expert RPL programmers, which I am not, can optimized it even more.

Thank you also for taking the time to concoct these S&SMCs. Even when we don’t participate it’s a joy just to read them, especially your final comments and solutions.

Best regards,

Gerson.