(* "Copyright 1999, Allan Swett. May be freely used, quoted, modified, and distributed, provided that reference to authorship, within these quotes, is included." A Mathematica program. This is Phase 2 of a software project to confirm that the "Erdos-Strauss Conjecture" holds to beyond 100 trillion ( = 10^14 ). The program is set to read a sequence of data files from a user - specified directory. The directory and data file numbers, in the function call at the end of the program, should be modified by end - users as required. Version date : 10/19/99. *) (* Part 1. A "toolbox" of procedures, to confirm ESC(n) for a particular prime n. *) setglobals[] := Module[{}, printlevel = 0; globalerr = 0; globaldepth = 0; biggestdepth = 0; bigone = 0; ]; setglobals[]; (* "keypair" returns a pair of divisors of the "denominator", dd*pp, whose sum is divisible by the "numerator", nn. If no such pair of divisors is found, keypair returns the empty list, {}. Usage : pp prime, and dd < pp. The algorithm is designed to avoid generating divisors of large integers, such as dd*pp. It is written under the assumption that pp is prime and relatively prime to dd, and is conclusive under these circumstances ... that is, if such a pair exists, it is returned. *) keypair[nn_, pp_, dd_] := Module[ {divs, divisorcount, modlist, ii, current, mm}, (* Return values quickly for easy or troublesome situations, where nn divides pp*dd + 1 : *) If[nn == 1, Return[{1, pp*dd}]]; If[Mod[pp*dd, nn] == (nn - 1), Return[{1, pp*dd}]]; (* Not so simple ... need to look through the list of divisors of pp*dd. Start with the divisors of dd : *) divs = Divisors[dd]; (* divs is a list of the divisors of dd. Since pp is prime and relatively prime to dd, this is half of the list of divisors of pp*dd. *) divisorcount = Length[divs]; (* Process the list of divisors, looking for a pair whose sum is congruent to 0 mod nn. "modlist" is used to mark the least residues, mod nn, of divisors already encountered; the marks also serve to retain the value of a divisor generating the least residue. The mark is 0 if not yet generated, and equals the most recent generating divisor if already produced. *) modlist = Table[0, {nn}]; ii = 1; While[ii <= divisorcount, current = divs[[ii]]; mm = Mod[current, nn]; modlist[[mm + 1]] = current; If[modlist[[nn - mm + 1]] != 0, Return[{modlist[[nn - mm + 1]], current}] ]; ii++; ]; (* Examine the remaining divisors ... divisors of dd*pp divisible by pp : *) ii = 1; While [ii <= divisorcount, current = divs[[ii]]*pp; mm = Mod[current, nn]; modlist[[mm + 1]] = current; If[modlist[[nn - mm + 1]] != 0, Return[{modlist[[nn - mm + 1]], current}] ]; ii++; ]; Return[{}]; ]; (* The module "trio" is the principal vehicle for numerically checking ESC. It produces three positive integers whose reciprocals have sum 4/givend, using a limited iteration based upon a "greedy" algorithm. For very large values of givend, a longer iteration may be needed ... the current iteration length works for all primes less than one billion. "trio" returns a triple of positive integers whose reciprocals have the claimed sum, or the empty list if it has failed. Usage : "givend" an odd prime. *) trio[givend_] := Module[ {rat, smd, ii, theden, nn, goodpair, f1, f2, tt, mm}, rat = 4/givend; smd = Floor[1/rat] + 1; (* "smd" is the smallest possible denominator *) ii = 0; While[ii < 30, (* 30 appears sufficient, for now, as the iteration depth. 4 does not suffice for 4/1201 nor 4/2521; 10 does not suffice for 4/118801; 20 does not suffice for 4/8803369, which requires 26 *) theden = smd + ii; nn = 4*theden - givend; goodpair = keypair[nn, givend, theden]; If[goodpair == {},(* then nothing *) , (* else *) {f1, f2} = goodpair; tt = nn/(givend*theden); (* target value *) mm = ((1/f1) + (1/f2))/tt; (* mm is the adjusting multiplier ... it is a positive integer, although this may not be immediately clear *) globaldepth = ii; If[globaldepth > biggestdepth, biggestdepth = globaldepth; bigone = givend]; Return[{theden, mm*f1, mm*f2}]; ]; (* endif *) ii++; ]; (* endwhile *) Return[{}]; ]; (* "toolbox", procedure "get3". This is the main procedure for confirming ESC(k) for k an odd prime. *) get3[givend_] := Module[ {ttt, f1, f2, f3, checkval}, ttt = trio[givend]; If[ttt == {}, Print[" "]; Print["No luck with ", givend]; globalerr++; Return[]; ]; (* play "dumb" ... make absolutely sure that we are dealing with three unit fractions ... *) {f1, f2, f3} = 1/Round[Abs[ttt]]; If[printlevel > 0, Print[" "]; Print[4/givend, " = ", f1, " + ", f2, " + ", f3, " (", globaldepth, ")"]; ]; (* Truly confirm ... check for possible mistake ... *) checkval = ((4/givend) == (f1 + f2 + f3)); If[checkval === True, , Print["Error!!! with ", givend]; globalerr++; ]; ]; (* Part 2. Procedures to check a list of exceptional cases. *) (* These procedures read and process a sequence of exception files. additional global variables ecases counts the cases computed primecases counts the prime cases rfile the file used for data input primeexceptions an optional list of prime cases uses the existing global variables globalerr, biggestdepth, bigone resets all toolbox globals with a call to setglobals *) (* Two procedures to open and close input files : *) startread[usedir_, nnn_] := Module[ {fname, flist}, SetDirectory[usedir]; fname = ToString[nnn] <> ".txt"; flist = FileNames[fname]; If[flist == {}, Print["ERROR: ", nnn, " th file missing!"]; globalerr++; ResetDirectory[]; Return[False]; , ClearAll[rfile]; rfile = OpenRead[fname]; ResetDirectory[]; Return[True]; ]; ]; endread[] := Module[{}, Close[rfile]; ]; (* Two procedures to comfirm ESC : confirm checks the primes in an input file (2 cannot be there!) confirmbatch calls confirm, checking a set of input files. *) confirm[usedir_, nnn_] := Module[ {ggg}, Print[nnn, "... "]; If[startread[usedir, nnn], ggg = Read[rfile, Number]; While[ggg != 0, ecases++; If [PrimeQ[ggg], (* AppendTo[primeexceptions, ggg]; *) get3[ggg]; primecases++; ]; ggg = Read[rfile, Number]; ]; endread[]; ]; (* endif *) ]; confirmbatch[usedir_, low_, high_] := Module[ {t1, t2, iii}, setglobals[]; ecases = 0; primecases = 0; (* primeexceptions = {} ; *) t1 = AbsoluteTime[]; Do[confirm[usedir, iii], {iii, low, high}]; t2 = AbsoluteTime[]; Print["Exceptional cases: ", ecases, ". ", primecases, " are prime."]; Print["Greatest depth = ", biggestdepth, ", at ", bigone, "."]; Print["Approximate time required: ", t2 - t1, " seconds."]; Print["Total errors: ", globalerr, "."]; ]; confirmbatch[ "C:\ESC Data", 1, 100];