將1到N^2表示成兩個長度為N的數列的正交和

我們可以找到兩個長度為N的數列，使得從兩個數列中各自挑選一個數相交得出的N^2個數正好是1到N^2之間的所有整數。

對於正樣的數列，如果我們將其中一個數列中所有數都加上一個常數，另外一個數列中所有數都減去一個常數，那麼必然也得到一個合格數列，所以我們可以認為這樣得到的數列對和原先的數列對 是等價的，不在我們考慮範圍之內。

問題是，相互不等價的數列對有多少個？

這個問題來源於CSDN中一道題目的討論：

http://community.csdn.net/Expert/TopicView.asp?id=4501259

而最終的計算程式由於使用大整數，將用到GMP庫。

最終結果挺有意思，同N的因子分解結果有關係。比如

N=p1^a1 *p2^a2*..*pt^at

那麼結果僅僅同a1,a2,...,at有關係，而同p1,p2,...,pt沒有關係。

   
  對於整數N,對於一個序列  
  n0=N,   n1,   n2,...,nk  
  其中n(i+1)|n(i)而且n(i+1)>n(i),   nk=1  
  我們稱這個序列為N一個長度為k的因子序列.  
  比如N長度為1的因子序列只有一個為N,1.  
  記N長度為s的因子序列的數目為L(N,s)  
  那麼L(N,0)=0,L(N,1)=1  
  上面計數問題的結果即  
    f(N)=Sum{L(N,k)*[L(N,k)+L(N,k-1)],   for   k=1,2,3,...}  
  比如對於N=6=2*3  
  長度為1的因子序列1個   {6,1}  
              2                     2個   {6,2,1},   {6,3,1}  
      長度大於3的因子序列沒有  
  所以f(6)=1*(1+0)+2*(2+1)=7  
  而對於N=p^k,   長度為t的因子序列相當於從p,p^2,...,p^(k-1)中選擇t-1個數,共有C(k-1,t-1)個  
  所以L(p^k,   t)=C(k-1,t-1)  
      f(p^k)=Sum{C(k-1,t-1)*[C(k-1,t-1)+C(k-1,t-2)],   t=1,2,...,k}  
                  =Sum{C(k-1,t-1)*C(k,t-1),   t=1,2,...,k}  
                  =Sum{C(k-1,s)*C(k,s),   s=0,1,...,k-1}  
                  =Sum{C(k-1,s)*C(k,k-s),   s=0,1,...,k-1}  
                  =C(2k-1,k)   (從2k-1個樣品中取k個的方案正好是從前k-1個中取s個,後k個中取k-s個,對所有可能的s進行求和)  
                  =(2k-1)!/(k!*(k-1)!),即gxqcn的結論.  
   
  我現在寫的程式就是對給定的N,先求出所有的L(N,s),然後計算出f(N).  
  還沒有找到更加好的辦法.  
   
  其實這個方法不僅僅對於N*N的方陣可用,對於將0,1,...,M*N-1填入一個M*N的矩陣的計數方案也可用.   
    
    
  其中L(N,s)也可以想象成下面的模型  
  假設  
  N=p1^r1*p2^r2*...*pt^rt  
  我們想象有一個袋子中總共裝了r1+r2+...+rt個球,這些球總共有t種顏色,r1個為第一種顏色,r2個為第二種顏色,...,rt個為第t種顏色.  
  現在要分s次將袋子種的球取光,請問有多少種方案?  
  這個方案數就是L(N,s).  
 

