Download or view Cipolla.frink in plain text format or (experimental) run in browser
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. <http://eudml.org/doc/123970>.
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[]
Download or view Cipolla.frink in plain text format or (experimental) run in browser
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen, eliasen@mindspring.com