求 n 使得 n! 以 2021 开头

wayne

mathe的这个公式是LogGamma[n+1]函数的在无穷远处的无限逼近
$LogGamma[n+1] =\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+1] =-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

我的高精度多核心版本，在 95s 完成 10^8 计算。
find 21440 sols, use time = 95.08599996566772

1. 
   
2. # set JULIA_NUM_THREADS=8
   
3. 
   
4. using Base.Threads
   
5. 
   
6. t1=time();
   
7. n=10^8; seg = 10000; ns=Int64(n/seg); u = log(big(2021.)); v = log(big(2022.)); w = log(big(10.0));
   
8. 
   
9. mylock = ReentrantLock()
   
10. result = Channel{Int64}(ns)
    
11. sols = Int64[]
    
12. 
    
13. function check(g)
    
14.   global sols, mylock
    
15.   g = g*seg
    
16.   for n=g+2:g+seg+1
    
17.     h = log(sqrt(2.0*pi*n)) + big(n)*log(n) - big(n)*(1-1/(0.4+12*n^2));
    
18.     #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))
    
19.     a = (h - u)/w; t = floor(Int, a);
    
20. 
    
21.     if h < v + t*w
    
22.       lock(mylock)
    
23.       sols=append!(sols, n)
    
24.       l = length(sols)
    
25.       if l % 1000 == 0 println("find $l sols ...") end
    
26.       unlock(mylock)
    
27.     end
    
28.   end
    
29.   put!(result, g)
    
30. end
    
31. 
    
32. @threads for i=0:ns-1
    
33.     check(i)
    
34. end
    
35. 
    
36. ct = 0; t = ns;
    
37. while t > 1
    
38.    take!(result)
    
39.    global t = t - 1
    
40. end
    
41. 
    
42. for x in sols
    
43.   println(x)
    
44.   global ct += 1
    
45. end
    
46. 
    
47. println("find $ct sols, use time = $(time() - t1)")
    

复制代码

uk702

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

1. 
   
2. # set JULIA_NUM_THREADS=8
   
3. using Base.Threads
   
4. 
   
5. t1=time();
   
6. 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))
   
7. 
   
8. mylock = ReentrantLock()
   
9. result = Channel{Int64}(ns)
   
10. sols = Int64[]
    
11. 
    
12. function check(g)
    
13.   global sols, mylock
    
14.   g = g*seg
    
15.   for n=g+2:g+seg+1
    
16.     h = c1 + log(n)*(n+0.5) - n;
    
17.     a = (h - u)/w; t = floor(Int, a);
    
18. 
    
19.     eps = v + t*w - h
    
20.     if eps > 0
    
21.       if eps > 0.001
    
22.         lock(mylock)
    
23.         sols=append!(sols, n)
    
24.         unlock(mylock)
    
25.       else
    
26.         h = c1 + log(n)*(big(n)+0.5) - big(n); a = (h - log(big(2021.)))/log(big(10.0)); t = floor(Int, a);
    
27.         if h < v + t*log(big(10.0))
    
28.           lock(mylock)
    
29.           sols=append!(sols, n)
    
30.           unlock(mylock)
    
31.         end
    
32.       end
    
33.       
    
34.       l = length(sols)
    
35.       if l % 1000 == 0 println("find $l sols ...") end
    
36.     end
    
37.   end
    
38.   put!(result, g)
    
39. end
    
40. 
    
41. @threads for i=0:ns-1
    
42.     check(i)
    
43. end
    
44. 
    
45. ct = 0; t = ns;
    
46. while t > 1
    
47.    take!(result)
    
48.    global t = t - 1
    
49. end
    
50. 
    
51. for x in sols
    
52.   println(x)
    
53.   global ct += 1
    
54. end
    
55. 
    
56. println("find $ct sols, use time = $(time() - t1)")
    

复制代码

无心人

我觉得某个帖子里我写过类似算法

https://bbs.emath.ac.cn/thread-969-1-1.html

uk702

uk702 发表于 2021-1-28 16:47
采用两个精度，先是 double 精度进行粗算，再对临界值进行精算，将用时降到 10s 左右。
find 21440 so ...

fix 多出一个解的 bug，16 核并行，耗时 10.6s。
find 21439 sols, use time = 10.594000101089478

1. # set JULIA_NUM_THREADS=16
   
2. using Base.Threads
   
3. 
   
4. t1=time();
   
5. 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))
   
6. 
   
7. mylock = ReentrantLock()
   
8. result = Channel{Int64}(ns)
   
9. sols = Int64[]
   
10. 
    
11. function check(g)
    
12.   global sols, mylock
    
13.   g = g*seg
    
14.   for n=g+2:g+seg+1
    
15.     h = c1 + log(n)*(n+0.5) - n;
    
16.     a = (h - u)/w; t = floor(Int, a);
    
17. 
    
18.     eps = v + t*w - h
    
19.     if eps > 0
    
20.       if eps > 0.001
    
21.         lock(mylock)
    
22.         sols=append!(sols, n)
    
23.         unlock(mylock)
    
24.       else
    
25.         h = c1 + log(big(n))*(big(n)+0.5) - big(n)*(big(1)-big(1)/(big(0.4)+big(12)*n^2));
    
26.         a = (h - log(big(2021.)))/log(big(10.0)); t = floor(Int64, a);
    
27.         if h < log(big(2022.)) + t*log(big(10.0))
    
28.           lock(mylock)
    
29.           sols=append!(sols, n)
    
30.           unlock(mylock)
    
31.         end
    
32.       end
    
33.       
    
34.       l = length(sols)
    
35.       if l % 1000 == 0 println("find $l sols ...") end
    
36.     end
    
37.   end
    
38.   put!(result, g)
    
39. end
    
40. 
    
41. @threads for i=0:ns-1
    
42.     check(i)
    
43. end
    
44. 
    
45. ct = 0; t = ns;
    
46. while t > 1
    
47.    take!(result)
    
48.    global t = t - 1
    
49. end
    
50. 
    
51. for x in sols
    
52.   println(x)
    
53.   global ct += 1
    
54. end
    
55. 
    
56. println("find $ct sols, use time = $(time() - t1)")
    

复制代码

wayne

这个问题其实是 一个方程求根的问题：定义$f(x)= x-\lfloor \frac{x}{\log (10)}\rfloor \log (10) $, 求方程 $\log (\frac{2021}{1000} ) <f(log \Gamma(n+1)) <\log (\frac{2022}{1000} )$的正整数根。