传送门:http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=5815&HPPROTID=f0b51c92
个人第一场ACM,虽然6题Ag滚粗但还是非常值得回味的。
虽然预料到会有高精度题,但是题目的疯狂程度还是超乎想象(当然没有队伍现场ac了)。
题目大意:定义p序列为{1,1,2,2,2,3,3,3,3……}即连续2个1,连续3个2,连续4个3……,再定义q序列为:
q1 = 1,q(n) = q(n - 1) + q(pi)
t(约1e4)组询问,给定n(<=1e40),求q(n)。
看到巨大的n,显然的想法是如何将其变成一个可以接受的范围。
第一步还是比较容易想到的:
q(n)=sigma i=1~n q(pi)
然后将q(pi)用同样的方式展开:
sigma i=1~n sigma j=1~pi q(pj)
分析p序列的性质,发现pi=最大的j使得j*(j+1)/2<=i。
所以,对于某一个q(pj),会在i=j*(j+1)/2~n被计算贡献
sigma j=1~pn q(pj) * (n - (j*(j+1)/2 - 1))
=n*sigma j=1~pn q(pj) - sigma j=1~pn q(pj)*(j*(j+1)/2-1)
到这一步为止,我们成功地把1~n变成了1~pn。
由于pn是O(sqrtn)级别,看起来只需要不停地递归迭代下去,迭代4次就能变成几百左右(实际上p(p(p(pn)))=605)。
但这就结束了吗?
观察上式,我们发现第一部分确实是可以直接递归的,然而第二部分的每一项还需要乘一个跟j有关的多项式。
怎么办?不急,我们继续推式子。
假设我们要求的是sigma i=1~n q(pi)*f(i),其中f(i)是某个多项式。
仍然沿用之前的推法,我们能够得到:
记sn = sigma i=1~n fi为f的前缀和。
sn * sigma i=1~pn q(pi) - sigma i=1~pn q(pi)*s(i*(i+1)/2-1)
这个式子与之前的式子形式很类似。实际上,如果把之前的式子看做是q(pi)*1的话,两者是相同的。
所以我们可以进行如下操作:
定义f(0,n)=1,f(i,n)=s(i-1,n*(n+1)/2-1),s(i,n)=sigma j=1~n f(j,i)
再定义ans(i,n)=sigma j=1~pn q(pj)*f(i,j),则ans(0,n)即为所求
且ans(i,n)=s(i,n)*ans(0,pn) - ans(i+1,pn)
我们只需要用插值等方法把i=0~3时的f(i,n)和s(i,n)算出来(直接保留多项式的形式,以便计算s(i,n)),再对j=1~605预处理出ans(i,j),然后每次询问时递归下去即可。
还有一些常数优化技巧,比如可以事先打出多项式的系数然后直接在程序里打表(这样多项式的系数可以直接开ll),以及每次预处理出所有可能用到的pn和s(i,n)值等。
当然写起来巨麻烦无比,需要手写高精度和多项式类,两者需要实现的东西也很多,因此诞生了我有史以来的第一个1000+行的程序!
#include<bits/stdc++.h>
using namespace std;
#define gc getchar()
#define pc putchar
#define li long long
inline li ll_in(){
li x = 0,y = 0,c = getchar();
while(!isdigit(c)) y = c,c = getchar();
while(c >= '0' && c <= '9') x = (x << 1) + (x << 3) + (c ^ '0'),c = getchar();
return y == '-' ? -x : x;
}
inline void ll_out(li q){
if(q < 0){
putchar('-');
q = -q;
}
if(q >= 10) ll_out(q / 10);
putchar(q % 10 + '0');
}
/*****************************以下为高精度***************************/
const li MAX_PER_UNIT = 100000000;
struct gjd{
vector<li> a;//存储高精度数据
int len;//该数据的位数
bool fh;//符号,0为非负1为负
gjd(){//初始化(清零)
len = fh = 0;a.clear();
}
inline gjd operator = (li x){
a.clear();len = fh = 0;
if(x < 0){
fh = 1;x = -x;
}
while(x){
a.push_back(x % MAX_PER_UNIT);
x /= MAX_PER_UNIT;
}
len = a.size();
return *this;
}
inline gjd operator = (string c){
a.clear();len = fh = 0;
if(c.size() == 1 && c[0] == '0') return *this;
int q = MAX_PER_UNIT,j = 0,d = 0;
while(j < c.size() && (c[j] < '0' || c[j] > '9')) d = c[j++];
while(j < c.size() && c[j] == '0') ++j;
if(j == c.size()) return *this;
for(int i = c.size() - 1;i >= j;--i){
if(q == MAX_PER_UNIT){
a.push_back(c[i] - '0');
++len;
q = 10;
}
else{
a[len - 1] += (c[i] - '0') * q;
q *= 10;
}
}
if(len && d == '-') fh = 1;
return *this;
}
inline gjd operator = (char *c){
a.clear();len = fh = 0;
int tp = strlen(c);
if(tp == 1 && c[0] == '0') return *this;
int q = MAX_PER_UNIT,j = 0,d = 0;
while(j < tp && (c[j] < '0' || c[j] > '9')) d = c[j++];
while(j < tp && c[j] == '0') ++j;
if(j == tp) return *this;
for(int i = tp - 1;i >= j;--i){
if(q == MAX_PER_UNIT){
a.push_back(c[i] - '0');
++len;
q = 10;
}
else{
a[len - 1] += (c[i] - '0') * q;
q *= 10;
}
}
if(len && d == '-') fh = 1;
return *this;
}
inline void read(){//读入
a.clear();len = fh = 0;
vector<char> c;
char d,e = 0;
int i = 0;
d = getchar();
while(d < '0' || d > '9') e = d,d = getchar();//屏蔽非数字字符
while(d == '0') d = getchar();//屏蔽前导0
if(d == -1) return;//处理读入0
while(d >= '0' && d <= '9') ++i,c.push_back(d),d = getchar();
int q = MAX_PER_UNIT;
for(--i;i >= 0;--i){
if(q == MAX_PER_UNIT){
a.push_back(c[i] - '0');
++len;
q = 10;
}
else{
a[len - 1] += (c[i] - '0') * q;
q *= 10;
}
}
if(len && e == '-') fh = 1;
}
inline void print(){//输出
int i = len;
if(i == 0){
putchar('0');
return;
}
if(fh) putchar('-');
ll_out(a[i - 1]);
for(i -= 2;i >= 0;--i){
if(a[i] < 10000000) putchar('0');
if(a[i] < 1000000) putchar('0');
if(a[i] < 100000) putchar('0');
if(a[i] < 10000) putchar('0');
if(a[i] < 1000) putchar('0');
if(a[i] < 100) putchar('0');
if(a[i] < 10) putchar('0');
ll_out(a[i]);
}
}
};
inline void clear(gjd &x){//初始化(清零)
x.a.clear();
x.len = 0;
}
inline gjd turn(li x){//将低精转化为高精
gjd y;
if(x < 0){
y.fh = 1;x = -x;
}
while(x){
y.a.push_back(x % MAX_PER_UNIT);
x /= MAX_PER_UNIT;
}
y.len = y.a.size();
return y;
}
inline li turn(gjd z){//将高精转化为低精(有爆ll风险,慎用)
li x = 0;
for(int i = z.len - 1;i >= 0;--i) x = x * MAX_PER_UNIT + z.a[i];
return z.fh ? -x : x;
}
inline void swap(gjd &x,gjd &y){//交换两个高精度数据
gjd t = x;
x = y;
y = t;
}
inline gjd abs(gjd x){//高精度数据取绝对值
x.fh = 0;
return x;
}
inline gjd operator - (gjd x){//相反数
if(x.len) x.fh ^= 1;
return x;
}
inline int cmp(gjd x,gjd y){//高精与高精比较大小
if(x.fh && !y.fh) return -1;
if(!x.fh && y.fh) return 1;
int tmp = 1;
if(x.fh && y.fh) tmp = -1;
int q = x.len,w = y.len;
if(q > w) return tmp;
if(w > q) return -tmp;
for(--q;q >= 0;--q){
if(x.a[q] > y.a[q]) return tmp;
if(x.a[q] < y.a[q]) return -tmp;
}
return 0;
}
inline bool operator == (gjd b,gjd c){
return cmp(b,c) == 0;
}
inline bool operator < (gjd b,gjd c){
return cmp(b,c) == -1;
}
inline bool operator > (gjd b,gjd c){
return cmp(b,c) == 1;
}
inline bool operator <= (gjd b,gjd c){
return cmp(b,c) != 1;
}
inline bool operator >= (gjd b,gjd c){
return cmp(b,c) != -1;
}
inline bool operator != (gjd b,gjd c){
return cmp(b,c) != 0;
}
inline int cmp(gjd x,li b){//高精与低精比较大小
return cmp(x,turn(b));
}
inline bool operator == (gjd b,li c){
return cmp(b,c) == 0;
}
inline bool operator < (gjd b,li c){
return cmp(b,c) == -1;
}
inline bool operator > (gjd b,li c){
return cmp(b,c) == 1;
}
inline bool operator <= (gjd b,li c){
return cmp(b,c) != 1;
}
inline bool operator >= (gjd b,li c){
return cmp(b,c) != -1;
}
inline bool operator != (gjd b,li c){
return cmp(b,c) != 0;
}
inline int cmp(li b,gjd c){//低精与高精比较大小
return cmp(c,b);
}
inline bool operator == (li b,gjd c){
return cmp(b,c) == 0;
}
inline bool operator < (li b,gjd c){
return cmp(b,c) == -1;
}
inline bool operator > (li b,gjd c){
return cmp(b,c) == 1;
}
inline bool operator <= (li b,gjd c){
return cmp(b,c) != 1;
}
inline bool operator >= (li b,gjd c){
return cmp(b,c) != -1;
}
inline bool operator != (li b,gjd c){
return cmp(b,c) != 0;
}
inline gjd max(gjd x,gjd y){//高精与高精取最大值
return x > y ? x : y;
}
inline gjd max(gjd x,li y){//高精与低精取最大值
return x >= y ? x : turn(y);
}
inline gjd max(li x,gjd y){//低精与高精取最大值
return x > y ? turn(x) : y;
}
inline gjd min(gjd x,gjd y){//高精与高精取最小值
return x < y ? x : y;
}
inline gjd min(gjd x,li y){//高精与低精取最小值
return x <= y ? x : turn(y);
}
inline gjd min(li x,gjd y){//低精与高精取最小值
return x < y ? turn(x) : y;
}
inline gjd jia(gjd b,gjd c) {//高精加高精的绝对值相加,要求b,c>=0
if(c.len == 0) return b;
if(b.len == 0) return c;
int z = 0;
gjd d;
int l1 = b.len,l2 = c.len,t1,t2;
d.a.resize(max(l1,l2) + 1);
while(z < l1 || z < l2){
if(l1 > z) t1 = b.a[z];
else t1 = 0;
if(l2 > z) t2 = c.a[z];
else t2 = 0;
d.a[z] += t1 + t2;
if(d.a[z] >= MAX_PER_UNIT) d.a[z + 1] = 1,d.a[z] -= MAX_PER_UNIT;
++z;
}
d.len = d.a.size();
while(d.len && !d.a[d.len - 1]) --d.len;
return d;
}
inline gjd jia(gjd b,li c){//高精加低精的绝对值相加,要求b,c>=0
if(c == 0) return b;
if(b.len == 0) return turn(c);
if(c > 9000000000000000000ll) return jia(b,turn(c));//不然会爆ll
int z = 0;
b.a[0] += c;
while(b.a[z] >= MAX_PER_UNIT){
if(b.a.size() == z + 1) b.a.push_back(b.a[z] / MAX_PER_UNIT);
else b.a[z + 1] += b.a[z] / MAX_PER_UNIT;
b.a[z] %= MAX_PER_UNIT;
++z;
}
b.len = max(b.len,z + 1);
return b;
}
inline gjd jian(gjd z,li x){//高精减低精的绝对值相减,要求z>=x>=0
if(x == 0) return z;
int q = 0,w = -1;
z.a[0] -= x;
while(z.a[q] < 0){
z.a[q + 1] += z.a[q] / MAX_PER_UNIT;
z.a[q] %= MAX_PER_UNIT;
if(z.a[q] < 0){
z.a[q] += MAX_PER_UNIT;
--z.a[q + 1];
}
if(z.a[q]) w = q;
++q;
}
if(z.len == q + 1 && z.a[q] == 0) z.len = w + 1;
return z;
}
inline gjd jian(gjd b,gjd c){//高精减高精的绝对值相减,要求b>=c>=0
if(c.len == 0) return b;
int z = 0,w = -1;
int l1 = b.len,l2 = c.len;
gjd d;
d.a.resize(b.len);
int tp;
while(z < l1){
if(l2 > z) tp = c.a[z];
else tp = 0;
d.a[z] += b.a[z] - tp;
if(d.a[z] < 0){
d.a[z] += MAX_PER_UNIT;
d.a[z + 1] = -1;
}
if(d.a[z]) w = z;
++z;
}
d.len = w + 1;
return d;
}
inline gjd operator + (gjd z,li x){//高精加低精
if(!z.fh && x >= 0) return jia(z,x);
gjd a;
if(z.fh && x < 0){
z.fh = 0;x = -x;
a = jia(z,x);
a.fh = 1;
return a;
}
if(!z.fh && x < 0){
x = -x;
if(z >= x){
a = jian(z,x);
a.fh = 0;
}
else{
a = jian(turn(x),z);
a.fh = 1;
}
return a;
}
if(z.fh && x >= 0){
z.fh = 0;
if(z > x){
a = jian(z,x);
a.fh = 1;
}
else{
a = jian(turn(x),z);
a.fh = 0;
}
return a;
}
return a;
}
inline gjd operator + (gjd z,gjd x){//高精加高精
if(!z.fh && !x.fh) return jia(z,x);
if(z > x) swap(z,x);
gjd a;
if(z.fh && x.fh){
z.fh = x.fh = 0;
a = jia(z,x);
a.fh = 1;
return a;
}
if(z.fh && !x.fh){
z.fh = 0;
if(z > x){
a = jian(z,x);
a.fh = 1;
}
else{
a = jian(x,z);
a.fh = 0;
}
return a;
}
return a;
}
inline gjd operator + (li z,gjd x){//低精加高精
return x + z;
}
inline gjd operator += (gjd &z,li x){
z = z + x;
return z;
}
inline gjd operator += (gjd &z,gjd x){
z = z + x;
return z;
}
inline li operator += (li &z,gjd x){//请注意爆long long风险
z = z + turn(x);
return z;
}
inline gjd operator ++ (gjd &z){//由于特殊原因,使用时只能写作++z,不能写作z++
z = z + 1;
return z;
}
inline gjd operator - (gjd z,li x){//高精减低精
return z + (-x);
}
inline gjd operator - (gjd z,gjd x){//高精减高精
if(x.len) x.fh ^= 1;
return z + x;
}
inline gjd operator - (li b,gjd c){//低精减高精
return turn(b) - c;
}
inline gjd operator -= (gjd &z,li x){
z = z - x;
return z;
}
inline gjd operator -= (gjd &z,gjd x){
z = z - x;
return z;
}
inline li operator -= (li &z,gjd x){//请注意爆long long风险
z = z - turn(x);
return z;
}
inline gjd operator -- (gjd &z){//由于特殊原因,使用时只能写作--z,不能写作z--
z = z - 1;
return z;
}
inline gjd operator * (gjd,gjd);
inline gjd operator * (gjd b,li c){//高精乘低精
if(b.len == 0 || c == 0) return turn(0);
if(c > 10000000000ll) return b * turn(c);//不然会爆long long
int l = b.len,i,j;
j = b.fh ^ (c < 0);
b.fh = 0;c = abs(c);
for(i = 0;i < l;++i) b.a[i] *= c;
for(i = 0;i < l - 1;++i){
b.a[i + 1] += b.a[i] / MAX_PER_UNIT;
b.a[i] %= MAX_PER_UNIT;
}
while(b.a[i] >= MAX_PER_UNIT){
if(b.a.size() == i + 1) b.a.push_back(b.a[i] / MAX_PER_UNIT);
else b.a[i + 1] = b.a[i] / MAX_PER_UNIT;
b.a[i] %= MAX_PER_UNIT;
++i;
}
b.len = i + 1;
b.fh = j;
return b;
}
inline gjd operator * (gjd b,gjd c){//高精乘高精
if(b.len == 0 || c.len == 0) return turn(0);
if(min(b.len,c.len) <= 200){
gjd d;
d.fh = b.fh ^ c.fh;
b.fh = c.fh = 0;
int i,j;
int l1 = b.len,l2 = c.len;
d.a.resize(l1 + l2 - 1);
for(i = 0;i < l1;++i){
for(j = 0;j < l2;++j){
d.a[i + j] += b.a[i] * c.a[j];
}
}
for(i = 0;i < l1 + l2 - 2;++i){
d.a[i + 1] += d.a[i] / MAX_PER_UNIT;
d.a[i] %= MAX_PER_UNIT;
}
while(d.a[i] >= MAX_PER_UNIT){
if(d.a.size() == i + 1) d.a.push_back(d.a[i] / MAX_PER_UNIT);
else d.a[i + 1] = d.a[i] / MAX_PER_UNIT;
d.a[i] %= MAX_PER_UNIT;
++i;
}
d.len = i + 1;
return d;
}
int l = max(b.len,c.len) / 2;
gjd ans,b1,b2,c1,c2,s1,s2,s3,s4,s5;
int tmp = b.fh ^ c.fh;
b.fh = c.fh = 0;
int i;
b1.len = min(l,b.len);
b1.a.resize(b1.len);
for(i = 0;i < l && i < b.len;++i) b1.a[i] = b.a[i];
c1.len = min(l,c.len);
c1.a.resize(c1.len);
for(i = 0;i < l && i < c.len;++i) c1.a[i] = c.a[i];
if(b.len > l){
b2.a.resize(b.len - 1);
for(i = l;i < b.len;++i) b2.a[i - l] = b.a[i];
b2.len = b.len - l;
}
if(c.len > l){
c2.a.resize(c.len - 1);
for(i = l;i < c.len;++i) c2.a[i - l] = c.a[i];
c2.len = c.len - l;
}
while(b1.len && !b1.a[b1.len - 1]) --b1.len;
while(c1.len && !c1.a[c1.len - 1]) --c1.len;
while(b2.len && !b2.a[b2.len - 1]) --b2.len;
while(c2.len && !c2.a[c2.len - 1]) --c2.len;
s1 = b1 * c1;
s2 = b2 * c2;
s3 = (b1 + b2) * (c1 + c2) - s1 - s2;
if(s2.len){
s4.a.resize((l << 1) + s2.len);
for(i = 0;i < s2.len;++i) s4.a[i + (l << 1)] = s2.a[i];
s4.len = (l << 1) + s2.len;
}
if(s3.len){
s5.a.resize(l + s3.len);
for(i = 0;i < s3.len;++i) s5.a[i + l] = s3.a[i];
s5.len = l + s3.len;
}
ans = s1 + s4 + s5;ans.fh = tmp;
return ans;
}
inline gjd operator * (li z,gjd x){//低精乘高精
return x * z;
}
inline gjd operator *= (gjd &z,li x){
z = z * x;
return z;
}
inline gjd operator *= (gjd &z,gjd x){
z = z * x;
return z;
}
inline li operator *= (li &z,gjd x){//请注意爆long long风险
z = z * turn(x);
return z;
}
inline gjd operator / (gjd b,li c){//高精除以低精
if(!c){
puts("error:divided by zero!");
exit(1);
}
if(b.len == 0) return turn(0);
gjd q;
q.fh = b.fh ^ (c < 0);
b.fh = 0;c = abs(c);
int l = b.len - 1;
int x = 0,y;
while(l){
y = b.a[l] / c;
if(!x && y){
x = l;
q.a.resize(x);
q.a.push_back(y);
}
if(x) q.a[l] = y;
b.a[l - 1] += b.a[l] % c * MAX_PER_UNIT;
--l;
}
if(!x){
y = b.a[0] / c;
if(!y) return q;
q.a.push_back(y);
q.len = x + 1;
return q;
}
q.a[0] = b.a[0] / c;
if(!q.a[0] && !x) q.len = 0;
else q.len = x + 1;
return q;
}
inline gjd operator / (gjd b,gjd c){//高精除以高精
if(c.len == 0){
puts("error:divided by zero!");
exit(1);
}
if(b.len == 0) return turn(0);
int tmp = b.fh ^ c.fh;
b.fh = c.fh = 0;
if(b < c) return turn(0);
gjd q,w;//q:商;w:余数
int i = b.len - 1,l = -1,al,ar,am;
for(;i >= 0;i--){
w = w * MAX_PER_UNIT + b.a[i];
if(w >= c){
if(l == -1) q.a.resize(i + 1);
l = max(l,i);
al = 1,ar = MAX_PER_UNIT - 1;
while(al <= ar){
am = al + ar >> 1;
if(w >= c * am){
q.a[i] = am;
al = am + 1;
}
else ar = am - 1;
}
w -= c * q.a[i];
}
}
q.len = l + 1;
q.fh = tmp;
return q;
}
inline gjd operator / (li x,gjd y){//低精除以高精
return turn(x) / y;
}
inline gjd operator /= (gjd &z,li x){
z = z / x;
return z;
}
inline gjd operator /= (gjd &z,gjd x){
z = z / x;
return z;
}
inline li operator /= (li &z,gjd x){//请注意爆long long风险
if(z < x){
z = 0;
return z;
}
z = z / turn(x);
return z;
}
inline gjd operator % (gjd b,li c){//高精模低精
if(!c){
puts("error:divided by zero!");
exit(1);
}
if(b.len == 0) return turn(0);
int l = b.len - 1,tmp = b.fh;
b.fh = 0;c = abs(c);
while(l){
b.a[l - 1] += b.a[l] % c * MAX_PER_UNIT;
--l;
}
b.a[0] %= c;
gjd y = turn(b.a[0]);
y.fh = tmp;
return y;
}
inline gjd operator % (gjd b,gjd c){//高精模高精
if(c.len == 0){
puts("error:divided by zero!");
exit(1);
}
if(b.len == 0) return turn(0);
gjd w,tp;//w:余数 tp:辅助计算
int tmp = b.fh;
b.fh = c.fh = 0;
if(b < c){
b.fh = tmp;
return b;
}
int i = b.len - 1,j;
for(;i >= 0;--i){
w = w * MAX_PER_UNIT + b.a[i];
j = 1 << 26;
while(j){
tp = c * j;
if(w >= tp) w -= tp;
j >>= 1;
}
}
w.fh = tmp;
return w;
}
inline gjd operator % (li x,gjd y){//低精模高精
return turn(x) % y;
}
inline gjd operator %= (gjd &z,li x){
z = z % x;
return z;
}
inline gjd operator %= (gjd &z,gjd x){
z = z % x;
return z;
}
inline li operator %= (li &z,gjd x) {//请注意爆long long风险
if(z < x) return z;
z = z % turn(x);
return z;
}
inline gjd spw(gjd q,li w){//高精度数据的单精度次幂运算,结果为高精度,要求指数是非负整数
if(w < 0){
puts("error:in function 'spw(gjd q,long long int w)',w is negative!");
exit(1);
}
if(!w) return turn(1);
if(w == 1) return q;
gjd z;
z = spw(q,w >> 1);
z = z * z;
if(w & 1) z = z * q;
return z;
}
inline gjd spw(li q,li w){//两个单精度数据的幂运算,结果为高精度,要求指数是非负整数
return spw(turn(q),w);
}
inline gjd operator ^ (gjd q,li w){
return spw(q,w);
}
inline gjd operator ^= (gjd &q,li w){
q = q ^ w;
return q;
}
//以下为指数是高精度的幂运算,小心结果可能很大!
inline gjd spw(gjd q,gjd w){
return spw(q,turn(w));
}
inline gjd spw(li q,gjd w){
return spw(turn(q),turn(w));
}
inline gjd operator ^ (gjd q,gjd w){
return spw(q,turn(w));
}
inline gjd operator ^= (gjd &q,gjd w){
q = q ^ w;
return q;
}
inline gjd operator ^ (li q,gjd w){
return spw(turn(q),turn(w));
}
inline li operator ^= (li &q,gjd w){//小心爆ll的风险
q = turn(q ^ w);
return q;
}
inline gjd sqrt(gjd b){//高精开平方(牛顿迭代法)
if(b.fh){
puts("error:in function 'sqrt(gjd b)',b is negative!");
exit(1);
}
if(b.len == 0) return turn(0);
gjd z[2],tmp;
bool i = 0;
z[1].a.resize((b.len + 1) / 2 - 1);
z[1].a.push_back(1);
z[1].len = (b.len + 1) / 2;
while(z[i] != z[!i]){
tmp = (b / z[!i] + z[!i]) / 2;
if(z[!i] > z[i] && z[!i] - z[i] == 1 && z[i] == tmp) break;
z[i] = tmp;
i = !i;
}
return z[i];
}
inline gjd sqr(gjd a,li b){//二分法高精开b次方根,要求次数为正整数
if(b <= 0){
puts("error:in function 'sqr(gjd a,li b)',b is not positive!");
exit(1);
}
if(a.fh && b % 2 == 0){
puts("error:in function 'sqr(gjd a,li b)',a is negative and b is even!");
exit(1);
}
if(a.len == 0) return turn(0);
if(b == 1) return a;
int tmp = a.fh;
a.fh = 0;
if(a <= b) return turn(1 * tmp);
gjd l,r,mid,ans;
l.a.resize((a.len - 1) / b);r.a.resize((a.len - 1) / b + 1);
l.a.push_back(1);l.len = (a.len - 1) / b + 1;r.a.push_back(1);r.len = l.len + 1;
if(b == 2){
while(l <= r){
mid = (l + r) / 2;
if(mid * mid <= a){
ans = mid;
l = mid + 1;
}
else r = mid - 1;
}
}
else{
while(l <= r){
mid = (l + r) / 2;
if(spw(mid,b) <= a){
ans = mid;
l = mid + 1;
}
else r = mid - 1;
}
}
ans.fh = tmp;
return ans;
}
inline gjd gcd(gjd x,gjd y){//高精度与高精度的最大公因数,返回值非负(下同)
x.fh = y.fh = 0;
if(!x.len) return y;
if(!y.len) return x;
gjd q = turn(1);
if(x < y) swap(x,y);
while(x != y){
if((x.a[0] & 1) && (y.a[0] & 1)) x -= y;
else if(x.a[0] & 1) y /= 2;
else if(y.a[0] & 1) x /= 2;
else{
x /= 2;
y /= 2;
q *= 2;
}
if(x < y) swap(x,y);
}
return q * x;
}
inline gjd gcd(gjd x,li y){//高精度与单精度的最大公因数
x.fh = 0;y = abs(y);
return y == 0 ? x : gcd(turn(y),x % y);
}
inline gjd gcd(li x,gjd y){//单精度与高精度的最大公因数
return gcd(y,x);
}
inline gjd gjd_gcd(li x,li y){//单精度与单精度的最大公因数,结果返回高精
x = abs(x);y = abs(y);
return y == 0 ? turn(x) : gjd_gcd(y,x % y);
}
inline li gcd(li x,li y){//单精度与单精度的最大公因数,结果返回单精
x = abs(x);y = abs(y);
return y == 0 ? x : gcd(y,x % y);
}
inline gjd lcm(gjd x,gjd y){//高精度与高精度的最小公倍数
x.fh = y.fh = 0;
if(!x.len || !y.len) return turn(0);
return x / gcd(x,y) * y;
}
inline gjd lcm(gjd x,li y){//高精度与单精度的最小公倍数
x.fh = 0;y = abs(y);
if(!x.len || !y) return turn(0);
return x / gcd(x,y) * y;
}
inline gjd lcm(li x,gjd y){//单精度与高精度的最小公倍数
return lcm(y,x);
}
inline gjd gjd_lcm(li x,li y){//单精度与单精度的最小公倍数,结果返回高精
x = abs(x);y = abs(y);
if(!x || !y) return turn(0);
return x / gjd_gcd(x,y) * y;
}
inline li lcm(li x,li y){//单精度与单精度的最小公倍数,结果返回单精
x = abs(x);y = abs(y);
if(!x || !y) return 0;
return x / gcd(x,y) * y;
}
/*****************************以上为高精度***************************/
/*****************************以下为多项式***************************/
const int sz = 17;
struct dxs{
gjd a[sz];
gjd div;
int len;
dxs(){
for(int i = 0;i < sz;++i) clear(a[i]);
div = 1;len = 0;
}
void clr(){
for(int i = 0;i < sz;++i) clear(a[i]);
div = 1;len = 0;
}
void out(){
ll_out(len);pc('\n');
for(int i = 0;i <= len;++i){
a[i].print();pc(' ');
}
pc('\n');div.print();pc('\n');
}
}a[6],s[6];
void yf(dxs &q){
if(q.div.fh){
q.div.fh ^= 1;
for(int i = 0;i <= q.len;++i) if(q.a[i].len) q.a[i].fh ^= 1;
}
gjd g = q.div;
for(int i = 0;i <= q.len;++i) g = gcd(g,q.a[i]);
q.div /= g;
for(int i = 0;i <= q.len;++i) q.a[i] /= g;
while(q.len && q.a[q.len] == 0) --q.len;
}
dxs operator + (dxs q,dxs w){
gjd a = lcm(q.div,w.div);
gjd b = a / q.div,c = a / w.div;
dxs ans;
ans.len = max(q.len,w.len);
ans.div = a;
for(int i = 0;i <= ans.len;++i) ans.a[i] = q.a[i] * b + w.a[i] * c;
yf(ans);
return ans;
}
dxs operator * (dxs q,dxs w){
dxs ans;
ans.len = q.len + w.len;
ans.div = q.div * w.div;
for(int i = 0;i <= q.len;++i){
for(int j = 0;j <= w.len;++j){
ans.a[i + j] += q.a[i] * w.a[j];
}
}
yf(ans);
return ans;
}
/*****************************以上为多项式***************************/
/*****************************以下为预处理部分***************************/
gjd tmp[10010],val[1010][6];
int p1[1010];
inline gjd getp(gjd q){
gjd ans = sqr(q * 2,2);
if((ans + 1) * ans > q * 2) --ans;
return ans;
}
inline gjd getx(dxs q,gjd x){
gjd ans = q.a[q.len];
for(int i = q.len - 1;i >= 0;--i) ans = ans * x + q.a[i];
return ans / q.div;
}
dxs tp2[35];
dxs cheng(int l,int r){
if(l == r) return tp2[l];
int mid = l + r >> 1;
return cheng(l,mid) * cheng(mid + 1,r);
}
inline void wk(dxs q,dxs &w,bool fg,int mx){
int i,j = 1;
for(i = 1;i <= mx;++i){
tmp[i] = tmp[i - 1];
for(;j <= (!fg ? i * (i + 1) / 2 - 1 : i);++j) tmp[i] += getx(q,turn(j));
}
w.clr();
for(i = 1;i <= mx;++i){
dxs tp;
gjd tp3;tp3 = 1;
tp.a[0] = 1;
for(j = 1;j <= mx;++j) if(i != j){
tp2[j].a[1] = 1;tp2[j].a[0] = -j;tp2[j].len = 1;
tp3 *= i - j;
}
tp2[i].clr();tp2[i].a[0] = 1;
tp = cheng(1,mx);
tp.div = tp3;
for(j = 0;j <= tp.len;++j) tp.a[j] *= tmp[i];
yf(tp);
w = w + tp;
}
}
const int mx = 605;
void init(){
register int i,j;
p1[1] = 1;
for(i = 2;i <= mx;++i) p1[i] = turn(getp(turn(i)));
a[0].a[0] = 1;
s[0].a[1] = 1;
s[0].len = 1;
val[1][0] = 1;
for(i = 2;i <= mx;++i) val[i][0] = val[i - 1][0] + val[p1[i]][0];
for(i = 1;i <= 3;++i){
wk(a[i - 1],a[i],0,pow(2,i + 1) - 1);
wk(a[i],s[i],1,pow(2,i + 1));
val[1][i] = getx(a[i],turn(1));
for(j = 2;j <= mx;++j) val[j][i] = val[j - 1][i] + val[p1[j]][0] * getx(a[i],turn(j));
}
//for(i = 0;i <= 3;++i) a[i].out();
//for(i = 0;i <= 3;++i) s[i].out();
for(i = 1;i <= mx;++i) val[i][4] = val[i - 1][4] + val[p1[i]][0] * getx(s[3],turn(i * (i + 1) / 2 - 1));
}
/*****************************以上为预处理部分***************************/
gjd pp[6],ss[6][6];
gjd work(int q,int k){
if(pp[q] <= mx) return val[turn(pp[q])][k];
return ss[q][k] * work(q + 1,0) - work(q + 1,k + 1);
}
inline void file(){
freopen("a.in","r",stdin);
freopen("a.out","w",stdout);
}
int main(){
//file();
init();
int t = ll_in();
while(t--){
pp[0].read();
ss[0][0] = pp[0];
for(int i = 1;i <= 4;++i){
pp[i] = getp(pp[i - 1]);
if(i != 4){
ss[i][0] = pp[i];
for(int j = 1;j <= i;++j) ss[i][j] = getx(s[j],pp[i]);
}
}
work(0,0).print();
pc('\n');
}
return 0;
}