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