【BZOJ3817/UOJ42】Sum（類歐）

題面

BZOJ
UOJ

題解

令\(x=\sqrt r\)，那麼要求的式子是\[\sum_{d=1}^n(-1)^{[dx]}\]
不難發現，對於每個\(d\)而言的取值只和\([dx]\)的奇偶性相關。
如果\(x\)是個整數，也就是\(r\)是完全平方數的時候，顯然是可以直接算答案的。
計算答案的時候顯然之和有幾個奇數或者幾個偶數相關（只要求一個另外一個就是補集）
比如說我們來求有幾個是偶數，那麼要滿足的條件就是\([dx]=2*[\frac{dx}{2}]\)。
先重新寫下式子，我們寫成這個樣子\[\sum_{d=1}^n(1-2*([dx]\ mod\ 2)\]

這個顯然成立，當\([dx]\)是偶數的時候貢獻是\(1\)，奇數的時候貢獻是\(0\)。
而\([dx]\ mod\ 2\)可以寫成減法的形式。所以原式可以寫成
\[\sum_{d=1}^n(1-2*([dx]-2*[\frac{dx}{2}]))\]
化簡之後得到了
\[n+\sum_{d=1}^n (4*[\frac{dx}{2}]-2*[dx])\]
現在把模型轉化一下，變得更加一般一點，那麼要求解的東西是\[\sum _{d=1}^n[\frac{ax+b}{c}d]\]
聽說這個玩意叫做類歐？類歐幾里得演算法。
令\(k=\frac{ax+b}{c}\)，那麼式子可以化簡為\(\sum[kd]\)
這裡進行分類討論
1.當\(k>=1\)的時候
\[\begin{aligned}\sum_{d=1}^n[kd]&=[d*(k-[\frac{ax+b}{c}])]+d*[\frac{ax+b}{c}]\end{aligned}\]
顯然後面那一半是可以直接計算的，而前面這一半令\(k'=k-[\frac{ax+b}{c}]\)
我們可以遞迴處理。
2.\(k<1\)
轉化到一個座標系上面去，我們要求的東西本質上就是有一條直線\(y=kx\)，要求解\(x\in[1,n]\)時，與\(x=n\)、\(x\)正半軸圍成的三角形內部整點的個數。
我們把這個三角形逆時針旋轉\(90\)度，再沿著\(y\)軸翻轉過來，讓長度為\(n\)的那條邊靠著\(y\)軸，這樣子翻轉過來之後，\(k\)變成了倒數，\(n\)變成了\([nk]\)，然後拿矩形減去多出來的部分，不就是矩形減去\(k>=1\)的情況了嗎？
再討論一下取倒數的結果，斜率\(k=\frac{ax+b}{c}\)，取倒數之後\(k=\frac{c}{ax+b}\)，然後分母有理化一下\(k=\frac{c(ax-b)}{a^2r-b^2}\)。
總的時間複雜度是\(log\)級別的。

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
inline int read()
{
    int x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
double x;ll n,r;
ll Calc(ll a,ll b,ll c,ll n)
{
    if(!n)return 0;
    ll d=__gcd(a,__gcd(b,c));a/=d;b/=d;c/=d;
    ll k=(a*x+b)/c;
    if(!k)
    {
        ll N=(a*x+b)/c*n;
        return n*N-Calc(a*c,-b*c,a*a*r-b*b,N);
    }
    else
        return k*n*(n+1)/2+Calc(a,b-c*k,c,n);
}
int main()
{
    int T=read();
    while(T--)
    {
        n=read();r=read();x=sqrt(r);
        ll k=x;
        if(k*k==r)
        {
            if(k&1)printf("%lld\n",n-2*((n+1)/2));
            else printf("%lld\n",n);
        }
        else printf("%lld\n",n+Calc(1,0,2,n)*4-Calc(1,0,1,n)*2);
    }
    return 0;
}