src/math/math_FunctionRoots.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math_DirectPolynomialRoots.hxx>
  23 #include <math_FunctionRoots.hxx>
  24 #include <math_FunctionWithDerivative.hxx>
  25 #include <math_BracketedRoot.hxx>
  26 #include <TColStd_Array1OfReal.hxx>
  27
  28 #define ITMAX  100
  29 #define EPS    1e-14
  30 #define EPSEPS 2e-14
  31 #define MAXBIS 100
  32
  33 #ifdef OCCT_DEBUG
  34 static Standard_Boolean myDebug = 0;
  35 static Standard_Integer nbsolve = 0;
  36 #endif
  37
  38 class DerivFunction: public math_Function
  39 {
  40   math_FunctionWithDerivative *myF;
  41
  42 public:
  43   DerivFunction(math_FunctionWithDerivative& theF)
  44     : myF (&theF)
  45   {}
  46
  47   virtual Standard_Boolean Value(const Standard_Real theX, Standard_Real& theFval)
  48   {
  49     return myF->Derivative(theX, theFval);
  50   }
  51 };
  52
  53 static void  AppendRoot(TColStd_SequenceOfReal& Sol,
  54                         TColStd_SequenceOfInteger& NbStateSol,
  55                         const Standard_Real X,
  56                         math_FunctionWithDerivative& F,
  57 //                      const Standard_Real K,
  58                         const Standard_Real ,
  59                         const Standard_Real dX) {
  60
  61   Standard_Integer n=Sol.Length();
  62   Standard_Real t;
  63 #ifdef OCCT_DEBUG
  64  if (myDebug) {
  65    std::cout << "   Ajout de la solution numero : " << n+1 << std::endl;
  66    std::cout << "   Valeur de la racine :" << X << std::endl;
  67  }
  68 #endif
  69
  70   if(n==0) {
  71     Sol.Append(X);
  72     F.Value(X,t);
  73     NbStateSol.Append(F.GetStateNumber());
  74   }
  75   else {
  76     Standard_Integer i=1;
  77     Standard_Integer pl=n+1;
  78     while(i<=n) {
  79       t=Sol.Value(i);
  80       if(t >= X) {
  81         pl=i;
  82         i=n;
  83       }
  84       if(Abs(X-t) <= dX) {
  85         pl=0;
  86         i=n;
  87       }
  88       i++;
  89     } //-- while
  90     if(pl>n) {
  91       Sol.Append(X);
  92       F.Value(X,t);
  93       NbStateSol.Append(F.GetStateNumber());
  94     }
  95     else if(pl>0) {
  96       Sol.InsertBefore(pl,X);
  97       F.Value(X,t);
  98       NbStateSol.InsertBefore(pl,F.GetStateNumber());
  99     }
 100   }
 101 }
 102
 103 static void  Solve(math_FunctionWithDerivative& F,
 104                    const Standard_Real K,
 105                    const Standard_Real x1,
 106                    const Standard_Real y1,
 107                    const Standard_Real x2,
 108                    const Standard_Real y2,
 109                    const Standard_Real tol,
 110                    const Standard_Real dX,
 111                    TColStd_SequenceOfReal& Sol,
 112                    TColStd_SequenceOfInteger& NbStateSol) {
 113 #ifdef OCCT_DEBUG
 114   if (myDebug) {
 115     std::cout <<"--> Resolution :" << ++nbsolve << std::endl;
 116     std::cout <<"   x1 =" << x1 << " y1 =" << y1 << std::endl;
 117     std::cout <<"   x2 =" << x2 << " y2 =" << y2 << std::endl;
 118   }
 119 #endif
 120
 121   Standard_Integer iter=0;
 122   Standard_Real tols2 = 0.5*tol;
 123   Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
 124   a=x1;b=c=x2;fa=y1; fb=fc=y2;
 125   for(iter=1;iter<=ITMAX;iter++) {
 126     if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
 127       c=a; fc=fa; e=d=b-a;
 128     }
 129     if(Abs(fc)<Abs(fb)) {
 130       a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
 131     }
 132     tol1 = EPSEPS*Abs(b)+tols2;
 133     xm=0.5*(c-b);
 134     if(Abs(xm)<tol1 || fb==0) {
 135       //-- On tente une iteration de newton
 136       Standard_Real Xp,Yp,Dp;
 137       Standard_Integer itern=5;
 138       Standard_Boolean Ok;
 139       Xp=b;
 140       do {
 141         Ok = F.Values(Xp,Yp,Dp);
 142         if(Ok) {
 143           Ok=Standard_False;
 144           if(Dp>1e-10 || Dp<-1e-10) {
 145             Xp = Xp-(Yp-K)/Dp;
 146           }
 147           if(Xp<=x2 && Xp>=x1) {
 148             F.Value(Xp,Yp); Yp-=K;
 149             if(Abs(Yp)<Abs(fb)) {
 150               b=Xp;
 151               fb=Yp;
 152               Ok=Standard_True;
 153             }
 154           }
 155         }
 156       }
 157       while(Ok && --itern>=0);
 158
 159       AppendRoot(Sol,NbStateSol,b,F,K,dX);
 160       return;
 161     }
 162     if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
 163       s=fb/fa;
 164       if(a==c) {
 165         p=xm*s; p+=p;
 166         q=1.0-s;
 167       }
 168       else {
 169         q=fa/fc;
 170         r=fb/fc;
 171         p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
 172         q=(q-1.0)*(r-1.0)*(s-1.0);
 173       }
 174       if(p>0.0) {
 175         q=-q;
 176       }
 177       p=Abs(p);
 178       min1=3.0*xm*q-Abs(tol1*q);
 179       min2=Abs(e*q);
 180       if((p+p)<( (min1<min2) ? min1 : min2)) {
 181         e=d;
 182         d=p/q;
 183       }
 184       else {
 185         d=xm;
 186         e=d;
 187       }
 188     }
 189     else {
 190       d=xm;
 191       e=d;
 192     }
 193     a=b;
 194     fa=fb;
 195     if(Abs(d)>tol1) {
 196       b+=d;
 197     }
 198     else {
 199       if(xm>=0) b+=Abs(tol1);
 200       else      b+=-Abs(tol1);
 201     }
 202     F.Value(b,fb);
 203     fb-=K;
 204   }
 205 #ifdef OCCT_DEBUG
 206   std::cout<<" Non Convergence dans math_FunctionRoots.cxx "<<std::endl;
 207 #endif
 208 }
 209
 210 #define NEWSEQ 1
 211
 212 #define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement
 213 //#define MATH_FUNCTIONROOTS_OLDCODE // Ancien
 214 //#define MATH_FUNCTIONROOTS_CHECK // Check
 215
 216 math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
 217                                        const Standard_Real A,
 218                                        const Standard_Real B,
 219                                        const Standard_Integer NbSample,
 220                                        const Standard_Real _EpsX,
 221                                        const Standard_Real EpsF,
 222                                        const Standard_Real EpsNull,
 223                                        const Standard_Real K )
 224 {
 225 #ifdef OCCT_DEBUG
 226   if (myDebug) {
 227     std::cout << "---- Debut de math_FunctionRoots ----" << std::endl;
 228     nbsolve = 0;
 229   }
 230 #endif
 231
 232 #if NEWSEQ
 233   TColStd_SequenceOfReal StaticSol;
 234 #endif
 235   Sol.Clear();
 236   NbStateSol.Clear();
 237   #ifdef MATH_FUNCTIONROOTS_NEWCODE
 238     {
 239     Done = Standard_True;
 240     Standard_Real X0=A;
 241     Standard_Real XN=B;
 242     Standard_Integer N=NbSample;
 243     //-- ------------------------------------------------------------
 244     //-- Verifications de bas niveau
 245     if(B<A) {
 246       X0=B;
 247       XN=A;
 248     }
 249     N*=2;
 250     if(N < 20) {
 251       N=20;
 252     }
 253     //--  On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
 254     Standard_Real EpsX   = _EpsX;
 255     Standard_Real DeltaU = Abs(X0)+Abs(XN);
 256     Standard_Real NEpsX  = 0.0000000001 * DeltaU;
 257     if(EpsX < NEpsX) {
 258       EpsX = NEpsX;
 259     }
 260
 261     //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
 262     //-- A .............................................................. B
 263     //-- X0   X1   X2 ........................................  Xn-1      Xn
 264     Standard_Integer i;
 265     Standard_Real X=X0;
 266     Standard_Boolean Ok;
 267     double dx = (XN-X0)/N;
 268     TColStd_Array1OfReal ptrval(0, N);
 269     Standard_Integer Nvalid = -1;
 270     Standard_Real aux = 0;
 271     for(i=0; i<=N ; i++,X+=dx) {
 272       if( X > XN) X=XN;
 273       Ok=F.Value(X,aux);
 274       if(Ok) ptrval(++Nvalid) = aux - K;
 275 //      ptrval(i)-=K;
 276     }
 277     //-- Toute la fonction est nulle ?
 278
 279     if( Nvalid < N ) {
 280       Done = Standard_False;
 281       return;
 282     }
 283
 284     AllNull=Standard_True;
 285 //    for(i=0;AllNull && i<=N;i++) {
 286     for(i=0;AllNull && i<=N;i++) {
 287       if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
 288         AllNull=Standard_False;
 289       }
 290     }
 291     if(AllNull) {
 292       //-- tous les points echantillons sont dans la tolerance
 293
 294     }
 295     else {
 296       //-- Il y a des points hors tolerance
 297       //-- on detecte les changements de signes STRICTS
 298       Standard_Integer ip1;
 299 //      Standard_Boolean chgtsign=Standard_False;
 300       Standard_Real tol = EpsX;
 301       Standard_Real X2;
 302       for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
 303         X2=X+dx;
 304         if(X2 > XN) X2 = XN;
 305         if(ptrval(i)<0.0) {
 306           if(ptrval(ip1)>0.0) {
 307             //-- --------------------------------------------------
 308             //-- changement de signe dans Xi Xi+1
 309             Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
 310           }
 311         }
 312         else {
 313           if(ptrval(ip1)<0.0) {
 314             //-- --------------------------------------------------
 315             //-- changement de signe dans Xi Xi+1
 316             Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
 317           }
 318         }
 319       }
 320       //-- On detecte les cas ou la fct s annule sur des Xi et est
 321       //-- non nulle au voisinage de Xi
 322       //--
 323       //-- On prend 2 points u0,u1 au voisinage de Xi
 324       //-- Si (F(u0)-K)*(F(u1)-K) <0   on lance une recherche
 325       //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
 326       for(i=0; i<=N; i++) {
 327         if(ptrval(i)==0) {
 328 //        Standard_Real Val,Deriv;
 329           X=X0+i*dx;
 330           if (X>XN) X=XN;
 331           Standard_Real u0,u1;
 332           u0=dx*0.5;      u1=X+u0;        u0+=X;
 333           if(u0<X0)  u0=X0;
 334           if(u0>XN)  u0=XN;
 335           if(u1<X0)  u1=X0;
 336           if(u1>XN)  u1=XN;
 337
 338           Standard_Real y0,y1;
 339           F.Value(u0,y0); y0-=K;
 340           F.Value(u1,y1); y1-=K;
 341           if(y0*y1 < 0.0) {
 342             Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
 343           }
 344           else {
 345             if(y0!=0.0 || y1!=0.0) {
 346               AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
 347             }
 348           }
 349         }
 350       }
 351       //-- --------------------------------------------------------------------------------
 352       //-- Il faut traiter differement le cas des points en bout :
 353       if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
 354         AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
 355       }
 356       if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
 357         AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
 358       }
 359
 360       //-- --------------------------------------------------------------------------------
 361       //-- --------------------------------------------------------------------------------
 362       //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
 363       //--                                                          un maximum avec f(x)<0
 364       //-- On reprend une discretisation plus fine au voisinage de ces extremums
 365       //--
 366       //-- Recherche d un minima positif
 367       Standard_Real xm,ym,dym,xm1,xp1;
 368       Standard_Real majdx = 5.0*dx;
 369       Standard_Boolean Rediscr;
 370 //      Standard_Real ptrvalbis[MAXBIS];
 371       Standard_Integer im1=0;
 372       ip1=2;
 373       for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
 374         Rediscr = Standard_False;
 375         if (xm > XN) xm=XN;
 376         if(ptrval(i)>0.0) {
 377           if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
 378             //-- Peut on traverser l axe Ox
 379             //-- -------------- Estimation a partir de Xim1
 380             xm1=xm-dx;
 381             if(xm1 < X0) xm1=X0;
 382             F.Values(xm1,ym,dym); ym-=K;
 383             if(dym<-1e-10 || dym>1e-10) {  // normalement dym < 0
 384               Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
 385               if(t<majdx && t > -majdx) {
 386                 Rediscr = Standard_True;
 387               }
 388             }
 389             //-- -------------- Estimation a partir de Xip1
 390             if(Rediscr==Standard_False) {
 391               xp1=xm+dx;
 392               if(xp1 > XN ) xp1=XN;
 393               F.Values(xp1,ym,dym); ym-=K;
 394               if(dym<-1e-10 || dym>1e-10) {  // normalement dym > 0
 395                 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
 396                 if(t<majdx && t > -majdx) {
 397                   Rediscr = Standard_True;
 398                 }
 399               }
 400             }
 401           }
 402         }
 403         else if(ptrval(i)<0.0) {
 404           if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
 405             //-- Peut on traverser l axe Ox
 406             //-- -------------- Estimation a partir de Xim1
 407             xm1=xm-dx;
 408             if(xm1 < X0) xm1=X0;
 409             F.Values(xm1,ym,dym); ym-=K;
 410             if(dym>1e-10 || dym<-1e-10) {  // normalement dym > 0
 411               Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
 412               if(t<majdx && t > -majdx) {
 413                 Rediscr = Standard_True;
 414               }
 415             }
 416             //-- -------------- Estimation a partir de Xim1
 417             if(Rediscr==Standard_False) {
 418                xm1=xm-dx;
 419                if(xm1 < X0) xm1=X0;
 420               F.Values(xm1,ym,dym); ym-=K;
 421               if(dym>1e-10 || dym<-1e-10) {  // normalement dym < 0
 422                 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
 423                 if(t<majdx && t > -majdx) {
 424                   Rediscr = Standard_True;
 425                 }
 426               }
 427             }
 428           }
 429         }
 430         if(Rediscr) {
 431           Standard_Real x0 = xm - dx;
 432           Standard_Real x3 = xm + dx;
 433           if (x0 < X0) x0 = X0;
 434           if (x3 > XN) x3 = XN;
 435           Standard_Real aSolX1 = 0., aSolX2 = 0.;
 436           Standard_Real aVal1 = 0., aVal2 = 0.;
 437           Standard_Real aDer1 = 0., aDer2 = 0.;
 438           Standard_Boolean isSol1 = Standard_False;
 439           Standard_Boolean isSol2 = Standard_False;
 440           //-- ----------------------------------------------------
 441           //-- Find minimum of the function |F| between x0 and x3
 442           //-- by searching for the zero of the function derivative
 443           DerivFunction aDerF(F);
 444           math_BracketedRoot aBR(aDerF, x0, x3, _EpsX);
 445           if (aBR.IsDone())
 446           {
 447             aSolX1 = aBR.Root();
 448             F.Value(aSolX1, aVal1);
 449             aVal1 = Abs(aVal1);
 450             if (aVal1 < EpsF)
 451             {
 452               isSol1 = Standard_True;
 453               aDer1 = aBR.Value();
 454             }
 455           }
 456
 457           //-- --------------------------------------------------
 458           //-- On recherche un extrema entre x0 et x3
 459           //-- x1 et x2 sont tels que x0<x1<x2<x3
 460           //-- et |f(x0)| > |f(x1)|   et |f(x3)| > |f(x2)|
 461           //--
 462           //-- En entree : a=xm-dx  b=xm c=xm+dx
 463           Standard_Real x1, x2, f0, f3;
 464           Standard_Real R = 0.61803399;
 465           Standard_Real C = 1.0 - R;
 466           Standard_Real tolCR = NEpsX*10.0;
 467           f0 = ptrval(im1);
 468           f3 = ptrval(ip1);
 469           Standard_Boolean recherche_minimum = (f0 > 0.0);
 470
 471           if (Abs(x3 - xm) > Abs(x0 - xm)) { x1 = xm; x2 = xm + C*(x3 - xm); }
 472           else                        { x2 = xm; x1 = xm - C*(xm - x0); }
 473           Standard_Real f1, f2;
 474           F.Value(x1, f1); f1 -= K;
 475           F.Value(x2, f2); f2 -= K;
 476           //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
 477           Standard_Real tolX = 0.001 * NEpsX;
 478           while (Abs(x3 - x0) > tolCR*(Abs(x1) + Abs(x2)) && (Abs(x1 - x2) > tolX)) {
 479             //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
 480             //--    x0,f0,x1,f1,x2,f2,x3,f3);
 481             if (recherche_minimum) {
 482               if (f2 < f1) {
 483                 x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
 484                 f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
 485               }
 486               else {
 487                 x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
 488                 f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
 489               }
 490             }
 491             else {
 492               if (f2 > f1) {
 493                 x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
 494                 f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
 495               }
 496               else {
 497                 x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
 498                 f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
 499               }
 500             }
 501             //-- On ne fait pas que chercher des extremas. Il faut verifier
 502             //-- si on ne tombe pas sur une racine
 503             if (f1*f0 < 0.0) {
 504               //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
 505               Solve(F, K, x0, f0, x1, f1, tol, NEpsX, Sol, NbStateSol);
 506             }
 507             if (f2*f3 < 0.0)  {
 508               //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
 509               Solve(F, K, x2, f2, x3, f3, tol, NEpsX, Sol, NbStateSol);
 510             }
 511           }
 512           if ((recherche_minimum && f1<f2) || (!recherche_minimum && f1>f2)) {
 513             //-- x1,f(x1) minimum
 514             if (Abs(f1) < EpsF) {
 515               isSol2 = Standard_True;
 516               aSolX2 = x1;
 517               aVal2 = Abs(f1);
 518             }
 519           }
 520           else {
 521             //-- x2.f(x2) minimum
 522             if (Abs(f2) < EpsF) {
 523               isSol2 = Standard_True;
 524               aSolX2 = x2;
 525               aVal2 = Abs(f2);
 526             }
 527           }
 528           // Choose the best solution between aSolX1, aSolX2
 529           if (isSol1 && isSol2)
 530           {
 531             if (aVal2 - aVal1 > EpsF)
 532               AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
 533             else if (aVal1 - aVal2 > EpsF)
 534               AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
 535             else
 536             {
 537               aDer1 = Abs(aDer1);
 538               F.Derivative(aSolX2, aDer2);
 539               aDer2 = Abs(aDer2);
 540               if (aDer1 < aDer2)
 541                 AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
 542               else
 543                 AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
 544             }
 545           }
 546           else if (isSol1)
 547             AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
 548           else if(isSol2)
 549             AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
 550         } //-- Recherche d un extrema
 551       } //-- for
 552     }
 553
 554 #if NEWSEQ
 555     #ifdef MATH_FUNCTIONROOTS_CHECK
 556     {
 557       StaticSol.Clear();
 558       Standard_Integer n=Sol.Length();
 559       for(Standard_Integer ii=1;ii<=n;ii++) {
 560         StaticSol.Append(Sol.Value(ii));
 561       }
 562       Sol.Clear();
 563       NbStateSol.Clear();
 564     }
 565 #endif
 566 #endif
 567 #endif
 568   }
 569 #ifdef MATH_FUNCTIONROOTS_OLDCODE
 570 {
 571     //-- ********************************************************************************
 572     //--                              ANCIEN TRAITEMENT
 573     //-- ********************************************************************************
 574
 575
 576   // calculate all the real roots of a function within the range
 577   // A..B. without condition on A and B
 578   // a solution X is found when
 579   //   abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
 580   // The function is considered as null between A and B if
 581   // abs(F-K) <= EpsNull within this range.
 582   Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
 583                               //-- Il ne faut pas EpsX = 0.000...001  car dans ce cas
 584                               //-- U + Nn*EpsX     ==     U
 585   Standard_Real Lowr,Upp;
 586   Standard_Real Increment;
 587   Standard_Real Null2;
 588   Standard_Real FLowr,FUpp,DFLowr,DFUpp;
 589   Standard_Real U,Xu;
 590   Standard_Real Fxu,DFxu,FFxu,DFFxu;
 591   Standard_Real Fyu,DFyu,FFyu,DFFyu;
 592   Standard_Boolean Finish;
 593   Standard_Real FFi;
 594   Standard_Integer Nbiter = 30;
 595   Standard_Integer Iter;
 596   Standard_Real Ambda,T;
 597   Standard_Real AA,BB,CC;
 598   Standard_Integer Nn;
 599   Standard_Real Alfa1=0,Alfa2;
 600   Standard_Real OldDF = RealLast();
 601   Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
 602   Standard_Boolean Ok;
 603
 604   Done = Standard_False;
 605
 606   StdFail_NotDone_Raise_if(NbSample <= 0, " ");
 607
 608   // initialisation
 609
 610   if (A > B) {
 611     Lowr=B;
 612     Upp=A;
 613   }
 614   else {
 615     Lowr=A;
 616     Upp=B;
 617   }
 618
 619   Increment = (Upp-Lowr)/NbSample;
 620   StdFail_NotDone_Raise_if(Increment < EpsX, " ");
 621   Done = Standard_True;
 622   //--  On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
 623   Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
 624   Standard_Real NEpsX  = 0.0000000001 * DeltaU;
 625   if(EpsX < NEpsX) {
 626     EpsX = NEpsX;
 627     //-- std::cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" )   EpsX = "<<EpsX<<std::endl;
 628   }
 629   //--
 630   Null2 = EpsNull*EpsNull;
 631
 632   Ok = F.Values(Lowr,FLowr,DFLowr);
 633
 634   if(!Ok) {
 635     Done = Standard_False;
 636     return;
 637   }
 638
 639   FLowr = FLowr - K;
 640
 641   Ok = F.Values(Upp,FUpp,DFUpp);
 642
 643   if(!Ok) {
 644     Done = Standard_False;
 645     return;
 646   }
 647
 648   FUpp = FUpp - K;
 649
 650   // Calcul sur U
 651
 652   U = Lowr-EpsX;
 653   Fyu = FLowr-EpsX*DFLowr;  // extrapolation lineaire
 654   DFyu = DFLowr;
 655   FFyu = Fyu*Fyu;
 656   DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
 657   AllNull = ( FFyu <= Null2 );
 658
 659   while ( U < Upp) {
 660
 661     Xu = U;
 662     Fxu = Fyu;
 663     DFxu = DFyu;
 664     FFxu = FFyu;
 665     DFFxu = DFFyu;
 666
 667     U = Xu + Increment;
 668     if (U <= Lowr) {
 669       Fyu = FLowr + (U-Lowr)*DFLowr;
 670       DFyu = DFLowr;
 671     }
 672     else if (U >= Upp) {
 673       Fyu = FUpp + (U-Upp)*DFUpp;
 674       DFyu = DFUpp;
 675     }
 676     else {
 677       Ok = F.Values(U,Fyu,DFyu);
 678
 679       if(!Ok) {
 680         Done = Standard_False;
 681         return;
 682       }
 683
 684       Fyu = Fyu - K;
 685     }
 686     FFyu = Fyu*Fyu;
 687     DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
 688
 689     if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
 690
 691       if (AllNull) {     //rechercher les vraix zeros depuis le debut
 692
 693         AllNull = Standard_False;
 694         Xu = Lowr-EpsX;
 695         Fxu = FLowr-EpsX*DFLowr;
 696         DFxu = DFLowr;
 697         FFxu = Fxu*Fxu;
 698         DFFxu = Fxu*DFxu; DFFxu+=DFFxu;  //-- DFFxu = 2.*Fxu*DFxu;
 699         U = Xu + Increment;
 700         Ok = F.Values(U,Fyu,DFyu);
 701
 702         if(!Ok) {
 703           Done = Standard_False;
 704           return;
 705         }
 706
 707         Fyu = Fyu - K;
 708         FFyu = Fyu*Fyu;
 709         DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
 710       }
 711       Standard_Real FxuFyu=Fxu*Fyu;
 712
 713       if (  (DFFyu > 0. && DFFxu <= 0.)
 714           || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
 715           || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
 716           || (FxuFyu <= 0.) ) {
 717         // recherche d 1 minimun possible
 718         Finish = Standard_False;
 719         Ambda = Increment;
 720         T = 0.;
 721         Iter=0;
 722         FFi=0.;
 723
 724         if (FxuFyu >0.) {
 725           // chercher si f peut s annuler pour eviter
 726           //  des iterations inutiles
 727           if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
 728               Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
 729
 730             Finish = Standard_True;
 731             FFi = Min ( FFxu , FFyu);  //pour ne pas recalculer yu
 732           }
 733           else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
 734                    (FFxu <= Standard_Underflow  && -FFxu  <= Standard_Underflow)) {
 735
 736             Finish = Standard_True;
 737             FFxu = 0.0;
 738             FFi = FFyu;   // pour recalculer yu
 739           }
 740           else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
 741                    (FFyu  <= Standard_Underflow && -FFyu  <= Standard_Underflow)) {
 742
 743             Finish = Standard_True;
 744             FFyu =0.0;
 745             FFi = FFxu;   // pour recalculer U
 746           }
 747         }
 748         else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
 749
 750           Finish = Standard_True;
 751           FFxu = 0.0;
 752           FFi = FFyu;
 753         }
 754         else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
 755
 756           Finish = Standard_True;
 757           FFyu =0.0;
 758           FFi = FFxu;
 759         }
 760         while (!Finish) {
 761
 762           // calcul des 2 solutions annulant la derivee de l interpolation cubique
 763           //    Ambda*t=(U-Xu)  F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
 764           //    df=aa*t*t+2*bb*t+cc
 765
 766           AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
 767           BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
 768           CC = Ambda*DFFxu;
 769
 770           if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
 771             AA=BB=CC=0;
 772           }
 773
 774
 775           math_DirectPolynomialRoots Solv(AA,BB,CC);
 776           if ( !Solv.InfiniteRoots() ) {
 777             Nn=Solv.NbSolutions();
 778             if (Nn <= 0) {
 779               if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
 780               else Finish = Standard_True;
 781             }
 782             else {
 783               Alfa1 = Solv.Value(1);
 784               if (Nn == 2) {
 785                 Alfa2 = Solv.Value(2);
 786                 if (Alfa1 > Alfa2){
 787                   AA = Alfa1;
 788                   Alfa1 = Alfa2;
 789                   Alfa2 = AA;
 790                 }
 791                 if (Alfa1 > 1. || Alfa2 < 0.){
 792                   // resolution par dichotomie
 793                   if (Fxu*Fyu < 0.) Alfa1 = 0.5;
 794                   else  Finish = Standard_True;
 795                 }
 796                 else if ( Alfa1 < 0. ||
 797                          ( DFFxu > 0. && DFFyu >= 0.) ) {
 798                   // si 2 derivees >0
 799                   //(cas changement de signe de la distance signee sans
 800                   // changement de signe de la derivee:
 801                   //cas de 'presque'tangence avec 2
 802                   // solutions proches) ,on prend la plus grane racine
 803                   if (Alfa2 > 1.) {
 804                     if (Fxu*Fyu < 0.) Alfa1 = 0.5;
 805                     else Finish = Standard_True;
 806                   }
 807                   else Alfa1 = Alfa2;
 808                 }
 809               }
 810             }
 811           }
 812           else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
 813           //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
 814           else  Finish = Standard_True;
 815
 816           if (!Finish) {
 817             // petits tests pour diminuer le nombre d iterations
 818             if (Alfa1 <= EpsX ) {
 819               Alfa1+=Alfa1;
 820             }
 821             else if (Alfa1 >= (1.-EpsX) ) {
 822               Alfa1 = Alfa1+Alfa1-1.;
 823             }
 824             Alfa1 = Ambda*(1. - Alfa1);
 825             U = U + T - Alfa1;
 826             if ( U <= Lowr ) {
 827               AA = FLowr + (U - Lowr)*DFLowr;
 828               BB = DFLowr;
 829             }
 830             else if ( U >= Upp ) {
 831               AA = FUpp + (U - Upp)*DFUpp;
 832               BB = DFUpp;
 833             }
 834             else {
 835               Ok = F.Values(U,AA,BB);
 836
 837               if(!Ok) {
 838                 Done = Standard_False;
 839                 return;
 840               }
 841
 842               AA = AA - K;
 843             }
 844             FFi = AA*AA;
 845             CC = AA*BB; CC+=CC;
 846             if (  ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
 847                 && AA*Fxu > 0. ) {
 848               FFxu = FFi;
 849               DFFxu = CC;
 850               Fxu = AA;
 851               DFxu = BB;
 852               T = Alfa1;
 853               if (Alfa1 > Ambda*0.5) {
 854                 // remarque (1)
 855                 // determination d 1 autre borne pour diviser
 856                 //le nouvel intervalle par 2 au -
 857                 Xu = U + Alfa1*0.5;
 858                 if ( Xu <= Lowr ) {
 859                   AA = FLowr + (Xu - Lowr)*DFLowr;
 860                   BB = DFLowr;
 861                 }
 862                 else if ( Xu >= Upp ) {
 863                   AA = FUpp + (Xu - Upp)*DFUpp;
 864                   BB = DFUpp;
 865                 }
 866                 else {
 867                   Ok = F.Values(Xu,AA,BB);
 868
 869                   if(!Ok) {
 870                     Done = Standard_False;
 871                     return;
 872                   }
 873
 874                   AA = AA -K;
 875                 }
 876                 FFi = AA*AA;
 877                 CC = AA*BB; CC+=CC;
 878                 if  ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
 879                      || Fxu * AA < 0. ) {
 880                   Fyu = AA;
 881                   DFyu = BB;
 882                   FFyu = FFi;
 883                   DFFyu = CC;
 884                   T = Alfa1*0.5;
 885                   Ambda = Alfa1*0.5;
 886                 }
 887                 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
 888                   // changement de signe sur l intervalle u,U
 889                   Fxu = AA;
 890                   DFxu = BB;
 891                   FFxu = FFi;
 892                   DFFxu = CC;
 893                   FFi = Min(FFxu,FFyu);
 894                   T = Alfa1*0.5;
 895                   Ambda = Alfa1*0.5;
 896                   U = Xu;
 897                 }
 898                 else { Ambda = Alfa1;}
 899               }
 900               else { Ambda = Alfa1;}
 901             }
 902             else {
 903               Fyu = AA;
 904               DFyu = BB;
 905               FFyu = FFi;
 906               DFFyu = CC;
 907               if ( (Ambda-Alfa1) > Ambda*0.5 ) {
 908                 // meme remarque (1)
 909                 Xu = U - (Ambda-Alfa1)*0.5;
 910                 if ( Xu <= Lowr ) {
 911                   AA = FLowr + (Xu - Lowr)*DFLowr;
 912                   BB = DFLowr;
 913                 }
 914                 else if ( Xu >= Upp ) {
 915                   AA = FUpp + (Xu - Upp)*DFUpp;
 916                   BB = DFUpp;
 917                 }
 918                 else {
 919                   Ok = F.Values(Xu,AA,BB);
 920
 921                   if(!Ok) {
 922                     Done = Standard_False;
 923                     return;
 924                   }
 925
 926                   AA = AA - K;
 927                 }
 928                 FFi = AA*AA;
 929                 CC = AA*BB; CC+=CC;
 930                 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
 931                   FFxu = FFi;
 932                   DFFxu = CC;
 933                   Fxu = AA;
 934                   DFxu = BB;
 935                   Ambda = ( Ambda - Alfa1 )*0.5;
 936                   T = 0.;
 937                 }
 938                 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
 939                   FFyu = FFi;
 940                   DFFyu = CC;
 941                   Fyu = AA;
 942                   DFyu = BB;
 943                   Ambda = ( Ambda - Alfa1 )*0.5;
 944                   T = 0.;
 945                   FFi = Min(FFxu , FFyu);
 946                   U = Xu;
 947                 }
 948                 else {
 949                   T =0.;
 950                   Ambda = Ambda - Alfa1;
 951                 }
 952               }
 953               else {
 954                 T =0.;
 955                 Ambda = Ambda - Alfa1;
 956               }
 957             }
 958             // tests d arrets
 959             if (Abs(FFxu) <= Standard_Underflow ||
 960                 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
 961                 ){
 962               Finish = Standard_True;
 963               if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
 964               FFi = FFyu;
 965             }
 966             else if (Abs(FFyu) <= Standard_Underflow ||
 967                      (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
 968                      ){
 969               Finish = Standard_True;
 970               if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
 971               FFi = FFxu;
 972             }
 973             else {
 974               Iter = Iter + 1;
 975               Finish = Iter >= Nbiter || (Ambda <= EpsX &&
 976                                           ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
 977             }
 978           }
 979         }  // fin interpolation cubique
 980
 981         // restitution du meilleur resultat
 982
 983         if ( FFxu < FFi ) { U = U + T -Ambda;}
 984         else if ( FFyu < FFi ) { U = U + T;}
 985
 986         if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
 987           U = Max( Lowr , Min( U , Upp));
 988           Ok = F.Value(U,FFi);
 989           FFi = FFi - K;
 990           if ( Abs(FFi) < EpsF ) {
 991             //coherence
 992             if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
 993             else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
 994             else if (Fxu*Fyu > 0.) { AA = 0.;}
 995             else { AA = Fyu-Fxu;}
 996             if (!Sol.IsEmpty()) {
 997               if (Abs(Sol.Last() - U) > 5.*EpsX
 998                   || (OldDF != RealLast() && OldDF*AA < 0.) ) {
 999                 Sol.Append(U);
1000                 NbStateSol.Append(F.GetStateNumber());
1001               }
1002             }
1003             else {
1004               Sol.Append(U);
1005               NbStateSol.Append(F.GetStateNumber());
1006             }
1007             OldDF = AA ;
1008           }
1009         }
1010         DFFyu = 0.;
1011         Fyu = 0.;
1012         Nn = 1;
1013         while ( Nn < 1000000 && DFFyu <= 0.) {
1014           // on repart d 1 pouyem plus loin
1015           U = U + Nn*EpsX;
1016           if ( U <= Lowr ) {
1017             AA = FLowr + (U - Lowr)*DFLowr;
1018             BB = DFLowr;
1019           }
1020           else if ( U >= Upp ) {
1021             AA = FUpp + (U - Upp)*DFUpp;
1022             BB = DFUpp;
1023           }
1024           else {
1025             Ok = F.Values(U,AA,BB);
1026             AA = AA - K;
1027           }
1028           if ( AA*Fyu < 0.) {
1029             U = U - Nn*EpsX;
1030             Nn = 1000001;
1031           }
1032           else {
1033             Fyu = AA;
1034             DFyu = BB;
1035             FFyu = AA*AA;
1036             DFFyu = 2.*DFyu*Fyu;
1037             Nn = 2*Nn;
1038           }
1039         }
1040       }
1041     }
1042   }
1043 #if NEWSEQ
1044   #ifdef MATH_FUNCTIONROOTS_CHECK
1045   {
1046     Standard_Integer n1=StaticSol.Length();
1047     Standard_Integer n2=Sol.Length();
1048     if(n1!=n2) {
1049       printf("\n mathFunctionRoots : n1=%d  n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
1050       for(Standard_Integer x1=1; x1<=n1;x1++) {
1051         Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
1052         printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
1053       }
1054       printf("\n");
1055       for(Standard_Integer x2=1; x2<=n2; x2++) {
1056         Standard_Real v; F.Value(Sol(x2),v); v-=K;
1057         printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
1058       }
1059       printf("\n");
1060     }
1061     Standard_Integer n=n1;
1062     if(n1>n2) n=n2;
1063     for(Standard_Integer i=1;i<=n;i++) {
1064       Standard_Real t = Sol(i)-StaticSol(i);
1065       if(Abs(t)>NEpsX) {
1066         printf("\n mathFunctionRoots : i:%d/%d  delta: %g",i,n,t);
1067       }
1068     }
1069   }
1070 #endif
1071 #endif
1072 }
1073 #endif
1074 }
1075
1076
1077 void math_FunctionRoots::Dump(Standard_OStream& o) const
1078 {
1079
1080   o << "math_FunctionRoots ";
1081   if(Done) {
1082     o << " Status = Done \n";
1083     o << " Number of solutions = " << Sol.Length() << std::endl;
1084     for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
1085       o << " Solution Number " << i << "= " << Sol.Value(i) << std::endl;
1086     }
1087   }
1088   else {
1089     o<< " Status = not Done \n";
1090   }
1091 }