求 n 使得 n! 以 2021 开头

uk702

某书证明了存在 n，使得 n! 可以以任意一串数字作为开头。现求正整数 n，使得 n! 以2021 开头。

mathe

16979
17791
24764
26619
27727
34990
42004
57216
71869
74494
77981
81253
87457
98508

uk702

使用 Mathematica 硬算，竟然也好使：

1. n = 11; prod = 10!;
   
2. While[n < 100000, prod = prod*n;
   
3. If[Take[RealDigits[prod][[1]], 4] == {2, 0, 2, 1}, Print[n]];
   
4. n = n + 1
   
5. ]
   

复制代码

wayne

1. {16979,17791,24764,26619,27727,34990,42004,57216,71869,74494,77981,81253,87457,98508,101646,106752,115082,117536,123759,124963,125097,128930,129496,143176,145430,145608,147954,150084,159119,160231,161163,168079,170851,175620,179179,186731,194395,197814,200531,202408,205634,220996,229949,232031,232165,232359,233538,233812,236569,239429,241831,244379,256844,258008,260086,266522,274849,279710,282887,283345,283429,287093,297038,313125,320693,324548,327849,339662,343693,344403,347492,347553,348061,355546,357943,358607,359862,366076,369742,375529,384520,385526,411111,414986,427070,428915,442399,446228,446642,446722,452782,453888,456909,459849,463837,465062,465878,470988,474609,478127,480111,490862,496341,499823,502028,506166,507946,510457,517036,528199,532074,543024,551055,554680,555384,556545,562137,562545,573616,574616,591316,597692,603304,606644,608358,615043,616439,621847,624892,630419,631890,641452,654517,661103,661844,665506,667020,675154,676942,678469,679987,681658,682597,682784,687780,687940,689781,699900,708254,709378,710487,722016,727589,728845,732825,736823,743496,749697,754686,765342,766274,770955,777918,783395,785148,786380,793419,797545,797657,799510,807757,819293,843222,845303,847777,849418,858218,867231,872432,878030,888901,889824,894485,906501,908758,908999,913891,916410,916542,920553,922326,922440,930970,937076,946271,949542,953722,972282,972778,981548,982177,988996}

复制代码

wayne

无聊的跑了10分钟，刷出了10^8的所有结果，突破了我就给代码，

