2117-一个区间内所有数乘积的缩写

Raphael Liu Lv10
  2023-09-13 00:47:26 2023-09-13 00:47:26 Created   2023-09-13 16:06:26 2023-09-13 16:06:26 Updated    

给你两个正整数 left 和 right ，满足 left <= right 。请你计算 闭区间 [left, right]
中所有整数的 乘积 。

由于乘积可能非常大，你需要将它按照以下步骤 缩写 ：

  1. 统计乘积中 后缀 0 的数目，并 移除 这些 0 ，将这个数目记为 C 。
    * 比方说，1000 中有 3 个后缀 0 ，546 中没有后缀 0 。
  2. 将乘积中剩余数字的位数记为 d 。如果 d > 10 ，那么将乘积表示为 <pre>...<suf> 的形式，其中 <pre> 表示乘积最 开始 的 5 个数位，<suf> 表示删除后缀 0 之后 结尾的 5 个数位。如果 d <= 10 ，我们不对它做修改。
    * 比方说，我们将 1234567654321 表示为 12345...54321 ，但是 1234567 仍然表示为 1234567 。
  3. 最后，将乘积表示为 字符串 "<pre>...<suf>eC" 。
    * 比方说，12345678987600000 被表示为 "12345...89876e5" 。

请你返回一个字符串，表示 闭区间 [left, right] 中所有整数 乘积 的 缩写 。

示例 1：

**输入：** left = 1, right = 4
**输出：** "24e0"
**解释：**
乘积为 1 × 2 × 3 × 4 = 24 。
由于没有后缀 0 ，所以 24 保持不变，缩写的结尾为 "e0" 。
因为乘积的结果是 2 位数，小于 10 ，所欲我们不进一步将它缩写。
所以，最终将乘积表示为 "24e0" 。

示例 2：

**输入：** left = 2, right = 11
**输出：** "399168e2"
**解释：** 乘积为 39916800 。
有 2 个后缀 0 ，删除后得到 399168 。缩写的结尾为 "e2" 。 
删除后缀 0 后是 6 位数，不需要进一步缩写。 
所以，最终将乘积表示为 "399168e2" 。

示例 3：

**输入：** left = 371, right = 375
**输出：** "7219856259e3"
**解释：** 乘积为 7219856259000 。

提示：

  • 1 <= left <= right <= 104

大家好，今天再次给大家带来一篇团灭系列的文章。这里我会分享一些构造极端数据的心得，介绍误差分析的一些知识，以及如何在合理的精度下通过本题。
要想成功叉掉标程，首先要有不错的直觉：这题很可能出现精度问题。LC的不少出题人似乎不会误差分析的样子，出的数据范围也比较随意(常见问题是只考虑时限，不考虑数据精度的限制)，所以这已经不是第一个被我叉掉的标程了…

误差分析

