1 2 /** 3 * @name CeL polynomial function 4 * @fileoverview 5 * 本檔案包含了數學多項式的 functions。 6 * @since 7 */ 8 9 10 if (typeof CeL === 'function') 11 CeL.setup_module('data.math.polynomial', 12 function(library_namespace, load_arguments) { 13 14 // no required 15 16 17 var 18 /** 19 * null module constructor 20 * @class 數學多項式相關之 function。 21 * @constructor 22 */ 23 CeL.data.math.polynomial 24 = function () { 25 // null module constructor 26 }; 27 28 /** 29 * for JSDT: 有 prototype 才會將之當作 Class 30 */ 31 CeL.data.math.polynomial 32 .prototype = {}; 33 34 35 36 37 38 // polynomial ----------------------------------- 39 40 /* 41 return [r1,r2,..[,餘式]] 42 ** 若有無法解的餘式，會附加在最後! 43 44 高次代數方程數值求根解法: http://www.journals.zju.edu.cn/sci/2003/200303/030305.pdf http://tsg.gxtvu.com.cn/eduwest/web_courseware/maths/0092/2/2-3.htm 45 修正牛頓法 1819年霍納法 伯努利法 勞思表格法 http://en.wikipedia.org/wiki/Ruffini%27s_rule 46 Newton's method牛頓法 x2=x1-f(x1)/f'(x1) http://zh.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95 47 四次方程Finding roots http://zh.wikipedia.org/wiki/%E5%9B%9B%E6%AC%A1%E6%96%B9%E7%A8%8B 48 一元三次方程的公式解 http://en.wikipedia.org/wiki/Cubic_equation http://math.xmu.edu.cn/jszg/ynLin/JX/jiaoxueKJ/5.ppt 49 50 var rootFindingFragment=1e-15; // 因為浮點乘除法而會產生的誤差 51 */ 52 rootFinding[generateCode.dLK]='rootFindingFragment'; 53 function rootFinding(polynomial){ 54 var r=[],a,q; 55 56 //alert(NewtonMethod(polynomial)); 57 58 while(a=polynomial.length,a>1){ 59 if(a<4){ 60 if(a==2)r.push(-polynomial[1]/polynomial[0]); 61 else{ 62 a=polynomial[1]*polynomial[1]-4*polynomial[0]*polynomial[2]; // b^2-4ac 63 q=2*polynomial[0]; 64 if(a<0)a=(Math.sqrt(-a)/Math.abs(q))+'i',q=-polynomial[1]/q,r.push(q+'+'+a,q+'-'+a); 65 else a=Math.sqrt(a)/q,q=-polynomial[1]/q,r.push(q+a,q-a); 66 } 67 polynomial=[];break; 68 }else if(a=NewtonMethod(polynomial),Math.abs(a[1])>rootFindingFragment){ 69 //alert('rootFinding: NewtonMethod 無法得出根!\n誤差:'+a[1]); 70 break; 71 } 72 a=qNum(a[0],1e6);//alert(a[0]+'/'+a[1]); 73 q=pLongDivision(polynomial,[a[1],-a[0]]); 74 if(Math.abs(q[1][0])>pLongDivisionFragment){alert('rootFinding error!\n誤差:'+q[1][0]);break;} 75 r.push(a[0]/a[1]),polynomial=q[0]; 76 //alert('get root: '+a[0]+'\n'+polynomial); 77 } 78 79 if(polynomial.length==5){ // 兩對共軛虛根四次方程 80 q=[],a=polynomial.length,i=0; 81 while(--a)q.push(polynomial[i++]*a); // 微分 82 if(q=rootFinding(q),q.length>1){ 83 //a=0;for(var i=0;i<polynomial.length;i++)a=a*q[0]+polynomial[i]; 84 // 將函數上下移動至原極值有根處，則會有二重根。原函數之根應為(-b +- (b^2-4ac)^.5)/2a，則此二重根即為-b/2a（？） 85 // 故可將原函數分解為(x^2-2*q[n]*x+&)(?x^2+?x+?) 86 // 以長除法解之可得&有三解:a*&^2+(-2*q[n]*(b+2*a*q[n])-c)*&+e=0 or .. 87 q=q[0],a=4*polynomial[0]*q+polynomial[1]; 88 if(a==0){a=rootFinding([polynomial[0],-2*q*(polynomial[1]+2*q*polynomial[0])-polynomial[2],polynomial[4]]);if(a.length<2)a=null;else a=a[0];} 89 else a=(2*polynomial[2]*q+polynomial[3]-2*polynomial[0]*q*(2*polynomial[0]*q+polynomial[1]))/a; 90 var o; 91 if(!isNaN(a)&&(q=pLongDivision(polynomial,o=[1,-2*q,a]),Math.abs(q[1][0])<pLongDivisionFragment&&Math.abs(q[1][1])<pLongDivisionFragment)) 92 a=rootFinding(q[0]),r.push(a[0],a[1]),a=rootFinding(o),r.push(a[0],a[1]),polynomial=[]; 93 } 94 } 95 96 if(polynomial.length>1){ 97 r.push(polynomial); 98 //if(polynomial.length%2==1)alert('rootFinding error!'); 99 } 100 return r; 101 } 102 //alert(rootFinding(getPbyR([1,4/3,5,2,6])).join('\n')); 103 //alert(NewtonMethod(getPbyR([1,4,5,2,6])).join('\n')); 104 //alert(rootFinding([1,4,11,14,10]).join('\n')); 105 //alert(rootFinding([1,2,3,2,1]).join('\n')); 106 107 /* 長除法 polynomial long division http://en.wikipedia.org/wiki/Polynomial_long_division 2005/3/4 18:48 108 dividend/divisor=quotient..remainder 109 110 input (dividend,divisor) 111 return [商,餘式] 112 113 var pLongDivisionFragment=1e-13; // 因為浮點乘除法而會產生的誤差 114 */ 115 pLongDivision[generateCode.dLK]='pLongDivisionFragment'; 116 function pLongDivision(dividend,divisor){ 117 if(typeof dividend!='object'||typeof divisor!='object')return; 118 while(!dividend[0])dividend.shift();while(!divisor[0])dividend.shift(); 119 if(!dividend.length||!divisor.length)return; 120 121 var quotient=[],remainder=[],r,r0=divisor[0],c=-1,l2=divisor.length,l=dividend.length-l2+1,i; 122 for(i=0;i<dividend.length;i++)remainder.push(dividend[i]); 123 while(++c<l) 124 for(quotient.push(r=remainder[c]/r0),i=1;i<l2;i++){ 125 remainder[c+i]-=r*divisor[i]; 126 //if(Math.abs(remainder[c+i])<Math.abs(.00001*divisor[i]*r))remainder[c+i]=0; 127 } 128 return [quotient,remainder.slice(l)]; 129 } 130 //alert(pLongDivision([4,-5,3,1/3+2/27-1],[3,-1]).join('\n')); 131 132 /* 133 // polynomial multiplication乘法 134 function polynomialMultiplication(pol1,pol2){ 135 //for() 136 } 137 */ 138 139 /* Newton Iteration Function 2005/2/26 1:4 140 return [root,誤差] 141 */ 142 function NewtonMethod(polynomial,init,diff,count){ 143 var x=0,d,i,t,l,o,dp=[]; 144 if(!polynomial||!(d=l=polynomial.length))return; 145 while(--d)dp.push(polynomial[x++]*d); // dp:微分derivative 146 if(!diff)diff=rootFindingFragment;diff=Math.abs(diff); 147 if(!count)count=15; 148 x=init||0,o=diff+1,l--; 149 //alert(polynomial+'\n'+dp+'\n'+diff+',l:'+l); 150 while(o>diff&&count--){ 151 //alert(count+':'+x+','+d); 152 for(d=t=i=0;i<l;i++)d=d*x+polynomial[i],t=t*x+dp[i]; 153 d=d*x+polynomial[l]; 154 //alert(d+'/'+t); 155 if(t)d/=t;else d=1;//alert(); 156 t=Math.abs(d); 157 if(o<=t)if(o<rootFindingFragment)break;else x++; // test 158 o=t,x-=d; 159 } 160 return [x,d]; 161 } 162 163 // 從roots得到多項式 2005/2/26 0:45 164 function getPbyR(roots){ 165 var p,r,i,c=0,l; 166 if(!roots||!(l=roots.length))return; 167 p=[1,-roots.pop()]; 168 while(++c<l) 169 if(r=roots.pop()){p.push(-r*p[i=c]);while(i)p[i]-=p[--i]*r;} 170 else p.push(0); 171 return p; 172 } 173 174 //alert(getPbyR([1,2,3])); 175 //document.write(Newton1(getPbyR([2,32,5,3]))); 176 177 // ↑polynomial ----------------------------------- 178 179 180 181 182 183 184 return ( 185 CeL.data.math.polynomial 186 ); 187 } 188 189 190 ); 191 192