洛谷 P3636 曲面
题目背景
xht喜欢研究数学函数，他特别喜欢反比例函数。
题目描述
我们知道，反比例函数xy=a的图象是双曲线。
xht于是想：把它推广到三维是什么样的呢？
定义曲面C(k)为方程xyz=k所确定的曲面。
又定义曲面的美观程度P(k)为曲面C(k)上所有整点（x,y,z坐标均为整数）到原点的曼哈顿距离的平方之和。
（点（x,y,z）到原点的曼哈顿距离为|x|+|y|+|z|）。
现在，xht把一些曲面{C(a),C(a+1)...C(b)}排成一列，你要求出它们美观程度之和对10007取模的结果。
输入输出格式
输入格式：
一行两个正整数数a，b
输出格式：
一行一个数
输入输出样例
说明
样例1的解释：
在曲面xyz=3上共有12个整点(1,1,3),(1,3,1),(3,1,1),(-1-1,3),(-1,-3,1),(-3,-1,1),(1,-1,-3),(1,-3,-1),(3,-1,-1),(-1,1,-3),(-1,3,-1),(-3,1,-1)。它们到原点的曼哈顿距离的平方之和为5^2*12=300。
对于20%的数据，a=b<=100
对于另外40%的数据，a,b<=3*10^5
对于100%的数据，1<=a,b<=3*10^8
记得是洛谷的某次上古比赛的题，，当时too young too naive的我只骗了60分hhhh，现在一看貌似也不是那么难了。
但这确实是一道好题。
首先可以只求坐标都是正的的和，然后*4,因为三维乘积要是正的话，必须满足坐标是负数的维数是偶数，所以C(3,0)+C(3,2)=4。
设solve(x)=ΣΣΣ(i+j+k)^2 ,其中i*j*k<=x。
那么答案显然就是（solve(b)-solve(a-1)）*4
求solve(x)的时候，我们把括号拆开，转化成 i^2 + 2*i * (j+k) + j^2+k^2+2*j*k。 我们固定i的时候左式就可以拆成三部分求。
考虑到x/i一定的时候，j*k的上界是一定的，而x/i的取值不是很多，而且x/i相等的i都是连续的，这就告诉我们可以数列分块。
内层求j和k的相关值时再套一层分块即可。
时间复杂度 O(玄学)
#include<bits/stdc++.h>
#define ll long long
#define ha 10007
using namespace std;
struct node{
int a,b,c;
//a表示i^2要乘的系数
//b表示2*i要乘的系数
//c表示∑j^2+k^2+2*k*j
};
int l,r;
/*
inline int add(int x,int y){
x+=y;
while(x>=ha) x-=ha;
return x;
}
*/
inline int c1(int x){
//一次前缀和
return (x*(ll)(x+1)>>1ll)%ha;
}
inline int c2(int x){
//二次前缀和
ll now=x*(ll)(x+1)>>1;
if(!(now%3)) return now/3%ha*(ll)(2*x+1)%ha;
else return (ll)(2*x+1)/3*(now%ha)%ha;
}
inline node work(int x){
node an=(node){0,0,0};
for(int i=1,j;i<=x;i=j+1){
int now=x/i;
j=x/now;
an.a=((ll)an.a+(j-i+1)*(ll)now)%ha;
an.b=((ll)an.b+(ll)(j-i+1)*(ll)c1(now)+(ll)now*(ll)(c1(j)-c1(i-1)+ha))%ha;
an.c=((ll)an.c+(ll)(c1(now)*(ll)(c1(j)-c1(i-1)+ha)<<1ll)+(ll)(c2(j)-c2(i-1)+ha)*(ll)now+(ll)c2(now)*(ll)(j-i+1))%ha;
}
return an;
}
inline int solve(int x){
int an=0;
for(int i=1,j;i<=x;i=j+1){
int now=x/i;
j=x/now;
node w=work(now);
an=((ll)an+(ll)(j-i+1)*(ll)w.c+(ll)((c1(j)-c1(i-1)+ha)<<1ll)*(ll)w.b+(ll)(c2(j)-c2(i-1)+ha)*(ll)w.a)%ha;
}
return an;
}
int main(){
scanf("%d%d",&l,&r);
int ans=solve(r)-solve(l-1);
if(ans<0) ans+=ha;
ans<<=2;
while(ans>=ha) ans-=ha;
printf("%d\n",ans);
return 0;
}
