抛开题目不管,先来认识一下类欧几里得算法
类欧几里得
就我所知(我自然是不懂什么的啦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
若纵坐标为
感性认知一下
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 如果是完全平方数显然可以直接得出答案
考虑
怎么办?
显然精度误差不是那么能承受的
考虑用一个三元组(a,b,c)来存表示斜率
斜率为
操作中斜率有两个变化
- 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
注意递归中要约分啊
代码
#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;
}