强网杯2021 - ezmath

发表于2021-06-15更新于2021-06-15阅读次数

这题做了我两天，最后看到别人全是猜结论（还不知道有多少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=e−n∗In−1I_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，这个结论其实是由积分的极限保证的

lim⁡n→∞n∫01xnexdx=lim⁡n→∞n∫01xn∑i=0∞xii!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

=lim⁡n→∞n∫01∑i=0∞xi+ni!dx=lim⁡n→∞n∑i=0∞∫01xi+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!}

=lim⁡n→∞n∑i=0∞1(i+n+1)i!=lim⁡n→∞∑i=0∞1(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=0∞1i!=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;
}