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