// Magma program to compute terms of OEIS sequence A363245. // // Program runs in about 40 seconds on Magma Online Calculator at // http://magma.maths.usyd.edu.au/calc/ // // Jon E. Schoenfield, Jun 10 2023 a := []; // terms of the sequence nMax := 36; // number of terms to compute ssMaxStore := 46000000; // size of Boolean array for storing the // status of integers as subset sums of terms // (must be at least Sum_{n=1..nMax-1} a(n)) // Initialize Boolean array: NSS[1 + s] = true iff s is NOT a subset sum of // the terms found thus far. (Note on use of "1 +": array values are offset // by 1 so as to allow the array to include an element for s = 0; I'm told // that Magma does not allow arrays to have an element with index 0.) NSS := [true : sPlus1 in [1 .. ssMaxStore + 1]]; NSS[1 + 0] := false; // the empty sum is a subset sum sumOfTerms := 0; // sum of known terms squarefreeProduct := 1; // product of distinct prime factors of known terms for n in [1 .. nMax] do if n eq 1 then aTestMin := 3; // 1 = 2^0 and 2 = 2^1 are excluded else aTestMin := a[n - 1] + 1; end if; for aTest in [aTestMin .. ssMaxStore] do if Gcd(aTest, squarefreeProduct) eq 1 then bOK := true; for e2 in [Ilog(2, aTest - 1) + 1 .. Ilog(2, sumOfTerms + aTest)] do if not NSS[1 + 2^e2 - aTest] then bOK := false; break; end if; end for; if bOK then a[n] := aTest; if n eq nMax then a; break n; end if; for f in Factorization(aTest) do squarefreeProduct *:= f[1]; end for; for sPlus1 in [1 + sumOfTerms + aTest .. 1 + aTest by -1] do if NSS[sPlus1] then if not NSS[sPlus1 - aTest] then NSS[sPlus1] := false; end if; end if; end for; sumOfTerms +:= aTest; break; end if; end if; end for; end for;