1. 程式人生 > >BZOJ3817 清華集訓2014 Sum 類歐幾里得

BZOJ3817 清華集訓2014 Sum 類歐幾里得

傳送門


\(\sqrt r = x\)

考慮將\(-1^{\lfloor d \sqrt r \rfloor}\)魔改一下

顯然它等於\(1-2 \times (\lfloor dx \rfloor \mod 2)\),也就等於\(1 - 2 \times \lfloor dx \rfloor + 4 \times \lfloor \frac{dx}{2} \rfloor\)

那麼我們現在就要求\(\sum\limits_{i=1}^n \lfloor ix \rfloor\)的值,求\(\sum\limits_{i=1}^n \lfloor \frac{ix}{2} \rfloor\)

方法顯然一致

先說自己的一種low到炸精度的做法

首先\(x \geq 1\)時可以直接提出整數項,所以只要考慮\(x<1\)的情況

\(\sum\limits_{i=1}^n \lfloor ix \rfloor=\sum\limits_{i=1}^n \sum\limits_{j=1}^{\lfloor nx \rfloor} [ix > j] = \sum\limits_{j=1}^{\lfloor nx \rfloor}n-\lfloor \frac{j}{x} \rfloor=\lfloor nx \rfloor n - \sum\limits_{j=1}^{\lfloor nx \rfloor}\lfloor \frac{j}{x} \rfloor\)

然後因為超多次的實數除法精度直接爆掉

所以需要一個靠譜一點的做法,將上面的\(i\)轉換一下

考慮求解\(\sum\limits_{i=1}^n \lfloor \frac{ax+b}{c}i \rfloor\),其中\(a,b,c\)為整數

首先避免爆longlong,對\(a,b,c\)同除\(gcd(a,b,c)\)

然後將\(\frac{ax+b}{c}\)代入上面長式子中的\(x\),我們可以得到

\(\sum\limits_{i=1}^n \lfloor \frac{ax+b}{c}i \rfloor=\lfloor \frac{ax+b}{c} n \rfloor n - \sum\limits_{j=1}^{\lfloor nx \rfloor}\lfloor \frac{cj}{ax + b} \rfloor\)

,發現分母裡有根號,有理化一下

就得到\(\sum\limits_{i=1}^n \lfloor \frac{ax+b}{c}i \rfloor = \lfloor \frac{ax+b}{c} n \rfloor n - \sum\limits_{j=1}^{\lfloor nx \rfloor}\lfloor \frac{c(ax-b)}{a^2r - b^2} j\rfloor\),這樣我們的係數都在整數域內,就不會出現太大的精度誤差了。

注意一點:當\(r\)為完全平方數的時候,這樣做是不可行的,因為上式中\(\sum\limits_{i=1}^n \lfloor ix \rfloor=\sum\limits_{i=1}^n \sum\limits_{j=1}^{\lfloor nx \rfloor} [ix > j]\)默認了\(x\)為無理數,若\(x\)為整數應當為\(\sum\limits_{i=1}^n \lfloor ix \rfloor=\sum\limits_{i=1}^n \sum\limits_{j=1}^{\lfloor nx \rfloor} [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) << '\n';
    }
}

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