编程示例:计算组合数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