我们设总共有n个节点，AB相连的方案数为g[n]。如图，图中从1到k共k各节点是连续相连的，1与1左边那个不相连，k与k右边那个不相连，AB是这一段边上相邻的两个点。因为AB包含在1...k中，所以这条边的位置有k-1种选择；同样，这条边中有且只有一个与中心相连，从而本身有k种选择；除去1...k剩下的那一部分有n-k各节点，本身有f[n-k]种选择。

           所以，图形中A、B两点相连，并且与A、B连续相连的节点有k个时，排列的方案数为(k-1)*k*f[n-k]，枚举k,得到：

           g[n]=(2-1)*2*f[n-2]+(3-1)*3*f[n-3]+......((n-1)-1)*(n-1)*f[1];

          为了得出g[n]与f[n]的关系，计算g[n]-g[n-1]：

           g[n]-g[n-1]=2*f[n-2]+4*f[n-3]+6*f[n-4]+....2*(n-2)*f[1]

                           =2*f[n-1]

           从而：

           g[n]=2*(f[1]+f[2]+f[3]+....f[n-1]).

          f[n]是可以通过计算矩阵的幂得到的，设这个矩阵是A，A={{2,1},{1,1}},计算A+A^2+A^3+.....A^(n-1)就可以得到 f[1]+f[2]+f[3]+....f[n-1] 。如果你不清楚怎样快速计算A+A^2+A^3+.....A^(n-1)，去看一下这个题：http://acm.pku.edu.cn/JudgeOnline/problem?id=3233

          现在已经有办法快速计算如果不允许旋转，合法的方案总数了。题目中说了旋转之后算相同的方案，怎么办呢？当然是使用Burnside引理，总共n个置换，统计每个置换下保持不变的方案总数之和，再除以n就可以了。假设置换c[i]表示将整个图形顺时针转过i个节点，那么置换i的循环节总数就是gcd(i,n)，每个循环的长度是n/gcd(i,n)。对于循环节总个数是k,是每个循环节都相等的连接方案总数 与 总共有k个节点，连接不受限制的连接方案总数是相同的。因为你可以假想在整个图形中截取一小段长为gcd(i,n)，然后把他们首尾相连就构成了一个新的有gcd(i,n)个节点的图形，它的连接方案总数与原来的整个图形把每个循环节都置成相同的连接方案总数是一致的。

          n<=10^9，我们不能逐一统计每个置换的gcd(i,n)值，好在循环节长度为gcd(i,n)的置换总数量，就是（n/gcd(i,n)）的欧拉函数。按照这个方法枚举gcd(i,n)就可以快速统计所有置换中保持不变的方案总数之和了。这一部分的具体的细节可以参考POJ2154，就是刚才说了楼天成的color那个题，color到处都有解题报告（但是Toy的解题报告就不好找了，所以我决定写一份）。

         基本思路就是这样的，最后在取模的环境下除N，且不保证MN互素，我的办法是先把一开始就把M乘上N，这样最后的答案就能直接除N了，附上代码：

