Wed, 26 Jun 2013 00:40:13 +0200
8017510: Add a regression test for 8005956
Summary: Regression test for 8005956
Reviewed-by: kvn, twisti
adlertz@5318 | 1 | //package com.polytechnik.utils; |
adlertz@5318 | 2 | /* |
adlertz@5318 | 3 | * (C) Vladislav Malyshkin 2010 |
adlertz@5318 | 4 | * This file is under GPL version 3. |
adlertz@5318 | 5 | * |
adlertz@5318 | 6 | */ |
adlertz@5318 | 7 | |
adlertz@5318 | 8 | /** Polynomial root. |
adlertz@5318 | 9 | * @version $Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $ |
adlertz@5318 | 10 | * @author Vladislav Malyshkin mal@gromco.com |
adlertz@5318 | 11 | */ |
adlertz@5318 | 12 | |
adlertz@5318 | 13 | /** |
adlertz@5318 | 14 | * @test |
adlertz@5318 | 15 | * @bug 8005956 |
adlertz@5318 | 16 | * @summary C2: assert(!def_outside->member(r)) failed: Use of external LRG overlaps the same LRG defined in this block |
adlertz@5318 | 17 | * |
adlertz@5318 | 18 | * @run main PolynomialRoot |
adlertz@5318 | 19 | */ |
adlertz@5318 | 20 | |
adlertz@5318 | 21 | public class PolynomialRoot { |
adlertz@5318 | 22 | |
adlertz@5318 | 23 | |
adlertz@5318 | 24 | public static int findPolynomialRoots(final int n, |
adlertz@5318 | 25 | final double [] p, |
adlertz@5318 | 26 | final double [] re_root, |
adlertz@5318 | 27 | final double [] im_root) |
adlertz@5318 | 28 | { |
adlertz@5318 | 29 | if(n==4) |
adlertz@5318 | 30 | { |
adlertz@5318 | 31 | return root4(p,re_root,im_root); |
adlertz@5318 | 32 | } |
adlertz@5318 | 33 | else if(n==3) |
adlertz@5318 | 34 | { |
adlertz@5318 | 35 | return root3(p,re_root,im_root); |
adlertz@5318 | 36 | } |
adlertz@5318 | 37 | else if(n==2) |
adlertz@5318 | 38 | { |
adlertz@5318 | 39 | return root2(p,re_root,im_root); |
adlertz@5318 | 40 | } |
adlertz@5318 | 41 | else if(n==1) |
adlertz@5318 | 42 | { |
adlertz@5318 | 43 | return root1(p,re_root,im_root); |
adlertz@5318 | 44 | } |
adlertz@5318 | 45 | else |
adlertz@5318 | 46 | { |
adlertz@5318 | 47 | throw new RuntimeException("n="+n+" is not supported yet"); |
adlertz@5318 | 48 | } |
adlertz@5318 | 49 | } |
adlertz@5318 | 50 | |
adlertz@5318 | 51 | |
adlertz@5318 | 52 | |
adlertz@5318 | 53 | static final double SQRT3=Math.sqrt(3.0),SQRT2=Math.sqrt(2.0); |
adlertz@5318 | 54 | |
adlertz@5318 | 55 | |
adlertz@5318 | 56 | private static final boolean PRINT_DEBUG=false; |
adlertz@5318 | 57 | |
adlertz@5318 | 58 | public static int root4(final double [] p,final double [] re_root,final double [] im_root) |
adlertz@5318 | 59 | { |
adlertz@5318 | 60 | if(PRINT_DEBUG) System.err.println("=====================root4:p="+java.util.Arrays.toString(p)); |
adlertz@5318 | 61 | final double vs=p[4]; |
adlertz@5318 | 62 | if(PRINT_DEBUG) System.err.println("p[4]="+p[4]); |
adlertz@5318 | 63 | if(!(Math.abs(vs)>EPS)) |
adlertz@5318 | 64 | { |
adlertz@5318 | 65 | re_root[0]=re_root[1]=re_root[2]=re_root[3]= |
adlertz@5318 | 66 | im_root[0]=im_root[1]=im_root[2]=im_root[3]=Double.NaN; |
adlertz@5318 | 67 | return -1; |
adlertz@5318 | 68 | } |
adlertz@5318 | 69 | |
adlertz@5318 | 70 | /* zsolve_quartic.c - finds the complex roots of |
adlertz@5318 | 71 | * x^4 + a x^3 + b x^2 + c x + d = 0 |
adlertz@5318 | 72 | */ |
adlertz@5318 | 73 | final double a=p[3]/vs,b=p[2]/vs,c=p[1]/vs,d=p[0]/vs; |
adlertz@5318 | 74 | if(PRINT_DEBUG) System.err.println("input a="+a+" b="+b+" c="+c+" d="+d); |
adlertz@5318 | 75 | |
adlertz@5318 | 76 | |
adlertz@5318 | 77 | final double r4 = 1.0 / 4.0; |
adlertz@5318 | 78 | final double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0; |
adlertz@5318 | 79 | final double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0; |
adlertz@5318 | 80 | final int mt; |
adlertz@5318 | 81 | |
adlertz@5318 | 82 | /* Deal easily with the cases where the quartic is degenerate. The |
adlertz@5318 | 83 | * ordering of solutions is done explicitly. */ |
adlertz@5318 | 84 | if (0 == b && 0 == c) |
adlertz@5318 | 85 | { |
adlertz@5318 | 86 | if (0 == d) |
adlertz@5318 | 87 | { |
adlertz@5318 | 88 | re_root[0]=-a; |
adlertz@5318 | 89 | im_root[0]=im_root[1]=im_root[2]=im_root[3]=0; |
adlertz@5318 | 90 | re_root[1]=re_root[2]=re_root[3]=0; |
adlertz@5318 | 91 | return 4; |
adlertz@5318 | 92 | } |
adlertz@5318 | 93 | else if (0 == a) |
adlertz@5318 | 94 | { |
adlertz@5318 | 95 | if (d > 0) |
adlertz@5318 | 96 | { |
adlertz@5318 | 97 | final double sq4 = Math.sqrt(Math.sqrt(d)); |
adlertz@5318 | 98 | re_root[0]=sq4*SQRT2/2; |
adlertz@5318 | 99 | im_root[0]=re_root[0]; |
adlertz@5318 | 100 | re_root[1]=-re_root[0]; |
adlertz@5318 | 101 | im_root[1]=re_root[0]; |
adlertz@5318 | 102 | re_root[2]=-re_root[0]; |
adlertz@5318 | 103 | im_root[2]=-re_root[0]; |
adlertz@5318 | 104 | re_root[3]=re_root[0]; |
adlertz@5318 | 105 | im_root[3]=-re_root[0]; |
adlertz@5318 | 106 | if(PRINT_DEBUG) System.err.println("Path a=0 d>0"); |
adlertz@5318 | 107 | } |
adlertz@5318 | 108 | else |
adlertz@5318 | 109 | { |
adlertz@5318 | 110 | final double sq4 = Math.sqrt(Math.sqrt(-d)); |
adlertz@5318 | 111 | re_root[0]=sq4; |
adlertz@5318 | 112 | im_root[0]=0; |
adlertz@5318 | 113 | re_root[1]=0; |
adlertz@5318 | 114 | im_root[1]=sq4; |
adlertz@5318 | 115 | re_root[2]=0; |
adlertz@5318 | 116 | im_root[2]=-sq4; |
adlertz@5318 | 117 | re_root[3]=-sq4; |
adlertz@5318 | 118 | im_root[3]=0; |
adlertz@5318 | 119 | if(PRINT_DEBUG) System.err.println("Path a=0 d<0"); |
adlertz@5318 | 120 | } |
adlertz@5318 | 121 | return 4; |
adlertz@5318 | 122 | } |
adlertz@5318 | 123 | } |
adlertz@5318 | 124 | |
adlertz@5318 | 125 | if (0.0 == c && 0.0 == d) |
adlertz@5318 | 126 | { |
adlertz@5318 | 127 | root2(new double []{p[2],p[3],p[4]},re_root,im_root); |
adlertz@5318 | 128 | re_root[2]=im_root[2]=re_root[3]=im_root[3]=0; |
adlertz@5318 | 129 | return 4; |
adlertz@5318 | 130 | } |
adlertz@5318 | 131 | |
adlertz@5318 | 132 | if(PRINT_DEBUG) System.err.println("G Path c="+c+" d="+d); |
adlertz@5318 | 133 | final double [] u=new double[3]; |
adlertz@5318 | 134 | |
adlertz@5318 | 135 | if(PRINT_DEBUG) System.err.println("Generic Path"); |
adlertz@5318 | 136 | /* For non-degenerate solutions, proceed by constructing and |
adlertz@5318 | 137 | * solving the resolvent cubic */ |
adlertz@5318 | 138 | final double aa = a * a; |
adlertz@5318 | 139 | final double pp = b - q1 * aa; |
adlertz@5318 | 140 | final double qq = c - q2 * a * (b - q4 * aa); |
adlertz@5318 | 141 | final double rr = d - q4 * a * (c - q4 * a * (b - q3 * aa)); |
adlertz@5318 | 142 | final double rc = q2 * pp , rc3 = rc / 3; |
adlertz@5318 | 143 | final double sc = q4 * (q4 * pp * pp - rr); |
adlertz@5318 | 144 | final double tc = -(q8 * qq * q8 * qq); |
adlertz@5318 | 145 | if(PRINT_DEBUG) System.err.println("aa="+aa+" pp="+pp+" qq="+qq+" rr="+rr+" rc="+rc+" sc="+sc+" tc="+tc); |
adlertz@5318 | 146 | final boolean flag_realroots; |
adlertz@5318 | 147 | |
adlertz@5318 | 148 | /* This code solves the resolvent cubic in a convenient fashion |
adlertz@5318 | 149 | * for this implementation of the quartic. If there are three real |
adlertz@5318 | 150 | * roots, then they are placed directly into u[]. If two are |
adlertz@5318 | 151 | * complex, then the real root is put into u[0] and the real |
adlertz@5318 | 152 | * and imaginary part of the complex roots are placed into |
adlertz@5318 | 153 | * u[1] and u[2], respectively. */ |
adlertz@5318 | 154 | { |
adlertz@5318 | 155 | final double qcub = (rc * rc - 3 * sc); |
adlertz@5318 | 156 | final double rcub = (rc*(2 * rc * rc - 9 * sc) + 27 * tc); |
adlertz@5318 | 157 | |
adlertz@5318 | 158 | final double Q = qcub / 9; |
adlertz@5318 | 159 | final double R = rcub / 54; |
adlertz@5318 | 160 | |
adlertz@5318 | 161 | final double Q3 = Q * Q * Q; |
adlertz@5318 | 162 | final double R2 = R * R; |
adlertz@5318 | 163 | |
adlertz@5318 | 164 | final double CR2 = 729 * rcub * rcub; |
adlertz@5318 | 165 | final double CQ3 = 2916 * qcub * qcub * qcub; |
adlertz@5318 | 166 | |
adlertz@5318 | 167 | if(PRINT_DEBUG) System.err.println("CR2="+CR2+" CQ3="+CQ3+" R="+R+" Q="+Q); |
adlertz@5318 | 168 | |
adlertz@5318 | 169 | if (0 == R && 0 == Q) |
adlertz@5318 | 170 | { |
adlertz@5318 | 171 | flag_realroots=true; |
adlertz@5318 | 172 | u[0] = -rc3; |
adlertz@5318 | 173 | u[1] = -rc3; |
adlertz@5318 | 174 | u[2] = -rc3; |
adlertz@5318 | 175 | } |
adlertz@5318 | 176 | else if (CR2 == CQ3) |
adlertz@5318 | 177 | { |
adlertz@5318 | 178 | flag_realroots=true; |
adlertz@5318 | 179 | final double sqrtQ = Math.sqrt (Q); |
adlertz@5318 | 180 | if (R > 0) |
adlertz@5318 | 181 | { |
adlertz@5318 | 182 | u[0] = -2 * sqrtQ - rc3; |
adlertz@5318 | 183 | u[1] = sqrtQ - rc3; |
adlertz@5318 | 184 | u[2] = sqrtQ - rc3; |
adlertz@5318 | 185 | } |
adlertz@5318 | 186 | else |
adlertz@5318 | 187 | { |
adlertz@5318 | 188 | u[0] = -sqrtQ - rc3; |
adlertz@5318 | 189 | u[1] = -sqrtQ - rc3; |
adlertz@5318 | 190 | u[2] = 2 * sqrtQ - rc3; |
adlertz@5318 | 191 | } |
adlertz@5318 | 192 | } |
adlertz@5318 | 193 | else if (R2 < Q3) |
adlertz@5318 | 194 | { |
adlertz@5318 | 195 | flag_realroots=true; |
adlertz@5318 | 196 | final double ratio = (R >= 0?1:-1) * Math.sqrt (R2 / Q3); |
adlertz@5318 | 197 | final double theta = Math.acos (ratio); |
adlertz@5318 | 198 | final double norm = -2 * Math.sqrt (Q); |
adlertz@5318 | 199 | |
adlertz@5318 | 200 | u[0] = norm * Math.cos (theta / 3) - rc3; |
adlertz@5318 | 201 | u[1] = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - rc3; |
adlertz@5318 | 202 | u[2] = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - rc3; |
adlertz@5318 | 203 | } |
adlertz@5318 | 204 | else |
adlertz@5318 | 205 | { |
adlertz@5318 | 206 | flag_realroots=false; |
adlertz@5318 | 207 | final double A = -(R >= 0?1:-1)*Math.pow(Math.abs(R)+Math.sqrt(R2-Q3),1.0/3.0); |
adlertz@5318 | 208 | final double B = Q / A; |
adlertz@5318 | 209 | |
adlertz@5318 | 210 | u[0] = A + B - rc3; |
adlertz@5318 | 211 | u[1] = -0.5 * (A + B) - rc3; |
adlertz@5318 | 212 | u[2] = -(SQRT3*0.5) * Math.abs (A - B); |
adlertz@5318 | 213 | } |
adlertz@5318 | 214 | if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+" u[1]="+u[1]+" u[2]="+u[2]+" qq="+qq+" disc="+((CR2 - CQ3) / 2125764.0)); |
adlertz@5318 | 215 | } |
adlertz@5318 | 216 | /* End of solution to resolvent cubic */ |
adlertz@5318 | 217 | |
adlertz@5318 | 218 | /* Combine the square roots of the roots of the cubic |
adlertz@5318 | 219 | * resolvent appropriately. Also, calculate 'mt' which |
adlertz@5318 | 220 | * designates the nature of the roots: |
adlertz@5318 | 221 | * mt=1 : 4 real roots |
adlertz@5318 | 222 | * mt=2 : 0 real roots |
adlertz@5318 | 223 | * mt=3 : 2 real roots |
adlertz@5318 | 224 | */ |
adlertz@5318 | 225 | |
adlertz@5318 | 226 | |
adlertz@5318 | 227 | final double w1_re,w1_im,w2_re,w2_im,w3_re,w3_im,mod_w1w2,mod_w1w2_squared; |
adlertz@5318 | 228 | if (flag_realroots) |
adlertz@5318 | 229 | { |
adlertz@5318 | 230 | mod_w1w2=-1; |
adlertz@5318 | 231 | mt = 2; |
adlertz@5318 | 232 | int jmin=0; |
adlertz@5318 | 233 | double vmin=Math.abs(u[jmin]); |
adlertz@5318 | 234 | for(int j=1;j<3;j++) |
adlertz@5318 | 235 | { |
adlertz@5318 | 236 | final double vx=Math.abs(u[j]); |
adlertz@5318 | 237 | if(vx<vmin) |
adlertz@5318 | 238 | { |
adlertz@5318 | 239 | vmin=vx; |
adlertz@5318 | 240 | jmin=j; |
adlertz@5318 | 241 | } |
adlertz@5318 | 242 | } |
adlertz@5318 | 243 | final double u1=u[(jmin+1)%3],u2=u[(jmin+2)%3]; |
adlertz@5318 | 244 | mod_w1w2_squared=Math.abs(u1*u2); |
adlertz@5318 | 245 | if(u1>=0) |
adlertz@5318 | 246 | { |
adlertz@5318 | 247 | w1_re=Math.sqrt(u1); |
adlertz@5318 | 248 | w1_im=0; |
adlertz@5318 | 249 | } |
adlertz@5318 | 250 | else |
adlertz@5318 | 251 | { |
adlertz@5318 | 252 | w1_re=0; |
adlertz@5318 | 253 | w1_im=Math.sqrt(-u1); |
adlertz@5318 | 254 | } |
adlertz@5318 | 255 | if(u2>=0) |
adlertz@5318 | 256 | { |
adlertz@5318 | 257 | w2_re=Math.sqrt(u2); |
adlertz@5318 | 258 | w2_im=0; |
adlertz@5318 | 259 | } |
adlertz@5318 | 260 | else |
adlertz@5318 | 261 | { |
adlertz@5318 | 262 | w2_re=0; |
adlertz@5318 | 263 | w2_im=Math.sqrt(-u2); |
adlertz@5318 | 264 | } |
adlertz@5318 | 265 | if(PRINT_DEBUG) System.err.println("u1="+u1+" u2="+u2+" jmin="+jmin); |
adlertz@5318 | 266 | } |
adlertz@5318 | 267 | else |
adlertz@5318 | 268 | { |
adlertz@5318 | 269 | mt = 3; |
adlertz@5318 | 270 | final double w_mod2_sq=u[1]*u[1]+u[2]*u[2],w_mod2=Math.sqrt(w_mod2_sq),w_mod=Math.sqrt(w_mod2); |
adlertz@5318 | 271 | if(w_mod2_sq<=0) |
adlertz@5318 | 272 | { |
adlertz@5318 | 273 | w1_re=w1_im=0; |
adlertz@5318 | 274 | } |
adlertz@5318 | 275 | else |
adlertz@5318 | 276 | { |
adlertz@5318 | 277 | // calculate square root of a complex number (u[1],u[2]) |
adlertz@5318 | 278 | // the result is in the (w1_re,w1_im) |
adlertz@5318 | 279 | final double absu1=Math.abs(u[1]),absu2=Math.abs(u[2]),w; |
adlertz@5318 | 280 | if(absu1>=absu2) |
adlertz@5318 | 281 | { |
adlertz@5318 | 282 | final double t=absu2/absu1; |
adlertz@5318 | 283 | w=Math.sqrt(absu1*0.5 * (1.0 + Math.sqrt(1.0 + t * t))); |
adlertz@5318 | 284 | if(PRINT_DEBUG) System.err.println(" Path1 "); |
adlertz@5318 | 285 | } |
adlertz@5318 | 286 | else |
adlertz@5318 | 287 | { |
adlertz@5318 | 288 | final double t=absu1/absu2; |
adlertz@5318 | 289 | w=Math.sqrt(absu2*0.5 * (t + Math.sqrt(1.0 + t * t))); |
adlertz@5318 | 290 | if(PRINT_DEBUG) System.err.println(" Path1a "); |
adlertz@5318 | 291 | } |
adlertz@5318 | 292 | if(u[1]>=0) |
adlertz@5318 | 293 | { |
adlertz@5318 | 294 | w1_re=w; |
adlertz@5318 | 295 | w1_im=u[2]/(2*w); |
adlertz@5318 | 296 | if(PRINT_DEBUG) System.err.println(" Path2 "); |
adlertz@5318 | 297 | } |
adlertz@5318 | 298 | else |
adlertz@5318 | 299 | { |
adlertz@5318 | 300 | final double vi = (u[2] >= 0) ? w : -w; |
adlertz@5318 | 301 | w1_re=u[2]/(2*vi); |
adlertz@5318 | 302 | w1_im=vi; |
adlertz@5318 | 303 | if(PRINT_DEBUG) System.err.println(" Path2a "); |
adlertz@5318 | 304 | } |
adlertz@5318 | 305 | } |
adlertz@5318 | 306 | final double absu0=Math.abs(u[0]); |
adlertz@5318 | 307 | if(w_mod2>=absu0) |
adlertz@5318 | 308 | { |
adlertz@5318 | 309 | mod_w1w2=w_mod2; |
adlertz@5318 | 310 | mod_w1w2_squared=w_mod2_sq; |
adlertz@5318 | 311 | w2_re=w1_re; |
adlertz@5318 | 312 | w2_im=-w1_im; |
adlertz@5318 | 313 | } |
adlertz@5318 | 314 | else |
adlertz@5318 | 315 | { |
adlertz@5318 | 316 | mod_w1w2=-1; |
adlertz@5318 | 317 | mod_w1w2_squared=w_mod2*absu0; |
adlertz@5318 | 318 | if(u[0]>=0) |
adlertz@5318 | 319 | { |
adlertz@5318 | 320 | w2_re=Math.sqrt(absu0); |
adlertz@5318 | 321 | w2_im=0; |
adlertz@5318 | 322 | } |
adlertz@5318 | 323 | else |
adlertz@5318 | 324 | { |
adlertz@5318 | 325 | w2_re=0; |
adlertz@5318 | 326 | w2_im=Math.sqrt(absu0); |
adlertz@5318 | 327 | } |
adlertz@5318 | 328 | } |
adlertz@5318 | 329 | if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+"u[1]="+u[1]+" u[2]="+u[2]+" absu0="+absu0+" w_mod="+w_mod+" w_mod2="+w_mod2); |
adlertz@5318 | 330 | } |
adlertz@5318 | 331 | |
adlertz@5318 | 332 | /* Solve the quadratic in order to obtain the roots |
adlertz@5318 | 333 | * to the quartic */ |
adlertz@5318 | 334 | if(mod_w1w2>0) |
adlertz@5318 | 335 | { |
adlertz@5318 | 336 | // a shorcut to reduce rounding error |
adlertz@5318 | 337 | w3_re=qq/(-8)/mod_w1w2; |
adlertz@5318 | 338 | w3_im=0; |
adlertz@5318 | 339 | } |
adlertz@5318 | 340 | else if(mod_w1w2_squared>0) |
adlertz@5318 | 341 | { |
adlertz@5318 | 342 | // regular path |
adlertz@5318 | 343 | final double mqq8n=qq/(-8)/mod_w1w2_squared; |
adlertz@5318 | 344 | w3_re=mqq8n*(w1_re*w2_re-w1_im*w2_im); |
adlertz@5318 | 345 | w3_im=-mqq8n*(w1_re*w2_im+w2_re*w1_im); |
adlertz@5318 | 346 | } |
adlertz@5318 | 347 | else |
adlertz@5318 | 348 | { |
adlertz@5318 | 349 | // typically occur when qq==0 |
adlertz@5318 | 350 | w3_re=w3_im=0; |
adlertz@5318 | 351 | } |
adlertz@5318 | 352 | |
adlertz@5318 | 353 | final double h = r4 * a; |
adlertz@5318 | 354 | if(PRINT_DEBUG) System.err.println("w1_re="+w1_re+" w1_im="+w1_im+" w2_re="+w2_re+" w2_im="+w2_im+" w3_re="+w3_re+" w3_im="+w3_im+" h="+h); |
adlertz@5318 | 355 | |
adlertz@5318 | 356 | re_root[0]=w1_re+w2_re+w3_re-h; |
adlertz@5318 | 357 | im_root[0]=w1_im+w2_im+w3_im; |
adlertz@5318 | 358 | re_root[1]=-(w1_re+w2_re)+w3_re-h; |
adlertz@5318 | 359 | im_root[1]=-(w1_im+w2_im)+w3_im; |
adlertz@5318 | 360 | re_root[2]=w2_re-w1_re-w3_re-h; |
adlertz@5318 | 361 | im_root[2]=w2_im-w1_im-w3_im; |
adlertz@5318 | 362 | re_root[3]=w1_re-w2_re-w3_re-h; |
adlertz@5318 | 363 | im_root[3]=w1_im-w2_im-w3_im; |
adlertz@5318 | 364 | |
adlertz@5318 | 365 | return 4; |
adlertz@5318 | 366 | } |
adlertz@5318 | 367 | |
adlertz@5318 | 368 | |
adlertz@5318 | 369 | |
adlertz@5318 | 370 | static void setRandomP(final double [] p,final int n,java.util.Random r) |
adlertz@5318 | 371 | { |
adlertz@5318 | 372 | if(r.nextDouble()<0.1) |
adlertz@5318 | 373 | { |
adlertz@5318 | 374 | // integer coefficiens |
adlertz@5318 | 375 | for(int j=0;j<p.length;j++) |
adlertz@5318 | 376 | { |
adlertz@5318 | 377 | if(j<=n) |
adlertz@5318 | 378 | { |
adlertz@5318 | 379 | p[j]=(r.nextInt(2)<=0?-1:1)*r.nextInt(10); |
adlertz@5318 | 380 | } |
adlertz@5318 | 381 | else |
adlertz@5318 | 382 | { |
adlertz@5318 | 383 | p[j]=0; |
adlertz@5318 | 384 | } |
adlertz@5318 | 385 | } |
adlertz@5318 | 386 | } |
adlertz@5318 | 387 | else |
adlertz@5318 | 388 | { |
adlertz@5318 | 389 | // real coefficiens |
adlertz@5318 | 390 | for(int j=0;j<p.length;j++) |
adlertz@5318 | 391 | { |
adlertz@5318 | 392 | if(j<=n) |
adlertz@5318 | 393 | { |
adlertz@5318 | 394 | p[j]=-1+2*r.nextDouble(); |
adlertz@5318 | 395 | } |
adlertz@5318 | 396 | else |
adlertz@5318 | 397 | { |
adlertz@5318 | 398 | p[j]=0; |
adlertz@5318 | 399 | } |
adlertz@5318 | 400 | } |
adlertz@5318 | 401 | } |
adlertz@5318 | 402 | if(Math.abs(p[n])<1e-2) |
adlertz@5318 | 403 | { |
adlertz@5318 | 404 | p[n]=(r.nextInt(2)<=0?-1:1)*(0.1+r.nextDouble()); |
adlertz@5318 | 405 | } |
adlertz@5318 | 406 | } |
adlertz@5318 | 407 | |
adlertz@5318 | 408 | |
adlertz@5318 | 409 | static void checkValues(final double [] p, |
adlertz@5318 | 410 | final int n, |
adlertz@5318 | 411 | final double rex, |
adlertz@5318 | 412 | final double imx, |
adlertz@5318 | 413 | final double eps, |
adlertz@5318 | 414 | final String txt) |
adlertz@5318 | 415 | { |
adlertz@5318 | 416 | double res=0,ims=0,sabs=0; |
adlertz@5318 | 417 | final double xabs=Math.abs(rex)+Math.abs(imx); |
adlertz@5318 | 418 | for(int k=n;k>=0;k--) |
adlertz@5318 | 419 | { |
adlertz@5318 | 420 | final double res1=(res*rex-ims*imx)+p[k]; |
adlertz@5318 | 421 | final double ims1=(ims*rex+res*imx); |
adlertz@5318 | 422 | res=res1; |
adlertz@5318 | 423 | ims=ims1; |
adlertz@5318 | 424 | sabs+=xabs*sabs+p[k]; |
adlertz@5318 | 425 | } |
adlertz@5318 | 426 | sabs=Math.abs(sabs); |
adlertz@5318 | 427 | if(false && sabs>1/eps? |
adlertz@5318 | 428 | (!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps)) |
adlertz@5318 | 429 | : |
adlertz@5318 | 430 | (!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))) |
adlertz@5318 | 431 | { |
adlertz@5318 | 432 | throw new RuntimeException( |
adlertz@5318 | 433 | getPolinomTXT(p)+"\n"+ |
adlertz@5318 | 434 | "\t x.r="+rex+" x.i="+imx+"\n"+ |
adlertz@5318 | 435 | "res/sabs="+(res/sabs)+" ims/sabs="+(ims/sabs)+ |
adlertz@5318 | 436 | " sabs="+sabs+ |
adlertz@5318 | 437 | "\nres="+res+" ims="+ims+" n="+n+" eps="+eps+" "+ |
adlertz@5318 | 438 | " sabs>1/eps="+(sabs>1/eps)+ |
adlertz@5318 | 439 | " f1="+(!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps))+ |
adlertz@5318 | 440 | " f2="+(!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))+ |
adlertz@5318 | 441 | " "+txt); |
adlertz@5318 | 442 | } |
adlertz@5318 | 443 | } |
adlertz@5318 | 444 | |
adlertz@5318 | 445 | static String getPolinomTXT(final double [] p) |
adlertz@5318 | 446 | { |
adlertz@5318 | 447 | final StringBuilder buf=new StringBuilder(); |
adlertz@5318 | 448 | buf.append("order="+(p.length-1)+"\t"); |
adlertz@5318 | 449 | for(int k=0;k<p.length;k++) |
adlertz@5318 | 450 | { |
adlertz@5318 | 451 | buf.append("p["+k+"]="+p[k]+";"); |
adlertz@5318 | 452 | } |
adlertz@5318 | 453 | return buf.toString(); |
adlertz@5318 | 454 | } |
adlertz@5318 | 455 | |
adlertz@5318 | 456 | static String getRootsTXT(int nr,final double [] re,final double [] im) |
adlertz@5318 | 457 | { |
adlertz@5318 | 458 | final StringBuilder buf=new StringBuilder(); |
adlertz@5318 | 459 | for(int k=0;k<nr;k++) |
adlertz@5318 | 460 | { |
adlertz@5318 | 461 | buf.append("x."+k+"("+re[k]+","+im[k]+")\n"); |
adlertz@5318 | 462 | } |
adlertz@5318 | 463 | return buf.toString(); |
adlertz@5318 | 464 | } |
adlertz@5318 | 465 | |
adlertz@5318 | 466 | static void testRoots(final int n, |
adlertz@5318 | 467 | final int n_tests, |
adlertz@5318 | 468 | final java.util.Random rn, |
adlertz@5318 | 469 | final double eps) |
adlertz@5318 | 470 | { |
adlertz@5318 | 471 | final double [] p=new double [n+1]; |
adlertz@5318 | 472 | final double [] rex=new double [n],imx=new double [n]; |
adlertz@5318 | 473 | for(int i=0;i<n_tests;i++) |
adlertz@5318 | 474 | { |
adlertz@5318 | 475 | for(int dg=n;dg-->-1;) |
adlertz@5318 | 476 | { |
adlertz@5318 | 477 | for(int dr=3;dr-->0;) |
adlertz@5318 | 478 | { |
adlertz@5318 | 479 | setRandomP(p,n,rn); |
adlertz@5318 | 480 | for(int j=0;j<=dg;j++) |
adlertz@5318 | 481 | { |
adlertz@5318 | 482 | p[j]=0; |
adlertz@5318 | 483 | } |
adlertz@5318 | 484 | if(dr==0) |
adlertz@5318 | 485 | { |
adlertz@5318 | 486 | p[0]=-1+2.0*rn.nextDouble(); |
adlertz@5318 | 487 | } |
adlertz@5318 | 488 | else if(dr==1) |
adlertz@5318 | 489 | { |
adlertz@5318 | 490 | p[0]=p[1]=0; |
adlertz@5318 | 491 | } |
adlertz@5318 | 492 | |
adlertz@5318 | 493 | findPolynomialRoots(n,p,rex,imx); |
adlertz@5318 | 494 | |
adlertz@5318 | 495 | for(int j=0;j<n;j++) |
adlertz@5318 | 496 | { |
adlertz@5318 | 497 | //System.err.println("j="+j); |
adlertz@5318 | 498 | checkValues(p,n,rex[j],imx[j],eps," t="+i); |
adlertz@5318 | 499 | } |
adlertz@5318 | 500 | } |
adlertz@5318 | 501 | } |
adlertz@5318 | 502 | } |
adlertz@5318 | 503 | System.err.println("testRoots(): n_tests="+n_tests+" OK, dim="+n); |
adlertz@5318 | 504 | } |
adlertz@5318 | 505 | |
adlertz@5318 | 506 | |
adlertz@5318 | 507 | |
adlertz@5318 | 508 | |
adlertz@5318 | 509 | static final double EPS=0; |
adlertz@5318 | 510 | |
adlertz@5318 | 511 | public static int root1(final double [] p,final double [] re_root,final double [] im_root) |
adlertz@5318 | 512 | { |
adlertz@5318 | 513 | if(!(Math.abs(p[1])>EPS)) |
adlertz@5318 | 514 | { |
adlertz@5318 | 515 | re_root[0]=im_root[0]=Double.NaN; |
adlertz@5318 | 516 | return -1; |
adlertz@5318 | 517 | } |
adlertz@5318 | 518 | re_root[0]=-p[0]/p[1]; |
adlertz@5318 | 519 | im_root[0]=0; |
adlertz@5318 | 520 | return 1; |
adlertz@5318 | 521 | } |
adlertz@5318 | 522 | |
adlertz@5318 | 523 | public static int root2(final double [] p,final double [] re_root,final double [] im_root) |
adlertz@5318 | 524 | { |
adlertz@5318 | 525 | if(!(Math.abs(p[2])>EPS)) |
adlertz@5318 | 526 | { |
adlertz@5318 | 527 | re_root[0]=re_root[1]=im_root[0]=im_root[1]=Double.NaN; |
adlertz@5318 | 528 | return -1; |
adlertz@5318 | 529 | } |
adlertz@5318 | 530 | final double b2=0.5*(p[1]/p[2]),c=p[0]/p[2],d=b2*b2-c; |
adlertz@5318 | 531 | if(d>=0) |
adlertz@5318 | 532 | { |
adlertz@5318 | 533 | final double sq=Math.sqrt(d); |
adlertz@5318 | 534 | if(b2<0) |
adlertz@5318 | 535 | { |
adlertz@5318 | 536 | re_root[1]=-b2+sq; |
adlertz@5318 | 537 | re_root[0]=c/re_root[1]; |
adlertz@5318 | 538 | } |
adlertz@5318 | 539 | else if(b2>0) |
adlertz@5318 | 540 | { |
adlertz@5318 | 541 | re_root[0]=-b2-sq; |
adlertz@5318 | 542 | re_root[1]=c/re_root[0]; |
adlertz@5318 | 543 | } |
adlertz@5318 | 544 | else |
adlertz@5318 | 545 | { |
adlertz@5318 | 546 | re_root[0]=-b2-sq; |
adlertz@5318 | 547 | re_root[1]=-b2+sq; |
adlertz@5318 | 548 | } |
adlertz@5318 | 549 | im_root[0]=im_root[1]=0; |
adlertz@5318 | 550 | } |
adlertz@5318 | 551 | else |
adlertz@5318 | 552 | { |
adlertz@5318 | 553 | final double sq=Math.sqrt(-d); |
adlertz@5318 | 554 | re_root[0]=re_root[1]=-b2; |
adlertz@5318 | 555 | im_root[0]=sq; |
adlertz@5318 | 556 | im_root[1]=-sq; |
adlertz@5318 | 557 | } |
adlertz@5318 | 558 | return 2; |
adlertz@5318 | 559 | } |
adlertz@5318 | 560 | |
adlertz@5318 | 561 | public static int root3(final double [] p,final double [] re_root,final double [] im_root) |
adlertz@5318 | 562 | { |
adlertz@5318 | 563 | final double vs=p[3]; |
adlertz@5318 | 564 | if(!(Math.abs(vs)>EPS)) |
adlertz@5318 | 565 | { |
adlertz@5318 | 566 | re_root[0]=re_root[1]=re_root[2]= |
adlertz@5318 | 567 | im_root[0]=im_root[1]=im_root[2]=Double.NaN; |
adlertz@5318 | 568 | return -1; |
adlertz@5318 | 569 | } |
adlertz@5318 | 570 | final double a=p[2]/vs,b=p[1]/vs,c=p[0]/vs; |
adlertz@5318 | 571 | /* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0 |
adlertz@5318 | 572 | */ |
adlertz@5318 | 573 | final double q = (a * a - 3 * b); |
adlertz@5318 | 574 | final double r = (a*(2 * a * a - 9 * b) + 27 * c); |
adlertz@5318 | 575 | |
adlertz@5318 | 576 | final double Q = q / 9; |
adlertz@5318 | 577 | final double R = r / 54; |
adlertz@5318 | 578 | |
adlertz@5318 | 579 | final double Q3 = Q * Q * Q; |
adlertz@5318 | 580 | final double R2 = R * R; |
adlertz@5318 | 581 | |
adlertz@5318 | 582 | final double CR2 = 729 * r * r; |
adlertz@5318 | 583 | final double CQ3 = 2916 * q * q * q; |
adlertz@5318 | 584 | final double a3=a/3; |
adlertz@5318 | 585 | |
adlertz@5318 | 586 | if (R == 0 && Q == 0) |
adlertz@5318 | 587 | { |
adlertz@5318 | 588 | re_root[0]=re_root[1]=re_root[2]=-a3; |
adlertz@5318 | 589 | im_root[0]=im_root[1]=im_root[2]=0; |
adlertz@5318 | 590 | return 3; |
adlertz@5318 | 591 | } |
adlertz@5318 | 592 | else if (CR2 == CQ3) |
adlertz@5318 | 593 | { |
adlertz@5318 | 594 | /* this test is actually R2 == Q3, written in a form suitable |
adlertz@5318 | 595 | for exact computation with integers */ |
adlertz@5318 | 596 | |
adlertz@5318 | 597 | /* Due to finite precision some double roots may be missed, and |
adlertz@5318 | 598 | will be considered to be a pair of complex roots z = x +/- |
adlertz@5318 | 599 | epsilon i close to the real axis. */ |
adlertz@5318 | 600 | |
adlertz@5318 | 601 | final double sqrtQ = Math.sqrt (Q); |
adlertz@5318 | 602 | |
adlertz@5318 | 603 | if (R > 0) |
adlertz@5318 | 604 | { |
adlertz@5318 | 605 | re_root[0] = -2 * sqrtQ - a3; |
adlertz@5318 | 606 | re_root[1]=re_root[2]=sqrtQ - a3; |
adlertz@5318 | 607 | im_root[0]=im_root[1]=im_root[2]=0; |
adlertz@5318 | 608 | } |
adlertz@5318 | 609 | else |
adlertz@5318 | 610 | { |
adlertz@5318 | 611 | re_root[0]=re_root[1] = -sqrtQ - a3; |
adlertz@5318 | 612 | re_root[2]=2 * sqrtQ - a3; |
adlertz@5318 | 613 | im_root[0]=im_root[1]=im_root[2]=0; |
adlertz@5318 | 614 | } |
adlertz@5318 | 615 | return 3; |
adlertz@5318 | 616 | } |
adlertz@5318 | 617 | else if (R2 < Q3) |
adlertz@5318 | 618 | { |
adlertz@5318 | 619 | final double sgnR = (R >= 0 ? 1 : -1); |
adlertz@5318 | 620 | final double ratio = sgnR * Math.sqrt (R2 / Q3); |
adlertz@5318 | 621 | final double theta = Math.acos (ratio); |
adlertz@5318 | 622 | final double norm = -2 * Math.sqrt (Q); |
adlertz@5318 | 623 | final double r0 = norm * Math.cos (theta/3) - a3; |
adlertz@5318 | 624 | final double r1 = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - a3; |
adlertz@5318 | 625 | final double r2 = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - a3; |
adlertz@5318 | 626 | |
adlertz@5318 | 627 | re_root[0]=r0; |
adlertz@5318 | 628 | re_root[1]=r1; |
adlertz@5318 | 629 | re_root[2]=r2; |
adlertz@5318 | 630 | im_root[0]=im_root[1]=im_root[2]=0; |
adlertz@5318 | 631 | return 3; |
adlertz@5318 | 632 | } |
adlertz@5318 | 633 | else |
adlertz@5318 | 634 | { |
adlertz@5318 | 635 | final double sgnR = (R >= 0 ? 1 : -1); |
adlertz@5318 | 636 | final double A = -sgnR * Math.pow (Math.abs (R) + Math.sqrt (R2 - Q3), 1.0 / 3.0); |
adlertz@5318 | 637 | final double B = Q / A; |
adlertz@5318 | 638 | |
adlertz@5318 | 639 | re_root[0]=A + B - a3; |
adlertz@5318 | 640 | im_root[0]=0; |
adlertz@5318 | 641 | re_root[1]=-0.5 * (A + B) - a3; |
adlertz@5318 | 642 | im_root[1]=-(SQRT3*0.5) * Math.abs(A - B); |
adlertz@5318 | 643 | re_root[2]=re_root[1]; |
adlertz@5318 | 644 | im_root[2]=-im_root[1]; |
adlertz@5318 | 645 | return 3; |
adlertz@5318 | 646 | } |
adlertz@5318 | 647 | |
adlertz@5318 | 648 | } |
adlertz@5318 | 649 | |
adlertz@5318 | 650 | |
adlertz@5318 | 651 | static void root3a(final double [] p,final double [] re_root,final double [] im_root) |
adlertz@5318 | 652 | { |
adlertz@5318 | 653 | if(Math.abs(p[3])>EPS) |
adlertz@5318 | 654 | { |
adlertz@5318 | 655 | final double v=p[3], |
adlertz@5318 | 656 | a=p[2]/v,b=p[1]/v,c=p[0]/v, |
adlertz@5318 | 657 | a3=a/3,a3a=a3*a, |
adlertz@5318 | 658 | pd3=(b-a3a)/3, |
adlertz@5318 | 659 | qd2=a3*(a3a/3-0.5*b)+0.5*c, |
adlertz@5318 | 660 | Q=pd3*pd3*pd3+qd2*qd2; |
adlertz@5318 | 661 | if(Q<0) |
adlertz@5318 | 662 | { |
adlertz@5318 | 663 | // three real roots |
adlertz@5318 | 664 | final double SQ=Math.sqrt(-Q); |
adlertz@5318 | 665 | final double th=Math.atan2(SQ,-qd2); |
adlertz@5318 | 666 | im_root[0]=im_root[1]=im_root[2]=0; |
adlertz@5318 | 667 | final double f=2*Math.sqrt(-pd3); |
adlertz@5318 | 668 | re_root[0]=f*Math.cos(th/3)-a3; |
adlertz@5318 | 669 | re_root[1]=f*Math.cos((th+2*Math.PI)/3)-a3; |
adlertz@5318 | 670 | re_root[2]=f*Math.cos((th+4*Math.PI)/3)-a3; |
adlertz@5318 | 671 | //System.err.println("3r"); |
adlertz@5318 | 672 | } |
adlertz@5318 | 673 | else |
adlertz@5318 | 674 | { |
adlertz@5318 | 675 | // one real & two complex roots |
adlertz@5318 | 676 | final double SQ=Math.sqrt(Q); |
adlertz@5318 | 677 | final double r1=-qd2+SQ,r2=-qd2-SQ; |
adlertz@5318 | 678 | final double v1=Math.signum(r1)*Math.pow(Math.abs(r1),1.0/3), |
adlertz@5318 | 679 | v2=Math.signum(r2)*Math.pow(Math.abs(r2),1.0/3), |
adlertz@5318 | 680 | sv=v1+v2; |
adlertz@5318 | 681 | // real root |
adlertz@5318 | 682 | re_root[0]=sv-a3; |
adlertz@5318 | 683 | im_root[0]=0; |
adlertz@5318 | 684 | // complex roots |
adlertz@5318 | 685 | re_root[1]=re_root[2]=-0.5*sv-a3; |
adlertz@5318 | 686 | im_root[1]=(v1-v2)*(SQRT3*0.5); |
adlertz@5318 | 687 | im_root[2]=-im_root[1]; |
adlertz@5318 | 688 | //System.err.println("1r2c"); |
adlertz@5318 | 689 | } |
adlertz@5318 | 690 | } |
adlertz@5318 | 691 | else |
adlertz@5318 | 692 | { |
adlertz@5318 | 693 | re_root[0]=re_root[1]=re_root[2]=im_root[0]=im_root[1]=im_root[2]=Double.NaN; |
adlertz@5318 | 694 | } |
adlertz@5318 | 695 | } |
adlertz@5318 | 696 | |
adlertz@5318 | 697 | |
adlertz@5318 | 698 | static void printSpecialValues() |
adlertz@5318 | 699 | { |
adlertz@5318 | 700 | for(int st=0;st<6;st++) |
adlertz@5318 | 701 | { |
adlertz@5318 | 702 | //final double [] p=new double []{8,1,3,3.6,1}; |
adlertz@5318 | 703 | final double [] re_root=new double [4],im_root=new double [4]; |
adlertz@5318 | 704 | final double [] p; |
adlertz@5318 | 705 | final int n; |
adlertz@5318 | 706 | if(st<=3) |
adlertz@5318 | 707 | { |
adlertz@5318 | 708 | if(st<=0) |
adlertz@5318 | 709 | { |
adlertz@5318 | 710 | p=new double []{2,-4,6,-4,1}; |
adlertz@5318 | 711 | //p=new double []{-6,6,-6,8,-2}; |
adlertz@5318 | 712 | } |
adlertz@5318 | 713 | else if(st==1) |
adlertz@5318 | 714 | { |
adlertz@5318 | 715 | p=new double []{0,-4,8,3,-9}; |
adlertz@5318 | 716 | } |
adlertz@5318 | 717 | else if(st==2) |
adlertz@5318 | 718 | { |
adlertz@5318 | 719 | p=new double []{-1,0,2,0,-1}; |
adlertz@5318 | 720 | } |
adlertz@5318 | 721 | else |
adlertz@5318 | 722 | { |
adlertz@5318 | 723 | p=new double []{-5,2,8,-2,-3}; |
adlertz@5318 | 724 | } |
adlertz@5318 | 725 | root4(p,re_root,im_root); |
adlertz@5318 | 726 | n=4; |
adlertz@5318 | 727 | } |
adlertz@5318 | 728 | else |
adlertz@5318 | 729 | { |
adlertz@5318 | 730 | p=new double []{0,2,0,1}; |
adlertz@5318 | 731 | if(st==4) |
adlertz@5318 | 732 | { |
adlertz@5318 | 733 | p[1]=-p[1]; |
adlertz@5318 | 734 | } |
adlertz@5318 | 735 | root3(p,re_root,im_root); |
adlertz@5318 | 736 | n=3; |
adlertz@5318 | 737 | } |
adlertz@5318 | 738 | System.err.println("======== n="+n); |
adlertz@5318 | 739 | for(int i=0;i<=n;i++) |
adlertz@5318 | 740 | { |
adlertz@5318 | 741 | if(i<n) |
adlertz@5318 | 742 | { |
adlertz@5318 | 743 | System.err.println(String.valueOf(i)+"\t"+ |
adlertz@5318 | 744 | p[i]+"\t"+ |
adlertz@5318 | 745 | re_root[i]+"\t"+ |
adlertz@5318 | 746 | im_root[i]); |
adlertz@5318 | 747 | } |
adlertz@5318 | 748 | else |
adlertz@5318 | 749 | { |
adlertz@5318 | 750 | System.err.println(String.valueOf(i)+"\t"+p[i]+"\t"); |
adlertz@5318 | 751 | } |
adlertz@5318 | 752 | } |
adlertz@5318 | 753 | } |
adlertz@5318 | 754 | } |
adlertz@5318 | 755 | |
adlertz@5318 | 756 | |
adlertz@5318 | 757 | |
adlertz@5318 | 758 | public static void main(final String [] args) |
adlertz@5318 | 759 | { |
adlertz@5318 | 760 | final long t0=System.currentTimeMillis(); |
adlertz@5318 | 761 | final double eps=1e-6; |
adlertz@5318 | 762 | //checkRoots(); |
adlertz@5318 | 763 | final java.util.Random r=new java.util.Random(-1381923); |
adlertz@5318 | 764 | printSpecialValues(); |
adlertz@5318 | 765 | |
adlertz@5318 | 766 | final int n_tests=10000000; |
adlertz@5318 | 767 | //testRoots(2,n_tests,r,eps); |
adlertz@5318 | 768 | //testRoots(3,n_tests,r,eps); |
adlertz@5318 | 769 | testRoots(4,n_tests,r,eps); |
adlertz@5318 | 770 | final long t1=System.currentTimeMillis(); |
adlertz@5318 | 771 | System.err.println("PolynomialRoot.main: "+n_tests+" tests OK done in "+(t1-t0)+" milliseconds. ver=$Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $"); |
adlertz@5318 | 772 | } |
adlertz@5318 | 773 | |
adlertz@5318 | 774 | |
adlertz@5318 | 775 | |
adlertz@5318 | 776 | } |