求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...