HDU 5970 最大公约数（找规律）

最大公约数

Time Limit: 6000/3000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 446    Accepted Submission(s): 164

Problem Description

有这样一个有关最大公约数的函数:
函数 f(x, y):

{     c=0     当 y>0:     {          c +=1          t = x % y          x = y          y = t      }      返回 c * x * x}

给出三个正整数n,m,p，你需要计算:

∑i=1n∑j=1m⌊i∗jf(i,j)⌋

对p取模的结果。

Input

包含多组测试数据。
第一行有一个正整数，表示数据的组数。
接下来每行表示一组数据，这一行有三个空格隔开的正整数n,m,p。
保证 n <= 666,666,666, m <= 666, p <= 666,666,666。
最终的测试数据中共有66组数据，并且每一个n,m,p都是在上述范围内均匀随机生成的。

Output

对于每个输入数据输出一行，这一行只包含一个整数即答案。

Sample Input

310 5 23333100 10 233331000 20 23333

Sample Output

2712235910998

Source

2016年中国大学生程序设计竞赛（合肥）-重现赛（感谢安徽大学）

        最讨厌找规律的题目了，本来还想好好做一道数论题。

        一开始，我自己假设先不考虑c。那么就变成了ΣΣi/gcd*j/gcd=Σj/gcd*Σi/gcd，如此一来，由于m比较小，我就可以枚举j，然后对应求出j所有的因子作为gcd，gcd确定之后再根据容斥来统计i/gcd的和。具体统计方法和15年沈阳regional的frog那题类似，用n的因子来进行暴力的容斥。但是很显然这样子很难把c的影响带进来，而且这里的c还要向下取整，更加的麻烦。

        于是打表找规律，首先很容易知道f(i,j)=f(i+kj,j)。根据这个，在不考虑向下取整的情况下，对于同一个j，我们就可以列出一个等差数列，其中首项是i*j/f(i,j)，公差为j*j/f(i,j)。但是这里要考虑这个向下取整。我们设i为模9为7的数j为9，可以打出如下i*j/f(i,j)的表，括号内为相邻值的差：

可以发现，每c组i*j/f(i,j)是一个循环节，也就是说可以看作c个等差数列，然后对于每一个等差数列，它的首项我们可以暴力算出，而公差也很容易求出。利用等差数列求和公式可以很快速的计算出结果。复杂度的话，我们需要枚举j和在j剩余系下的i，然后还有c个等差数列，复杂度为O(m^2logN)在接受范围之内。具体见代码：

#include <bits/stdc++.h>#define LL long long#define N 200010using namespace std;LL n,m,mod;LL f(LL i,LL j,LL &c){    LL t;    while(j)    {        c++;        t=i%j;        i=j;        j=t;    }    return c*i*i;}int main(){    int T_T;    cin>>T_T;    while(T_T--)    {        LL ans=0;        scanf("%I64d%I64d%I64d",&n,&m,&mod);        for(int i=1;i<=m;i++)            for(int j=1;j<=i&&j<=n;j++)            {                LL c=0,ff=f(j,i,c);                for(int k=0;k<c;k++)                {                    if (j+k*i>n) break;                    LL a0=(j+k*i)*i/ff;                    LL d=c*i*i/ff;                    LL num=(n-j-k*i)/(c*i)+1;                    ans=(ans+num*a0%mod+num*(num-1)/2%mod*d%mod)%mod;                }            }        printf("%I64d\n",ans);    }    return 0;}