调和级数计算

求H(n).如果精度要求很高..

1. n<=1e8

调和级数是能直接apply Binary Splitting method的,这样能计算出全精度的H(n)(分数形式).

因为大家都知道好的Binary Splitting method是很快的(真的么..),用来计算Pi以复杂度的劣势能艹爆AGM..

2. n很大(monstrously large)

这个时候就需要Euler-Maclaurin展开..

大家都知道的这个式子:

H(n)=ln(n)+EulerGamma+1/(12n)-sum(i>=1,B_2i/(2i n^(2i)))

(我曾经以为这东西是发散的..事实上它就是发散的..不过似乎没什么毛关系)

现在考虑分别计算这些东西..

ln(n)计算肯定是不能直接Taylor一发的..会boom..

可以考虑先reduce到1+p,然后再对exp进行牛迭..

如果p特别小可以考虑taylor一下..不过要换个series..

EulerGamma没有直接好用的series,不过还是可以很快速的计算(quasi-linear!)的..比较崎岖..(还要吐槽一下它那个崎岖的做法里使用了调和级数..不过收敛得很快)

后面那个东西因为是Bernoulli不方便变成HGS Sum..不过n以下的Bernoulli数是可以计算的,选O(n)个>=n+3的小素数求出bernoulli mod p(这个可以转化为卷!积!mod p),然后CRT一下.Bernoulli它有个好的性质就是它的分母超级他妈好算..num(B2k)=prod(p-1|2k & prime(p),p)..

到wiki上翻一下bernoulli的asymptotic behavior..发现它B2n~4sqrt(pi n)(n/(pi e))^(2n)..那么它的位数大概是2.8n..

然后我们来看看那个Euler-Maclaurin Formula的渐进behavior..估计一下余项大概是B_(2n+2)/((2n+4) x^(2n+4))~4 sqrt(pi n)/(2n+4)*((x)/(n pi e))^(2n+4) 毛咕咕是2n位十进制精度..

好全篇民科hhh...