连分数因子分解法——C语言实现

参考网址：连分数分解法寻找整数的因子(Python)-优快云博客

大数运算：C语言实现 大数运算 加减乘除模运算 超详细_64编程 加减乘除取模 复杂运算-优快云博客

‌连分数因子分解法‌是一种用于大整数因子分解的算法，它是计算数论中的一个重要方法。连分数因子分解法通过寻找x2≡y2 (mod p)x2≡y2 (mod p)的形式来分解N。具体来说，这种方法涉及到计算N的简单连分数展开，并通过组合得到的同余式来寻找或组合得到因子。

连分数因子分解法的步骤大致如下：

1. ‌计算NN的简单连分数展开‌，同时计算PkPk​和QkQk​序列，其中P0=a0,P1=a0a1+1,Pk=akPk−1+Pk−2P0​=a0​,P1​=a0​a1​+1,Pk​=ak​Pk−1​+Pk−2​，Q0=1,Q1=a1,Qk=akQk−1Q0​=1,Q1​=a1​,Qk​=ak​Qk−1​。
2. ‌通过连分数展开寻找或组合得到同余式‌：x2≡y2 (mod p)x2≡y2 (mod p)，其中pp是素数，xx和yy是整数。
3. ‌利用找到的同余式‌，通过计算和组合这些同余式，可以进一步分解给定的整数NN。

这种方法在处理大整数时特别有用，因为它提供了一种有效的手段来分解那些无法通过简单试除法直接分解的大数。连分数因子分解法与二次筛法、数域筛法等方法一起，构成了处理大数因子分解的一系列工具‌1。

此外，连分数因子分解法在实际应用中也有其局限性，尤其是在处理非常大的数时，其计算复杂度和所需资源可能会成为限制因素。因此，尽管这种方法在理论和实践上都有其价值，但在实际应用中还需要考虑其他更高效的算法或技术的结合使用‌

我使用C语言实现的连分数因子分解法部分代码，由于是测试代码，目前还没有优化，存在一些问题：

/*
-----------------------连分数因子分解法----------------------
  实现时间：2024年9月5日
  作者：赵良军
*/
//求完全平方数
int is_square(unsigned char* seq_Q){
	  unsigned char* p = my_sqrt(seq_Q);
	  p = mul(p, p);
	  //printf( "s = %s\n",p );
	  if(maxer( p,seq_Q ) == 0){
	      return 1;	
	  }	  
	  return 0;
}
//求最大公约数
unsigned char* my_gcd(unsigned char* a,unsigned char* b){
	  unsigned char* c = NULL;
	  	  
    while (a){   	
        c = mol(b,a);   
        b = a;
        a = c;
                
        if( maxer( a, "0" ) == 0 || maxer( c, "0" ) == 0 ){
            return b;
        }     	
    } 
    return NULL;   
}
void cfrac(char* n){
	  unsigned char seq_P[3][BUFLEN] = {{'\0'},{'\0'},{'\0'}};
	  unsigned char seq_Q[3][BUFLEN] = {{'\0'},{'\0'},{'\0'}};
    unsigned char seq_a[3][BUFLEN] = {{'\0'},{'\0'},{'\0'}};
    unsigned char seq_p[3][BUFLEN] = {{'\0'},{'\0'},{'\0'}};
  
    strncpy(seq_P[0],"0",1);
    //strncpy(seq_P[1],"0",1);
    //strncpy(seq_P[2],"0",1);
    strncpy(seq_Q[0],"1",1);
    //strncpy(seq_Q[1],"1",1);
    //strncpy(seq_Q[2],"1",1);
      
    unsigned char* p = NULL;
    unsigned char buf[BUFLEN] = {'\0'};
    
    printf("N = %s\n",n);
      
    unsigned char* a0 = my_sqrt(n);
    
    strcpy(seq_a[0],a0);   
    strncpy(seq_p[0],"0",1);
    strcpy(seq_p[1],a0 ); 
          
    unsigned char i[BUFLEN] = {'\0'};
    strncpy(i,"1",1);
    
    unsigned char t[BUFLEN] = {'\0'};
	  while(1){
	  	  //P_{k},Q_{k}序列
	  	  printf("Start------------------------->\n");
	  	  //P序列的模数
	      p = mul(seq_a[0], seq_Q[0]);
	      strcpy(buf,p);
	  	
	  	  p = sub(buf,seq_P[0]);
	      strcpy( seq_P[0],p );
	      printf("P_M = %s\n",p );
	      
	      p = mul(seq_P[0], seq_P[0]);
	      strcpy(buf,p);
	      //printf("M0 = %s\n",p );
	      
	      p = sub(n,buf);
	      strcpy(buf,p);
	      //printf("M1 = %s\n",p );
	      
	      p = divi( buf,seq_Q[0]);
	      strcpy(seq_Q[0],p);
	      //printf("M2 = %s\n",p );
	      
	      //p = mol(buf,seq_Q[0]);
	      //strcpy(buf,p);
	      //printf("M3 = %s\n",p );
	      
	      p = add(seq_P[0],my_sqrt(n));
	      strcpy(buf,p);
	      //printf("M3 = %s\n",p );
	      
	      p = divi( buf,seq_Q[0]);
	      strcpy(seq_a[0],p);
	      //printf("M4 = %s\n",seq_a[0] );
	      //p_{k}序列
	      if( maxer(i,"1") == 0 ){
	          p = mul(a0, seq_a[0]);
	          strcpy(buf,p);
	          p = add(buf,"1");
	          strcpy(seq_p[2],p);
	      }else{
	      	  strcpy(seq_p[0],seq_p[1]);
	      	  strcpy(seq_p[1],seq_p[2]);
      	  
	      	  p = mul(seq_a[0], seq_p[1]);
	      	  strcpy(buf,p);
	      	  p = add( buf,seq_p[0] );
	      	  strcpy(seq_p[2],p);
	      }
	      
	      if( maxer( mol( i, "2" ), "0" ) == 0 && is_square(seq_Q[0]) ){
	      	  unsigned char* s = my_sqrt(seq_Q[0]);
	      	  unsigned char* p_p = NULL;
	      	  
	          printf("S = %s\n",s );
	          
	          p = add(seq_p[1],s);
	          p_p = my_gcd(p,n);     	  
	      	  printf("P+ = %s\n",p_p);
	      	  
	          p = sub(seq_p[1],s);
	          p_p = my_gcd(p,n);     	  
            printf("P- = %s\n",p_p);
            
	      	  if(p_p != NULL && maxer( p_p, "1" ) != 0){	      	  	  
	      	      break;
	      	  }
        } 
	       
	      if( maxer( mol( i, "1000" ), "0" ) == 0 ){
	      	  printf("Run to %s\n",i);
	      }
	      p = add(i,"1");
	      strcpy(i,p);
	  }
}

int main(){
	
	//unsigned char* n = "70191551";
	
	//2**67-1
	unsigned char* n = "147573952589676412927";
	
	cfrac(n);

}