$LogGamma =\frac{\log(2\pi n)}2+n\log(n)-n+\frac1{12n}-\frac{1}{360 n^3}+\frac{1}{1260 n^5}-\frac{1}{1680 n^7}+\frac{1}{1188 n^9}+O\left(\left(\frac{1}{n}\right)^{11}\right)$
要能同级别pk，只能搬出拉马努金了
$LogGamma =-n+n \log (n)+\frac{1}{6} \log (n (4 n (2 n+1)+1))+\frac{\log (\pi )}{2} +\frac{1}{1440 n^3}-\frac{1}{768 n^4}+\frac{17}{16128 n^5}-\frac{19}{30720 n^7}+\frac{1}{98304 n^8}+\frac{4085}{4866048 n^9}+O\left(\left(\frac{1}{n}\right)^{11}\right)$
哦，其实Wikipedia上更丰富。
https://en.wikipedia.org/wiki/Stirling%27s_approximation 本帖最后由 uk702 于 2021-1-28 11:03 编辑
我的高精度多核心版本，在 95s 完成 10^8 计算。
find 21440 sols, use time = 95.08599996566772
# set JULIA_NUM_THREADS=8
using Base.Threads
t1=time();
n=10^8; seg = 10000; ns=Int64(n/seg); u = log(big(2021.)); v = log(big(2022.)); w = log(big(10.0));
mylock = ReentrantLock()
result = Channel{Int64}(ns)
sols = Int64[]
function check(g)
global sols, mylock
g = g*seg
for n=g+2:g+seg+1
h = log(sqrt(2.0*pi*n)) + big(n)*log(n) - big(n)*(1-1/(0.4+12*n^2));
#h = log(sqrt(big(2.)*pi*n)) + log(n)*n - n + log(big(1.0)+1/(12*n)+1/(288*n^2)-139/(51840*n^3)-571/(2488320*n^4))
a = (h - u)/w; t = floor(Int, a);
if h < v + t*w
lock(mylock)
sols=append!(sols, n)
l = length(sols)
if l % 1000 == 0 println("find $l sols ...") end
unlock(mylock)
end
end
put!(result, g)
end
@threads for i=0:ns-1
check(i)
end
ct = 0; t = ns;
while t > 1
take!(result)
global t = t - 1
end
for x in sols
println(x)
global ct += 1
end
println("find $ct sols, use time = $(time() - t1)")
本帖最后由 uk702 于 2021-1-28 16:50 编辑
uk702 发表于 2021-1-28 11:01
我的高精度多核心版本，在 95s 完成 10^8 计算。
find 21440 sols, use time = 95.08599996566772
采用两个精度，先是 double 精度进行粗算，再对临界值进行精算，将用时降到 10s 左右。
find 21440 sols, use time = 10.601999998092651
# set JULIA_NUM_THREADS=8
using Base.Threads
t1=time();
n=10^8; seg = 10000; ns=Int64(n/seg); u = log(2021.); v = log(2022.); w = log(10.0); c1=log(sqrt(2.0*pi))
mylock = ReentrantLock()
result = Channel{Int64}(ns)
sols = Int64[]
function check(g)
global sols, mylock
g = g*seg
for n=g+2:g+seg+1
h = c1 + log(n)*(n+0.5) - n;
a = (h - u)/w; t = floor(Int, a);
eps = v + t*w - h
if eps > 0
if eps > 0.001
lock(mylock)
sols=append!(sols, n)
unlock(mylock)
else
h = c1 + log(n)*(big(n)+0.5) - big(n); a = (h - log(big(2021.)))/log(big(10.0)); t = floor(Int, a);
if h < v + t*log(big(10.0))
lock(mylock)
sols=append!(sols, n)
unlock(mylock)
end
end
l = length(sols)
if l % 1000 == 0 println("find $l sols ...") end
end
end
put!(result, g)
end
@threads for i=0:ns-1
check(i)
end
ct = 0; t = ns;
while t > 1
take!(result)
global t = t - 1
end
for x in sols
println(x)
global ct += 1
end
println("find $ct sols, use time = $(time() - t1)")
我觉得某个帖子里我写过类似算法
https://bbs.emath.ac.cn/thread-969-1-1.html uk702 发表于 2021-1-28 16:47
采用两个精度，先是 double 精度进行粗算，再对临界值进行精算，将用时降到 10s 左右。
find 21440 so ...
fix 多出一个解的 bug，16 核并行，耗时 10.6s。
find 21439 sols, use time = 10.594000101089478
# set JULIA_NUM_THREADS=16
using Base.Threads
t1=time();
n=10^8; seg = 10000; ns=Int64(n/seg); u = log(2021.); v = log(2022.); w = log(10.0); c1=log(sqrt(2.0*pi))
mylock = ReentrantLock()
result = Channel{Int64}(ns)
sols = Int64[]
function check(g)
global sols, mylock
g = g*seg
for n=g+2:g+seg+1
h = c1 + log(n)*(n+0.5) - n;
a = (h - u)/w; t = floor(Int, a);
eps = v + t*w - h
if eps > 0
if eps > 0.001
lock(mylock)
sols=append!(sols, n)
unlock(mylock)
else
h = c1 + log(big(n))*(big(n)+0.5) - big(n)*(big(1)-big(1)/(big(0.4)+big(12)*n^2));
a = (h - log(big(2021.)))/log(big(10.0)); t = floor(Int64, a);
if h < log(big(2022.)) + t*log(big(10.0))
lock(mylock)
sols=append!(sols, n)
unlock(mylock)
end
end
l = length(sols)
if l % 1000 == 0 println("find $l sols ...") end
end
end
put!(result, g)
end
@threads for i=0:ns-1
check(i)
end
ct = 0; t = ns;
while t > 1
take!(result)
global t = t - 1
end
for x in sols
println(x)
global ct += 1
end
println("find $ct sols, use time = $(time() - t1)")
这个问题其实是 一个方程求根的问题：定义$f(x)= x-\lfloor \frac{x}{\log (10)}\rfloor \log (10) $, 求方程 $\log (\frac{2021}{1000} ) <f(log \Gamma(n+1)) <\log (\frac{2022}{1000} )$的正整数根。
页:
1
[2]