猜到标程可能出问题之后，我们就可以动手寻找使程序出现精度问题的极端情况。首先估算一下这题极限情况下的精度要求，这需要一些误差分析的知识：如果我们把 n 个随机数乘起来，每乘一次后四舍五入截断产生 \varepsilon 的(相对)误差，那么最终期望的相对误差为 \sqrt{n}\cdot \varepsilon 的量级。这是因为根据随机游走的分析，每次四舍五入时可能有一半的概率高估，一半的概率低估，绝大部分的高估与低估会互相抵消。注意在最坏情况下，如果我们每次都低估(截断时舍去)，那么最终的相对误差会是 n\cdot \varepsilon 量级的，所以四舍五入会比手动舍去尾数的做法精度更好。注意我这里为了方便分析把 [l,r] 内的数看成随机的了。
如果我们用 double 进行计算的话，\varepsilon大概是 1e-16 的量级(2^{-52)。那么对于 n=10^6，期望的误差是\sqrt{n}\cdot \varepsilon\approx 1e-13 的量级。

极端数据的构造

下一步是构造一个极端情况下的数据，使得计算过程中微小的精度误差可以影响最终的输出结果。比如如果真实值是 xxxx40001… 的形式，而我们的(绝对)误差低估了 2… 的形式(尾部省略的位数相同)，那么我们实际输出的数会变成 xxxx39999… 的形式，就错了。也就是说我们要让构造的真实值的从高往低数第6位开始尽可能地接近 0000…或者9999… 的形式(具体是哪种需要根据想叉掉的程序是如何舍入的来确定)。
那么最坏情况下我们在第6位后可以构造出几个连续的9呢？如果我们把 O(n^2) 种可能输出的数也看作是随机的(这样简化处理是因为受到了分析技巧的限制)，那么对于最坏的一个输入我们大概可以构造出 12 个连续的 9，也就是说使得 10^{-16}\sim 10^{-17 量级的相对误差就可以产生错误。这远小于 1e-13，表明用double或者long long进行计算是不可靠的。有了量级估算之后，我写代码暴力枚举出了一些较坏的数据。下面就让我们来看看周赛前五名的代码能否顺利通过我们精心构造的数据吧！
第五名 @suspicious-keldyshlkl：

1
2
6
563035

输出 78782，实际值为 7878300000012…。

第四名 @tsreaper：
以上数据正确。加大难度。

1
2
18
237575

输出 76682，实际值为 76683000000061…。

第三名 @darrenhp：
以上数据正确。加大难度。

1
2
3940
836328

输出 24931，标程输出24930，实际值为 249310000000000016…。标程挂了，让我们恭喜这位选手。不过我们还是要继续加大难度：

1
2
2965
574229

输出 89071，实际值为 8907099999999951…。

第二名 @tonghuikang：
第二名用的是 python 啊，我很慌张。不过还是挂在叉掉第四名的数据上了。

第一名 @arignote：
第一名用的是 Java，虽然跟第三名的 C++ 代码计算方式是一样的(double 一直乘，超过 100000 就 /10)，但是换成 Java 就神奇地过了。以上数据全部正确。继续加大难度。

1
2
6148
215373

输出 36553，实际值为 36552999999999913…。至此前五名全部团灭。

如何通过本题

在平衡运算速度与精度后，可以使用C++的__int128或者__float128进行计算(对应的 \varepsilon 约为 1e-34 的量级)，或者用python自动进行高精度计算。C++的long double似乎也是可行的，我暂时没有叉掉，但是理论精度正好处在将要挂掉的边界上。Java 的大整数类我不确定够不够快，或者手写高精度也可以。我看到了 Go 的高精度还是挺快的。其他语言就自寻出路吧。
更新：我花了两天搜索了所有容易出错的数据，没有卡掉正确实现的 long double。其中最难的数据大概是

1
2
468620
595075

实际值为 1655700000000000135…，用 long double 相乘算出来的是 1655700000000000075…，偏差稍微差了那么一点，再翻一倍的话就有可能出错了。比赛中唯一通过的 C++ 选手的写法为

1
2
3
4
5
6
7
long double pre = 1.0;
for (int x = left; x <= right; x++) {
	pre = pre * x;
	while(pre > 100000) pre /= 10;
}
//注：换成 while(pre > 1000000000) pre /= 1000000000; 的话精度会更好，因为少做了几次除法。
//把数据范围放到 10^9 可以明显地看出区别。

我也试了用 long double 求 log 之和再 pow 恢复的方案，求完 log 之和精度就差很远了，算出的是 165569999999948…。尽力优化细节的话可以把精度提升两个数量级，但还是会差一些。


另外再更新一个 \log^{O(1)} n 复杂度的算法，思路是用 Stirling 近似公式来快速计算 \log (n!)。以下代码用的 long double，所以还是有些精度问题，但可以通过目前LC上的所有数据。我再手写了个 float128 版本，精度就有保障了，能算到三十几位。注意到 [l,r] 内的数乘积末尾的 L 位非零值也可以在 \mathrm{poly}(L,\log n) 的复杂度内算出(可以看[这里](https://emathgroup.github.io/blog/factorial-tail))，所以本题是有 \mathrm{polylog}~n 复杂度的做法的。
[code1-long double版]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#define double long double
const double PI=atan((double)1)*4,E=exp((double)1),
  B[]={1,1./2,1./6,0,-1./30,0,1./42,0,-1./30,0,5./66,0,-691./2730,0,7./6,0,-3617./510,0,43867./798,0,-174611./330};
const int M=10;
double log_fac(double n){
	double ans=n*log(n/E)+log(n)/2+log(sqrt(PI*2));
	for (int i=2;i<=M;++i)ans+=pow(-1,i)*B[i]/(i*(i-1)*pow(n,i-1));
	if (n<20){
		ans=0;
		for (int i=1;i<=n;++i)ans+=log(i);
	}
	return ans;
}
int leading(int l,int r){
	double t=log_fac(r)-log_fac(l-1);
	t/=log((double)10); t-=int(t)-4;
	return pow(10,t);
}
[code1-float128版]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
typedef __float128 D;
const D PI=3.1415926535897932384626433832795028841971Q,
	E=2.7182818284590452353602874713526624977572Q;
void print(__float128 x,int K=35){
	if (x<0){printf("-");x=-x;}
	int y=0,t;
	while (x>=1)x*=0.1Q,++y;
	if (!y)printf("0.");
	for (int i=1;i<=max(y,K);++i){
		x*=10; t=(int)x;
		if (x<t)--t;
		else if (x>=t+1)++t;
		printf("%d",t);
		if (i==y)printf(".");
		x-=t;
	}
	printf("\n");
}
inline __float128 sqrt(__float128 x){
	__float128 f=1.5Q,x1=x*0.5Q,y=1.Q/sqrt((double)x);
	for (int i=1;i<=4;++i)y=y*(f-x1*y*y);
	return x*y;
}
inline __float128 exp(__float128 x){
	if (x<0)return 1/exp(-x);
	int t=0; for (;x>1e-5;++t)x/=2;
	__float128 ans=0,y=x;
	for (int i=1;i<=7;++i,y*=x/i)ans+=y;
	while (t--)ans=ans*2+ans*ans;
	return ans+1;
}
inline __float128 log(__float128 x){
	if (x<1)return -log(1/x);
	int t=0; for (x-=1;x>1e-5;++t)x=x/(sqrt(x+1)+1);
	__float128 ans=0,y=x;
	for (int i=1;i<=7;++i,y*=-x)ans+=y/i;
	ans*=1<<t;
	return ans;
}
inline __float128 pow(__float128 x,__float128 y){return exp(y*log(x));}
const D B[]={1,D(1)/2,D(1)/6,0,-D(1)/30,0,D(1)/42,0,-D(1)/30,0,D(5)/66,0,
	-D(691)/2730,0,D(7)/6,0,-D(3617)/510,0,D(43867)/798,0,-D(174611)/330};
const int M=10;
D log_fac(D n){
	D ans=0;
	if (n>20){
		ans=n*log(n/E)+log(n)/2+log(sqrt(PI*2));
		for (int i=2;i<=M;++i)ans+=(i%2?-1:1)*B[i]/(D(i)*(i-1)*pow(n,D(i-1)));
	}
	else for (int i=1;i<=n;++i)ans+=log(D(i));
	return ans;
}
int leading(int l,int r){
	D t=log_fac(r)-log_fac(l-1);
	t/=log(10.Q); t-=int(t)-4;
	return pow(10.Q,t);
}

三更：现在数据范围被改小到 n\leq 10000 了，不过还是得小心处理细节才能用 double 通过本题。比如我测试了一下比赛中前 150 名中的所有 82 位 C++ 选手，其中只有 24 位通过了下面这个数据： 
1
2
4838
6186
实际值为 36088000000027...。 另外再提供一些数据： 
1
2
3
4
5
6
7
8
9
10
11
12
13
//格式为left, right, 实际值的前5位。
5451,7344,"36274..."
1368,1577,"16375..."
1385,5673,"46222..."
5387,9507,"11343..."
445,7931,"15669..."
2230,5489,"34409..."
3175,8949,"31926..."
3317,5963,"19446..."
3372,5890,"10907..."
8229,8924,"31684..."
3543,4907,"45166..."
5975,9674,"36300..."

最后贴一下完整版的 O(\mathrm{poly}\log n) 复杂度的代码，可以跑到 0ms： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#define double long double
const double PI=atan((double)1)*4,E=exp((double)1),
  B[]={1,1./2,1./6,0,-1./30,0,1./42,0,-1./30,0,5./66,0,-691./2730,0,7./6,0,-3617./510,0,43867./798,0,-174611./330};
const int M=10;
double log_fac(double n){
	double ans=n*log(n/E)+log(n)/2+log(sqrt(PI*2));
	for (int i=2;i<=M;++i)ans+=pow(-1,i)*B[i]/(i*(i-1)*pow(n,i-1));
	if (n<20){
		ans=0;
		for (int i=1;i<=n;++i)ans+=log(i);
	}
	return ans;
}
int leading(int l,int r){
	double t=log_fac(r)-log_fac(l-1);
	t/=log((double)10); t-=int(t)-4;
	return pow(10,t);
}

const int P=3125;
template<class T> T extend_gcd(T a,T b,T &x,T &y){
	if (!b){x=1;y=0;return a;}
	T res=extend_gcd(b,a%b,y,x); y-=x*(a/b);
	return res;
}
template<class T>
inline T inv(T a,T M){
	T x,y; extend_gcd(a,M,x,y);
	return (x%M+M)%M;
}
const int inv2=inv(2,P);
int exp(int x,int y,int z){
	int ans=1;
	while (y){
		if (y&1)ans=ans*x%z;
		x=x*x%z;y>>=1;
	}
	return ans;
}
int Chinese_Remainder_Theorem(int A[],int M[],int n){
	for (int i=1;i<n;++i){
		int x,y,P=M[0]*M[i];
		extend_gcd(M[0],M[i],x,y);
		x=(x*(A[i]-A[0])%P+P)%P;
		A[0]=(x*M[0]+A[0])%P; M[0]=P;
	}
	return A[0];
}
// PolynomialMod[Product[If[Mod[i,5]>0,5^4*x+i,1],{i,1,5^4}],3125]
int calc(int x,int d){
	switch(d){
		case 0: return x%5?x%P:1;
		case 1: return (24+250*x+875*x%P*x+1250*x%P*x%P*x+625*x%P*x%P*x%P*x)%P;
		case 2: return 124;
		case 3: return 624;
		case 4:
		case 5: return 3124;
	}
	return 0;
}
int get(int l,int r,int d){
	if (l>r||d<0)return 1;
	int mod=pow(5,d),l0=(l-1)/mod+1,r0=(r+1)/mod,res=1;
	for (int i=l0;i<r0;++i)res=res*calc(i,d)%P;
	if (l0<r0)return res*get(l,l0*mod-1,d-1)%P*get(r0*mod,r,d-1)%P;
	else return get(l,r,d-1);
}
string abbreviateProduct1(int left, int right) {  // brute force. 随便抄了一个
    long long suff = 1, c = 0, total = 0, max_suff = 100000000000LL;
    double pref = 1.0;
    for (int i = left; i <= right; ++i) {
        pref *= i;
        suff *= i;
        while (pref >= 100000) {
            pref /= 10;
            total = total == 0 ? 6 : total + 1; 
        }
        while (suff % 10 == 0) {
            suff /= 10;
            ++c;
        }
        suff %= max_suff;
    }
    string s = to_string(suff);
    return to_string((int)pref) + (total - c <= 10 ? "" : "...") 
    	+ (total - c < 5 ? "" : s.substr(s.size() - min(5LL, total - c - 5))) + "e" + to_string(c);
}
class Solution {
public:
	string abbreviateProduct(int l, int r) {
		if (r-l<20)return abbreviateProduct1(l,r);
		int last=1,l0=l,r0=r,zeros=0;
		while (l0<=r0){
			last=last*get(l0,r0,5)%P;
			l0=(l0-1)/5+1; r0/=5;
			zeros+=r0-l0+1;
		}
		last=last*exp(inv2,zeros,P)%P;
		int A[]={0,last},M[]={32,3125};
		last=Chinese_Remainder_Theorem(A,M,2);
		char str1[15]; sprintf(str1,"%05d",last);
		return to_string(leading(l,r))+"..."+string(str1)+"e"+to_string(zeros);
	}
};
![2117_0ms_cn.png](https://pic.leetcode-cn.com/1656929244-XDorEn-2117_0ms_cn.png)
 Comments
On this page
2117-一个区间内所有数乘积的缩写
  1. 误差分析
  2. 极端数据的构造
  3. 如何通过本题