bzoj|Bzoj3817:Sum bzoj|类欧几里得

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

∑i=1n?d?i?
我们先写一个正比例函数，把d看作斜率
y=d?x
把它放进平面直角坐标系中观察
y=OA

文章图片

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

文章图片

又转化成了原先的形式
n′=?d?n?
d′=1d
可以递归求解
复杂度是 logn的（并不会证明）
回到题目
r=R??√
(?1)?i?r?=1?2(?i?r?mod2)=1?2?(?i?r??2???i?r?2?)
∑i=1n(?1)i??r?=n?2?∑i=1n?i?r?+4?∑i=1n??i?r?2?
R如果是完全平方数显然可以直接得出答案
考虑 R不是完全平方数，这时 r为无理数
怎么办？
~~直接用实数存斜率，然后迭代~~
显然精度误差不是那么能承受的
考虑用一个三元组（a,b,c）来存表示斜率
斜率为 b?r+ca
操作中斜率有两个变化

d=d??d?
d=1d

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

d??d?=b?r+ca??d?=b?r+(c?a??d?)a
第二个怎么办呢？
ab?r+c=a?(b?r?c)(b?r+c)?(b?r?c)=a?b?r+(?a?c)b2?R?c2
a′=b2?R?c2
b′=a?b
c′=?a?c
注意递归中要约分啊 【bzoj|Bzoj3817:Sum】
代码

#include #include #include #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; }