001package cnslab.cnsmath;
002import java.math.BigDecimal;
003import java.math.RoundingMode;
004public class PolynomialLongDivision
005{
006//      public static final double  ACC = 1e-15; // accurace for polynomial division
007        public static final double  ACC = 0.0; // accurace for polynomial division
008
009        //returns the remainder of the polynomial long division:
010        public static Polynomial getRemainder(Polynomial numerator, Polynomial denominator)
011        {
012                double[][] table=new double[numerator.coefficients.length-denominator.coefficients.length+3][numerator.coefficients.length];//row,col
013                //fill in first row:
014                for(int a=0;a<denominator.coefficients.length;a++)
015                        table[0][a]=denominator.coefficients[denominator.coefficients.length-a-1];
016                //fill in second row:
017                for(int a=0;a<numerator.coefficients.length;a++)
018                        table[1][a]=numerator.coefficients[numerator.coefficients.length-a-1];
019                
020                //fill in the rest of the rows:
021                for(int row=2;row<table.length;row++)
022                {
023                        for(int col=0;col<table[0].length-row+1;col++)
024                        {
025//                              table[row][col]=(table[0][0]*table[row-1][col+1]-table[0][1+col]*table[row-1][0])/table[0][0];
026                                table[row][col]=(table[row-1][col+1]-table[0][1+col]*table[row-1][0]/table[0][0]);
027//                              if(Math.abs(table[row][col])<1e-15)table[row][col]=0.0;
028                        }
029                }
030                
031                //last row of table makes the coefficients, in reverse order!
032                int highestDegree=denominator.coefficients.length-2;
033                for(int a=0;a<denominator.coefficients.length-1;a++)
034                {
035//                      if(table[table.length-1][a]==0.0)
036                        if(Math.abs(table[table.length-1][a])<=ACC)
037                                highestDegree--;
038                        else
039                                break;
040                }
041                if(highestDegree<0)
042                        highestDegree=0;
043                //System.out.println("highestDegree: " +highestDegree);
044                double[] remainder=new double[highestDegree+1];
045                for(int a=0;a<remainder.length;a++)
046                        remainder[a]=(Math.abs(table[table.length-1][denominator.coefficients.length-2-a])<=ACC?0:table[table.length-1][denominator.coefficients.length-2-a]);
047//                      remainder[a]=table[table.length-1][denominator.coefficients.length-2-a];
048                
049                //print table:
050                //System.out.println("Table:");
051                //for(int a=0;a<table.length;a++)
052                //{
053                //      for(int b=0;b<table[0].length;b++)
054                //              System.out.print(table[a][b] + " ");
055                //      System.out.println("");
056                //}
057                //table=null;
058                return new Polynomial(remainder);
059        }
060        //old implementation
061/*
062        public static Polynomial getNegRemainder(Polynomial numerator, Polynomial denominator)
063        {
064                double[][] table=new double[numerator.coefficients.length-denominator.coefficients.length+3][numerator.coefficients.length];//row,col
065                //fill in first row:
066                for(int a=0;a<denominator.coefficients.length;a++)
067                        table[0][a]=denominator.coefficients[denominator.coefficients.length-a-1]/denominator.coefficients[denominator.coefficients.length-1];
068                //fill in second row:
069                for(int a=0;a<numerator.coefficients.length;a++)
070                        table[1][a]=numerator.coefficients[numerator.coefficients.length-a-1];
071                
072                //fill in the rest of the rows:
073                for(int row=2;row<table.length;row++)
074                {
075                        if(table[row-1][0]==0) //copy the line row-1 to row, left shift by 1
076                        {
077                                System.arraycopy(table[row-1],1,table[row],0,table[0].length-row+1);
078                        }
079                        else
080                        {
081                                for(int col=0;col<table[0].length-row+1;col++)
082                                {
083                                        //                      table[row][col]=(table[0][0]*table[row-1][col+1]-table[0][1+col]*table[row-1][0])/table[0][0];
084                                        //                              table[row][col]=(table[row-1][col+1]-table[0][1+col]*table[row-1][0]/table[0][0]);
085                                        table[row][col]=(table[row-1][col+1]-table[0][1+col]*table[row-1][0]);
086                                        //                              if(Math.abs(table[row][col])<1e-15)table[row][col]=0.0;
087                                }
088                        }
089                }
090                
091                //last row of table makes the coefficients, in reverse order!
092                int highestDegree=denominator.coefficients.length-2;
093                for(int a=0;a<denominator.coefficients.length-1;a++)
094                {
095//                      if(table[table.length-1][a]==0.0)
096                        if(Math.abs(table[table.length-1][a])<=ACC)
097                                highestDegree--;
098                        else
099                                break;
100                }
101                if(highestDegree<0)
102                        highestDegree=0;
103                //System.out.println("highestDegree: " +highestDegree);
104                double[] remainder=new double[highestDegree+1];
105                for(int a=0;a<remainder.length;a++)
106                        remainder[a]=(Math.abs(table[table.length-1][denominator.coefficients.length-2-a])<=ACC?0:-table[table.length-1][denominator.coefficients.length-2-a]);
107//                      remainder[a]=table[table.length-1][denominator.coefficients.length-2-a];
108                
109                //print table:
110                //System.out.println("Table:");
111                //for(int a=0;a<table.length;a++)
112                //{
113                //      for(int b=0;b<table[0].length;b++)
114                //              System.out.print(table[a][b] + " ");
115                //      System.out.println("");
116                //}
117                //table=null;
118                return new Polynomial(remainder);
119        }
120*/
121
122        //new one from nr
123        /*
124        public static Polynomial getNegRemainder(Polynomial numerator, Polynomial denominator)
125        {
126                int k,j;
127                int n=numerator.coefficients.length-1;
128                int nv=denominator.coefficients.length-1;
129                BigDecimal [] r = new BigDecimal [n+1];
130                BigDecimal [] q = new BigDecimal [n+1];
131                double [] out;
132
133                for (j=0;j<=n;j++) {
134                        r[j]=new BigDecimal(numerator.coefficients[j]);
135                        q[j]=new BigDecimal(0.0);
136                }
137                for (k=n-nv;k>=0;k--) {
138                        q[k]=new BigDecimal(r[nv+k].doubleValue()/denominator.coefficients[nv]);
139                        for (j=nv+k-1;j>=k;j--) r[j] = r[j].subtract(q[k].multiply(new BigDecimal(denominator.coefficients[j-k])));
140                }
141        // we need negative Remainder
142                for (j=0;j<nv;j++) r[j]=r[j].negate();
143                for(int a=nv-1;a>=0;a--)
144                {
145                        if(Math.abs(r[a].doubleValue())!=0)
146                        {
147                                out=new double[a+1];
148                                for (int e=0;e<a+1;e++) out[e]=r[a].doubleValue();
149                                return new Polynomial(out);
150                        }
151                }
152                return new Polynomial(new double [] {0.0});
153        }
154        */
155        //new one from nr for doesfire
156        public static Polynomial getNegRemainder(Polynomial numerator, Polynomial denominator)
157        {
158                int k,j;
159                int n=numerator.coefficients.length-1;
160                int nv=denominator.coefficients.length-1;
161                double [] r = new double [n+1];
162                double [] q = new double [n+1];
163                double [] out;
164
165                for (j=0;j<=n;j++) {
166                        r[j]=numerator.coefficients[j];
167                        q[j]=0.0;
168                }
169                for (k=n-nv;k>=0;k--) {
170                        q[k]=r[nv+k]/denominator.coefficients[nv];
171                        for (j=nv+k-1;j>=k;j--) r[j] -= q[k]*denominator.coefficients[j-k];
172                }
173        // we need negative Remainder
174                for (j=0;j<nv;j++) r[j]=-r[j];
175                for(int a=nv-1;a>=0;a--)
176                {
177//                      if(Math.abs(r[a])!=0)
178                        if(Math.abs(r[a])>1e-9)
179                        {
180                                out=new double[a+1];
181                                System.arraycopy(r,0,out,0,a+1);
182                                return new Polynomial(out);
183                        }
184                }
185                return new Polynomial(new double [] {0.0});
186        }
187
188        //new one from nr for time of fire
189        public static Polynomial getNegRemainderA(Polynomial numerator, Polynomial denominator)
190        {
191                int k,j;
192                int n=numerator.coefficients.length-1;
193                int nv=denominator.coefficients.length-1;
194                double [] r = new double [n+1];
195                double [] q = new double [n+1];
196                double [] out;
197
198                for (j=0;j<=n;j++) {
199                        r[j]=numerator.coefficients[j];
200                        q[j]=0.0;
201                }
202                for (k=n-nv;k>=0;k--) {
203                        q[k]=r[nv+k]/denominator.coefficients[nv];
204                        for (j=nv+k-1;j>=k;j--) r[j] -= q[k]*denominator.coefficients[j-k];
205                }
206        // we need negative Remainder
207                for (j=0;j<nv;j++) r[j]=-r[j];
208                for(int a=nv-1;a>=0;a--)
209                {
210                        if(Math.abs(r[a])!=0)
211//                      if(Math.abs(r[a])>1e-13)
212                        {
213                                out=new double[a+1];
214                                System.arraycopy(r,0,out,0,a+1);
215                                return new Polynomial(out);
216                        }
217                }
218                return new Polynomial(new double [] {0.0});
219        }
220}