找到一種時間複雜度為O(L^2)次乘法的演算法.  
  其中L=r1+r2+...+rt  
   
  定義C(m,n)=m!/(n!*(m-n)!)是組合數.  
   
  i)計算群組合數C(m,n),其中m=1,2,...,L+max(r)-1.   0<=n<=m   (max(r)=max(r1,r2,...,rt))  
      計算方法  
          C(1,0)=C(1,1)=1  
      for(m=2;m<L+max(r);m++){  
                C(m,0)=1;  
                for(n=1;n<=m;n++)C(m,n)=C(m-1,n-1)+C(m-1,n);  
      }  
      總共花費O(L^2)次加法.  
   
  然後我們定義S(N,1)=1,S(0,k)=1   S(N,k)=S(N,k-1)+S(N-1,k),  
        其中S(N,k)就是k個非負整數和為N的表示方法數目(k個整數有序,也就是0+1和1+0不同)  
  ii)計算S(n,m),其中n=0,1,2,...,L;   1<=m<=L-n+1  
          使用上面遞推式計算,時間複雜度也是O(L^2)次加法  
   
  最後一步直接計算L(N,k)  
      其中L(N,1)=1  
  iii)   L(N,k)=C(r1+k-1,k-1)*C(r2+k-1,k-1)*...*C(rt+k-1,k-1)  
                  -Sum{   S(j+1,k-j)*L(N,k-j-1),   for   j=0   to   k-1}  
    從L(N,1)開始遞推計算到L(N,L),其中每步max(k,t)<=L次乘法,共L步,所以時間複雜度為O(L^2)次乘法.  
   
  有了L(N,k)後,我們就可以使用O(L)次乘法計算出f(N)了.

