强网杯2021 - ezmath

这题做了我两天,最后看到别人全是猜结论(还不知道有多少py),就我是正解,心情还是比较复杂的盒盒

函数1:

v3是i!,v4是x^i,v5是∑v4/v3,这实际上是一个8225阶麦克劳林公式,近似等于 exe^x

最后又把v4和v5乘起来,所以这个函数返回值为 x8225exx^{8225}e^x ,设它为f(x)

函数2:

上界 n1n_1 下届 n0n_0,实际上这个函数只调了一次,n1和n0一定是1和0

v9=(n1-n0)/1000,这是划分为1000段,v9是段长

v7=∑f(v9/2+i*v9+n0);

v8=∑f(i*v9+n0)

v7有第0项,v8没有

最后返回(4*v7+2*v8+f(n0)+f(n1))*v9/6;

注意到4+2=6,因此很可能是平衡v7和v8用的,又最后乘了个v9,所以这很可能是一个近似计算积分的公式

事实上,这是2阶拉格朗日插值,也叫辛普森公式

函数的结果是从0到1积分f(x),设它为g(x)

函数3:

搞了个smc,把比较函数的初值设置为g(x)

函数4:

比较函数,输入n,然后从8225开始枚举到n,每次让f=e-i*f,其中i=8225->n

把结果跟一个数组enc中的值比较,如果全部匹配成功则得到correct

注意到,g(x)=01x8225exdxg(x)=\int_0^1x^{8225}e^xdx,尝试计算此函数的解析解可以发现

积分In=01xnexdxI_n=\int_0^1x^ne^xdx是一个递推式,In=enIn1I_n=e-n*I_{n-1}

因此事实上,此函数的结果其实是计算In=01xnexdxI_n=\int_0^1x^ne^xdx

但是题目程序的计算方式不合理,由于精度问题,很快结果就会变成inf

因此本题的目标应该是已知InI_n的值,找出对应的n

易证InI_n为递减数列,因此可以用二分查找快速找n,计算InI_n可以直接抄题目中的辛普森公式

这道题还有一个比较有意思的地方,就是喜欢乱搞的选手可能会发现enc[i]=e/n,这个结论其实是由积分的极限保证的

limnn01xnexdx=limnn01xni=0xii!dx\lim_{n\rightarrow \infty}n\int_0^1x^ne^xdx=\lim_{n\rightarrow \infty}n\int_0^1x^n\sum_{i=0}^{\infty}\frac{x^i}{i!}dx

=limnn01i=0xi+ni!dx=limnni=001xi+ndxi!=\lim_{n\rightarrow \infty}n\int_0^1\sum_{i=0}^{\infty}\frac{x^{i+n}}{i!}dx=\lim_{n\rightarrow \infty}n\sum_{i=0}^{\infty}\frac{\int_0^1x^{i+n}dx}{i!}

=limnni=01(i+n+1)i!=limni=01(i+1n+1)i!=\lim_{n\rightarrow \infty}n\sum_{i=0}^{\infty}\frac{1}{(i+n+1)i!}=\lim_{n\rightarrow \infty}\sum_{i=0}^{\infty}\frac{1}{(\frac{i+1}{n}+1)i!}

=i=01i!=e=\sum_{i=0}^{\infty}\frac{1}{i!}=e

因此当n比较大的时候(比如达到0x2020这种级别),就会出现n近似等于e/enc[i]的情况

从这个角度来看,这道题出的还是非常有意思的,只是可惜我出了两次失误,卡了两天才做出来

脚本:

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
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
long double enc[19] = {
0.00009794904266317233, 0.00010270456917442, 0.00009194256152777895,
0.0001090322021913372, 0.0001112636336217534, 0.0001007442677411854,
0.0001112636336217534, 0.0001047063607908828, 0.0001112818534005219,
0.0001046861985862495, 0.0001112818534005219, 0.000108992856167966,
0.0001112636336217534, 0.0001090234561758122, 0.0001113183108652088,
0.0001006882924839248, 0.0001112590796092291, 0.0001089841164633298,
0.00008468431512187874
};
const long double e=exp(1);
//long double ans[110000];
long double ex[11000000];
int bc=1000000;
long double p(long double x,int n){
long double z=1;
for(;n;n>>=1){
if(n&1) z=z*x;
x=x*x;
}
return z;
}
long double g(long double x,int n){
return p(x,n)*exp(x);
}
long double f(int n){
/*
long double b=1.0/bc;
long double s=0;
long double x=0.5*b;
for(int i=0;i<bc;++i){
s+=p(x,n)*ex[i]*b;
x+=b;
}
return s;
*/
long double v9=(long double)1/bc;
long double v7=g(v9/2,n);
long double v8=0;
for(int i=1;i<bc;++i){
v7=g(v9/2+i*v9,n)+v7;
v8=g(i*v9,n)+v8;
}
return (4*v7+2*v8+e)*v9/6;
}
void ini(){
long double b=1.0/bc;
long double x=0.5*b;
for(int i=0;i<bc;++i){
ex[i]=exp(x);
x+=b;
}
}
int bsc(long double x){
int l=15935,r=40865,md;
while(l+1<r){
md=(l+r)>>1;
long double fm=f(md);
(fm>x ? l : r)=md;
//printf("%d %d\n",l,r);
}
return fabs(f(l)-x)<1e-12 ? l : r;
}
int main(){
//printf("%.20LF\n",f(27750));
for(int k=0;k<19;++k){
int ans=bsc(enc[k]);
printf("%c%c",ans%256,ans/256);
fflush(stdout);
}
printf("\n");
return 0;
}