src/math/math_FunctionRoots.cxx

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