抛开题目不管,先来认识一下类欧几里得算法
类欧几里得 就我所知(我自然是不懂什么的啦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
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;
}
推荐阅读
- BZOJ|BZOJ2763[JLOI2011]飞行路线【分层图最短路】
- BZOJ3817(Sum(类欧几里得))
- 类欧几里得算法总结
- 类欧几里得|bzoj2987 Earthquake 类欧几里得
- bzoj2712 -- 类欧几里得算法
- Bzoj|[BZOJ2187][fraction][类欧几里得算法]
- 类欧几里得算法|[BZOJ2712][[Violet 2]棒球][类欧几里得算法]
- 2017|[BZOJ3817][Sum][类欧几里得算法 数论]
- 凸包|bzoj5317: [Jsoi2018]部落战争【凸包/Minkowski sum】