1. {16979,17791,24764,26619,27727,34990,42004,57216,71869,74494,77981,81253,87457,98508,101646,106752,115082,117536,123759,124963,125097,128930,129496,143176,145430,145608,147954,150084,159119,160231,161163,168079,170851,175620,179179,186731,194395,197814,200531,202408,205634,220996,229949,232031,232165,232359,233538,233812,236569,239429,241831,244379,256844,258008,260086,266522,274849,279710,282887,283345,283429,287093,297038,313125,320693,324548,327849,339662,343693,344403,347492,347553,348061,355546,357943,358607,359862,366076,369742,375529,384520,385526,411111,414986,427070,428915,442399,446228,446642,446722,452782,453888,456909,459849,463837,465062,465878,470988,474609,478127,480111,490862,496341,499823,502028,506166,507946,510457,517036,528199,532074,543024,551055,554680,555384,556545,562137,562545,573616,574616,591316,597692,603304,606644,608358,615043,616439,621847,624892,630419,631890,641452,654517,661103,661844,665506,667020,675154,676942,678469,679987,681658,682597,682784,687780,687940,689781,699900,708254,709378,710487,722016,727589,728845,732825,736823,743496,749697,754686,765342,766274,770955,777918,783395,785148,786380,793419,797545,797657,799510,807757,819293,843222,845303,847777,849418,858218,867231,872432,878030,888901,889824,894485,906501,908758,908999,913891,916410,916542,920553,922326,922440,930970,937076,946271,949542,953722,972282,972778,981548,982177,988996,1010616,1013325,1031953,1033245,1035162,1043571,1044638,1045835,1046346,1050392,1050579,1051136,1055948,1058351,1083839,1087892,1093823,1101817,1115838,1115922,1121485,1127931,1128313,1128808,1128865,1130864,1130995,1134176,1135142,1136426,1136498,1140001,1144434,1144792,1145302,1148876,1162926,1168305,1170702,1179688,1183182,1186279,1188371,1210638,1220499,1221272,1225837,1227651,1230031,1231528,1232411,1240899,1249064,1260279,1278026,1282555,1286628,1287048,1292226,1294330,1298675,1301238,1303339,1305928,1306066,1308378,1310036,1310760,1312754,1321407,1322257,1324025,1328826,1337609,1340299,1340739,1341547,1341876,1354060,1363894,1364539,1370147,1373860,1376015,1376923,1388870,1395426,1395944,1398237,1400512,1406591,1407508,1410128,1427201,1429410,1435955,1437645,1458822,1466927,1469535,1473386,1479761,1483737,1484227,1501671,1504186,1508340,1512522,1526343,1526490,1526958,1529786,1531090,1542711,1550132,1553888,1555125,1557504,1559678,1561127,1568230,1572536,1578871,1579118,1594172,1606070,1620370,1624220,1632212,1632865,1632973,1633372,1636627,1648284,1650580,1652558,1664419,1665196,1666279,1667230,1677012,1680541,1692051,1697981,1702549,1705120,1707107,1717127,1734340,1736366,1743935,1746321,1746544,1752082,1753436,1763439,1768286,1770465,1770731,1785561,1788888,1800624,1809649,1816780,1823025,1823669,1824033,1827842,1828487,1830275,1831284,1839497,1839565,1844469,1849187,1854251,1854512,1857855,1860380,1879112,1881730,1890020,1909653,1912142,1914704,1917052,1921280,1926200,1934752,1936969,1938257,1944328,1946615,1949719,1955414,1962540,1967418,1973681,1977543,1989783,1997106,2007688,2013520,2016028,2022836,2029244,2039315,2042981,2045200,2047602,2052106,2054436,2055312,2056443,2060675,2061385,2070363,2077371,2077582,2078432,2079530,2089561,2094605,2095539,2102805,2105603,2106858,2114324,2116681,2117338,2121657,2133299,2143013,2150166,2152198,2155736,2156671,2162049,2163980,2179251,2185073,2186089,2201412,2202780,2229301,2233627,2234452,2242992,2243505,2261649,2264581,2266034,2267978,2275084,2277925,2279172,2279390,2280407,2284450,2289805,2290333,2294524,2295447,2298739,2300235,2303436,2308617,2310724,2311642,2318442,2318672,2322239,2323015,2334488,2338181,2348019,2352029,2353528,2355808,2356117,2359824,2363294,2363551,2378075,2386796,2394091,2402438,2404452,2404906,2410315,2417816,2423827,2426421,2437562,2439662,2439915,2440888,2447054,2450029,2451981,2455016,2455521,2457842,2458825,2468282,2470857,2470969,2472338,2486814,2492161,2493119,2494064,2495978,2496197,2517703,2519281,2519737,2525555,2532127,2532915,2536856,2537363,2540257,2542200,2544849,2548689,2553196,2555450,2561372,2566449,2571569,2575941,2579501,2581534,2582221,2583017,2584841,2586259,2587991,2589337,2592445,2594320,2595129,2595636,2602037,2605925,2609214,2615348,2615985,2618706,2622296,2626768,2631056,2633554,2636565,2636736,2639748,2641031,2641645,2647233,2656613,2660115,2660468,2667234,2669709,2672266,2672875,2673917,2680137,2690638,2691780,2694713,2697485,2699918,2703333,2704583,2705534,2709762,2727072,2730210,2732978,2739847,2741021,2744077,2745557,2749175,2755059,2758583,2760205,2762302,2766437,2769341,2772732,2778443,2790100,2790735,2794672,2797375,2804194,2804424,2813057,2814991,2825757,2832364,2835388,2836265,2839627,2862173,2870980,2873385,2881496,2884597,2889707,2895794,2896803,2902513,2910235,2911653,2912603,2912631,2912659,2916669,2929883,2936621,2937478,2939533,2940560,2944460,2952452,2958616,2971266,2973684,2986650,2989515,2994444,2999810,3001566,3007002,3024162,3033342,3033425,3035518,3039087,3046809,3055398,3056171,3058215,3058767,3068618,3069503,3074226,3078275,3085472,3088628,3092544,3095337,3099907,3101785,3106436,3107086,3118387,3126771,3143066,3147316,3148063,3149125,3152632,3154749,3159880,3172660,3176254,3179098,3179737,3184022,3185485,3192103,3201348,3209581,3211853,3220441,3220756,3222061,3229090,3229255,3229911,3230074,3232736,3235604,3239614,3241915,3244825,3256753,3263989,3271431,3274905,3280493,3283804,3284048,3284989,3286013,3294796,3297572,3299374,3302380,3306258,3309986,3318593,3321993,3322739,3323389,3326055,3327327,3333243,3333396,3335826,3339389,3353935,3358385,3367663,3385637,3388687,3389053,3396289,3400874,3416494,3416747,3420660,3422588,3422908,3427163,3438536,3444751,3447630,3449609,3450323,3469782,3470476,3472463,3479366,3481755,3485443,3487425,3495068,3495436,3497967,3499561,3505621,3513652,3519517。。。。。。。。。。。,99352598,99354373,99360445,99369117,99370574,99371669,99372400,99379388,99379758,99383470,99384588,99385708,99393985,99395882,99397023,99399694,99405072,99405458,99409332,99415584,99415977,99421110,99423893,99427491,99435569,99453762,99455868,99458406,99460956,99464375,99466523,99472584,99476079,99480481,99482696,99486260,99489399,99489849,99491653,99492105,99493917,99500771,99501692,99504002,99513352,99525299,99526753,99531142,99534091,99539048,99542550,99546079,99548108,99549126,99552707,99553221,99556317,99565738,99576452,99581363,99591925,99595886,99596455,99598167,99600461,99601613,99603927,99611543,99612728,99615109,99620521,99629717,99637876,99648832,99651460,99657447,99659466,99662858,99676072,99683975,99686165,99689109,99698876,99700407,99704269,99705828,99711351,99716982,99723558,99724391,99732005,99741616,99744300,99745201,99746105,99747923,99748837,99750675,99764918,99774914,99781131,99783243,99787530,99790803,99794128,99814091,99817841,99821670,99824269,99832308,99842204,99849671,99852764,99854335,99860797,99862460,99864143,99867574,99869323,99871096,99884265,99886271,99903812,99908723,99916630,99919438,99922348,99925371,99931818,99935282,99938942,99942836,99947016,99951555,99956566,99956567,99962237,99968926,99977524,99977525,99977526,99993316,99993317,99993318,99993319,99993320,99993321,99993322}

