UOJ42/BZOJ3817 清华集训2014 Sum 类欧几里得

传送门

---

令(sqrt r = x)

考虑将(-1^{lfloor d sqrt r floor})魔改一下

它等于(1-2 imes (lfloor dx floor mod 2))，也就等于(1 - 2 imes lfloor dx floor + 4 imes lfloor frac{dx}{2} floor)

那么我们现在就要求(sumlimits_{i=1}^n lfloor ix floor)的值，求(sumlimits_{i=1}^n lfloor frac{ix}{2} floor)方法一致

先说自己的一种low到炸精度的做法

首先(x geq 1)时可以直接提出整数项，所以只要考虑(x<1)的情况

(sumlimits_{i=1}^n lfloor ix floor=sumlimits_{i=1}^n sumlimits_{j=1}^{lfloor nx floor} [ix > j] = sumlimits_{j=1}^{lfloor nx floor}n-lfloor frac{j}{x} floor=lfloor nx floor n - sumlimits_{j=1}^{lfloor nx floor}lfloor frac{j}{x} floor)

然后因为超多次的实数除法精度直接爆掉

所以需要一个靠谱一点的做法，将上面的(x)转换一下

考虑求解(sumlimits_{i=1}^n lfloor frac{ax+b}{c}i floor)，其中(a,b,c)为整数

首先避免爆longlong，对(a,b,c)同除(gcd(a,b,c))

然后将(frac{ax+b}{c})代入上面长式子中的(x)，我们可以得到

(sumlimits_{i=1}^n lfloor frac{ax+b}{c}i floor=lfloor frac{ax+b}{c} n floor n - sumlimits_{j=1}^{lfloor nx floor}lfloor frac{cj}{ax + b} floor)，发现分母里有根号，有理化一下

就得到(sumlimits_{i=1}^n lfloor frac{ax+b}{c}i floor = lfloor frac{ax+b}{c} n floor n - sumlimits_{j=1}^{lfloor nx floor}lfloor frac{c(ax-b)}{a^2r - b^2} j floor)，这样我们的系数都在整数域内，就不会出现太大的精度误差了。

注意一点：当(r)为完全平方数的时候，这样做是不可行的，因为上式中(sumlimits_{i=1}^n lfloor ix floor=sumlimits_{i=1}^n sumlimits_{j=1}^{lfloor nx floor} [ix > j])默认了(x)为无理数，若(x)为整数应当为(sumlimits_{i=1}^n lfloor ix floor=sumlimits_{i=1}^n sumlimits_{j=1}^{lfloor nx floor} [ix geq j])。特判其实比较方便

#include<bits/stdc++.h>
#define int long long
#define ld long double
//This code is written by Itst
using namespace std;

inline int read(){
    int a = 0;
    char c = getchar();
    bool f = 0;
    while(!isdigit(c) && c != EOF){
        if(c == '-')
            f = 1;
        c = getchar();
    }
    if(c == EOF)
        exit(0);
    while(isdigit(c)){
        a = (a << 3) + (a << 1) + (c ^ '0');
        c = getchar();
    }
    return f ? -a : a;
}

int N , R;
ld P;

inline int gcd(int a , int b){
    if(!b)
        return a;
    int r = a % b;
    while(r){
        a = b;
        b = r;
        r = a % b;
    }
    return b;
}

int solve(int a , int b , int c , int rg){
    if(rg <= 0)
        return 0;
    int t = gcd(a , gcd(b , c));
    a /= t;
    b /= t;
    c /= t;
    int cur = (a * P + b) / c;
    if(!cur)
        return (int)((a * P + b) / c * rg) * rg - solve(a * c , -b * c , a * a * R - b * b , (a * P + b) / c * rg);
    else
        return cur * (rg * (rg + 1) / 2) + solve(a , b - c * cur , c , rg);
}

void work(){
    for(int T = read() ; T ; --T){
        N = read();
        R = read();
        P = sqrt(R);
        if((int)P * (int)P == R)
            if((int)P & 1)
                cout << (N & 1 ? -1 : 0) << endl;
            else
                cout << N << endl;
        else
            cout << N - 2 * solve(1 , 0 , 1 , N) + 4 * solve(1 , 0 , 2 , N) << '
';
    }
}

signed main(){
#ifndef ONLINE_JUDGE
    freopen("in" , "r" , stdin);
    freopen("out" , "w" , stdout);
#endif
    work();
    return 0;
}