發現S(n,m)其實就是C(n+m-1,m-1),所以只要事先計算好組合數就夠了

 

  #include   <stdio.h>  
  #include   <stdlib.h>  
  #include   <memory.h>  
  #include   <gmp.h>  
   
  int   *factor_list;  
  int   factor_count;  
  int   L,   LAR;  
  mpz_t   *triangle;  
  mpz_t   *f;  
  mpz_t   N;  
   
  #define   TRIANGLE(x,y)     triangle[(x)*LAR+(y)]  
   
  void   set_triangle()  
  {  
  int   i,j;  
  triangle   =   (mpz_t   *)malloc(sizeof(mpz_t)*LAR*LAR);  
  for(i=0;i<LAR;i++)for(j=0;j<=i;j++)mpz_init(TRIANGLE(i,j));  
  mpz_set_ui(TRIANGLE(0,0),1);  
          mpz_set_ui(TRIANGLE(1,0),1);  
          mpz_set_ui(TRIANGLE(1,1),1);  
          for(i=2;i<LAR;i++){  
  mpz_set_ui(TRIANGLE(i,0),1);  
  mpz_set_ui(TRIANGLE(i,i),1);  
  for(j=1;j<i;j++){  
  mpz_add(TRIANGLE(i,j),TRIANGLE(i-1,j-1),TRIANGLE(i-1,j));  
  }  
          }  
  }  
   
   
  void   calc_f(){  
  int   i,k;  
  mpz_t   M,R;  
  mpz_init(M);  
  mpz_init(R);  
  f=   (mpz_t   *)malloc(sizeof(mpz_t)*(L+1));  
  for(i=0;i<=L;i++)mpz_init(f[i]);  
  mpz_set_ui(f[1],1);  
  printf("L[1]=1/n");  
  for(k=2;k<=L;k++){  
  mpz_set_ui(M,1);  
  for(i=0;i<factor_count;i++){  
  mpz_mul(M,M,TRIANGLE(factor_list[i]+k-1,k-1));  
  }  
  for(i=0;i<k-1;i++){  
  mpz_set(R,TRIANGLE(k,i+1));  
  mpz_mul(R,R,f[k-i-1]);  
  mpz_sub(M,M,R);  
  }  
  mpz_set(f[k],M);  
  printf("L[%d]=",k);  
  mpz_out_str(stdout,10,f[k]);  
  printf("/n");  
  }  
  mpz_clear(R);  
  mpz_clear(M);  
  }  
   
  void   output_result()  
  {  
  mpz_t   M,R;  
  int   i;  
  mpz_init(M);  
  mpz_init(R);  
  mpz_set_ui(M,1);  
  for(i=2;i<=L;i++){  
  mpz_set(R,f[i]);  
  mpz_add(R,R,f[i-1]);  
  mpz_mul(R,R,f[i]);  
  mpz_add(M,M,R);  
  }  
  printf("The   result   is:/n");  
  mpz_out_str(stdout,10,M);  
  printf("/n");  
  mpz_clear(R);  
  mpz_clear(M);  
  }  
   
  #define   PRIME_LIMIT   1048576  
  #define   SQR_LIMIT       1024  
  int   prime_flag[PRIME_LIMIT];  
  int   *prime_list;  
  int   prime_count;  
  void   init_prime(){  
  int   i,count;  
  memset(prime_flag,-1,sizeof(prime_flag));  
  prime_flag[0]=prime_flag[1]=0;  
  for(i=2;i<=SQR_LIMIT;i++){  
  if(prime_flag[i]){  
  int   j;  
  for(j=i*i;j<PRIME_LIMIT;j+=i)  
  prime_flag[j]=0;  
  }  
  }  
  for(i=2,count=0;i<PRIME_LIMIT;i++){  
  if(prime_flag[i])count++;  
  }  
  prime_count=count;  
  prime_list=(int   *)malloc(sizeof(int)*count);  
  for(i=2,count=0;i<PRIME_LIMIT;i++){  
  if(prime_flag[i]){  
  prime_list[count++]=i;  
  }  
  }  
  }  
   
  void   defactor()  
  {  
  int   i,pc,maxr;  
  mpz_t   R,M;  
  mpz_init(R);  
  mpz_init(M);  
  mpz_set(M,N);  
  for(i=0,pc=0;i<prime_count;i++){  
  if(!mpz_mod_ui(R,M,prime_list[i])){  
  pc++;  
  do{  
  mpz_fdiv_q_ui(M,M,prime_list[i]);  
  }while(!mpz_mod_ui(R,M,prime_list[i]));  
  if(!mpz_cmp_si(M,1))break;  
  }  
  }  
  if(mpz_cmp_si(M,1)){  
  int   r;  
  if(r=mpz_probab_prime_p(M,5)){//Using   prime   test   to   verify   the   left   is   a   prime  
  pc++;  
  if(r==1){  
  fprintf(stderr,"WARNING:The   left   factor   could   not   be   verified   to   be   prime   but   looks   like   to   a   prime.   It   is   treated   as   a   prime/n");  
  }  
  }else{  
  fprintf(stderr,"Cannot   find   all   prime   factors   of   input/n");  
  exit(-1);  
  }  
  }  
  factor_list=(int   *)malloc(sizeof(int)*pc);  
  factor_count=pc;L=0,maxr=0;  
  mpz_set(M,N);  
  for(i=0,pc=0;i<prime_count;i++){  
  if(!mpz_mod_ui(R,M,prime_list[i])){  
  factor_list[pc]=0;  
  do{  
  mpz_fdiv_q_ui(M,M,prime_list[i]);  
  factor_list[pc]++;  
  }while(!mpz_mod_ui(R,M,prime_list[i]));  
  L+=factor_list[pc];  
  if(maxr<factor_list[pc])maxr=factor_list[pc];  
  pc++;  
  if(!mpz_cmp_si(M,1))break;  
  }  
  }  
  if(mpz_cmp_si(M,1)){  
  factor_list[pc++]=1;  
  L++;  
  if(maxr<1)maxr=1;  
  }  
  LAR=L+maxr+1;  
  mpz_clear(R);  
  mpz_clear(M);  
  }  
   
   
  int  
  main(void)  
  {  
  mpz_init(N);  
          printf("Please   input   the   large   number:/n");  
          mpz_inp_str(N,stdin,10);  
  init_prime();  
  defactor();  
  set_triangle();  
  calc_f();  
  output_result();  
          mpz_clear(N);  
          return   0;  
  }