Threaded Mode | Linear Mode

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.