复制代码

wayne

单核版本的代码：

1. n=1;Monitor[Reap[While[n<10^5,t=LogGamma[n+1.];If[Log[2021]<t-(Floor[t/Log[10]]-3)Log[10]<Log[2022],Print[Sow[n]]];n++]],n]

复制代码

多核版本的代码：

1. CloseKernels[];LaunchKernels[12];
   
2. Timing[ans=ParallelCombine[Select[#,Log[2021]<Last[#]-(Floor[Last[#]/Log[10]]-3)Log[10]<Log[2022]&]&,Table[{n,LogGamma[n+1.]},{n,10^8}]][[All,1]];]

复制代码

uk702

wayne 发表于 2021-1-27 21:46
单核版本的代码：

多核版本的代码：

我将你的多核版本的结果打印出来，似乎各有各的不准。
31419078!=2022000269906979489591697...，你的多了这个。
39600422!=2021000126106549851567511...，你的少了这个。
48400100!=2020999793431651954566836...，我的多了这个。
48844860!=2021000416575056021475393...，你的少了这个。
80495570!=2020999980052118200212907...，我的多了这个。

我实际给出的解是 21439 个（显示时多算了一个），多了 2 个，你给出了 21436 个，少了 2 个并多了 1 个。故检查后共同的解是 21437 个，暂不排除其中有我们都算错的。

wayne

这是有可能的，因为在做比较的时候用的是double精度,我稍微做了修改，用30位高精度计算，重新跑了下程序，速度慢了不少，算出正确的结果应该有21439 个。

1. CloseKernels[];LaunchKernels[12];
   
2. Timing[ans=ParallelCombine[Select[#,Log[2021]<Last[#]-(Floor[Last[#]/Log[10]]-3)Log[10]<Log[2022]&]&,Table[{n,LogGamma[n+1`30]},{n,10^8}]][[All,1]];]

复制代码

发现前面的double精度的版本产生了错误解31419078，60459888，还少了39600422，48844860，61543437，88289946，90711389。

uk702

我的单核版本，julia 代码，double 精度，使用的某位同学论文中的斯特林加强公式（https://wenku.baidu.com/view/77419c44336c1eb91a375d3d.html），直接计算 log(n) 的累加和在 double 精度下累积误差似乎相当大。

1. t1=time();
   
2. h(n) = (log(sqrt(2.0*pi*n)) + n*log(n)) - n * (1-1/(0.4+12*n^2));
   
3. n=10; ct=1; u = log(2021.); v = log(2022.)
   
4. while n < 10^8
   
5.   global n, ct
   
6.   a = (h(n) - u)/log(10.); t = floor(Int, a);
   
7.   if h(n) < v + t*log(10) ct = ct+1;println(n) end
   
8.   n=n+1
   
9. end
   
10. 
    
11. # 打印出来的 ct 比实际多 1 个
    
12. println("find $ct sols, use time = $(time() - t1)")
    

复制代码

mathe

斯泰林公式可以采用
\(\log(n!) \approx \frac{\log(2\pi n)}2+n\log(n)-n+\frac1{12n}+O(\frac{1}{n^3})\)