模积和 洛谷2260 bzoj2956 整除分块、逆元

2 篇文章 0 订阅
1 篇文章 0 订阅

题意:

i=1nj=1m(n mod i)(m mod j),ij ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) , i ≠ j

mod 19940417 m o d   19940417 的值
首先我们直接求是 O(nm) O ( n ∗ m ) 的,那么我们考虑推一下这个式子。我们发现如果没有 ij i ≠ j 的限制,那么我们只需要分别求出 ni=1n mod i ∑ i = 1 n n   m o d   i mj=1m mod j ∑ j = 1 m m   m o d   j ,因为 n mod i=nnii n   m o d   i = n − ⌊ n i ⌋ ∗ i ,所以可以整数分块来算每一个的结果,可以用然后做乘积即可。但是有了这个限制我们应该怎么办呢?我们考虑用包含 i=j i = j 的总和减去 i=j i = j 的情况。
i=1nj=1m(n mod i)(m mod j),ij ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) , i ≠ j
i=1nj=1m(n mod i)(m mod j)i=1min(n,m)(n mod i)(m mod i) ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) − ∑ i = 1 m i n ( n , m ) ( n   m o d   i ) ∗ ( m   m o d   i )
i=1n(n mod i)j=1m(m mod j)i=1min(n,m)(nnii)(mmii) ∑ i = 1 n ( n   m o d   i ) ∗ ∑ j = 1 m ( m   m o d   j ) − ∑ i = 1 m i n ( n , m ) ( n − ⌊ n i ⌋ ∗ i ) ∗ ( m − ⌊ m i ⌋ ∗ i )
i=1n(nnii)j=1m(mmii)i=1min(n,m)nmmniinmii+nimii2 ∑ i = 1 n ( n − ⌊ n i ⌋ ∗ i ) ∗ ∑ j = 1 m ( m − ⌊ m i ⌋ ∗ i ) − ∑ i = 1 m i n ( n , m ) n ∗ m − m ∗ ⌊ n i ⌋ ∗ i − n ∗ ⌊ m i ⌋ ∗ i + ⌊ n i ⌋ ∗ ⌊ m i ⌋ ∗ i 2

好了,推导完毕,我们发现这个式子的前半部分已经讲过可以数论分块来用 O(n+m) O ( n + m ) 的时间求了,后半部分同理也可以。
这里需要补充一个公式:
i=1ni2=16n(n+1)(2n+1) ∑ i = 1 n i 2 = 1 6 n ( n + 1 ) ( 2 n + 1 )

这里我就不给出证明了,有兴趣的话可以自行搜索。
最后吐槽一句,这种取模是真的烦,式子本来就长,再加上不停地取模操作,还不能忘了 (x+mod)%mod ( x + m o d ) % m o d ,不然还会出现负数,真的很讨厌。
代码:

#include <bits/stdc++.h>
using namespace std;

long long n,m,ans,mod=19940417,ans1,ans2,ni;
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
    }
    else
    {
        exgcd(b,a%b,y,x);
        y-=a/b*x;
    }
}
long long calc(long long x)
{
    return (x*(x+1)%mod*(2*x+1)%mod*ni%mod)%mod;
}
int main()
{
    scanf("%lld%lld",&n,&m);
    long long r;
    long long x,y;
    exgcd(6,mod,x,y);
    x=(x%mod+mod)%mod;
    if(!x)
    x+=mod;
    ni=x;
    ans1=(n*n)%mod;
    ans2=(m*m)%mod;
    for(long long l=1;l<=n;l=r+1)
    {
        if(n/l!=0)
        r=min(n/(n/l),n);
        else
        r=n;
        ans1=(ans1-(((r-l+1)*(l+r)/2%mod*(n/l))%mod+mod)%mod)%mod;
        ans1=(ans1+mod)%mod;
    }
    for(long long l=1;l<=m;l=r+1)
    {
        if(m/l!=0)
        r=min(m/(m/l),m);
        else
        r=m;
        ans2=(ans2-(((r-l+1)*(l+r)/2%mod*(m/l))%mod+mod)%mod)%mod;
        ans2=(ans2+mod)%mod;
    }
    ans=(ans1*ans2)%mod;
    for(int l=1;l<=min(n,m);l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        ans=(ans-((r-l+1)*((m*n)%mod)%mod+(calc(r)-calc(l-1))%mod*((n/l)*(m/l)%mod)%mod-((r-l+1)*(l+r)/2)%mod*m%mod*(n/l)%mod-((r-l+1)*(l+r)/2)%mod*n%mod*(m/l)%mod))%mod;
        ans=(ans+mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}
  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值