// OEIS A092671: // "Numbers n such that there exists a solution to the equation // 1 = 1/x_1 + ... + 1/x_k (for any k), 0 < x_1 < ... < x_k = n." bOutputForbFile:=true; bShowSolns:=false; itos:=IntegerToString; nMax:=1810; // maximum n to test // (running this program on the Magma Calculator at // http://magma.maths.usyd.edu.au/calc/ // yields 1000 terms in ~0.3 sec) //nMax:=15934; // maximum n to test // // (running this program on the Magma Calculator at // // http://magma.maths.usyd.edu.au/calc/ // // yields 10000 terms in ~4 sec) mStoreMax:=14; // 23 causes the Magma Calculator to run out of memory // Minimum n at which mStoreMax is insufficient to prove whether n // is or is not in the sequence, for mStoreMax = 1..22: // [ 10, 21, 44, 78, 152, 259, 423, 1199, 2453, 4188, 8359, 8359, 15472, 37095, // 52144, 102014, 133220, 180622, 427161, 427161, 678696, 979238] // Build a table of all the subset sums of reciprocals of integers 1..mStoreMax; // the m-th row will contain all the subset sums that include 1/m in the sum. s:=[0/1]; // empty sum, stored as a fraction for m in [1..mStoreMax] do r:=1/m; d:=2^(m-1); // offset for new row of table for i in [1..d] do s[i+d]:=s[i]+r; end for; end for; provedInCount:=0; provedOutCount:=0; failureCount:=0; // Special case: n=1 n:=1; bIsInSeq:=[true]; // 1 is in sequence A092671 provedInCount+:=1; termCount:=1; // 1st term in the sequence if bOutputForbFile then termCount, n; // output to b-file end if; for n in [2..nMax] do if IsPrimePower(n) then bIsInSeq[n]:=false; provedOutCount+:=1; else // find largest prime power that divides n F:=Factorization(n); pAtPPmax:=F[1][1]; PPmax:=pAtPPmax^F[1][2]; for i in [2..#F] do p:=F[i][1]; PP:=p^F[i][2]; if PP gt PPmax then PPmax:=PP; pAtPPmax:=p; end if; end for; m:=n div PPmax; if m le mStoreMax then bSolnProbable:=false; // init bSolnImpossible:=true; // init for i in [2^(m-1)+1..2^m] do if Numerator(s[i]) mod pAtPPmax eq 0 then i1stGoodSumInS:=i; bSolnProbable:=true; bSolnImpossible:=false; break; end if; end for; if bSolnImpossible then provedOutCount+:=1; end if; else bSolnProbable:=false; // init bSolnImpossible:=false; // can't rule it out using the table s[] rm:=1/m; for i in [1..#s] do t:=s[i]+rm; if Numerator(t) mod pAtPPmax eq 0 then i1stGoodSumInS:=i; bSolnProbable:=true; break; end if; end for; if not bSolnProbable then ""; "First n at which mStoreMax is too small: ", n; break; end if; end if; if bSolnImpossible then bIsInSeq[n]:=false; else // search for a solution (my conjecture: a solution exists) bSolnFound:=false; sum:=1/n; xSet:={n}; // set of integers in denominators of unit fractions summed bSuccess:=false; if PPmax^2 gt n then // (empirically, this threshold seems to work okay) // use the first good result found in table s[] // to get a set of integers (including n) // whose sum of reciprocals does not have PPmax // as a factor in the denominator sumBest:=0/1; xSetBest:={}; if m gt mStoreMax then // m isn't encoded in i1stGoodSumInS-1 xAdd:=m*PPmax; xSetBest:=Include(xSetBest, xAdd); sumBest+:=1/xAdd; end if; t:=i1stGoodSumInS-1; mm:=0; while t gt 0 do mm+:=1; if t mod 2 eq 1 then xAdd:=mm*PPmax; xSetBest:=Include(xSetBest,xAdd); sumBest+:=1/xAdd; end if; t:=t div 2; end while; if Numerator(sumBest) eq 1 then DBest:=Denominator(sumBest); if bIsInSeq[DBest] then bSuccess:=true; elif DBest mod 2 eq 0 then // denom is even if IsPrimePower(DBest) then // denom is pwr of 2 bSuccess:=true; end if; end if; end if; if not bSuccess then sum:=sumBest; xSet:=xSetBest; end if; end if; if not bSuccess then while sum lt 1 do xSet0:=xSet; D0:=Denominator(sum); F0:=Factorization(D0); FF0:=F0[#F0]; PMax0:=FF0[1]; PMaxPwr0:=PMax0^FF0[2]; mMax:=n div PMaxPwr0; FBest:=F0; xSetBest:=xSet0; sumBest:=sum; bSuccess:=false; // It looks as though (for n <= 90000, anyway) adding more // than one reciprocal is rarely needed, but when it is, // three are needed, and more than three are never needed. // (But this depends on the empirical "if PPmax^2 gt n" // threshold above.) // try adding one reciprocal for m1 in [1..mMax] do xAdd1:=m1*PMaxPwr0; if not xAdd1 in xSet0 then sumTest:=sum+1/xAdd1; DTest:=Denominator(sumTest); if DTest mod PMaxPwr0 ne 0 then if DTest eq 1 then bImproved:=true; else bImproved:=false; // init FTest:=Factorization(DTest); for kTest in [#FTest..1 by -1] do kBest:=kTest+#FBest-#FTest; if kBest eq 0 then break; elif FTest[kTest][1] lt FBest[kBest][1] then bImproved:=true; break; elif FTest[kTest][1] gt FBest[kBest][1] then break; else // i.e., FTest[kTest][1] = FBest[kBest][1] if FTest[kTest][2] lt FBest[kBest][2] then bImproved:=true; break; end if; end if; end for; end if; if bImproved then xSetBest:=Include(xSet0,xAdd1); FBest:=FTest; sumBest:=sumTest; if Numerator(sumBest) eq 1 then DBest:=Denominator(sumBest); if bIsInSeq[DBest] then bSuccess:=true; elif DBest mod 2 eq 0 then // denom is even if IsPrimePower(DBest) then // dnm is pwr of 2 bSuccess:=true; end if; end if; if bSuccess then break; end if; end if; end if; end if; end if; end for; if bSuccess then break; end if; if sumBest eq sum then // no improvement // try adding two reciprocals for m1 in [1..mMax-1] do xAdd1:=m1*PMaxPwr0; if not xAdd1 in xSet0 then for m2 in [m1+1..mMax] do xAdd2:=m2*PMaxPwr0; if not xAdd2 in xSet0 then sumTest:=sum+1/xAdd1+1/xAdd2; DTest:=Denominator(sumTest); if DTest mod PMaxPwr0 ne 0 then if DTest eq 1 then bImproved:=true; else bImproved:=false; // init FTest:=Factorization(DTest); for kTest in [#FTest..1 by -1] do kBest:=kTest+#FBest-#FTest; if kBest eq 0 then break; elif FTest[kTest][1] lt FBest[kBest][1] then bImproved:=true; break; elif FTest[kTest][1] gt FBest[kBest][1] then break; else // i.e., FTest[kTest][1] = FBest[kBest][1] if FTest[kTest][2] lt FBest[kBest][2] then bImproved:=true; break; end if; end if; end for; end if; if bImproved then xSetBest:=Include(xSet0,xAdd1); xSetBest:=Include(xSetBest,xAdd2); FBest:=FTest; sumBest:=sumTest; if Numerator(sumBest) eq 1 then DBest:=Denominator(sumBest); if bIsInSeq[DBest] then bSuccess:=true; elif DBest mod 2 eq 0 then // denom is even if IsPrimePower(DBest) then // denom is pwr of 2 bSuccess:=true; end if; end if; if bSuccess then break; end if; end if; end if; end if; end if; end for; if bSuccess then break; end if; end if; end for; if bSuccess then break; end if; end if; if sumBest eq sum then // no improvement // try adding three reciprocals for m1 in [1..mMax-2] do xAdd1:=m1*PMaxPwr0; if not xAdd1 in xSet0 then for m2 in [m1+1..mMax-1] do xAdd2:=m2*PMaxPwr0; if not xAdd2 in xSet0 then for m3 in [m2+1..mMax] do xAdd3:=m3*PMaxPwr0; if not xAdd3 in xSet0 then sumTest:=sum+1/xAdd1+1/xAdd2+1/xAdd3; DTest:=Denominator(sumTest); if DTest mod PMaxPwr0 ne 0 then if DTest eq 1 then bImproved:=true; else bImproved:=false; // init FTest:=Factorization(DTest); for kTest in [#FTest..1 by -1] do kBest:=kTest+#FBest-#FTest; if kBest eq 0 then break; elif FTest[kTest][1] lt FBest[kBest][1] then bImproved:=true; break; elif FTest[kTest][1] gt FBest[kBest][1] then break; else // i.e., FTest[kTest][1] = FBest[kBest][1] if FTest[kTest][2] lt FBest[kBest][2] then bImproved:=true; break; end if; end if; end for; end if; if bImproved then xSetBest:=Include(xSet0,xAdd1); xSetBest:=Include(xSetBest,xAdd2); xSetBest:=Include(xSetBest,xAdd3); FBest:=FTest; sumBest:=sumTest; if Numerator(sumBest) eq 1 then DBest:=Denominator(sumBest); if bIsInSeq[DBest] then bSuccess:=true; elif DBest mod 2 eq 0 then // denom is even if IsPrimePower(DBest) then // denom is pwr of 2 bSuccess:=true; end if; end if; if bSuccess then break; end if; end if; end if; end if; end if; end for; if bSuccess then break; end if; end if; end for; if bSuccess then break; end if; end if; end for; if bSuccess then break; end if; end if; if sumBest eq sum then // no improvement break; // failure! end if; xSet:=xSetBest; sum:=sumBest; end while; end if; bIsInSeq[n]:=true; // (but this is proved only if no failures have occurred) if bSuccess then provedInCount+:=1; if bShowSolns then itos(n) * ":", sumBest, xSetBest; end if; termCount+:=1; if bOutputForbFile then termCount, n; end if; else // failure failureCount+:=1; n , "=", itos(n div PPmax) * "*" * itos(PPmax), s[i1stGoodSumInS], Intseq(i1stGoodSumInS-1, 2); if sum ne 1/n then "*** FAILED AFTER IMPROVEMENT ON AT LEAST 1 PASS ***"; sum, "=", xSet; break; end if; end if; end if; end if; end for; ""; "Number of integers n in [1.." * itos(nMax) * "] for which"; " it was proved that n is in A092671:", provedInCount; " it was proved that n is NOT in A092671:", provedOutCount; " the program failed to prove either result:", failureCount;