有沒有自然數 n，使得 n! 的前四個數字是 2020？
有，其中最小的自然數為 n=384。
384!
=
20205002129876021521366889509458521857103131384424164703865633675390
18051160495716743506601574409795917922928563661299169269617068846773
01229576788637571823309033713973447373497759772077557498674190878950
85555858041382225358897077307262143988184652395676694547987597481082
42045386512721839333046204349301762565136352275077524224647524747714
86396742989664230239433214649824117145479984513359446662469851311138
75121097651996768775622775133017838802100903333902268129172961937277
25084255908023383476149641735328778812131475694211759651482790988260
72666361443374874234852854408882112301367626134918709620006585891596
43594934505150789217346436153635132813435806722718188852357091873258
00311071868796569251630544786139037035591559504566681600000000000000
00000000000000000000000000000000000000000000000000000000000000000000
000000000000
其實只要把滿足 的所有
找到就行了，所以不需要大整數運算即可在很快的時間內得出
之內滿足要求的數恰好有
個。
另外，我們可以把 的十進位表示看成一個各位均為隨機的數（但是這一點未經證明），所以根據著名的本福特定律，滿足條件的數在全體自然數中所佔比例並不是
，而是
，可以看到逼近於此比例的速度還是比較慢的。
UPD:
其實壓根就不需要用 去算，這樣不僅引入精度誤差，而且慢。只需要在算階乘的過程之中，如果大於10000了馬上除個10，就可以保證一直在精度範圍之內了。如果說您非要使用
去算，需要注意一些數值穩定性的小trick，然後精度開double就夠了。
在我的電腦上，使用 計算的代碼需要8.6s左右才能跑出
以內的答案，只需要浮點數乘除法的代碼只需要0.33秒。（Ubuntu 20.04.2, g++ 10.2.0, 編譯開關-Ofast -march=native）
並且有人提到浮點誤差的問題，我用GMP庫的高精度驗算了一下，保留200位小數的意義下也還是這個答案，可以基本排除浮點誤差的可能。
大家都喜歡貼代碼，那麼我也貼一個C++代碼吧：
#include&
long double x=1,b=1;
int ans=0;
int main(){
for(int i=1,a=10;i&<=1e8;i++){
if(i==a)a=i*10,b/=10;
x*=i*b;
if(x&>1e4)x/=10;
if(x&>=2020 x&<2021)ans++;
}
printf("%d
",ans);
return 0;
}
先給個不那麼嚴謹的證明
題目等價於存在一對正整數 滿足：
簡單變形一下：
令 ，
，上式就可以寫成：
故只要存在正整數 滿足
就可以了（這裡
表示
的小數部分）
取一個足夠大的正整數 ，使得區間
內的正整數個數大於
，對於其中任一正整數
,容易知道：
假設該區間內最小的正整數為 ，最大的正整數為
，那麼
，故：
所以在 中必然存在正整數
滿足
，故
。
直觀上理解就是 「增加」的過程中「超過」了
，而每次的「步長」又小於
，所以一定存在一個臨界的
滿足
。
其實編程簡單地枚舉一下就知道 是滿足條件的一個正整數，但是精確地計算階乘對於計算機來說開銷不少。 @SarvaTathagata 的回答裏說可以利用
這個式子來避免大整數運算，可是在實際操作中也會遇到精度相關的問題
我們用C++來簡單測試一下，找到所有 以內滿足條件的正整數
// test.cpp
#include &
#include &
template&
void count_factorial(FILE* output) {
int count = 0;
F log_sum = 0;
for (int i = 2; i &< 100000000; ++i) {
log_sum += std::log10((F)i);
F tmp = std::pow((F)10, (F)3 + log_sum - std::floor(log_sum));
if (tmp &>= 2020 tmp &< 2021) {
++count;
fprintf(output, "%d
", i);
}
}
printf("%d
", count);
}
int main() {
count_factorial&(fopen("float.txt", "w"));
count_factorial&(fopen("double.txt", "w"));
count_factorial&(fopen("ldouble.txt", "w"));
return 0;
}
編譯運行，查看輸出，結果令人驚訝：使用float類型(4byte)時一個也沒找著，使用double(8byte)時有21456個，而使用long double類型(16byte)時則有21488個。
這表明在實際運算的過程中，各種近似誤差的疊加最終超過閾值，影響到了最終結果。
利用GMP庫，我們來做個測試。計算 的前16位數字
// gmp_check.c
#include &
int main() {
mpz_t fact, base;
mpz_init_set_ui(base, 10);
mpz_fac_ui(fact, 99999999);
size_t digits = mpz_sizeinbase(fact, 10);
mpz_pow_ui(base, base, digits - 20);
mpz_t result;
mpz_div(result, fact, base);
gmp_printf("%Zd
", result);
return 0;
}
// 輸出16172037949214623863
使用C++標準庫，用三種浮點類型分別計算這個結果
// test2.cpp
#include &
#include &
#include &
template&
F test() {
F log_sum = 0;
for (int i = 2; i &< 100000000; ++i) {
log_sum += std::log10((F)i);
}
return std::pow((F)10, (F)3 + log_sum - std::floor(log_sum));
}
int main() {
using namespace std;
auto s1 = test&();
auto s2 = test&();
auto s3 = test&();
cout &
得到
1
1617.089863098135
1617.203354540494
（float, float你在幹什麼啊float）
重載float類型的庫函數可能實現略有不同，兩次測試都顯現出巨大的差異。（事實上，這裡float類型計算出的對數之和僅在1.3e8左右，與實際值7.6e8差距巨大）
對比正確數值1617.2037949214623863，可以發現double類型在第5位出現不同，而long double類型在第8位出現差異。
BTW，雖然在第一個測試中使用long double類型僅比double類型多找到32個數，但是對比輸出文件可知並情況沒有看上去這麼簡單
comm -3 & diff.txt
這個diff.txt僅有2048行，非常的巧
在dalao @SarvaTathagata 的提點下做了一些改進，取得了顯著的效果。
// test_opt.cpp
#include &
#include &
const long double SUP = std::log10(2.021L);
const long double INF = std::log10(2.020L);
template&
void count_factorial(FILE* output) {
int count = 0;
F log_sum = 0;
for (int i = 2; i &< 100000000; ++i) {
log_sum += std::log10((F)i);
log_sum -= std::floor(log_sum);
if (log_sum &>= INF log_sum &< SUP) {
++count;
fprintf(output, "%d
", i);
}
}
printf("%d
", count);
}
int main() {
count_factorial&(fopen("float2.txt", "w"));
count_factorial&(fopen("double2.txt", "w"));
count_factorial&(fopen("ldouble2.txt", "w"));
return 0;
}
// 個數分別為21658, 21486 和 21486double與long double所得數字完全相同
（我確實太菜了
環境：x86_64，gcc10.2.0 on archlinux
顯然
384! = 20205002129876021521366889509458521857103131384424164703......
8297! = 20207983458785287090983863080884101922314442223694268578......
8795! = 20206439534140183506320813566892771922507666440925556515......
16624! = 20204447515212337619006649257632217013568552466164768543......
17775! = 20201706710997231611432820358540894429946023160331657707......
故至少存在自然數n，使得n！的前四個數字是2020
證明過程如下
#! /usr/bin/python3
sum = 1
math = 1
while True:
sum *= math
sum_str = str(sum)
if (sum_str[0:4] == "2020"):
print("{}! = {}".format(math, sum))
math += 1
#EOF
384! = 202050021298760215213668895094585218571031313844241647038656336753901805116049571674350660157440979591792292856366129916926961706884677301229576788637571823309033713973447373497759772077557498674190878950855558580413822253588970773072621439881846523956766945479875974810824204538651272183933304620434930176256513635227507752422464752474771486396742989664230239433214649824117145479984513359446662469851311138751210976519967687756227751330178388021009033339022681291729619372772508425590802338347614964173532877881213147569421175965148279098826072666361443374874234852854408882112301367626134918709620006585891596435949345051507892173464361536351328134358067227181888523570918732580031107186879656925163054478613903703559155950456668160000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
之前看到過字元串乘法的題目，今天趁著無聊寫著玩~不過好久不寫python了，有點生疏~
def multiply(a, b):
"""字元串乘法，不限制長度
by [email protected]"""
if len(a)&= len(b)
result = "0"
for i in range(len(b)):
result = multiply_a(result, multiply_m(a, b[-i-1], i))
return result
def multiply_m(a, bi, n=0): # """單個數與a相乘，n為bi後面的0
if len(bi)&>1:
return
add = 0
result = ""
bi = int(bi)
for i in a[::-1]:
result = str((bi*int(i)+add) % 10) + result
add = (bi*int(i)+add) // 10
if add:
return str(add)+result + "0"*n
return result + "0"*n
def multiply_a(a, b): #字元串相加
if len(a)&= len(b)
b = "0"*(len(a)-len(b)) + b
add = 0
result = ""
for i in range(-1, -1*len(a)-1, -1):
result = str((int(a[i])+int(b[i])+add) % 10) + result
add = (int(a[i])+int(b[i]) + add) // 10
return str(add)+result if add&>0 else result
if __name__ == __main__:
s = "1" #階乘結果
i = "1" #i的階乘
while 1:
i = multiply_a(i, "1") #i+1
s = multiply(s, i) #計算新的階乘
if s.startswith("2020"):
print(i+"! = ", s, "
")
break
2021年2月8日 15:27:02
沒想到還有人看。我看了一下自己前兩天寫的代碼，什麼玩意，莫名其妙的命名，連個注釋都沒加~
經評論區提醒，python3的整形沒有長度限制，所以不需要用字元串去做，直接算就好了。
result = 1 #階乘結果
i = 1 #i的階乘
n = 2 #結果數量
while n:
i += 1
result *= i #計算階乘
if str(result).startswith("2020"):
print("%d! = %d
"%(i,result))
n -= 1
推薦閱讀：
