// 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;