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