我们学到了求格点三角形的皮克公式：
(S=a+frac{b}{2}+1)（(a)为格点三角形内部点数，(b)为格点三角形边上点数）
嗯，然后不会了
但是我找到了(code) ，如果有好心人做出来了请教我(qaq)
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define mid ((l+r)>>1)
#define lowbit(i) ((i)&(-i))
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=3010,mod=1004535809,inf=0x3f3f3f3f;
int n,m,ret;
int mu[N];
inline int add(int x,int y)
{
return x+y>=mod?x+y-mod:x+y;
}
inline void dec(int &x,int y)
{
if((x-=y)<0) x+=mod;
}
inline int inv(int x,int k=mod-2)
{
int ret=1;
while(k)
{
if(k&1) ret=ret*x%mod;
x=x*x%mod;
k>>=1;
}
return ret;
}
inline void main()
{
n=read(),m=read();
if(n>m) swap(n,m);
for(int i=1;i<=n;++i) mu[i]=i*i;
for(int i=1;i<=n;++i)
{
for(int j=(i<<1);j<=n;j+=i)
{
mu[j]-=mu[i];
}
}
ret=n*(n-1)%mod*m%mod*(m-1)%mod;
ret=ret*add(11ll*n*(n+1)%mod*m%mod*(m+1)%mod,6ll*add(1ll*n*(n+1)%mod,1ll*m*(m+1)%mod)%mod)%mod*inv(144)%mod;
for(int i=1;i<=n;++i) dec(ret,1ll*mu[i]*(n-i+n%i)%mod*(n/i)%mod*inv(2)%mod*(m-i+m%i)%mod*(m/i)%mod*inv(2)%mod);
printf("%lld
",ret*inv(3)%mod);
}
}
signed main()
{
red::main();
return 0;
}
设(m=sumlimits_{i=1}^{n}a_i)
首先总方案数是一个可重集排列
[frac{m!}{prodlimits_{i=1}^{n}a_i}
]
满足要求的排列方式是一个第(i)行有(a_i)个元素的杨表
根据钩子公式，满足条件的方案数为
[frac{m!}{prod h(i,j)}
]
其中(h(i,j))为每个格子的钩子长度(+1)
考虑下怎么快速求(prod h(i,j))
对于第(i)行的第一列，钩长为(a_i+n-i)
那么第(i)行所有格子的乘积为(frac{(a_i+n-i)!}{prodlimits_{j>i}(j-i)+(a_i-a_j)})
（手玩一行就比较好理解了）
得到
[prod h(i,j)=frac{prodlimits_{i=1}^{n}(a_i+n-i)!}{prodlimits_{1le i<jle n}(a_i-i)-(a_j-j)}
]
合并一下得到
[ret=frac{frac{m!}{prod h(i,j)}}{frac{m!}{prodlimits_{i=1}^{n}a_i}}
]
[=frac{frac{m!prodlimits_{1le i<jle n}(a_i-i)-(a_j-j)}{prodlimits_{i=1}^{n}(a_i+n-i)}}{frac{m!}{prodlimits_{i=1}^{n}a_i}}
]
[=frac{prodlimits_{i=1}^{n}a_i}{prodlimits_{i=1}^{n}(a_i+n-i)}×prodlimits_{1le i<jle n}(a_i-i)-(a_j-j)
]
左边这个可以(O(n))求出来，右边暴力算需要(n^2)
根据(c_k=sumlimits_{i=0}^{k}f_i*g_{k-i})
我们可以设多项式(f_k=sumlimits_{i=1}^{n}[a_i-i==k],g_k=sumlimits_{i=1}^{n}[-(a_i-i)==k])
然后求(c_k)，那么右边就等于(prod i^{c_i})
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define lowbit(i) ((i)&(-i))
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=4e6+10,p=1004535809,G=3,Gi=334845270;
int n,tn,tm,limit=1,len,ans=1;
int a[N],f[N],g[N],h[N],w[N],pos[N],fac[N],inv[N];
inline int fast(int x,int k)
{
int ret=1;
while(k)
{
if(k&1) ret=ret*x%p;
x=x*x%p;
k>>=1;
}
return ret;
}
inline void setup(int *a,int &l,int n)
{
++a[n],l=max(l,n);
}
inline void init_fac(int n)
{
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<=n;++i) fac[i]=fac[i-1]*i%p;
for(int i=2;i<=n;++i) inv[i]=(p-p/i)*inv[p%i]%p;
for(int i=2;i<=n;++i) inv[i]=inv[i-1]*inv[i]%p;
}
inline void ntt(int *a,int inv)
{
for(int i=0;i<limit;++i)
if(i<pos[i]) swap(a[i],a[pos[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
int Wn=fast(inv?G:Gi,(p-1)/(mid<<1));
for(int r=mid<<1,j=0;j<limit;j+=r)
{
int w=1;
for(int k=0;k<mid;++k,w=w*Wn%p)
{
int x=a[j+k],y=w*a[j+k+mid]%p;
a[j+k]=(x+y)%p;
a[j+k+mid]=(x-y)%p;
if(a[j+k+mid]<0) a[j+k+mid]+=p;
}
}
}
}
inline void main()
{
n=read();
for(int i=1;i<=n;++i) a[i]=read();
init_fac(n+a[1]);
for(int i=1;i<=n;++i)
{
ans=ans*(fac[a[i]]*inv[a[i]+n-i]%p)%p;
setup(f,tn,a[i]-i+n);
setup(g,tm,-a[i]+i+a[1]);
}
while(limit<tn+tm+1) limit<<=1,++len;
for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
ntt(f,1),ntt(g,1);
for(int i=0;i<limit;++i) f[i]=f[i]*g[i]%p;
ntt(f,0);
int invs=fast(limit,p-2);
for(int i=0;i<limit;++i) f[i]=f[i]*invs%p;
for(int i=1,d=n+a[1];i+d<limit;++i)
{
ans=ans*fast(i,f[i+d])%p;
}
printf("%lld
",ans);
}
}
signed main()
{
red::main();
return 0;
}
求((a+bi)(c+di)=x),则
[left{
egin{array}{**lr**}
ac-bd=x \
ad+bc=0 \
end{array}
ight.]
可得到(frac{a}{b}=-frac{c}{d})
设(gcd(p,q)==1)，且(frac{a}{b}=frac{p}{q})
那么有
[x=s(p+qi)·t(p-qi)
]
[x=s·t·(p^2+q^2)
]
一个数对答案有贡献，需要满足既是((p^2+q^2))的倍数，又是(x)的约数
所以一个((p^2+q^2))的贡献为
[sumlimits_{w|frac{x}{p^2+q^2}}(p·w)=p*σ(frac{x}{p^2+q^2})
]
枚举(y=p^2+q^2)
[ans=sumlimits_{y}sumlimits_{gcd(p,q)==1&(p^2+q^2==y)}psumlimits_{y|x}σ(k)
]
设(D(i)=sumlimits_{j=1}^{i}σ(j)，f(i)=sumlimits_{y}sumlimits_{gcd(p,q)==1&(p^2+q^2==y)}p,F(i)=sumlimits_{j=1}^{i}f(i))
[ans=sumlimits_{y}f(y)·D(lfloorfrac{n}{y}
floor)
]
其中
[D(n)=sumlimits_{i=1}^{n}i*lfloorfrac{n}{i}
floor
]
求(F)考虑杜教筛
设(G(n)=sumlimits_{p^2+q^2le n}p=sumlimits_{p=1}^{sqrt{n}}p·lfloorsqrt{n-p^2} floor)
枚举(gcd(p,q)==d)有
[G(n)=sumlimits_{dge 1}d·F(lfloorfrac{n}{d^2}
floor)
]
[F(n)=G(n)-sumlimits_{dge 2}d·F(lfloorfrac{n}{d^2}
floor)
]
上面讨论的是(b>0&d>0)的情况，(b<0&d<0)时同理
当(b==0& d==0)时
[sumlimits_{a}sumlimits_{c}[acle n]=D(n)
]
所以最终答案为
[2*ans+D(n)
]
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define LL long long
#define CI const int&
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define lowbit(i) ((i)&(-i))
inline LL read()
{
LL x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=5e6+10,mod=1004535809,inv2=502267905;
LL n,ans;
bool vis[N];
LL prime[N],ds[N],fs[N],num;
map<LL,int> _ds,_fs;
inline int gcd(int a,int b)
{
return b?gcd(b,a%b):a;
}
inline void inc(LL& x,const LL y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void inc(int& x,CI y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void dec(int& x,CI y)
{
if ((x-=y)<0) x+=mod;
}
inline int sum(CI x,CI y)
{
int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
int t=x-y; return t<0?t+mod:t;
}
inline int Sum(const LL& l,const LL& r)
{
return (l+r)%mod*((r-l+1)%mod)%mod*inv2%mod;
}
inline void oula(const int n=5e6)
{
vis[1]=ds[1]=1;
for(int i=2;i<=n;++i)
{
if(!vis[i]) prime[++num]=i,ds[i]=i+1;
for(int j=1;j<=num&&i*prime[j]<=n;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
ds[i*prime[j]]=ds[i]*(prime[j]+1)-prime[j]*ds[i/prime[j]];
break;
}
ds[i*prime[j]]=ds[i]*(prime[j]+1);
}
}
for(int i=1;i*i<=n;++i)
{
int t=i*i;
for(int j=1;j*j+t<=n;++j)
if(gcd(i,j)==1) fs[j*j+t]+=i;
}
for(int i=1;i<=n;++i)
{
ds[i]=(ds[i]%mod+ds[i-1])%mod;
fs[i]=(fs[i]%mod+fs[i-1])%mod;
}
}
inline int Ds(const LL &x)
{
if(x<=5e6) return ds[x];
if(_ds.count(x)) return _ds[x];
int ret=0;
for(LL l=1,r;l<=x;l=r+1)
r=x/(x/l),inc(ret,1LL*Sum(l,r)*((x/r)%mod)%mod);
return _ds[x]=ret;
}
inline int Fs(const LL &x)
{
if(x<=5e6) return fs[x];
if(_fs.count(x)) return _fs[x];
int ret=0;
for(LL i=1;i*i<=x;++i)
inc(ret,i*((LL)floor(sqrt(x-i*i))%mod)%mod);
for(LL i=2;i*i<=x;++i)
dec(ret,i*Fs(x/(i*i))%mod);
return _fs[x]=ret;
}
inline void main()
{
n=read();oula();
for(LL l=1,r;l<=n;l=r+1)
r=n/(n/l),inc(ans,1LL*sub(Fs(r),Fs(l-1))*Ds(n/l)%mod);
printf("%d
",sum(sum(ans,ans),Ds(n)));
}
}
signed main()
{
red::main();
return 0;
}