use functionUtils.frink /** This gives an asymptotic estimate for nthPrime[n]. Asymptotic means that it gets better as n gets larger, but you should be warned that *if you use too many terms of the series, it gets worse!* It typically guesses the first half of the digits of nthPrime correctly if the number of terms is chosen correctly. This gives a very good first guess to a Newton's method solver (which essentially doubles the number of correct digits with each iteration.) It uses a 1902 paper by Michele Cipolla: M. Cipolla: La determinazione assintotica dell nimo numero primo, Rend. Acad. Sci. Fis. Mat. Napoli 8(3) (1902) pp. 132-166. See, for a newer (but possibly broken) version Mincu, Gabriel. "An asymptotic expansion.." JIPAM. Journal of Inequalities in Pure & Applied Mathematics [electronic only] 4.2 (2003): Paper No. 30, 7 p., electronic only-Paper No. 30, 7 p., electronic only. . or https://web.archive.org/web/20240901000000/https://www.emis.de/journals/JIPAM/article268.html However NOTE THAT THE ABOVE PAPER SEEMS TO BE MISSING FACTORS of / j! and may have other errors in the recurrence relation. See corrections below: https://mathoverflow.net/questions/279520/the-nth-prime-in-terms-of-lambert-w https://arxiv.org/abs/1011.1667 The coefficients can be found (except for sign which alternate) at https://oeis.org/A200265 We need to figure out how many terms to use for arguments of each size. (This is now done via the numTerms function and automatically applied by estimateNthPrime.) */ class Cipolla { // This is https://oeis.org/A200265 if you want more terms. // Array of functions for successive terms of approximating function. // These are all added together. class var P = [{|x| 1}, {|x| x - 2}, {|x| x^2 - 6x + 11}, {|x| 2x^3 - 21x^2 + 84x - 131}, {|x| 6x^4 - 92x^3 + 588x^2 - 1908x + 2666}, {|x| 24x^5 - 490x^4 + 4380x^3 - 22020x^2 + 62860x - 81534}, {|x| 120x^6 - 3084x^5 + 35790x^4 - 246480x^3 + 1075020x^2 - 2823180x + 3478014}, {|x| 720x^7 - 22428x^6 + 322224x^5 - 2838570x^4 + 16775640x^3 - 66811920x^2 + 165838848x - 196993194}, {|x| 5040x^8 - 185184x^7 + 3186848x^6 - 34369776x^5 + 257567520x^4 - 1381360960x^3 + 5177983104x^2 - 12358329648x + 14297456816}, {|x| 40320x^9 - 1712016x^8 + 34445664x^7 - 441118944x^6 + 4018300272x^5 - 27182604960x^4 + 136703843136x^3 - 491501715264x^2 + 1140641857584x - 1293975782352}] /** Estimate the nth prime. */ class estimateNthPrime[n] := { terms = numTerms[n] return estimateNthPrime[n, terms] } /** Estimate the nth prime. */ class estimateNthPrime[n, terms] := { if n < 1 { println["estimateNthPrime: invalid argument $n"] return undef } if n <= 1 return 2 if n <= 2 return 3 // TODO: Better estimate for small n. Maybe lookup in sieved primes or // return nth1[primes[], n] or put this into a local cache. // UPDATE: Now we can do something like: // if (n <= Frink.function.Factor.getHighestIndexInSieve()) // return Frink.function.Factor.getOneBasedNthPrimeFromSieve(n); ln = ln[n] lln = ln[ln] part = 0 for t = terms to 1 step -1 part = part + (-1)^(t-1) apply[P@t, lln] / (t! ln^t) return round[n (ln + lln - 1 + part)] } /** This is a private method that returns the number of terms of P to optimally calculate the nth prime number n. These numbers were found empirically from the test[] function. */ class numTerms[n] := { if n > 3ee22 return 9 if n > 3ee20 return 8 if n > 3ee18 return 7 if n > 3ee16 return 6 if n > 6ee13 return 5 if n > 6ee10 return 4 if n > 6ee8 return 3 if n > 2ee6 return 2 if n > 2ee3 return 1 return 0 } class test[n, correct, terms] := { nth = estimateNthPrime[n, terms] if nth <= 1 nth = 2 err = (correct-nth)/nth // println["Value is " + round[nth]] // println["should be " + correct + ", err is " + formatSci[err, 1, 3]] return [nth, err] } class test[n, correct] := { minerr = 1e100 bestTerms = 0 bestVal = 0 for terms = 0 to length[P]-1 { [nth, err] = test[n, correct, terms] if abs[err] < abs[minerr] { minerr = err bestVal = nth bestTerms = terms } } println["$n\t" + formatSci[n,1,1]] println["best value is " + round[bestVal] + ", err is " + formatSci[minerr, 1, 3] + " with $bestTerms"] println["should be " + correct + "\n"] return [minerr, bestTerms] } class test[] := { Cipolla.test[1ee0, 2] Cipolla.test[1ee1, 29] Cipolla.test[3ee1, 113] Cipolla.test[1ee2, 541] Cipolla.test[3ee2, 1987] Cipolla.test[1ee3, 7919] Cipolla.test[3ee3, 27449] Cipolla.test[1ee4, 104729] Cipolla.test[3ee4, 350377] Cipolla.test[1ee5, 1299709] Cipolla.test[3ee5, 4256233] Cipolla.test[1ee6, 15485863] Cipolla.test[3ee6, 49979687] Cipolla.test[1ee7, 179424673] Cipolla.test[3ee7, 573259391] Cipolla.test[1ee8, 2038074743] Cipolla.test[3ee8, 6461335109] Cipolla.test[1ee9, 22801763489] Cipolla.test[3ee9, 71856445751] Cipolla.test[1ee10, 252097800623] Cipolla.test[3ee10, 790645490053] Cipolla.test[1ee11, 2760727302517] Cipolla.test[3ee11, 8624419641811] Cipolla.test[1ee12, 29996224275833] Cipolla.test[3ee12, 93400375993241] Cipolla.test[1ee13, 323780508946331] Cipolla.test[3ee13, 1005368767096627] Cipolla.test[1ee14, 3475385758524527] Cipolla.test[3ee14, 10765662794071351] Cipolla.test[1ee15, 37124508045065437] Cipolla.test[1ee16, 394906913903735329] Cipolla.test[1ee17, 4185296581467695669] Cipolla.test[1ee18, 44211790234832169331] Cipolla.test[1ee19, 465675465116607065549] Cipolla.test[1ee20, 4892055594575155744537] Cipolla.test[1ee21, 51271091498016403471853] Cipolla.test[1ee22, 536193870744162118627429] Cipolla.test[1ee23, 5596564467986980643073683] Cipolla.test[1ee24, 58310039994836584070534263] } class exportBodies[] := { inits = new set exprs = new array var neg var num for f = Cipolla.P { body = inputForm[exportBody[f]] orig = inputForm[functionBody[f]] for [sign, num] = body =~ %r/(\-?)(\d+)/g { neg = "" if sign == "-" neg = "_NEG" inits.put[" private static final FrinkInteger N${neg}_$num = FrinkInteger.construct(\"$sign\$num\", 10);"]; } body =~ %s/(\-?)(\d+)/N=$2;$1=="-" ? "N_NEG_$N" : "N_$N"/eg println[body] body =~ %s/\[/\(/g body =~ %s/\]/\)/g bstr = """ // $orig new CipollaEquation() { public Numeric apply(FrinkInteger x, MathContext mc) throws NumericException { return $body; } },\n""" exprs.push[bstr] } println[joinln[exprs]] println[joinln[lexicalSort[toArray[inits]]]] } class exportBody[func] := { body = functionBody[func] body = transformExpression[body] // body = substituteExpression[body, noEval[_a ^ _b], noEval[NumericMath.power[_a, _b, mc]]] // body = substituteExpression[body, noEval[_a * _b], noEval[NumericMath.multiply[_a, _b, mc]]] // body = substituteExpression[body, noEval[_a + _b], noEval[NumericMath.add[_a, _b, mc]]] return body } } transformations FI { _a + _b <-> NumericMath.add[_a, _b, mc] _a * _b <-> NumericMath.multiply[_a, _b, mc] _a ^ _b <-> NumericMath.power[_a, _b, mc] } Cipolla.test[] Cipolla.exportBodies[]