CoreFlows/Models/src/utilitaire_algebre.cxx

   1 #include "utilitaire_algebre.h"
   2 #include <iostream>
   3 using namespace std;
   4
   5 roots_polynoms::roots_polynoms(): smalno(1.2e-38), infin(3.4e38),
   6                                   eta(2.22e-16), base(2)
   7 {
   8 }
   9
  10 int roots_polynoms::rpoly(double *op, int degree, double *zeror, double *zeroi)
  11 {
  12   double t,aa,bb,cc,*temp,factor,rot;
  13   double *pt;
  14   double lo,max,min,xx,yy,cosr,sinr,xxx,x,sc,bnd;
  15   double xm,ff,df,dx;
  16   int cnt,nz,i,j,jj,l,nm1,zerok;
  17
  18   /*  The following statements set machine constants. */
  19   are = eta;
  20   mre = eta;
  21   lo = smalno/eta;
  22   /*  Initialization of constants for shift rotation. */
  23   xx = sqrt(0.5);
  24   yy = -xx;
  25   rot = 94.0;
  26   rot *= 0.017453293;
  27   cosr = cos(rot);
  28   sinr = sin(rot);
  29   n = degree;
  30   /*  Algorithm fails of the leading coefficient is zero. */
  31   if (fabs(op[0]) <=  smalno) return -1;
  32   /*  Remove the zeros at the origin, if any. */
  33   while (op[n] == 0.0) {
  34     j = degree - n;
  35     zeror[j] = 0.0;
  36     zeroi[j] = 0.0;
  37     n--;
  38   }
  39   if (n < 1) return -1;
  40
  41   /*
  42    *  Allocate memory here
  43    */
  44   temp = new double [degree+1];
  45   pt = new double [degree+1];
  46   p = new double [degree+1];
  47   qp = new double [degree+1];
  48   k = new double [degree+1];
  49   qk = new double [degree+1];
  50   svk = new double [degree+1];
  51   /*  Make a copy of the coefficients. */
  52   for (i=0;i<=n;i++)
  53     p[i] = op[i];
  54   /*  Start the algorithm for one zero. */
  55  _40:
  56   if (n == 1) {
  57     zeror[degree-1] = -p[1]/p[0];
  58     zeroi[degree-1] = 0.0;
  59     n -= 1;
  60     goto _99;
  61   }
  62   /*  Calculate the final zero or pair of zeros. */
  63   if (n == 2) {
  64     quad(p[0],p[1],p[2],&zeror[degree-2],&zeroi[degree-2],
  65          &zeror[degree-1],&zeroi[degree-1]);
  66     n -= 2;
  67     goto _99;
  68   }
  69   /*  Find largest and smallest moduli of coefficients. */
  70   max = 0.0;
  71   min = infin;
  72   for (i=0;i<=n;i++) {
  73     x = fabs(p[i]);
  74     if (x > max) max = x;
  75     if (x != 0.0 && x < min){ min = x;
  76     //cerr << " x= " << x << " smalno " << smalno << endl;
  77     }
  78   }
  79   /*  Scale if there are large or very small coefficients.
  80    *  Computes a scale factor to multiply the coefficients of the
  81    *  polynomial. The scaling si done to avoid overflow and to
  82    *  avoid undetected underflow interfering with the convergence
  83    *  criterion. The factor is a power of the base.
  84    */
  85   sc = lo/min;
  86   if (sc > 1.0 && infin/sc < max) goto _110;
  87   if (sc <= 1.0) {
  88     if (max < 10.0) goto _110;
  89     if (sc == 0.0)
  90       sc = smalno;
  91   }
  92   l = (int)(log(sc)/log(base) + 0.5);
  93   factor = pow(base*1.0,l);
  94   if (factor != 1.0) {
  95     for (i=0;i<=n;i++)
  96       p[i] = factor*p[i];     /* Scale polynomial. */
  97   }
  98  _110:
  99   /*  Compute lower bound on moduli of roots. */
 100   for (i=0;i<=n;i++) {
 101     pt[i] = (fabs(p[i]));
 102   }
 103   pt[n] = - pt[n];
 104   /*  Compute upper estimate of bound. */
 105   x = exp((log(-pt[n])-log(pt[0])) / (double)n);
 106   /*  If Newton step at the origin is better, use it. */
 107   //if (pt[n-1]  !=  0.0) {
 108   if (fabs(pt[n-1] ) >  smalno*fabs(pt[n])) {
 109     xm = -pt[n]/pt[n-1];
 110     if (xm < x)  x = xm;
 111   }
 112   /*  Chop the interval (0,x) until ff <= 0 */
 113   while (1) {
 114     xm = x*0.1;
 115     ff = pt[0];
 116     for (i=1;i<=n;i++)
 117       ff = ff*xm + pt[i];
 118     if (ff <= 0.0) break;
 119     x = xm;
 120   }
 121   dx = x;
 122   /*  Do Newton interation until x converges to two
 123    *  decimal places.
 124    */
 125   while (fabs(dx/x) > 0.005) {
 126     ff = pt[0];
 127     df = ff;
 128     for (i=1;i<n;i++) {
 129       ff = ff*x + pt[i];
 130       df = df*x + ff;
 131     }
 132     ff = ff*x + pt[n];
 133     dx = ff/df;
 134     x -= dx;
 135   }
 136   bnd = x;
 137   /*  Compute the derivative as the initial k polynomial
 138    *  and do 5 steps with no shift.
 139    */
 140   nm1 = n - 1;
 141   for (i=1;i<n;i++)
 142     k[i] = (double)(n-i)*p[i]/(double)n;
 143   k[0] = p[0];
 144   aa = p[n];
 145   bb = p[n-1];
 146   zerok = (fabs(k[n-1])< smalno);
 147   for(jj=0;jj<5;jj++) {
 148     cc = k[n-1];
 149     if (!zerok) {//cerr << "cc= " << cc << endl;
 150       /*  Use a scaled form of recurrence if value of k at 0 is nonzero. */
 151       t = -aa/cc;
 152       for (i=0;i<nm1;i++) {
 153         j = n-i-1;
 154         k[j] = t*k[j-1]+p[j];
 155       }
 156       k[0] = p[0];
 157       zerok = (fabs(k[n-1]) <= fabs(bb)*eta*10.0 || (fabs(k[n-1])< smalno));
 158     }
 159     else {
 160       /*  Use unscaled form of recurrence. */
 161       for (i=0;i<nm1;i++) {
 162         j = n-i-1;
 163         k[j] = k[j-1];
 164       }
 165       k[0] = 0.0;
 166       zerok = (fabs(k[n-1])< smalno);
 167     }
 168   }
 169   /*  Save k for restarts with new shifts. */
 170   for (i=0;i<n;i++)
 171     temp[i] = k[i];
 172   /*  Loop to select the quadratic corresponding to each new shift. */
 173   for (cnt = 0;cnt < 20;cnt++) {
 174     /*  Quadratic corresponds to a double shift to a
 175      *  non-real point and its complex conjugate. The point
 176      *  has modulus bnd and amplitude rotated by 94 degrees
 177      *  from the previous shift.
 178      */
 179     xxx = cosr*xx - sinr*yy;
 180     yy = sinr*xx + cosr*yy;
 181     xx = xxx;
 182     sr = bnd*xx;
 183     si = bnd*yy;
 184     u = -2.0 * sr;
 185     v = bnd;
 186     fxshfr(20*(cnt+1),&nz);
 187     if (nz != 0) {
 188       /*  The second stage jumps directly to one of the third
 189        *  stage iterations and returns here if successful.
 190        *  Deflate the polynomial, store the zero or zeros and
 191        *  return to the main algorithm.
 192        */
 193       j = degree - n;
 194       zeror[j] = szr;
 195       zeroi[j] = szi;
 196       n -= nz;
 197       for (i=0;i<=n;i++)
 198         p[i] = qp[i];
 199       if (nz != 1) {
 200         zeror[j+1] = lzr;
 201         zeroi[j+1] = lzi;
 202       }
 203       goto _40;
 204     }
 205     /*  If the iteration is unsuccessful another quadratic
 206      *  is chosen after restoring k.
 207      */
 208     for (i=0;i<n;i++) {
 209       k[i] = temp[i];
 210     }
 211   }
 212   /*  Return with failure if no convergence with 20 shifts. */
 213  _99:
 214   delete [] svk;
 215   delete [] qk;
 216   delete [] k;
 217   delete [] qp;
 218   delete [] p;
 219   delete [] pt;
 220   delete [] temp;
 221
 222 //   cerr << " rpoly roots " << endl;
 223 //   for (int ii=0; ii<degree;ii++) {
 224 //     cerr << " (" << zeror[ii] << "," << zeroi[ii] << ")";
 225 //   }
 226 //   cerr << endl;
 227 //   cerr << " rpoly roots point " << &zeror[0] << " " <<  &zeroi[0] << endl;
 228 //   cerr << " rpoly roots  point 1 " << &zeror[1] << " " <<  &zeroi[1] << endl;
 229
 230
 231   return degree - n;
 232 }
 233
 234 /*  Computes up to L2 fixed shift k-polynomials,
 235  *  testing for convergence in the linear or quadratic
 236  *  case. Initiates one of the variable shift
 237  *  iterations and returns with the number of zeros
 238  *  found.
 239  */
 240  void roots_polynoms::fxshfr(int l2,int *nz)
 241 {
 242   double svu,svv,ui,vi,s;
 243   double betas,betav,oss,ovv,ss,vv,ts,tv;
 244   double ots,otv,tvv,tss;
 245   int type, i,j,iflag,vpass,spass,vtry,stry;
 246
 247   *nz = 0;
 248   betav = 0.25;
 249   betas = 0.25;
 250   oss = sr;
 251   ovv = v;
 252   /*  Evaluate polynomial by synthetic division. */
 253   quadsd(n,&u,&v,p,qp,&a,&b);
 254   calcsc(&type);
 255   for (j=0;j<l2;j++) {
 256     /*  Calculate next k polynomial and estimate v. */
 257     nextk(&type);
 258     calcsc(&type);
 259     newest(type,&ui,&vi);
 260     vv = vi;
 261     /*  Estimate s. */
 262     ss = 0.0;
 263     if (k[n-1] != 0.0) ss = -p[n]/k[n-1];
 264     tv = 1.0;
 265     ts = 1.0;
 266     if (j == 0 || type == 3) goto _70;
 267     /*  Compute relative measures of convergence of s and v sequences. */
 268     if (vv != 0.0) tv = fabs((vv-ovv)/vv);
 269     if (ss != 0.0) ts = fabs((ss-oss)/ss);
 270     /*  If decreasing, multiply two most recent convergence measures. */
 271     tvv = 1.0;
 272     if (tv < otv) tvv = tv*otv;
 273     tss = 1.0;
 274     if (ts < ots) tss = ts*ots;
 275     /*  Compare with convergence criteria. */
 276     vpass = (tvv < betav);
 277     spass = (tss < betas);
 278     if (!(spass || vpass)) goto _70;
 279     /*  At least one sequence has passed the convergence test.
 280      *  Store variables before iterating.
 281      */
 282     svu = u;
 283     svv = v;
 284     for (i=0;i<n;i++) {
 285       svk[i] = k[i];
 286     }
 287     s = ss;
 288     /*  Choose iteration according to the fastest converging
 289      *  sequence.
 290      */
 291     vtry = 0;
 292     stry = 0;
 293     if (spass && (!vpass) || tss < tvv) goto _40;
 294   _20:
 295     quadit(&ui,&vi,nz);
 296     if (*nz > 0) return;
 297     /*  Quadratic iteration has failed. Flag that it has
 298      *  been tried and decrease the convergence criterion.
 299      */
 300     vtry = 1;
 301     betav *= 0.25;
 302     /*  Try linear iteration if it has not been tried and
 303      *  the S sequence is converging.
 304      */
 305     if (stry || !spass) goto _50;
 306     for (i=0;i<n;i++) {
 307       k[i] = svk[i];
 308     }
 309   _40:
 310     realit(&s,nz,&iflag);
 311     if (*nz > 0) return;
 312     /*  Linear iteration has failed. Flag that it has been
 313      *  tried and decrease the convergence criterion.
 314      */
 315     stry = 1;
 316     betas *=0.25;
 317     if (iflag == 0) goto _50;
 318     /*  If linear iteration signals an almost double real
 319      *  zero attempt quadratic iteration.
 320      */
 321     ui = -(s+s);
 322     vi = s*s;
 323     goto _20;
 324     /*  Restore variables. */
 325   _50:
 326     u = svu;
 327     v = svv;
 328     for (i=0;i<n;i++) {
 329       k[i] = svk[i];
 330     }
 331     /*  Try quadratic iteration if it has not been tried
 332      *  and the V sequence is convergin.
 333      */
 334     if (vpass && !vtry) goto _20;
 335     /*  Recompute QP and scalar values to continue the
 336      *  second stage.
 337      */
 338     quadsd(n,&u,&v,p,qp,&a,&b);
 339     calcsc(&type);
 340   _70:
 341     ovv = vv;
 342     oss = ss;
 343     otv = tv;
 344     ots = ts;
 345   }
 346 }
 347 /*  Variable-shift k-polynomial iteration for a
 348  *  quadratic factor converges only if the zeros are
 349  *  equimodular or nearly so.
 350  *  uu, vv - coefficients of starting quadratic.
 351  *  nz - number of zeros found.
 352  */
 353  void roots_polynoms::quadit(double *uu,double *vv,int *nz)
 354 {
 355   double ui,vi;
 356   double mp,omp,ee,relstp,t,zm;
 357   int type,i,j,tried;
 358
 359   *nz = 0;
 360   tried = 0;
 361   u = *uu;
 362   v = *vv;
 363   j = 0;
 364   /*  Main loop. */
 365  _10:
 366   quad(1.0,u,v,&szr,&szi,&lzr,&lzi);
 367   /*  Return if roots of the quadratic are real and not
 368    *  close to multiple or nearly equal and of opposite
 369    *  sign.
 370    */
 371   if (fabs(fabs(szr)-fabs(lzr)) > 0.01 * fabs(lzr)) return;
 372   /*  Evaluate polynomial by quadratic synthetic division. */
 373   quadsd(n,&u,&v,p,qp,&a,&b);
 374   mp = fabs(a-szr*b) + fabs(szi*b);
 375   /*  Compute a rigorous bound on the rounding error in
 376    *  evaluating p.
 377    */
 378   zm = sqrt(fabs(v));
 379   ee = 2.0*fabs(qp[0]);
 380   t = -szr*b;
 381   for (i=1;i<n;i++) {
 382     ee = ee*zm + fabs(qp[i]);
 383   }
 384   ee = ee*zm + fabs(a+t);
 385   ee *= (5.0 *mre + 4.0*are);
 386   ee = ee - (5.0*mre+2.0*are)*(fabs(a+t)+fabs(b)*zm)+2.0*are*fabs(t);
 387   /*  Iteration has converged sufficiently if the
 388    *  polynomial value is less than 20 times this bound.
 389    */
 390   if (mp <= 20.0*ee) {
 391     *nz = 2;
 392     return;
 393   }
 394   j++;
 395   /*  Stop iteration after 20 steps. */
 396   if (j > 20) return;
 397   if (j < 2) goto _50;
 398   if (relstp > 0.01 || mp < omp || tried) goto _50;
 399   /*  A cluster appears to be stalling the convergence.
 400    *  Five fixed shift steps are taken with a u,v close
 401    *  to the cluster.
 402    */
 403   if (relstp < eta) relstp = eta;
 404   relstp = sqrt(relstp);
 405   u = u - u*relstp;
 406   v = v + v*relstp;
 407   quadsd(n,&u,&v,p,qp,&a,&b);
 408   for (i=0;i<5;i++) {
 409     calcsc(&type);
 410     nextk(&type);
 411   }
 412   tried = 1;
 413   j = 0;
 414  _50:
 415   omp = mp;
 416   /*  Calculate next k polynomial and new u and v. */
 417   calcsc(&type);
 418   nextk(&type);
 419   calcsc(&type);
 420   newest(type,&ui,&vi);
 421   /*  If vi is zero the iteration is not converging. */
 422   if (vi == 0.0) return;
 423   relstp = fabs((vi-v)/vi);
 424   u = ui;
 425   v = vi;
 426   goto _10;
 427 }
 428 /*  Variable-shift H polynomial iteration for a real zero.
 429  *  sss - starting iterate
 430  *  nz  - number of zeros found
 431  *  iflag - flag to indicate a pair of zeros near real axis.
 432  */
 433  void roots_polynoms::realit(double *sss, int *nz, int *iflag)
 434 {
 435   double pv,kv,t,s;
 436   double ms,mp,omp,ee;
 437   int i,j;
 438
 439   *nz = 0;
 440   s = *sss;
 441   *iflag = 0;
 442   j = 0;
 443   /*  Main loop */
 444   while (1) {
 445     pv = p[0];
 446     /*  Evaluate p at s. */
 447     qp[0] = pv;
 448     for (i=1;i<=n;i++) {
 449       pv = pv*s + p[i];
 450       qp[i] = pv;
 451     }
 452     mp = fabs(pv);
 453     /*  Compute a rigorous bound on the error in evaluating p. */
 454     ms = fabs(s);
 455     ee = (mre/(are+mre))*fabs(qp[0]);
 456     for (i=1;i<=n;i++) {
 457       ee = ee*ms + fabs(qp[i]);
 458     }
 459     /*  Iteration has converged sufficiently if the polynomial
 460      *  value is less than 20 times this bound.
 461      */
 462     if (mp <= 20.0*((are+mre)*ee-mre*mp)) {
 463       *nz = 1;
 464       szr = s;
 465       szi = 0.0;
 466       return;
 467     }
 468     j++;
 469     /*  Stop iteration after 10 steps. */
 470     if (j > 10) return;
 471     if (j < 2) goto _50;
 472     if (fabs(t) > 0.001*fabs(s-t) || mp < omp) goto _50;
 473     /*  A cluster of zeros near the real axis has been
 474      *  encountered. Return with iflag set to initiate a
 475      *  quadratic iteration.
 476      */
 477     *iflag = 1;
 478     *sss = s;
 479     return;
 480     /*  Return if the polynomial value has increased significantly. */
 481   _50:
 482     omp = mp;
 483     /*  Compute t, the next polynomial, and the new iterate. */
 484     kv = k[0];
 485     qk[0] = kv;
 486     for (i=1;i<n;i++) {
 487       kv = kv*s + k[i];
 488       qk[i] = kv;
 489     }
 490     if (fabs(kv) <= fabs(k[n-1])*10.0*eta) {
 491       /*  Use unscaled form. */
 492       k[0] = 0.0;
 493       for (i=1;i<n;i++) {
 494         k[i] = qk[i-1];
 495       }
 496     }
 497     else {
 498       /*  Use the scaled form of the recurrence if the value
 499        *  of k at s is nonzero.
 500        */
 501       t = -pv/kv;
 502       k[0] = qp[0];
 503       for (i=1;i<n;i++) {
 504         k[i] = t*qk[i-1] + qp[i];
 505       }
 506     }
 507     kv = k[0];
 508     for (i=1;i<n;i++) {
 509       kv = kv*s + k[i];
 510     }
 511     t = 0.0;
 512     if (fabs(kv) > fabs(k[n-1]*10.0*eta)) t = -pv/kv;
 513     s += t;
 514   }
 515 }
 516
 517 /*  This routine calculates scalar quantities used to
 518  *  compute the next k polynomial and new estimates of
 519  *  the quadratic coefficients.
 520  *  type - integer variable set here indicating how the
 521  *  calculations are normalized to avoid overflow.
 522  */
 523  void roots_polynoms::calcsc(int *type)
 524 {
 525   /*  Synthetic division of k by the quadratic 1,u,v */
 526   quadsd(n-1,&u,&v,k,qk,&c,&d);
 527   if (fabs(c) > fabs(k[n-1]*100.0*eta)) goto _10;
 528   if (fabs(d) > fabs(k[n-2]*100.0*eta)) goto _10;
 529   *type = 3;
 530   /*  Type=3 indicates the quadratic is almost a factor of k. */
 531   return;
 532  _10:
 533   if (fabs(d) < fabs(c)) {
 534     *type = 1;
 535     /*  Type=1 indicates that all formulas are divided by c. */
 536     e = a/c;
 537     f = d/c;
 538     g = u*e;
 539     h = v*b;
 540     a3 = a*e + (h/c+g)*b;
 541     a1 = b - a*(d/c);
 542     a7 = a + g*d + h*f;
 543     return;
 544   }
 545   *type = 2;
 546   /*  Type=2 indicates that all formulas are divided by d. */
 547   e = a/d;
 548   f = c/d;
 549   g = u*b;
 550   h = v*b;
 551   a3 = (a+g)*e + h*(b/d);
 552   a1 = b*f-a;
 553   a7 = (f+u)*a + h;
 554 }
 555 /*  Computes the next k polynomials using scalars
 556  *  computed in calcsc.
 557  */
 558  void roots_polynoms::nextk(int *type)
 559 {
 560   double temp;
 561   int i;
 562
 563   if (*type == 3) {
 564     /*  Use unscaled form of the recurrence if type is 3. */
 565     k[0] = 0.0;
 566     k[1] = 0.0;
 567     for (i=2;i<n;i++) {
 568       k[i] = qk[i-2];
 569     }
 570     return;
 571   }
 572   temp = a;
 573   if (*type == 1) temp = b;
 574   if (fabs(a1) <= fabs(temp)*eta*10.0) {
 575     /*  If a1 is nearly zero then use a special form of the
 576      *  recurrence.
 577      */
 578     k[0] = 0.0;
 579     k[1] = -a7*qp[0];
 580     for(i=2;i<n;i++) {
 581       k[i] = a3*qk[i-2] - a7*qp[i-1];
 582     }
 583     return;
 584   }
 585   /*  Use scaled form of the recurrence. */
 586   a7 /= a1;
 587   a3 /= a1;
 588   k[0] = qp[0];
 589   k[1] = qp[1] - a7*qp[0];
 590   for (i=2;i<n;i++) {
 591     k[i] = a3*qk[i-2] - a7*qp[i-1] + qp[i];
 592   }
 593 }
 594 /*  Compute new estimates of the quadratic coefficients
 595  *  using the scalars computed in calcsc.
 596  */
 597  void roots_polynoms::newest(int type,double *uu,double *vv)
 598 {
 599   double a4,a5,b1,b2,c1,c2,c3,c4,temp;
 600   //  double smalno = 1.2e-38;
 601   //  double smalno = 0.0;
 602
 603   /* Use formulas appropriate to setting of type. */
 604   if (type == 3) {
 605     /*  If type=3 the quadratic is zeroed. */
 606     *uu = 0.0;
 607     *vv = 0.0;
 608     return;
 609   }
 610   if (type == 2) {
 611     a4 = (a+g)*f + h;
 612     a5 = (f+u)*c + v*d;
 613   }
 614   else {
 615     a4 = a + u*b +h*f;
 616     a5 = c + (u+v*f)*d;
 617   }
 618   /*  Evaluate new quadratic coefficients. */
 619   b1 = -k[n-1]/p[n];
 620   b2 = -(k[n-2]+b1*p[n-1])/p[n];
 621   c1 = v*b2*a1;
 622   c2 = b1*a7;
 623   c3 = b1*b1*a3;
 624   c4 = c1 - c2 - c3;
 625   temp = a5 + b1*a4 - c4;
 626   if (fabs(temp) <= smalno) {
 627     *uu = 0.0;
 628     *vv = 0.0;
 629     return;
 630   }
 631   *uu = u - (u*(c3+c2)+v*(b1*a1+b2*a7))/temp;
 632   *vv = v*(1.0+c4/temp);
 633   return;
 634 }
 635
 636 /*  Divides p by the quadratic 1,u,v placing the quotient
 637  *  in q and the remainder in a,b.
 638  */
 639  void roots_polynoms::quadsd(int my_n,double *my_u,double *my_v,
 640                                          double *my_p,double *my_q,
 641                                          double *my_a,double *my_b)
 642 {
 643   double my_c;
 644   int i;
 645   *my_b = my_p[0];
 646   my_q[0] = *my_b;
 647   *my_a = my_p[1] - (*my_b)*(*my_u);
 648   my_q[1] = *my_a;
 649   for (i=2;i<=my_n;i++) {
 650     my_c = my_p[i] - (*my_a)*(*my_u) - (*my_b)*(*my_v);
 651     my_q[i] = my_c;
 652     *my_b = *my_a;
 653     *my_a = my_c;
 654   }
 655 }
 656 /*  Calculate the zeros of the quadratic a*z^2 + b1*z + c.
 657  *  The quadratic formula, modified to avoid overflow, is used
 658  *  to find the larger zero if the zeros are real and both
 659  *  are complex. The smaller real zero is found directly from
 660  *  the product of the zeros c/a.
 661  */
 662  void roots_polynoms::quad(double my_a,double my_b1,double my_c,double *my_sr,double *my_si,
 663                                        double *my_lr,double *my_li)
 664 {
 665   double my_b,my_d,my_e;
 666   //  double smalno = 1.2e-38;
 667
 668   if (fabs(my_a) <= smalno) {         /* less than two roots */
 669     if (fabs(my_b1) >= smalno)
 670       *my_sr = -my_c/my_b1;
 671     else
 672     *my_sr = 0.0;
 673     *my_lr = 0.0;
 674     *my_si = 0.0;
 675     *my_li = 0.0;
 676     return;
 677   }
 678   if (fabs(my_c) <= smalno) {         /* one real root, one zero root */
 679     *my_sr = 0.0;
 680     *my_lr = -my_b1/my_a;
 681     *my_si = 0.0;
 682     *my_li = 0.0;
 683     return;
 684   }
 685   /* Compute discriminant avoiding overflow. */
 686   my_b = my_b1/2.0;
 687   if (fabs(my_b) < fabs(my_c)) {
 688     if (my_c < 0.0)
 689       my_e = -my_a;
 690     else
 691       my_e = my_a;
 692     my_e = my_b*(my_b/fabs(my_c)) - my_e;
 693     my_d = sqrt(fabs(my_e))*sqrt(fabs(my_c));
 694   }
 695   else {
 696     my_e = 1.0 - (my_a/my_b)*(my_c/my_b);
 697     my_d = sqrt(fabs(my_e))*fabs(my_b);
 698   }
 699   if (my_e < 0.0) {      /* complex conjugate zeros */
 700     *my_sr = -my_b/my_a;
 701     *my_lr = *my_sr;
 702     *my_si = fabs(my_d/my_a);
 703     *my_li = -(*my_si);
 704   }
 705   else {
 706     if (my_b >= 0.0)   /* real zeros. */
 707       my_d = -my_d;
 708     *my_lr = (-my_b+my_d)/my_a;
 709     *my_sr = 0.0;
 710     if (*my_lr != 0.0)
 711       *my_sr = (my_c/ *my_lr)/my_a;
 712     *my_si = 0.0;
 713     *my_li = 0.0;
 714   }
 715 }
 716
 717  //qques outils d'algebre lineaire
 718
 719  /** \fn add
 720   * \brief Calcule la somme de 2vecteurs ( U=U+V)
 721   * @param *U,*V les 2 vecteurs
 722   * @param N taille des 2 vecteurs
 723   * @return la somme de U et V dans U*/
 724  void Polynoms::add
 725  (
 726                 double *U,              //vecteur auquel on ajoute
 727                 int N,                  //taille du vecteur
 728                 const double *V //ajout
 729  )
 730  {
 731         for(int i=0;i<N;i++) U[i] += V[i];
 732  };
 733
 734  /** \fn tensor
 735   * \brief Calcule le tenseur de 2 vecteurs
 736   * @param *a le 1er vecteur (*double)
 737   * @param N la taille du premier vecteur (int)
 738   * @param *b le 2ieme vecteur (*double)
 739   * @param M la taille du 2ieme vecteur (int)
 740   * @param *ab le resultat
 741   * @return */
 742
 743  void Polynoms::tensor
 744  (      const double *a,        //vecteur gauche
 745                 int N,                  //taille du vecteur gauche
 746                 const double *b,        //vecteur droit
 747                 int M,                  //taille du vecteur droit
 748                 double *ab              //conteneur de sortie
 749  )
 750  {
 751         for(int i=0;i<N;i++)
 752                 for(int j=0;j<M;j++)
 753                         ab[i*M+j] = a[i]*b[j];
 754  };
 755
 756  /** \fn matrixProduct
 757   * \brief calcule le produit de 2 matrices .
 758   * @param *Matrix1 la 1ere matrice
 759   * @param R1,C1 la taille de la matrice1 Nbr de ligne et le nbr de colonne
 760   * @param Matrix2 la 2ieme matrice
 761   * @param R2,C2 la taille de la matrice2 Nbr de ligne et le nbr de colonne
 762   * @param out le resultat de Matrix1*Matrix2
 763   * @return */
 764  void Polynoms::matrixProduct
 765  (
 766                 const double *Matrix1,
 767                 const int &R1,
 768                 const int &C1,
 769                 const double *Matrix2,
 770                 const int &R2,
 771                 const int &C2,
 772                 double *out
 773  )
 774  {
 775         for(int i=0; i<R1; i++)
 776         {
 777                 for(int j=0; j<C2; j++)
 778                 {
 779                         out[i*C2+j] = 0;
 780                         for(int k=0; k<C1; k++)
 781                         {
 782                                 out[i*C2+j] += Matrix1[i*C1+k]*Matrix2[k*C2+j];
 783                         }
 784                 }
 785         }
 786  };
 787  /** \fn matrixProdVec
 788   * \brief Calcule le produit Matrice vecteur
 789   * @param *Matrix correspond à la matrice du produit
 790   * @param R1,C1 la taille de la matrice( Nbr de ligne et le nbr de colonne)
 791   * @param Vector correspond au vecteur du produit
 792   * @param out le resultat de Matrix*Vector
 793   * @return   */
 794  void Polynoms::matrixProdVec
 795  (      const double *Matrix,
 796                 const int &R1,
 797                 const int &C1,
 798                 const double *Vector,
 799                 double *out
 800  )
 801  {
 802         for(int i=0; i<R1; i++)
 803         {
 804                 out[i] = 0;
 805                 for(int j=0; j<C1; j++)
 806                 {
 807                         out[i] += Matrix[i*C1+j]*Vector[j];
 808                 }
 809         }
 810  };
 811
 812  /** \fn shift_diagonal
 813   * \brief Calcule la trace d'une matrice carrée
 814   * @param *Matrix correspond  à la matrice
 815   * @param N taille de la matrice
 816   * @param shift résultat
 817   * @return renvoie la trace de la matrice Matrix dans shift */
 818  void Polynoms::shift_diagonal( double *Matrix, const int N, double shift)
 819  {
 820         for (int i=0; i<N; i++)
 821                 Matrix[i*N + i] += shift;
 822  };
 823  //qques outils d'algebre des polynomes
 824
 825  /** \fn remplir_trinome
 826   * \brief remplie un trinome
 827   */
 828  void Polynoms::remplir_trinome(double coeff, double u, double alpha, vector< double >& trinome)
 829  {
 830         trinome[2] = coeff;
 831         coeff *= u;
 832         trinome[1] = -2*coeff;
 833         coeff *= u;
 834         trinome[0] = coeff - alpha;
 835  };
 836
 837
 838  /** \fn additionne
 839   * \brief Additionne 2 vecteurs de tailles differentes
 840   * @param P,Q les 2 vecteurs à additionner
 841   * @param n,m les 2 tailles respectives des vecteurs P et Q
 842   * @return un vecteur qui correspond a P+Q*/
 843  vector< double > Polynoms::additionne(const vector< double >& P, const vector< double >& Q, int n, int m)
 844  {
 845         if(m < n) {
 846                 vector< double > P_plus_Q = P;
 847                 for(int j = 0; j<m+1; j++){
 848                         P_plus_Q[j] += Q[j];
 849                 }
 850                 return P_plus_Q;
 851         }
 852         else{
 853                 vector< double > P_plus_Q = Q;
 854                 for(int j = 0; j<n+1; j++){
 855                         P_plus_Q[j] += P[j];
 856                 }
 857                 return P_plus_Q;
 858         }
 859  };
 860
 861
 862  /** \fn multiplie
 863   * \brief Calcule le produit scalaire de 2 vecteurs de tailles differentes
 864   * @param P,Q les 2 vecteurs
 865   * @param n,m les 2 tailles respectives des 2 vecteurs
 866   * @return un vecteur correspond à (P,Q) */
 867  vector< double > Polynoms::multiplie(const vector< double >& P, const vector< double >& Q, int n, int m)
 868  {
 869         vector< double > P_fois_Q (n+m+1,0.);
 870         for(int i = 0; i<n+1; i++){
 871                 for(int j = 0; j<m+1; j++){
 872                         P_fois_Q[i+j] += P[i]*Q[j];
 873                 }
 874         }
 875         return P_fois_Q;
 876  }
 877
 878
 879  /** \fn polynome_caracteristique
 880   * \brief Calcule le polynome caracteristique
 881   * @return */
 882  vector< double > Polynoms::polynome_caracteristique(double alpha_v, double alpha_l, double Uv, double Ul, double rhov, double rhol, double   cv_2, double cl_2, double dPiv, double dPil)//c= inverse vitesse du son!
 883  {
 884         vector< double > trinome_u(3);
 885         vector< double > trinome_c(3);
 886         vector< double > pol_char;
 887
 888         remplir_trinome(alpha_v*cv_2,Uv,alpha_v,trinome_c);
 889         remplir_trinome(rhol,Ul,dPil,trinome_u);
 890
 891         pol_char = multiplie(trinome_u, trinome_c, 2, 2);
 892
 893         remplir_trinome(alpha_l*cl_2,Ul,alpha_l,trinome_c);
 894         remplir_trinome(rhov,Uv,dPiv,trinome_u);
 895
 896         pol_char = additionne(pol_char,multiplie(trinome_u, trinome_c, 2, 2),4,4);
 897
 898         return pol_char;
 899  }
 900
 901  /** \fn polynome_caracteristique
 902   * \brief
 903   * @param alpha1
 904   * @param alpha2
 905   * @param u1
 906   * @param u2
 907   * @param rho1
 908   * @param rho2
 909   * @param invc1sq
 910   * @param invc2sq
 911   * @param dpi1
 912   * @param dpi2
 913   * @param g2press
 914   * @param g2alpha
 915   * @param g2
 916   * @param epsilon
 917   * @return */
 918  vector< double > Polynoms::polynome_caracteristique(double alpha1, double alpha2, double u1, double u2, double rho1, double rho2, double invc1sq, double invc2sq, double dpi1, double dpi2, double g2press, double g2alpha, double g2, double epsilon)
 919  {
 920         vector< double > trinome1(3);
 921         vector< double > trinome2(3);
 922         vector< double > pol_char;
 923         double K1,K2,K3;
 924
 925         /* (x-u1)^2(x-u2)^2-K1(x-u2)^2-K2(x-u1)^2+K3 */
 926         K1=alpha1*rho2*g2press + dpi1*alpha2*invc2sq*g2alpha;
 927         K2=alpha2*rho1*g2press + dpi2*alpha1*invc1sq*g2alpha;
 928         K3=(alpha1*dpi2+alpha2*dpi1)*g2alpha*g2press/g2;
 929         //cout<<"K1= " <<K1<<", K2= " <<K2<<", K3= " <<K3<<endl;
 930
 931         if(fabs(K1)>epsilon && fabs(K2)>epsilon )
 932         {
 933                 remplir_trinome(1/K1,u1,1,trinome1);
 934                 //cout<<"coeff constant trinome1= " << trinome1[0]<<endl;
 935                 remplir_trinome(1/K2,u2,1,trinome2);
 936                 //cout<<"coeff constant trinome2= " << trinome2[0]<<endl;
 937
 938                 pol_char = multiplie(trinome1, trinome2, 2, 2);
 939                 //cout<<"coeff constant produit= " << pol_char[0]<<endl;
 940
 941                 pol_char[0]+=K3/K2/K1-1;
 942                 //cout<<"coeff constant apres soustraction= " << pol_char[0]<<endl;
 943         }
 944         else
 945         {
 946                 remplir_trinome(1,u1,K1,trinome1);
 947                 remplir_trinome(1,u2,K2,trinome2);
 948
 949                 pol_char = multiplie(trinome1, trinome2, 2, 2);
 950
 951                 pol_char[0]+=K3-K2*K1;
 952         }
 953
 954         /*   for(int ct =0; ct<pol_char.size()-1; ct++)  */
 955         /*     if(fabs(pol_char[ct])< epsilon)  */
 956         /*       pol_char[ct]=0.; */
 957
 958         return pol_char;
 959  }
 960
 961  //Pour le tri des valeurs propres
 962
 963  /** \fn module
 964   * \brief calcule le module d'un nombre complexe
 965   * @param z est un nombre complexe
 966   * @return  calcule le module de z*/
 967  double Polynoms::module(complex< double > z)
 968  {
 969         return abs(z);
 970  }
 971
 972  /** \fn modulec calcule le module² de z */
 973  double Polynoms::modulec(complex< double > z)
 974  {
 975         return norm(z);
 976  }
 977
 978  /** \fn abs_generalise
 979   * \brief calcule la valeur absolue d'un nombre complexe
 980   * \Details calcule la valeur absolue d'un nombre complexe en prenant en compte
 981   * que la partie réelle c-a-d si z= a+ib abs_generalize() renvoie |a|+ib
 982   * @param z
 983   * @return si z = a+ib elle renvoie |a|+ib */
 984  complex< double > Polynoms::abs_generalise(complex< double > z)
 985  {
 986         if(z.real()>=0)
 987                 return z ;
 988         else
 989                 return -z ;
 990  }
 991  void Polynoms::ordre_croissant_abs(vector< complex< double > > &L, int n)
 992  {
 993         vector< complex< double > > copieL = L;
 994         int i_max;
 995         double Lmaxc;
 996         double Lmax_iterc;
 997
 998         for(int j=0; j<n; j++){
 999                 i_max = 0;
1000                 Lmaxc = copieL[0].real();
1001                 for (int i=1; i<n; i++){
1002                         Lmax_iterc = copieL[i].real();
1003                         if( Lmax_iterc > Lmaxc ){
1004                                 Lmaxc = Lmax_iterc;
1005                                 i_max = i;
1006                         }
1007                 }
1008                 L[n-1-j] = copieL[i_max];
1009                 copieL[i_max] = -INFINITY;
1010         }
1011  }
1012
1013  int Polynoms::new_tri_selectif(vector< complex< double > > &L, int n, double epsilon)
1014  {
1015         int size=1;
1016         /*
1017         cout<<"avant ordre croissant"<<endl;
1018         for(int i=0;i<n;i++)
1019                 cout<<L[i]<<", "<<endl;
1020          */
1021         ordre_croissant_abs(L,n);
1022         /*
1023         cout<<"après ordre croissant"<<endl;
1024         for(int i=0;i<n;i++)
1025                 cout<<L[i]<<", "<<endl;
1026          */
1027         vector< complex< double > >result(1,L[0]);
1028         for(int i=1;i<n;i++)
1029                 if(fabs(L[i].real()-result[size-1].real())>epsilon || fabs(L[i].imag())>epsilon)
1030                 {
1031                         result.push_back(L[i]);
1032                         size++;
1033                         /*
1034                         cout<<"tri selectif step "<< i<<endl;
1035                         for(int j=0;j<size;j++)
1036                                 cout<<result[j]<<", "<<endl;
1037                          */
1038                 }
1039         L=result;
1040         return size;
1041  }
1042
1043  /** \fn tri_selectif
1044   * \brief
1045   * @param L
1046   * @param n
1047   * @param epsilon
1048   * @return */
1049  template< typename T >
1050  int Polynoms::tri_selectif(T& L, int n, double epsilon)
1051  {
1052         int imin;
1053         int size_vp =n;
1054         int i=0;
1055         int j;
1056         complex< double > Lmin;
1057
1058         for(i = 0 ; i<size_vp ; i++){//On travaille sur le sous vecteur des indices de  i � size_vp-1
1059                 imin = i;
1060                 Lmin = L[i];
1061
1062                 for(j = i+1 ; j < size_vp ; j++)
1063                         if( Lmin.real() > L[j].real() )
1064                         {Lmin = L[j]; imin = j;}
1065
1066                 //on met la vp � zero si elle est trop petite
1067                 if( fabs(Lmin.real()) < epsilon) Lmin.real()=0;
1068                 if( fabs(Lmin.imag()) < epsilon) Lmin.imag()=0;
1069
1070                 switch (i) {
1071                 case 0: {
1072                         if (imin != i) {
1073                                 L[imin]=L[i];
1074                                 L[i]=Lmin;
1075                         }
1076                         break;
1077                 }
1078                 case 1: {
1079                         if ( module(L[i-1] - Lmin) > epsilon )  {
1080                                 if (imin != i) {
1081                                         L[imin]=L[i];
1082                                         L[i]=Lmin;
1083                                 }
1084
1085                         } else {
1086                                 L[imin]=L[size_vp - 1];
1087                                 L[size_vp - 1]=Lmin;
1088
1089                                 size_vp--;i--;
1090                         }
1091                         break;
1092                 }
1093                 default : {
1094                         if ( (module(L[i-1] - Lmin) > epsilon) &&  (module(L[i-2] - Lmin) > epsilon) )  {
1095                                 if (imin != i) {
1096                                         L[imin]=L[i];
1097                                         L[i]=Lmin;
1098                                 }
1099
1100                         } else {
1101                                 L[imin]=L[size_vp - 1];
1102                                 L[size_vp - 1]=Lmin;
1103
1104                                 size_vp--;i--;
1105                         }
1106                 }
1107                 }
1108         }
1109
1110         return size_vp;
1111  }
1112
1113  //Calcul des coefficients du polynome d'interpolation x->y dans la base de Newton
1114  template<class T>
1115
1116  /** \fn dif_div
1117   * \brief
1118   * @param n
1119   * @param x
1120   * @param y
1121   * @param epsilon
1122   * @return */
1123  vector< complex< double > > Polynoms::dif_div(int n, const vector< complex< double > >& x, const vector< T >& y, double epsilon)//n=nb valeurs à interpoler
1124  {
1125         complex< double > tab_dif_div[n*n];//La matrice de taille n*n
1126         vector< complex< double > > dif_div (n);
1127         int i, j;
1128         complex< double > delta_x;
1129
1130         for( i=0 ; i<n ; i++) {
1131                 tab_dif_div[i*n]=y[i];
1132                 //   cerr << tab_dif_div[i*n+0] << endl;
1133         }
1134         for(j=1; j<n ; j++){
1135                 for(i=0;i < n-j;i++) {
1136                         delta_x = x[i+j]-x[i];
1137                         if(module(delta_x) < epsilon) {
1138                                 cout << " attention, dif_div, delta_x <epsilon, " << " delta_x= " << delta_x <<" epsilon= "<<epsilon<<endl;
1139                                 cout<<" vp= " <<endl;
1140                                 for(int k=0; k<n; k++)
1141                                         cout << x[k]<<", " << endl;
1142                         }
1143                         tab_dif_div[i*n+j]=(tab_dif_div[(i+1)*n+j-1]-tab_dif_div[i*n+j-1])/delta_x;
1144                 }
1145         }
1146         for(i=0 ; i<n ; i++)
1147                 dif_div[i] = tab_dif_div[i*n+n-1-i];
1148         return dif_div;
1149  }
1150  //attention, Les vp complexes conjugu�es doivent se suivre dans la liste x
1151  void Polynoms::appliquer_dif_div(int n, const vector< complex< double > >& dif_div, const vector< complex< double > >& x, const double* A, const int sizeA, double epsilon, double *p)//p=result
1152  {
1153         int i, debut = 1;
1154         double Id[sizeA*sizeA], aux[sizeA*sizeA], matProd[sizeA*sizeA];
1155
1156         for(int k=0; k<sizeA*sizeA; k++)
1157                 Id[k]=0;
1158         for(int k=0; k<sizeA; k++)
1159                 Id[k*sizeA+k]=1;
1160         /*   for(int k=0; k<sizeA*sizeA; k++) */
1161         /*     aux[k] = A[k] ; */
1162         //cerr << x << endl;
1163         if( fabs(x[0].imag()) > epsilon){
1164                 for(int k=0; k<sizeA*sizeA; k++)
1165                         p[k] = A[k] * dif_div[0].real() ;
1166                 shift_diagonal(p, sizeA,(dif_div[1] - dif_div[0]*x[1]).real());
1167                 debut = 2;
1168         }else{
1169                 for(int k=0; k<sizeA*sizeA; k++)
1170                         p[k] = Id[k] * dif_div[0].real() ;
1171         }
1172         for(i=debut ; i<n ; i++){
1173                 for(int k=0; k<sizeA*sizeA; k++)
1174                         aux[k] = A[k];
1175                 //cerr << " on traite la racine reele " << x[i] << endl;
1176                 if ( fabs(x[i].imag()) < epsilon ) {
1177                         shift_diagonal(aux, sizeA, (- x[i].real()));
1178                         matrixProduct( p, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1179                         for(int k=0; k<sizeA*sizeA; k++)
1180                                 p[k] = matProd[k] ;
1181                         shift_diagonal(p, sizeA, dif_div[i].real());
1182                 }
1183                 else{ //couple de conjugees
1184                         if(fabs(x[i].imag() + x[i+1].imag()) < epsilon){
1185                                 matrixProduct( aux, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1186                                 for(int k=0; k<sizeA*sizeA; k++)
1187                                         aux[k] = matProd[k];
1188                                 for(int k=0; k<sizeA*sizeA; k++)
1189                                         aux[k]+=(-2*x[i].real()) * A[k];
1190                                 shift_diagonal(aux, sizeA, modulec(x[i]));
1191
1192                                 matrixProduct( p, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1193                                 for(int k=0; k<sizeA*sizeA; k++)
1194                                         p[k] = matProd[k] ;
1195                                 for(int k=0; k<sizeA*sizeA; k++)
1196                                         p[k]+= dif_div[i].real() * A[k];
1197                                 shift_diagonal(p, sizeA, (dif_div[i+1] - dif_div[i]*x[i+1]).real());
1198                                 i++;
1199                         }
1200                         else{
1201                                 cout << " appliquer_dif_div : les racines complexes conjuguees ne se suivent pas " <<  endl;
1202                                 for(int k=0; k<n; k++)
1203                                         cout << x[k]<<", " << endl;
1204                                 throw(" appliquer_dif_div : les racines complexes conjuguees ne se suivent pas ");
1205                         }
1206                 }
1207         }
1208  }
1209
1210  void Polynoms::abs_par_interp_directe(int n,  vector< complex< double > > vp,   double * A, int sizeA, double epsilon, double *result,vector< complex< double > > y)
1211  {
1212         vector< complex< double > > differences_divisees = dif_div(n, vp, y, epsilon);
1213         /*
1214         cout<<"valeurs propres"<<endl;
1215         for( int i=0 ; i<n ; i++)
1216                 cout<<vp[i] <<", "<<endl;
1217         cout<<"valeurs à atteindre"<<endl;
1218         for( int i=0 ; i<n ; i++)
1219                 cout<<y[i] <<", "<<endl;
1220         cout<<"differences_divisees= "<<endl;
1221         for(int i=0 ; i<n ; i++)
1222                 cout<<differences_divisees[i]<<", ";
1223         cout<<endl;
1224          */
1225         appliquer_dif_div(n, differences_divisees, vp, A, sizeA, epsilon, result);
1226  }
1227
1228  bool Polynoms::belongTo(complex< double > a , vector< complex <double > > v, double epsilon)
1229  {
1230         bool result = false;
1231         int size = v.size();
1232         int i = 0;
1233         double epsilon2 = epsilon*epsilon;
1234
1235         while(i<size && !result)
1236         {
1237                 result = modulec(a-v[i]) < epsilon2;
1238                 i++;
1239         }
1240         return result;
1241  }
1242  double Polynoms::avr(double a, double b)
1243  {
1244          return 0.5*(a+b);
1245  }