#include<iostream>
using namespace std;
__int64 pow[40][2][2]={{{2,1},{1,1}}};   //pow[i]表示A^(2^i);
__int64 sum[40][2][2]={{{2,1},{1,1}}};   //sum[i]表示A+A^2+...+A^(2^i);
__int64 M,N;
int facts[100],p;
int fnum[100];
int gcds[10000],g;
__int64 mul_mod(__int64 a,__int64 b)
{
      __int64 res=0;
      while(b!=0)
      {
            if(b&1)
            {
                  res+=a;
                  if(res>=M)
                        res-=M;
            }
            a=a<<1;
            if(a>=M)
                  a-=M;
            b>>=1;
      }
      return res;
}
inline void add(__int64 a[2][2],__int64 b[2][2])
{
      a[0][0]+=b[0][0];
      if(a[0][0]>=M)
            a[0][0]-=M;
      a[0][1]+=b[0][1];
      if(a[0][1]>=M)
            a[0][1]-=M;
      a[1][0]+=b[1][0];
      if(a[1][0]>=M)
            a[1][0]-=M;
      a[1][0]+=b[1][0];
      if(a[1][1]>=M)
            a[1][1]-=M;
}
inline void mov(__int64 a[2][2],__int64 b[2][2])
{
      a[0][0]=b[0][0];
      a[0][1]=b[0][1];
      a[1][0]=b[1][0];
      a[1][0]=b[1][0];
}
void mul(__int64 c[2][2],__int64 a[2][2],__int64 b[2][2])
{
      c[0][0]=mul_mod(a[0][0],b[0][0])+mul_mod(a[0][1],b[1][0]);
      if(c[0][0]>=M)
            c[0][0]-=M;
      c[0][1]=mul_mod(a[0][0],b[0][1])+mul_mod(a[0][1],b[1][1]);
      if(c[0][1]>=M)
            c[0][1]-=M;
      c[1][0]=mul_mod(a[1][0],b[0][0])+mul_mod(a[1][1],b[1][0]);
      if(c[1][0]>=M)
            c[1][0]-=M;
      c[1][1]=mul_mod(a[1][0],b[0][1])+mul_mod(a[1][1],b[1][1]);
      if(c[1][1]>=M)
            c[1][1]-=M;
}
void mat_init()
{
      for(int i=1;i<35;i++)
      {
            mul(pow[i],pow[i-1],pow[i-1]);
            mul(sum[i],sum[i-1],pow[i-1]);
            add(sum[i],sum[i-1]);
      }
}
__int64 fib(int n)
{
      int i=0;
      __int64 tmp[2][2],res[2][2]={{1,0},{0,1}};
      while(n!=0)
      {
            if(n&1)
            {
                  mul(tmp,res,pow[i]);
                  mov(res,tmp);
            }
            n>>=1;
            i++;
      }
      return res[0][1];
}
__int64 sigma_fib(int n)
{
      int i=0;
      __int64 tmp[2][2],res[2][2]={0};
      while(n!=0)
      {
            if(n&1)
            {
                  mul(tmp,res,pow[i]);
                  mov(res,tmp);
                  add(res,sum[i]);
            }
            n>>=1;
            i++;
      }
      return res[0][1];
}
__int64 stable_stat(int n)
{
      __int64 res=fib(n)+2*sigma_fib(n-1);
      res%=M;
      return res;
}
void ana()
{
    p=0;
    int i,tn=N;
    for(i=2;i*i<=tn;i++)
        if(tn%i==0)
        {
            facts[p]=i;
            fnum[p]=0;
            do
            {
                 fnum[p]++;
                 tn/=i;
            }
            while(tn%i==0);
            p++;
        }
    if(tn!=1)
    {
          fnum[p]=1;
          facts[p++]=tn;
    }
}
void allgcds(int dp,int now)
{
      if(dp==p)
            gcds[g++]=now;
      else
            for(int i=0;i<=fnum[dp];i++)
            {
                  allgcds(dp+1,now);
                  now*=facts[dp];
            }
}
__int64 phi(int n)
{
    int i,res=n;
    if(n==1)
       return 1;
    for(i=0;i<p;i++)
    {
        if(n%facts[i]==0)
        {
            res/=facts[i];
            res*=(facts[i]-1);
        }
    }
    return res;
}
__int64 work()
{
    if(N==1)
       return 1;
    __int64 t,res=phi(N)%M;
    int i;
    for(i=1;i<g;i++)
    {
          t=mul_mod(stable_stat(gcds[i]),phi(N/gcds[i]));
          res=(res+t)%M;
    }
    return res/N;
}
int main()
{
      while(scanf("%I64d%I64d",&N,&M)==2)
      {
            M*=N;
            mat_init();
            ana();
            g=0;
            allgcds(0,1);
            printf("%I64d\n",work());
      }
      return 0;
}

         顺便说一下，这份代码在HDU上跑得特别快，A掉此题的另外三个人中有两个分别跑了7000多ms和6000多ms，最快的acrush也跑了5000多ms，而这份代码只要400多ms，快了10倍多。我不知道别人用的什么算法，我觉得统计不旋转的方案总数可能有更简单的算法，我这种分析真是太麻烦了。

        在现场赛中清华它们应该用了更好的办法，它们讲的难点似乎主要是枚举循环节的一部分。如果在时间紧迫的的现场赛，统计不旋转的方案总数用我这种方法，那真是太恐怖了。希望清华他们也发一份解题报告，分享一下他们更为巧妙的方法