隱藏在素數規律中的Pi -- BZOJ1041解題報告

退役狗在刷程書的過程中看到了一個有趣的視頻, 講解了一個有趣的問題.

在網上隨便搜索了一下居然還真的找到了一道以它爲背景的OI題目, BZOJ1041.

下面的內容會首先回顧一下視頻所討論的知識, 有了這些知識, 自然就明白應該怎麼去做這道題了.

另外, 這道題在網路上的解法大多是暴力求解, 這篇博客也提出了本題一個複雜度爲(mathcal{O}(分解質因數))的解決方案.

隱藏在素數規律中的(pi)

首先介紹收斂於(pi)的一個無窮級數:

[pi = lim_{N ightarrow infty} frac 4N sum_{i=1}^N sum_{d|i}chi(d) ]

怎麼證明它呢? 一種常見方法是使用Calculus, 另外一種就是使用數論方法.

Consider有一個網格, 我們畫出一個半徑爲(r)的圓:

[x^2+y^2=r^2 ]

那麼, 如果我們能夠數出來這裏面有多少個點, 我們就可以間接計算出(pi)的值.

我們把這個圓分成若干個小圓, 假設每個圓的方程是:

[x^2+y^2=n ]

我們可以枚舉所有可能的(n)來確定圓中有多少個格點.

考慮到格點的(x, y)都是integer, 那麼我們只需要去枚舉所有小於(r^2)的整數(n)就可以了.

那麼現在問題就轉化成:

給定一個正整數(n), 有多少個pair((x, y))使得(x^2+y^2=n)成立?

考察(n)的任意一個分解(n = x^2+y^2), 若存在這樣一個分解, 這樣的分解也可以寫成:

[n = (x+ymathrm{i})(x-ymathrm{i}) ]

考慮(n)的質因數分解:

(n = prod p_i^{a_i})

對於其中的一個素數, 我們不加證明地給出幾個結論:

• (p_i)能夠被共軛分解的衝要條件是(p_iequiv1 mod 4), 或者(p_i = 2).
• 上面的分解被視爲"Almost Unique"的, 因爲如果我們對於一個分解中的偶數個因子都乘以(1, -1, i, -i), 那麼新的分解仍然是成立的.

---

然後, 我們可以對幾種因子作分類討論:

1. (p_i mod 4 = 1)

若出現這樣的情形, 在不考慮乘以(1, -1, i, -i)的情況下, 我們可以將這分解出來的(2a_i)個複數隨意組合相乘, 可以證明得到的兩個共軛複數不同的情形有(a_i+1)種.

2. (p_i mod 4 = 3)

這樣的素數叫做高斯素數. 若出現這樣的情形, 如果(a_i)是偶數, 我們可以在之前得到了的兩個共軛複數上都乘以(p_i^{a_i/2})而不影響結果. 否則, 方案數就是0.

3. (p_i = 2)

若出現這樣的情形, 理論上來說2可以分解成(2 = (1+i)(1-i)), 但是有((1-i)*i = (1+i)), 從而不管怎麼分組, 都會與最終的乘以(1, -1, i, -i)相重複, 所以對於方案數目的貢獻只有( imes 1)

---

經過上面的討論, 我們就得到了計數拆分(n = x^2+y^2)方案的一個方法. 具體的方案與代碼在下面的題解中給出.

到現在, 這個結果已經蠻不錯, 但是我們還可以做的更優美.

定義函數:

[egin{eqnarray}chi(x)= egin{cases} 1, &x mod 4 = 1cr 0, &x mod 2 = 0 cr -1, &x mod 4 = 3end{cases} end{eqnarray} ]

這個函數是一個積性函數.

通過(chi(x))定義中的分類討論, 在求解的時候, 我們便不必再分類討論了.

這樣, 我們要求得的答案就可以寫成:

[mathcal{ans} = 4 imes prod sum_{0 leqslant i leqslant a_i} chi(p_i^i) ]

考慮到上述函數是一個積性函數, 我們有:

[ans = 4 imes sum_{i=1}^N sum_{d|i}chi(d) ]

對上述函數稍作整理, 並令(lim_{n ightarrow infty})近似, 就得到了要證的式子.

[mathcal{Q.E.D} ]

回到本題

本題是一個上面結論的退化形式, 令(n = r^2), 直接應用算法即可.

下面是代碼.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll xjb=10;
ll mmul(ll a, ll b, ll m){
    ll d=((long double)a/m*b+1e-8);
    ll r=a*b-d*m;
    return r<0?r+m:r;
}
ll mpow(ll a, ll b, ll m){ll r=1;for(;b;b>>=1,a=mmul(a,a,m))if(b&1)r=mmul(r,a,m);return r;}
ll gcd(ll a, ll b){return a?gcd(b%a,a):b;}
ll prime(ll n){
    if(n==1) return 0;
    if(n==2||n==3||n==5) return 1;
    if(!(n&1)||(n%3==0)||(n%5==0)) return 0;
    ll m=n-1; ll k=0;
    while(!(m&1)) m>>=1, k++;
    for(ll tt=0; tt<xjb; ++tt){
	ll x=mpow(rand()%(n-2)+2,m,n), y=x;
	for(ll i=0; i<k; ++i){
	    x=mmul(x,x,n);
	    if(x==1&&y!=1&&y!=n-1) return 0;
	    y=x;
	}
	if(x!=1) return 0;
    }
    return 1;
}
ll f[105]; ll M;
ll rho(ll n, ll c){
    ll x=rand()%n, y=x, t=1;
    for(ll i=1, k=2; t==1; ++i){
	x=(mmul(x,x,n)+c)%n;
	t=gcd(x>y?x-y:y-x, n);
	if(i==k) y=x, k<<=1;
    }
    return t;
}
void work(ll n){
    if(n==1) return;
    if(prime(n)){f[M++]=n; return;}
    ll t=n;
    while(t==n) t=rho(n, rand()%5+1);
    work(t); work(n/t);
}
int main(){
    srand(19260817);
    ll n; 
    scanf("%lld", &n);
    n*=n;
    M=0;
    work(n);
    sort(f, f+M);
    ll ans = 1;
    for(ll i=0, c=1; i<M; ++i){
	if(f[i]!=f[i+1]) {
//	    cout << f[i] << ' ' << c << endl;
	    if(f[i] == 2) { } 
	    else if(f[i] % 4 == 1) {
		ans *= c+1;
	    } else {
		if(c % 2 == 1) ans = 0;
	    }
	    c = 1;
	}
	else c++;
    }
    printf("%lld
", ans*4);
    return 0;
}

---

1. 退役老年選手碼力下降太多了啊....告別OI以後AC的第一道題...感覺還是OI比較美好啊...可惜再也回不去了.

2. 樸素的根號分解過不去這個題, 從網上找了一個Pollard-Rho模板改了改AC了這個題.

3. 又頹廢了一個下午....還有100天就高考了....這樣下去是要gg啊...