bzoj2154||洛谷P1829 Crash的数字表格&&JZPTAB && bzoj3309 DZY Loves Math

bzoj2154||洛谷P1829

https://www.lydsy.com/JudgeOnline/problem.php?id=2154

https://www.luogu.org/problemnew/show/P1829

不妨设n<=m

就是求$ans=\sum_{k=1}^m{\frac{1}{k}}\sum_{i=1}^n\sum_{j=1}^m{ij[(i,j)=k]}$

把1/k后面的那一部分提出来，设为f(k)，

然后莫比乌斯反演得到f(k)较简易的计算式，代回ans，并化简（过程略）

最后化简出来是$\sum_{k=1}^m{k}\sum_{i=1}^{{\lfloor}\frac{m}{k}{\rfloor}}\mu(i){i^2}g({{\lfloor}\frac{{\lfloor}\frac{n}{k}{\rfloor}}{i}{\rfloor}},{{\lfloor}\frac{{\lfloor}\frac{m}{k}{\rfloor}}{i}{\rfloor}})$

其中$g(x,y)=\frac{(x+1)x(y+1)y}{4}$

可以通过两重的整除分块以O(n)的复杂度计算（第一重枚举n/k,m/k，第二重枚举(n/k)/i,(m/k)/i）

 #pragma GCC optimize(3)

 #include<cstdio>

 #include<algorithm>

 #include<cstring>

 #include<vector>

 using namespace std;

 #define fi first

 #define se second

 #define mp make_pair

 #define pb push_back

 typedef long long ll;

 typedef unsigned long long ull;

 typedef pair<int,int> pii;

 #define md 20101009

 #define N 10000000

 ll prime[N+],len,mu[N+],dd[N+];

 bool nprime[N+];

 ll Mod(ll a,ll b)

 {

     if(a>=)    return a%b;

     else if(a%b==)    return ;

     else    return b+a%b;

 }

 ll G(ll x,ll y)

 {

     return (x+)*x%md*(y+)%md*y%md*%md;

 }

 ll calc(ll n,ll m)

 {

     ll i,j,ans=;

     for(i=;i<=m;i=j+)

     {

         if(n<i)    j=min(m,m/(m/i));

         else    j=min(m,min(n/(n/i),m/(m/i)));

         ans=Mod(ans+(dd[j]-dd[i-])*G(n/i,m/i),md);

     }

     return ans;

 }

 ll n,m;

 int main()

 {

     ll i,j,ans=;

     mu[]=;

     for(i=;i<=N;i++)

     {

         if(!nprime[i])    prime[++len]=i,mu[i]=-;

         for(j=;j<=len&&i*prime[j]<=N;j++)

         {

             nprime[i*prime[j]]=;

             if(i%prime[j]==)    {mu[i*prime[j]]=;break;}

             else    mu[i*prime[j]]=-mu[i];

         }

     }

     for(i=;i<=N;i++)    dd[i]=dd[i-]+i*i*mu[i];

     scanf("%lld%lld",&n,&m);

     //n=10000000;m=1000;

     if(n>m)    swap(n,m);

     for(i=;i<=m;i=j+)

     {

         if(n<i)    j=min(m,m/(m/i));

         else    j=min(m,min(n/(n/i),m/(m/i)));

         ans=(ans+(i+j)*(j-i+)/*calc(n/i,m/i))%md;

     }

     printf("%lld",ans);

     return ;

 }

额，好像还有一道题面基本一样的题（然而bzoj权限题，不能交），然而是很多组数据，不能每次O(n)算

对于这个式子$\sum_{k=1}^m{k}\sum_{i=1}^{{\lfloor}\frac{m}{k}{\rfloor}}\mu(i){i^2}g({{\lfloor}\frac{n}{ik}{\rfloor}},{{\lfloor}\frac{m}{ik}{\rfloor}})$

令T=ik

原式=$\sum_{k=1}^m{k}\sum_{k|T}\mu(\frac{T}{k})(\frac{T}{k})^2g({{\lfloor}\frac{n}{T}{\rfloor}},{{\lfloor}\frac{m}{T}{\rfloor}})$

=$\sum_{T=1}^mg({{\lfloor}\frac{n}{T}{\rfloor}},{{\lfloor}\frac{m}{T}{\rfloor}})\sum_{k|T}k\mu(\frac{T}{k})(\frac{T}{k})^2$

最后从sigma k|T开始那一部分，等于$(id*(\mu\cdot id\cdot id))(T)$

是个积性函数，且是可以线性筛出来的（然而我不会，又去查了题解），将其设为函数h

对于1，h[1]=1

对于一个质数p，h[p]=$1*\mu(p)*p^2+p*\mu(1)*1^2$

对于i%p!=0，h(i*p)=h(i)*h(p)（由于积性）

对于i%p==0，由于h(i*p)拆开后式子中比h(i)多的项的$\mu$值均为0，因此只有对于h(i)中已有的项增加贡献；每一项由$\mu(d)*d^2*(i/d)$变到$\mu(d)*d^2*((i*p)/d)$，因此h(i*p)=h(i)*p

筛出来之后做一下前缀和，这样子每一次询问就可以数论分块根号解决了

 #include<cstdio>

 #include<algorithm>

 #include<cstring>

 #include<vector>

 using namespace std;

 #define fi first

 #define se second

 #define mp make_pair

 #define pb push_back

 typedef long long ll;

 typedef unsigned long long ull;

 typedef pair<int,int> pii;

 #define md 20101009

 #define N 10000000

 ll prime[N+],len,dd[N+];

 bool nprime[N+];

 ll h[N+];

 ll Mod(ll a,ll b)

 {

     if(a>=)    return a%b;

     else if(a%b==)    return ;

     else    return b+a%b;

 }

 ll G(ll x,ll y)

 {

     return (x+)*x%md*(y+)%md*y%md*%md;

 }

 ll n,m;

 int main()

 {

     ll i,j,ans=;

     h[]=;

     for(i=;i<=N;i++)

     {

         if(!nprime[i])    prime[++len]=i,h[i]=Mod(-i*i+i,md);

         for(j=;j<=len&&i*prime[j]<=N;j++)

         {

             nprime[i*prime[j]]=;

             if(i%prime[j]==)    {h[i*prime[j]]=h[i]*prime[j]%md;break;}

             else    h[i*prime[j]]=h[i]*h[prime[j]]%md;

         }

     }

     for(i=;i<=N;i++)    dd[i]=(dd[i-]+h[i])%md;

     scanf("%lld%lld",&n,&m);

     //n=10000000;m=1000;

     if(n>m)    swap(n,m);

     for(i=;i<=m;i=j+)

     {

         if(n<i)    j=min(m,m/(m/i));

         else    j=min(m,min(n/(n/i),m/(m/i)));

         ans=(ans+Mod(dd[j]-dd[i-],md)*G(n/i,m/i))%md;

     }

     printf("%lld",ans);

     return ;

 }

upd181010:

$\sum_{k|T}k\mu(\frac{T}{k})(\frac{T}{k})^2=T\sum_{k|T}\mu(\frac{T}{k})\frac{T}{k}=T\sum_{k|T}\mu(k)k$

这个式子似乎并没有什么用...但是可以记一下，说不定下次就给出最后的那个呢

bzoj3309

https://www.lydsy.com/JudgeOnline/problem.php?id=3309

不妨设a<=b

显然原式=$\sum_{k=1}^bf(k)\sum_{i=1}^{{\lfloor}\frac{b}{k}{\rfloor}}\mu(i){{\lfloor}\frac{a}{ik}{\rfloor}}{{\lfloor}\frac{b}{ik}{\rfloor}}$

用跟上面类似的方法变到$\sum_{T=1}^b{\lfloor}\frac{a}{T}{\rfloor}{\lfloor}\frac{b}{T}{\rfloor}\sum_{k|T}f(k)\mu(\frac{T}{k})$

后面那个不太好筛啊，根本就不会，我又去查题解了。。。

（弃用）如果i%p==0，可以丢掉那些新产生因子的$\mu$，每一项是由$\mu(d)f(i/d)$变到$\mu(d)f(i/d*p)$，显然会加(i的(i/d)中"最大幂指数等于最小质因子的幂指数"的数的$\mu(d)$之和)，那个求时依赖的值也可以递推出

（弃用）如果i%p!=0，那么首先会多一些产生贡献的项，这些项是由原来i的任意因子乘上一个p得到的，是$\mu(i/d)f(d)$变到$\mu(i*p/(d*p))f(d*p)=\mu(i/d)f(dp)$，产生的贡献是原来的总和加上一个值（i的因子中"最大幂指数等于最小质因子的幂指数"的数(设为d)的$\mu(i/d)$之和）；

筛法的解释：https://blog.csdn.net/phenix_2015/article/details/50799021

这条结论，简单来讲，就是那个函数，只有各个质因子的指数相等的有值((-1)^{质因子数+1})，其他都为0

直接想好像很难想到的样子。。以后记得做这类题要打打表分解分解质因数找找规律

算是A了吧。。。

 #pragma GCC optimize(3)

 #include<cstdio>

 #include<algorithm>

 #include<cstring>

 #include<vector>

 using namespace std;

 #define fi first

 #define se second

 #define mp make_pair

 #define pb push_back

 typedef long long ll;

 typedef unsigned long long ull;

 typedef pair<int,int> pii;

 #define N 10000000

 ll prime[N+],len,mu[N+];

 ll n1[N+],n2[N+],n3[N+],h[N+];

 //分别为最小质因子次数，其他质因子次数，质因子数，函数值

 //n2为-1:第一个质因子;n2为-2:已经失败

 bool nprime[N+];

 ll a,b,ans,T;

 template<class T>

 inline void read(T &x) {

     int f=;x=;char ch=getchar();

     while(ch>''||ch<''){if(ch=='-')f=-;ch=getchar();}

     while(ch>=''&&ch<=''){x*=;x+=(ch-'');ch=getchar();}

     x*=f;

 }

 template<class T>

 inline void write(T x) {

     if(x<) putchar('-'),x=-x;

     if(x>) write(x/);

     putchar(x%+'');

 }

 int main()

 {

     ll i,j;

     mu[]=;

     for(i=;i<=N;i++)

     {

         if(!nprime[i])

         {

             prime[++len]=i;

             mu[i]=-;

             n2[i]=-;

             n3[i]=;

             n1[i]=;

             h[i]=;

         }

         for(j=;j<=len&&i*prime[j]<=N;j++)

         {

             nprime[i*prime[j]]=;

             if(i%prime[j]==)

             {

                 mu[i*prime[j]]=;

                 n1[i*prime[j]]=n1[i]+;

                 n2[i*prime[j]]=n2[i];

                 n3[i*prime[j]]=n3[i];

                 if(n2[i*prime[j]]==-||n1[i*prime[j]]==n2[i*prime[j]])

                     h[i*prime[j]]=(n3[i*prime[j]]%==)?-:;

                 else

                     h[i*prime[j]]=;

                 break;

             }

             else

             {

                 mu[i*prime[j]]=-mu[i];

                 n1[i*prime[j]]=;

                 n2[i*prime[j]]=(n2[i]==-||n2[i]==n1[i])?n1[i]:-;

                 n3[i*prime[j]]=n3[i]+;

                 if(n2[i*prime[j]]==-||n1[i*prime[j]]==n2[i*prime[j]])

                     h[i*prime[j]]=(n3[i*prime[j]]%==)?-:;

                 else

                     h[i*prime[j]]=;

             }

         }

     }

     //for(i=1;i<=1000;i++)    printf("%lld %lld %lld %lld\n",i,h[i],n1[i],n2[i]);

     //return 0;

     for(i=;i<=N;i++)    h[i]+=h[i-];

     read(T);

     //T=1000;

     while(T--)

     {

         read(a);read(b);

         //a=7558588;b=9653114;

         if(a==||b==)

         {

             puts("");

             continue;

         }

         if(a>b)    swap(a,b);

         ans=;

         for(i=;i<=a;i=j+)

         {

             j=min(a/(a/i),b/(b/i));

             ans+=(a/i)*(b/i)*(h[j]-h[i-]);

         }

         write(ans);putchar('\n');

     }

     return ;

 }

巴特西

bzoj2154||洛谷P1829 Crash的数字表格&&JZPTAB && bzoj3309 DZY Loves Math

最新文章

热门文章