编程示例:计算组合数c(2000,1000)

    科技2022-07-21  163

    编程示例:计算组合数c(2000,1000)

    组合数的计算公式是 c(m,n)=m!/n!/(m-n)! 所以它的计算是依赖于阶乘的计算,当n>21时n!的结果值超出了编程语言的整数表示范围的最大值。 需要使用数组来保存计算结果。这就让阶乘计算依赖于大整数的乘法计算。大整数的乘法,使用了傅立叶变换法进行加速。 这个公式中涉及到了巨大的整数的除法计算,这是十分麻烦的。如果巨大的整数能够表示成素数的幂的乘积形式,那么计算它们的除法,就是以各个素数因子的幂为元素的向量,两个向量相对应的元素相减幂次即可。虽然巨大的整数的素数的因子分解很困难,但是对于几百,几千的阶乘结果,这样的特殊的巨大的整数,求素因子的分解形式,可以用组合数学的方式高效求解。 在对阶乘的结果进行素因子分解时,它依赖于一个素数的有序数组,这个数组是用素数的筛法生成的。 在把一个巨大的整数,从素因子的分解形式,转换成由几千位,几万位组成的十进制的精确值时,使用了指数函数的大整数版本的子函数。 编程系统提供的原始版的Math.pow(x,y) 无法计算形如 Math.pow(2,y)

    当m,n都小于60时,程序是非常短小的。如下所示: //二项式系数   Binomial(8,6) ----->28 function  Binomial(n,k) {   var result=0;      if(k>=n/2)         {  var r1=range([k+1,n]);            var pr=array_product(r1);             result=pr/Factorial(n-k);         }      else         {    var r2=range([n-k+1,n]);              var pr2=array_product(r2);                  result=pr2/Factorial(k);         }   return result; }

    //阶乘   Factorial(5)---->120    Factorial(17)---->355687428096000 function Factorial(n) {   var result=1;     for(var i=1;i<=n;i++)      {result=result*i;}     return result; }

    //求数组中的数的积 array_product([2,3,5]) function  array_product(arr) {         var result=reduce(arr,function(a,b){if (a=='') {return b;}    else {return a*b;} });        return result;   }

    //生成一个数组  range([9])---->[1,2,3,4,5,6,7,8,9]     range([0,5])---->[0,1,2,3,4,5]   range([1,5,2])--->[1,3,5] function range(arr)  {  var result=[];     if(arr.length==1)     {           for(var i=1;i<=arr[0];i++)          { result.push(i); }     }     else if (arr.length==2)     {         for(var i=arr[0];i<=arr[1];i++)          { result.push(i); }     }     else if (arr.length==3)     {         for(var i=arr[0];i<=arr[1];i=i+arr[2])          { result.push(i); }     }    return result; }

    function reduce(arr,functionname) {      var result=0;   var temp=0;   var t='';     if (isArray(arr))        {            for(var i=0;i<arr.length;i++)              {                                               if(isArray(arr[i]))                                               {                                                    t=  reduce(arr[i],functionname);                                               }                                              else {t=arr[i];}                                             temp=functionname(temp,t);                                                                         }              result=temp;        }                   else if (isObject(arr))                       { result='error this is object';}        else        {        result='error that is primative';        }        return result; }

    function isArray(obj){ return (typeof obj=='object')&&obj.constructor==Array; }

    function isObject(A) {   if( (typeof A === "object") && (A !== null) )   {     return 1;   }  else return 0; }

    当m,n都大于60时,程序是有非常多的依赖的子函数。

    代码如下: //计算 组合数 c(8,4)=70 function combination() { var m=parseInt(document.getElementById("txt0").innerText,10);   var n=parseInt(document.getElementById("txt1").innerText,10);  // var res=Binomial(m,n);  var m_arr=factor(m)  var n_arr=factor(n);  var mn_arr=factor(m-n);  var res=subtract_for_prime_factor(m_arr,n_arr);  res=subtract_for_prime_factor(res,mn_arr);    var ans=[1];  if(res.length==1)   {ans=pow_for_big_data(res[0][0],res[0][1]);}   else if (res.length==0) {ans=[];}   else {         ans=pow_for_big_data(res[0][0],res[0][1]);         for(var i=1;i<res.length;i++)         {  var temp=pow_for_big_data(res[i][0],res[i][1]);            ans=fft_mul(ans,temp);         }       }   document.getElementById("txt2").innerText=ans;//res; } // 7^14=  [72849,678223]  即678223072849 function pow_for_big_data(b,n) { var len_index=Math.floor(Math.log(n)/Math.log(2));     var bit=[1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768];      var res=[b];      if(n>1){     for(var i=len_index-1;i>=0;i--)      {             if((n&bit[i])==bit[i])               {res=fft_mul(res,res);                res=fft_mul(res,[b]);               }              else              {res=fft_mul(res,res);}      }      }      else if (n==1) {res=[b];}      else if (n==0) {res=[1];}   return res; } //素因子分解形式的幂次相减 function subtract_for_prime_factor(m_arr,n_arr) {   var min_len=Math.min(m_arr.length,n_arr.length);    for(var i=0;i<min_len;i++)    {m_arr[i][1]=m_arr[i][1]-n_arr[i][1];}    return m_arr; }

    //阶乘结果的素因子分解,利用于现有的素数的有序数组prime_array[] // 2^10*3^5*5^2*7^1*11^1*1=12! function factor(n) {     var result=[];     var str='';     for(var i=0;prime_array[i]<=n;i++)      {       var len_index=Math.floor(Math.log(n)/Math.log(prime_array[i]));       var temp=0;         for(var j=1;j<=len_index;j++)         {temp+=Math.floor(n/Math.pow(prime_array[i],j));         }         /*if (prime_array[i+1]>d.innerText)         str+=prime_array[i]+'^'+temp;         else str+=prime_array[i]+'^'+temp+'*';*/       result.push([prime_array[i],temp]);      }      return result; }

    //蝴蝶变换N function change(y, len) {     // var y=[0,1,2,3,4,5,6,7]; // var y=[0,4,2,6,1,5,3,7];      var rev=initArray([len],function(a){return 0;}); //[0,0,0,0,0,0,0,0];   // var len=8;   for (var i = 0; i < len; ++i) {     rev[i] = rev[i >> 1] >> 1;     if (i & 1) {  // 如果最后一位是 1,则翻转成 len/2       rev[i] |= len >> 1;     }   }   for (var i = 0; i < len; ++i) {     if (i < rev[i]) {  // 保证每对数只翻转一次       y=swap(y,i,rev[i]);//swap(y[i], y[rev[i]]);     }   }   return y; }

    //傅立叶变换 function fft(yy, len, on) { var PI=3.14159265358979; var  y=change(yy, len);   for (var h = 2; h <= len; h <<= 1) {                  // 模拟合并过程     var wn=[Math.cos(2 * PI / h), Math.sin(on * 2 * PI / h)];  // 计算当前单位复根     for (var j = 0; j < len; j += h) {       var w=[1, 0]; //Complex w(1, 0);  // 计算当前单位复根  var cur=[1, 0];       for (var k = j; k < j + h / 2; k++) {         var u = y[k];         var t = complex.mul(w,y[k + h / 2]);//w * y[k + h / 2];         y[k] = complex.add(u,t);//u + t;  // 这就是吧两部分分治的结果加起来         y[k + h / 2] = complex.subtract(u,t);//u - t;         // 后半个 “step” 中的ω一定和 “前半个” 中的成相反数         // “红圈”上的点转一整圈“转回来”,转半圈正好转成相反数         // 一个数相反数的平方与这个数自身的平方相等         w = complex.mul(w,wn);//w * wn;       }     }   }   if (on == -1) {     for (var i = 0; i < len; i++) {       y[i][0] /= len;     }   }   return y; }

     

    //利用傅立叶变换,计算大整数乘法 function fft_mul(a,b) { var a_arr=a; var b_arr=b; var a_plus_b=a.length+b.length; var len_index=Math.ceil(Math.log(a_plus_b)*Math.LOG2E);     var len=Math.pow(2,len_index);//8;     var a_mat= initArray([len,2],function(a,b){ if(b==1) return 0;      else if(b==0&&a>=a_arr.length) return 0;      else {return a_arr[a];}      });     var b_mat= initArray([len,2],function(a,b){ if(b==1) return 0;      else if(b==0&&a>=b_arr.length) return 0;      else {return b_arr[a];}      });           //[[2,0],[8,0],[2,0],[8,0],[0,0],[0,0],[0,0],[0,0]]      //[[8,0],[2,0],[2,0],[8,0],[0,0],[0,0],[0,0],[0,0]]     var list_str= fft(a_mat,len,1);     var list_str2= fft(b_mat,len,1);      for (var i = 0; i < len; ++i)      {      list_str[i]=complex.mul(list_str[i],list_str2[i]);      }      var  list_str3= fft(list_str,len,-1);      var ans=initArray([a_plus_b+1],function(a){return 0;});//[0,0,0,0,0,0,0,0,0,0];       for(var i=0;i<a_plus_b;i++){             ans[i]+=parseInt(list_str3[i][0]+0.5,10);             ans[i+1]+=parseInt(ans[i]/1000000,10);             ans[i]=ans[i]00000;         }         return ans; }

    //筛法求素数   function prime() {var d=document.getElementById("txt1");  var array_size=70000;  var rev=initArray([array_size],function(a){return 0;});  //var prime=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97];  var result=[];  var stop_value=Math.sqrt(d.innerText);  var low_limit=d.innerText-array_size;  for(var i=0;prime_array[i]<stop_value;i++)    {        for(var j=2*prime_array[i];j<d.innerText;j+=prime_array[i])         { if(j>=low_limit)         rev[j-low_limit]=1;}    }       for(var k=0;k<array_size;k++)     if(rev[k]==0)     {result.push(k+low_limit);}   return result;     }

    结果如下:从低位到高位显示的数组如下,每个元素代表6位的十进制数。

    c(2000,1000)= 149120,963991,976733,911108,293155,342502,906270,96311,306110,904511,653674,941966,658887,253534,181135,499101,602525,674585,96685,232980,18187,30331,749238,496946,146807,884755,465651,470377,177859,414308,408435,837648,166198,914170,274627,637248,100851,443708,837964,550307,814317,438594,890497,322553,846829,637730,762536,731662,285511,145518,687185,462063,798887,595144,104559,713065,896300,738145,738556,137700,336176,544445,684641,73786,561479,395462,886115,277685,516512,523952,148511,320647,896634,563012,261578,549769,18047,524792,798004,415481,23692,319998,464766,206201,182610,932994,755828,202083,748186,637671,820382,397033,887981,396424,825044,502980,335162,489714,626989,48151,2

     

    Processed: 0.013, SQL: 8