Source code

001package cnslab.cnsmath;
002import java.util.TreeSet;
003//this could be implemented as a sparse array
004//which would make evaluation faster, but PLD slower!
005public class Polynomial
006{
007        //the ith order term coefficient is in the ith spot:
008        public double[] coefficients;
009        public TreeSet<KCelement> poolKC;
010        public Polynomial Q;
011        public double minAccuracy=1e-8;
012
013        public Polynomial(double[] coefficients)
014        {
015                this.coefficients=coefficients;
016        }
017
018        public int getDegree()
019        {
020                return coefficients.length-1;
021        }
022/*
023        public double evaluate(double x)
024        {
025                double sum=0.0;
026                for(int exponent=0;exponent<coefficients.length;exponent++)
027                        sum+=coefficients[exponent]*Math.pow(x,exponent);
028                return sum;
029        }
030*/
031
032        public double evaluate(double x)
033        {
034                double sum=0.0;
035                sum =  coefficients[coefficients.length-1];
036                for(int exponent=coefficients.length-2 ;exponent >= 0;exponent--) 
037                {
038                        sum = x*sum+coefficients[exponent];
039                }
040                return sum;
041        }
042
043        public double evaluateAtZero()
044        {
045                return coefficients[0];
046        }
047
048        public double evaluateAtOne()
049        {
050                double sum=0.0;
051                for(int exponent=0;exponent<coefficients.length;exponent++)
052                        sum+=coefficients[exponent];
053                return sum;
054        }
055
056        public Polynomial talorShiftOne() //g(x)=f(x+1)
057        {
058                double [] out = new double[coefficients.length];
059                out[0] = coefficients[coefficients.length-1];
060                for(int i=1; i< coefficients.length; i++)
061                {
062                        out[i] = out[i-1];
063                        double tmp= out[0];
064                        double tmp2= out[1];
065                        for(int j=1; j< i; j++)
066                        {
067                                out[j] = out[j]+tmp;
068                                tmp=tmp2;
069                                tmp2=out[j+1];
070                        }
071                        out[0] = out[0]+coefficients[coefficients.length-i-1];
072                        /*
073                        for(int j=0;j<i+1;j++)
074                                System.out.print(out[j]+" ");
075                        System.out.println();
076                        */
077                }
078                return new Polynomial(out);
079        }
080
081
082        public Polynomial talorShift(double n) //g(x)=f(x+n)
083        {
084                return this.H(n).talorShiftOne().H(1.0/n);
085
086        }
087
088        public Polynomial H(double n) //g(x)=f(nx)
089        {
090                double [] out = new double[coefficients.length];
091                out[0] = coefficients[coefficients.length-1];
092                for(int i=1; i< coefficients.length; i++)
093                {
094                        for(int j=i; j>= 1; j--)
095                        {
096                                out[j] = out[j-1]*n;
097                        }
098                        out[0] = coefficients[coefficients.length-i-1];
099                        /*
100                        for(int j=0;j<i+1;j++)
101                                System.out.print(out[j]+" ");
102                        System.out.println();
103                        */
104                }
105                return new Polynomial(out);
106        }
107
108        public Polynomial reverse() {
109                double [] out = new double[coefficients.length];
110                System.arraycopy(coefficients,0,out,0,coefficients.length);
111                int i = 0;
112                int j = coefficients.length - 1;
113                double tmp;
114                while (j > i) {
115                        tmp = out[j];
116                        out[j] = out[i];
117                        out[i] = tmp;
118                        j--;
119                        i++;
120                }
121                return new Polynomial(out);
122        }
123
124        public int count() {
125                int count=0;
126                boolean lastNeg=false;
127                boolean lastPos=false;
128
129                if(coefficients[0]>0)
130                        lastPos=true;
131                else if(coefficients[0]<0)
132                        lastNeg=true;
133
134
135                for(int i=1; i< coefficients.length; i++)
136                {
137                        if(coefficients[i]>0&&lastNeg){
138                                count++;
139                                lastPos=true;
140                                lastNeg=false;
141                        }
142                        else if(coefficients[i]<0&&lastPos){
143                                count++;
144                                lastPos=false;
145                                lastNeg=true;
146                        }
147                        else if(coefficients[i]<0&&!lastPos){
148                                lastNeg=true;
149                        }
150                        else if(coefficients[i]>0&&!lastNeg){
151                                lastPos=true;
152                        }
153                }
154                return count;
155        }
156
157        public int desBound() {
158                return this.reverse().talorShiftOne().count();
159        }
160
161        public Polynomial scale(double n)
162        {
163                double [] out = new double[coefficients.length];
164                System.arraycopy(coefficients,0,out,0,coefficients.length);
165                for(int i=0;i<out.length;i++)out[i]*=n;
166                return new Polynomial(out);
167        }
168
169        public double largestRoot() // return the largest Root in interval (0,1) or -1.0 for no root
170        {
171                if(this.count()==0)
172                {
173                        return -1.0;
174                }
175                else
176                {
177                        return 1.0-this.H(-1.0).talorShift(-1.0).smallestRoot();
178                }
179        }
180
181        public double pow2(int expo)
182        {
183                if(expo>0)
184                {
185                        return (double)(1<<expo);
186                }
187                else if(expo<0)
188                {
189                        return 1.0/((double)(1<<(-expo)));
190                }
191                else
192                {
193                        return 1.0;
194                }
195        }
196
197        public double smallestRoot() { // return the smallerst Root in interval (0,1) or -1.0 for no root
198                if(this.count()==0)return -1.0;
199//              System.out.println("root of "+ this);
200                poolKC  = new TreeSet<KCelement>();
201                poolKC.add(new KCelement(0,0));
202                double [] a = new double [coefficients.length];
203                System.arraycopy(coefficients,0,a,0,coefficients.length);
204                Q= new Polynomial(a);
205
206                KCelement first,old;
207                old = new KCelement(0,0);
208
209                while((first=poolKC.pollFirst())!=null)
210                {
211                        if(first.k==32){
212                                throw new RuntimeException("k:"+first.k+"c:"+first.c+" poly:"+this.toString());
213                                //continue; //it is a mutiply root
214                        }
215//                      System.out.println(first);
216//                      System.out.println(Q);
217                        if((int)(pow2(old.k-first.k)*(double)first.c-old.c)==1)
218                        {
219                                Q = Q.talorShiftOne().H(pow2(old.k-first.k)).scale(pow2(((coefficients.length-1)*(first.k-old.k))));
220                        }
221                        else if((int)(pow2(old.k-first.k)*(double)first.c-old.c)==0)
222                        {
223                                Q = Q.H(pow2(old.k-first.k)).scale(pow2(((coefficients.length-1)*(first.k-old.k))));
224                        }
225                        else
226                        {
227                                throw new RuntimeException("no shift ("+first.c+","+first.k+") old kc ("+old.c+","+old.k+")");
228                        }
229
230                        int countRoot = Q.desBound();
231//                      System.out.println("desBound count:"+countRoot);
232                        if(countRoot ==1)
233                        {
234//                              System.out.println("left:"+(double)first.c/(double)(1<<first.k)+" right:"+(double)(first.c+1)/(double)(1<<first.k));
235//
236                                return safeRootFind(this, (double)first.c/(double)(1<<first.k) ,(double)(first.c+1)/(double)(1<<first.k));
237                        }
238                        else if (countRoot >1)
239                        {
240                                if(Math.abs(Q.evaluate(0.5))<minAccuracy)return (2*(double)first.c+1.0)/(double)(1<<first.k+1);
241                                poolKC.add(new KCelement(2*first.c,first.k+1));
242                                poolKC.add(new KCelement(2*first.c+1,first.k+1));
243//                              System.out.println("add two size "+poolKC.size());
244
245                        }
246                        old = first;
247                }
248//              System.out.println("not found");
249                return -1.0;
250        }
251
252        public Polynomial df()
253        {
254                double[] dPrimeCoefs=new double[coefficients.length-1];
255                for(int a=1;a<coefficients.length;a++)
256                {
257                                dPrimeCoefs[a-1]=coefficients[a]*a;
258                }
259                return (new Polynomial(dPrimeCoefs));
260        }
261
262        public double safeRootFind(Polynomial Q, double guess, double end)
263        {
264                Polynomial dQ=Q.df();
265                int MAXIT=100;
266                double f1,fh;
267                double xl,xh;
268                f1=Q.evaluate(guess);
269                fh=Q.evaluate(end);
270                if( (f1>0 && fh >0) || (f1<0 && fh <0) )throw new RuntimeException("not bracketed");
271                if(Math.abs(f1)<minAccuracy)return guess;
272                //if(Math.abs(fh)<minAccuracy)return end;
273                if(f1 <0.0)
274                {
275                        xl=guess;
276                        xh=end;
277                }
278                else
279                {
280                        xh=guess;
281                        xl=end;
282                }
283                double rts = 0.5*(xl+xh); 
284                double fvalue=Q.evaluate(rts);
285                double dfvalue=dQ.evaluate(rts);
286                double dxold = Math.abs(xh-xl);
287                double dx=dxold;
288                double temp;
289                for( int j=0; j< MAXIT; j++)
290                {
291                        if((((rts-xh)*dfvalue-fvalue)*((rts-xl)*dfvalue-fvalue)>0.0) || (Math.abs(2.0*fvalue) > Math.abs(dxold*dfvalue))) //out of range or slow convergent then bisection
292                        {
293                                dxold = dx;
294                                dx=0.5*(xh-xl);
295                                rts=xl+dx;
296                                if(xl == rts) return rts;
297                        }
298                        else //Newton method
299                        {
300                                dxold=dx;
301                                dx = fvalue/dfvalue;
302                                temp = rts;
303                                rts -= dx;
304                                if(temp == rts) return rts;
305                        }
306                        if( Math.abs(dx)<minAccuracy) return rts;
307                        fvalue=Q.evaluate(rts);
308                        dfvalue=dQ.evaluate(rts);
309                        if( fvalue<0)
310                        {
311                                xl=rts;
312                        }
313                        else
314                        {
315                                xh=rts;
316                        }
317                }
318                throw new RuntimeException("max iteration");
319        //      return 0.0;
320        }
321
322
323        public String toString()
324        {
325                String ans="";
326                if(coefficients.length>0)
327                        ans+=coefficients[0];
328                
329                for(int a=1;a<coefficients.length;a++)
330                {
331                        if(coefficients[a]!=0.0)
332                                ans+=" + " + coefficients[a] + "*x^" + a;
333                }
334                return ans;
335        }
336
337        
338}