『HDU 5970』最大公约数

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=5970

题意:(中文题,不用我说了吧.)

分析:
当时在比赛的时候并没有做出来,主要是时间不够,有点分心在其他题上,然后就水了.其实看到这种题,第一想法不用想太复杂,可能有暴力的方法,找规律!
看看数据范围
n <= 666,666,666, m <= 666, p <= 666,666,666
突破口必然在m上,那怎么用呢?
不可能暴力gcd和T的,这样必然超时,我们知道gcd（a,b）=gcd(b,a%b)
那又怎么样呢?那么给你一个提示,例如gcd(5,9)=gcd(16,9)
手写一下,你就发现确实是这样啊~！
不仅如此,f(5,9)=f(16,9)
那么我们就知道,如果以m为基数的话,那么对于所有的n,莫非就余0,余1,余2…余m-1
我们就可以这样做,然后就可以统计同余的有多少个,然后看一下规律,然后一大片的求,
就能在m^2的时间求出答案了…
别住,这个规律还是有问题啊.
它有向下取整,不可能把个i什么的拆出来,就必须一个个才行.就是**(i*j)/f(i,j)**.
那就会出现多1,少1的情况…因为取整问题…

那么我们先暴力一段序列.看看数据有什么问题.
例如 9余5的所有数的f(i,j)

这明显有规律啊,循环节啊…

这下就懵逼了,循环节该怎么搞啊?
他们的差值又不一样…,这是个好问题,但是既然是循环节,肯定等差..

我也写了很久,然后我灵感来的猜出来的.
除去第一个f（y,m）=11;
他的差值循环节就是c
然后 d=c（f（y，m）=11）（代表c倍的,余数y和m的f）.*
以上是我找到的规律,为什么是这样我就不太懂了,我也想了很久.
只要我们求出 **（31+51+72+92）**这段的初值,就可以用等差数列求循环节一整段
然后就是细节的处理了.

我的程序是卡过这题的,我也不知道为什么那么慢
算法复杂度在O（m^2*c） 我试过**On(m^2*2c)**都不行,因为这个循环节规律就是逼我在这种情况下找出来的.

/* Author:GavinjouElephant
 * Title:
 * Number:
 * main meanning：
 *
 *
 *
 */


#include <iostream>
using namespace std;
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <sstream>
#include <cctype>
#include <vector>
#include <set>
#include <cstdlib>
#include <map>
#include <queue>
//#include<initializer_list>
//#include <windows.h>
//#include <fstream>
//#include <conio.h>
#define MaxN 0x7fffffff
#define MinN -0x7fffffff
#define lson 2*k
#define rson 2*k+1
typedef long long ll;
const ll INF=0x3f3f3f3f;
const ll maxn=700;
ll Scan()//读入整数外挂.
{
    ll res = 0, ch, flag = 0;

    if((ch = getchar()) == '-')             //判断正负
        flag = 1;

    else if(ch >= '0' && ch <= '9')           //得到完整的数
        res = ch - '0';
    while((ch = getchar()) >= '0' && ch <= '9' )
        res = res * 10 + ch - '0';
    return flag ? -res : res;
}
void Out(ll a)    //输出外挂
{
    if(a>9)
        Out(a/10);
    putchar(a%10+'0');
}
ll Fc[maxn][maxn];
ll Fcx[maxn][maxn];


void  f(ll x,ll y)
{
    ll tx=x,ty=y;
    ll c=0;
    while(y>0)
    {
        c+=1;
        ll t=x%y;
        x=y;
        y=t;
    }

    Fc[tx][ty]=c;
    Fcx[tx][ty]=c*x*x;

}

void init()
{

    for(ll j=1; j<=666; j++)
    {
        for(ll i=0; i<=j; i++)
        {
            f(i,j);
        }
    }
}


ll cas(ll x,ll y)
{
    return (x*y)/Fcx[x%y][y];
}

int T;
ll n,m,p;
ll top;
ll xhj[5000];
ll Sum[5000];

int main()
{
#ifndef ONLINE_JUDGE
    freopen("coco.txt","r",stdin);
    freopen("lala.txt","w",stdout);
#endif

    init();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d%I64d%I64d",&n,&m,&p);
        ll ans=0;
        ll d;

        for(ll j=1; j<=m; j++)
        {
            for(ll y=0; y<j; y++)
            {
                ll  l=1-j-y;
                ll  r=n-y;
                ll num=(r/j - l/j);

                if(y==0)
                {
                    ans=(ans+num%p+(num*(num-1)/2)%p)%p;
                }
                else
                {
                    //-->找循环节
                    if(num==0) continue;

                    ll t=(num-1)/Fc[y][j];
                    ll tt=(num-1)%Fc[y][j];
                    ll s=cas(y,j);

                    ans=(ans+s)%p;

                    ll ss=0;
                    ll sum=0;
                    top=0;

                    int i=y+j;
                    for(top=0;top<Fc[y][j];top++,i+=j)
                    {
                        ll k=cas(i,j)-cas(i-j,j);
                        sum=(sum+k)%p;
                        xhj[top]=k;

                        if(top==0){Sum[top]=k;}
                        else {Sum[top]=Sum[top-1]+sum;}
                    }

                    if(t)
                    {   ss=(Sum[top-1]+s*top)%p;
                        d=(Fc[y][j]*sum)%p;
                        ans=(ans+(ss*t)%p+(t*(t-1)/2*d)%p)%p;
                        s=(s+(sum*t))%p;
                    }


                    if(tt){ans=(ans+Sum[tt-1]+(s*tt)%p)%p;}


                }

            }

        }
        printf("%I64d\n",(ans+p)%p);

    }



    return 0;
}