Bzoj3817:Sum

抛开题目不管，先来认识一下类欧几里得算法

类欧几里得

就我所知(我自然是不懂什么的啦TAT)，类欧几里得算法大致是用来求解一类问题形如

$$\sum_{i=1}^n \lfloor d*i\rfloor$$
我们先写一个正比例函数，把d看作斜率
$$y=d*x$$
把它放进平面直角坐标系中观察
$$y=OA$$

感性认知一下，是不是像求一个三角形（上图中的 $OAB$）内的整点个数
这里的整点指横纵坐标皆为正整数的点，后文也都是这个意思
若d>=1那么式子可以变成
$$\sum_{i=1}^n \lfloor d\rfloor*i+(d-\lfloor d\rfloor)*i$$
显然的 $\sum_{i=1}^n\lfloor d\rfloor*i$ 可以直接得出
$(d-\lfloor d\rfloor)$<1
下面我们只考虑d<1的情况
由于d<1，
所以上图由 $y=OA$ 变成了 $y=OC$
所以纵坐标范围 $d*n$小于横坐标范围 $n$
若纵坐标为$j$，横坐标即为 $\frac{j}{d}$
感性认知一下
$y=j$ 与三角形$OCD$相交的那一部分的长度为$\frac{j}{d}$ 吧
三角形$OBC$的整点数=四边形$OBCD$的整点数-三角形$OCD$ 的整点数+$OC$上的整点数
四边形$OBCD$的整点数显然为$\lfloor d*n\rfloor*n$ 个吧
若d为有理数，可以表示成 $\frac{a}{b}$（自然还是要约分的）的形式
$OC$上的整点数= $\lfloor \frac{n}{b} \rfloor$
无理数就不会有交点了是吧
把三角形$OCD$ 放进另一个平面直角坐标系中

又转化成了原先的形式
$$n'=\lfloor d*n\rfloor$$
$$d'=\frac {1}{d}$$
可以递归求解
复杂度是 $\log n$ 的（并不会证明）

回到题目

$$r=\sqrt{R}$$
$$(-1)^{\lfloor i*r \rfloor}=1-2(\lfloor i*r \rfloor\mod 2)=1-2*(\lfloor i*r \rfloor-2*\lfloor \frac {\lfloor i*r \rfloor}{2}\rfloor)$$
$$\sum_{i=1}^n (-1)^{i*\lfloor r \rfloor}= n-2*\sum_{i=1}^n\lfloor i*r \rfloor +4*\sum_{i=1}^n\lfloor \frac {\lfloor i*r \rfloor}{2}\rfloor$$
$R$ 如果是完全平方数显然可以直接得出答案
考虑$R$ 不是完全平方数，这时 $r$ 为无理数
怎么办？
直接用实数存斜率，然后迭代
显然精度误差不是那么能承受的
考虑用一个三元组（a,b,c）来存表示斜率
斜率为$\frac{b*r+c}{a}$
操作中斜率有两个变化

1. $d=d-\lfloor d\rfloor$
2. $d=\frac{1}{d}$

第一个可以直接在c上面开刀

$$d-\lfloor d\rfloor=\frac{b*r+c}{a}-\lfloor d\rfloor=\frac{b*r+(c-a*\lfloor d\rfloor)}{a}$$
第二个怎么办呢？
$$\frac{a}{b*r+c}=\frac{a*(b*r-c)}{(b*r+c)*(b*r-c)}=\frac{a*b*r+(-a*c)}{b^2*R-c^2}$$
$$a'=b^2*R-c^2$$
$$b'=a*b$$
$$c'=-a*c$$
注意递归中要约分啊

代码

#include <cstdio>
#include <cmath>
#include <algorithm>
#define gec getchar
using namespace std;
typedef long long ll;
char ch;
bool fl;
inline void read(int &a){
    for(fl=0,ch=gec();ch<'0'||ch>'9';ch=gec()) fl|=(ch=='-');
    for(a=0;ch>='0'&&ch<='9';ch=gec())a=(a<<3)+(a<<1)+(ch^'0');
    if(fl)a=-a;
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int n,R,T;
double x;
int Solve(int n,int a,int b,int c){
    if(!n)return 0;
    int tmp=gcd(gcd(a,b),c);
    a/=tmp;b/=tmp;c/=tmp;
    tmp=(b*x+c)/a;c-=tmp*a;
    int Sum=((ll)n*(n+1)>>1)*tmp;
    int P=(b*x+c)/a*n;
    return Sum-Solve(P,b*b*R-c*c,a*b,-a*c)+(ll)P*n;
}
int main()
{
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&R);
        x=sqrt(R);
        int w=x;
        if(w*w==R) printf("%d\n",(w&1)?((n&1)?-1:0):n);
        else printf("%d\n",n+(Solve(n,2,1,0)<<2)-(Solve(n,1,1,0)<<1));
    }
    return 0;
}