1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
22 #include <math_DirectPolynomialRoots.hxx>
23 #include <math_FunctionRoots.hxx>
24 #include <math_FunctionWithDerivative.hxx>
25 #include <Standard_RangeError.hxx>
26 #include <StdFail_NotDone.hxx>
27 #include <TColStd_Array1OfReal.hxx>
36 static Standard_Boolean myDebug = 0;
37 static Standard_Integer nbsolve = 0;
40 static void AppendRoot(TColStd_SequenceOfReal& Sol,
41 TColStd_SequenceOfInteger& NbStateSol,
42 const Standard_Real X,
43 math_FunctionWithDerivative& F,
44 // const Standard_Real K,
46 const Standard_Real dX) {
48 Standard_Integer n=Sol.Length();
52 cout << " Ajout de la solution numero : " << n+1 << endl;
53 cout << " Valeur de la racine :" << X << endl;
60 NbStateSol.Append(F.GetStateNumber());
64 Standard_Integer pl=n+1;
80 NbStateSol.Append(F.GetStateNumber());
83 Sol.InsertBefore(pl,X);
85 NbStateSol.InsertBefore(pl,F.GetStateNumber());
90 static void Solve(math_FunctionWithDerivative& F,
91 const Standard_Real K,
92 const Standard_Real x1,
93 const Standard_Real y1,
94 const Standard_Real x2,
95 const Standard_Real y2,
96 const Standard_Real tol,
97 const Standard_Real dX,
98 TColStd_SequenceOfReal& Sol,
99 TColStd_SequenceOfInteger& NbStateSol) {
102 cout <<"--> Resolution :" << ++nbsolve << endl;
103 cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
104 cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
108 Standard_Integer iter=0;
109 Standard_Real tols2 = 0.5*tol;
110 Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
111 a=x1;b=c=x2;fa=y1; fb=fc=y2;
112 for(iter=1;iter<=ITMAX;iter++) {
113 if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
116 if(Abs(fc)<Abs(fb)) {
117 a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
119 tol1 = EPSEPS*Abs(b)+tols2;
121 if(Abs(xm)<tol1 || fb==0) {
122 //-- On tente une iteration de newton
123 Standard_Real Xp,Yp,Dp;
124 Standard_Integer itern=5;
128 Ok = F.Values(Xp,Yp,Dp);
131 if(Dp>1e-10 || Dp<-1e-10) {
134 if(Xp<=x2 && Xp>=x1) {
135 F.Value(Xp,Yp); Yp-=K;
136 if(Abs(Yp)<Abs(fb)) {
144 while(Ok && --itern>=0);
146 AppendRoot(Sol,NbStateSol,b,F,K,dX);
149 if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
158 p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
159 q=(q-1.0)*(r-1.0)*(s-1.0);
165 min1=3.0*xm*q-Abs(tol1*q);
167 if((p+p)<( (min1<min2) ? min1 : min2)) {
186 if(xm>=0) b+=Abs(tol1);
193 cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
199 #define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement
200 //#define MATH_FUNCTIONROOTS_OLDCODE // Ancien
201 //#define MATH_FUNCTIONROOTS_CHECK // Check
203 math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
204 const Standard_Real A,
205 const Standard_Real B,
206 const Standard_Integer NbSample,
207 const Standard_Real _EpsX,
208 const Standard_Real EpsF,
209 const Standard_Real EpsNull,
210 const Standard_Real K )
214 cout << "---- Debut de math_FunctionRoots ----" << endl;
220 TColStd_SequenceOfReal StaticSol;
224 #ifdef MATH_FUNCTIONROOTS_NEWCODE
226 Done = Standard_True;
229 Standard_Integer N=NbSample;
230 //-- ------------------------------------------------------------
231 //-- Verifications de bas niveau
240 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
241 Standard_Real EpsX = _EpsX;
242 Standard_Real DeltaU = Abs(X0)+Abs(XN);
243 Standard_Real NEpsX = 0.0000000001 * DeltaU;
248 //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
249 //-- A .............................................................. B
250 //-- X0 X1 X2 ........................................ Xn-1 Xn
254 double dx = (XN-X0)/N;
255 TColStd_Array1OfReal ptrval(0, N);
256 Standard_Integer Nvalid = -1;
257 Standard_Real aux = 0;
258 for(i=0; i<=N ; i++,X+=dx) {
261 if(Ok) ptrval(++Nvalid) = aux - K;
264 //-- Toute la fonction est nulle ?
267 Done = Standard_False;
271 AllNull=Standard_True;
272 // for(i=0;AllNull && i<=N;i++) {
273 for(i=0;AllNull && i<=N;i++) {
274 if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
275 AllNull=Standard_False;
279 //-- tous les points echantillons sont dans la tolerance
283 //-- Il y a des points hors tolerance
284 //-- on detecte les changements de signes STRICTS
285 Standard_Integer ip1;
286 // Standard_Boolean chgtsign=Standard_False;
287 Standard_Real tol = EpsX;
289 for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
293 if(ptrval(ip1)>0.0) {
294 //-- --------------------------------------------------
295 //-- changement de signe dans Xi Xi+1
296 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
300 if(ptrval(ip1)<0.0) {
301 //-- --------------------------------------------------
302 //-- changement de signe dans Xi Xi+1
303 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
307 //-- On detecte les cas ou la fct s annule sur des Xi et est
308 //-- non nulle au voisinage de Xi
310 //-- On prend 2 points u0,u1 au voisinage de Xi
311 //-- Si (F(u0)-K)*(F(u1)-K) <0 on lance une recherche
312 //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
313 for(i=0; i<=N; i++) {
315 // Standard_Real Val,Deriv;
319 u0=dx*0.5; u1=X+u0; u0+=X;
326 F.Value(u0,y0); y0-=K;
327 F.Value(u1,y1); y1-=K;
329 Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
332 if(y0!=0.0 || y1!=0.0) {
333 AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
338 //-- --------------------------------------------------------------------------------
339 //-- Il faut traiter differement le cas des points en bout :
340 if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
341 AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
343 if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
344 AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
347 //-- --------------------------------------------------------------------------------
348 //-- --------------------------------------------------------------------------------
349 //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
350 //-- un maximum avec f(x)<0
351 //-- On reprend une discretisation plus fine au voisinage de ces extremums
353 //-- Recherche d un minima positif
354 Standard_Real xm,ym,dym,xm1,xp1;
355 Standard_Real majdx = 5.0*dx;
356 Standard_Boolean Rediscr;
357 // Standard_Real ptrvalbis[MAXBIS];
358 Standard_Integer im1=0;
360 for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
361 Rediscr = Standard_False;
364 if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
365 //-- Peut on traverser l axe Ox
366 //-- -------------- Estimation a partir de Xim1
369 F.Values(xm1,ym,dym); ym-=K;
370 if(dym<-1e-10 || dym>1e-10) { // normalement dym < 0
371 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
372 if(t<majdx && t > -majdx) {
373 Rediscr = Standard_True;
376 //-- -------------- Estimation a partir de Xip1
377 if(Rediscr==Standard_False) {
379 if(xp1 > XN ) xp1=XN;
380 F.Values(xp1,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;
390 else if(ptrval(i)<0.0) {
391 if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
392 //-- Peut on traverser l axe Ox
393 //-- -------------- Estimation a partir de Xim1
396 F.Values(xm1,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;
403 //-- -------------- Estimation a partir de Xim1
404 if(Rediscr==Standard_False) {
407 F.Values(xm1,ym,dym); ym-=K;
408 if(dym>1e-10 || dym<-1e-10) { // normalement dym < 0
409 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
410 if(t<majdx && t > -majdx) {
411 Rediscr = Standard_True;
418 //-- --------------------------------------------------
419 //-- On recherche un extrema entre x0 et x3
420 //-- x1 et x2 sont tels que x0<x1<x2<x3
421 //-- et |f(x0)| > |f(x1)| et |f(x3)| > |f(x2)|
423 //-- En entree : a=xm-dx b=xm c=xm+dx
424 Standard_Real x0,x1,x2,x3,f0,f3;
425 Standard_Real R=0.61803399;
426 Standard_Real C=1.0-R;
427 Standard_Real tolCR=NEpsX*10.0;
434 Standard_Boolean recherche_minimum = (f0>0.0);
436 if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); }
437 else { x2=xm; x1=xm-C*(xm-x0); }
439 F.Value(x1,f1); f1-=K;
440 F.Value(x2,f2); f2-=K;
441 //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
442 Standard_Real tolX = 0.001 * NEpsX;
443 while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > tolX)) {
444 //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
445 //-- x0,f0,x1,f1,x2,f2,x3,f3);
446 if(recherche_minimum) {
448 x0=x1; x1=x2; x2=R*x1+C*x3;
449 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
452 x3=x2; x2=x1; x1=R*x2+C*x0;
453 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
458 x0=x1; x1=x2; x2=R*x1+C*x3;
459 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
462 x3=x2; x2=x1; x1=R*x2+C*x0;
463 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
466 //-- On ne fait pas que chercher des extremas. Il faut verifier
467 //-- si on ne tombe pas sur une racine
469 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
470 Solve(F,K,x0,f0,x1,f1,tol,NEpsX,Sol,NbStateSol);
473 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
474 Solve(F,K,x2,f2,x3,f3,tol,NEpsX,Sol,NbStateSol);
478 //-- x1,f(x1) minimum
480 AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
484 //-- x2.f(x2) minimum
486 AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
489 } //-- Recherche d un extrema
494 #ifdef MATH_FUNCTIONROOTS_CHECK
497 Standard_Integer n=Sol.Length();
498 for(Standard_Integer ii=1;ii<=n;ii++) {
499 StaticSol.Append(Sol.Value(ii));
508 #ifdef MATH_FUNCTIONROOTS_OLDCODE
510 //-- ********************************************************************************
511 //-- ANCIEN TRAITEMENT
512 //-- ********************************************************************************
515 // calculate all the real roots of a function within the range
516 // A..B. whitout condition on A and B
517 // a solution X is found when
518 // abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
519 // The function is considered as null between A and B if
520 // abs(F-K) <= EpsNull within this range.
521 Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
522 //-- Il ne faut pas EpsX = 0.000...001 car dans ce cas
523 //-- U + Nn*EpsX == U
524 Standard_Real Lowr,Upp;
525 Standard_Real Increment;
527 Standard_Real FLowr,FUpp,DFLowr,DFUpp;
529 Standard_Real Fxu,DFxu,FFxu,DFFxu;
530 Standard_Real Fyu,DFyu,FFyu,DFFyu;
531 Standard_Boolean Finish;
533 Standard_Integer Nbiter = 30;
534 Standard_Integer Iter;
535 Standard_Real Ambda,T;
536 Standard_Real AA,BB,CC;
538 Standard_Real Alfa1=0,Alfa2;
539 Standard_Real OldDF = RealLast();
540 Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
543 Done = Standard_False;
545 StdFail_NotDone_Raise_if(NbSample <= 0, " ");
558 Increment = (Upp-Lowr)/NbSample;
559 StdFail_NotDone_Raise_if(Increment < EpsX, " ");
560 Done = Standard_True;
561 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
562 Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
563 Standard_Real NEpsX = 0.0000000001 * DeltaU;
566 //-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
569 Null2 = EpsNull*EpsNull;
571 Ok = F.Values(Lowr,FLowr,DFLowr);
574 Done = Standard_False;
580 Ok = F.Values(Upp,FUpp,DFUpp);
583 Done = Standard_False;
592 Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
595 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
596 AllNull = ( FFyu <= Null2 );
608 Fyu = FLowr + (U-Lowr)*DFLowr;
612 Fyu = FUpp + (U-Upp)*DFUpp;
616 Ok = F.Values(U,Fyu,DFyu);
619 Done = Standard_False;
626 DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
628 if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
630 if (AllNull) { //rechercher les vraix zeros depuis le debut
632 AllNull = Standard_False;
634 Fxu = FLowr-EpsX*DFLowr;
637 DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
639 Ok = F.Values(U,Fyu,DFyu);
642 Done = Standard_False;
648 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
650 Standard_Real FxuFyu=Fxu*Fyu;
652 if ( (DFFyu > 0. && DFFxu <= 0.)
653 || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
654 || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
655 || (FxuFyu <= 0.) ) {
656 // recherche d 1 minimun possible
657 Finish = Standard_False;
664 // chercher si f peut s annuler pour eviter
665 // des iterations inutiles
666 if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
667 Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
669 Finish = Standard_True;
670 FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
672 else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
673 (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
675 Finish = Standard_True;
677 FFi = FFyu; // pour recalculer yu
679 else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
680 (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
682 Finish = Standard_True;
684 FFi = FFxu; // pour recalculer U
687 else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
689 Finish = Standard_True;
693 else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
695 Finish = Standard_True;
701 // calcul des 2 solutions annulant la derivee de l interpolation cubique
702 // Ambda*t=(U-Xu) F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
703 // df=aa*t*t+2*bb*t+cc
705 AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
706 BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
709 if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
714 math_DirectPolynomialRoots Solv(AA,BB,CC);
715 if ( !Solv.InfiniteRoots() ) {
716 Nn=Solv.NbSolutions();
718 if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
719 else Finish = Standard_True;
722 Alfa1 = Solv.Value(1);
724 Alfa2 = Solv.Value(2);
730 if (Alfa1 > 1. || Alfa2 < 0.){
731 // resolution par dichotomie
732 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
733 else Finish = Standard_True;
735 else if ( Alfa1 < 0. ||
736 ( DFFxu > 0. && DFFyu >= 0.) ) {
738 //(cas changement de signe de la distance signee sans
739 // changement de signe de la derivee:
740 //cas de 'presque'tangence avec 2
741 // solutions proches) ,on prend la plus grane racine
743 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
744 else Finish = Standard_True;
751 else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
752 //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
753 else Finish = Standard_True;
756 // petits tests pour diminuer le nombre d iterations
757 if (Alfa1 <= EpsX ) {
760 else if (Alfa1 >= (1.-EpsX) ) {
761 Alfa1 = Alfa1+Alfa1-1.;
763 Alfa1 = Ambda*(1. - Alfa1);
766 AA = FLowr + (U - Lowr)*DFLowr;
769 else if ( U >= Upp ) {
770 AA = FUpp + (U - Upp)*DFUpp;
774 Ok = F.Values(U,AA,BB);
777 Done = Standard_False;
785 if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
792 if (Alfa1 > Ambda*0.5) {
794 // determination d 1 autre borne pour diviser
795 //le nouvel intervalle par 2 au -
798 AA = FLowr + (Xu - Lowr)*DFLowr;
801 else if ( Xu >= Upp ) {
802 AA = FUpp + (Xu - Upp)*DFUpp;
806 Ok = F.Values(Xu,AA,BB);
809 Done = Standard_False;
817 if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
826 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
827 // changement de signe sur l intervalle u,U
832 FFi = Min(FFxu,FFyu);
837 else { Ambda = Alfa1;}
839 else { Ambda = Alfa1;}
846 if ( (Ambda-Alfa1) > Ambda*0.5 ) {
848 Xu = U - (Ambda-Alfa1)*0.5;
850 AA = FLowr + (Xu - Lowr)*DFLowr;
853 else if ( Xu >= Upp ) {
854 AA = FUpp + (Xu - Upp)*DFUpp;
858 Ok = F.Values(Xu,AA,BB);
861 Done = Standard_False;
869 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
874 Ambda = ( Ambda - Alfa1 )*0.5;
877 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
882 Ambda = ( Ambda - Alfa1 )*0.5;
884 FFi = Min(FFxu , FFyu);
889 Ambda = Ambda - Alfa1;
894 Ambda = Ambda - Alfa1;
898 if (Abs(FFxu) <= Standard_Underflow ||
899 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
901 Finish = Standard_True;
902 if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
905 else if (Abs(FFyu) <= Standard_Underflow ||
906 (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
908 Finish = Standard_True;
909 if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
914 Finish = Iter >= Nbiter || (Ambda <= EpsX &&
915 ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
918 } // fin interpolation cubique
920 // restitution du meilleur resultat
922 if ( FFxu < FFi ) { U = U + T -Ambda;}
923 else if ( FFyu < FFi ) { U = U + T;}
925 if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
926 U = Max( Lowr , Min( U , Upp));
929 if ( Abs(FFi) < EpsF ) {
931 if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
932 else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
933 else if (Fxu*Fyu > 0.) { AA = 0.;}
934 else { AA = Fyu-Fxu;}
935 if (!Sol.IsEmpty()) {
936 if (Abs(Sol.Last() - U) > 5.*EpsX
937 || (OldDF != RealLast() && OldDF*AA < 0.) ) {
939 NbStateSol.Append(F.GetStateNumber());
944 NbStateSol.Append(F.GetStateNumber());
952 while ( Nn < 1000000 && DFFyu <= 0.) {
953 // on repart d 1 pouyem plus loin
956 AA = FLowr + (U - Lowr)*DFLowr;
959 else if ( U >= Upp ) {
960 AA = FUpp + (U - Upp)*DFUpp;
964 Ok = F.Values(U,AA,BB);
983 #ifdef MATH_FUNCTIONROOTS_CHECK
985 Standard_Integer n1=StaticSol.Length();
986 Standard_Integer n2=Sol.Length();
988 printf("\n mathFunctionRoots : n1=%d n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
989 for(Standard_Integer x1=1; x1<=n1;x1++) {
990 Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
991 printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
994 for(Standard_Integer x2=1; x2<=n2; x2++) {
995 Standard_Real v; F.Value(Sol(x2),v); v-=K;
996 printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
1000 Standard_Integer n=n1;
1002 for(Standard_Integer i=1;i<=n;i++) {
1003 Standard_Real t = Sol(i)-StaticSol(i);
1005 printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
1016 void math_FunctionRoots::Dump(Standard_OStream& o) const
1019 o << "math_FunctionRoots ";
1021 o << " Status = Done \n";
1022 o << " Number of solutions = " << Sol.Length() << endl;
1023 for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
1024 o << " Solution Number " << i << "= " << Sol.Value(i) << endl;
1028 o<< " Status = not Done \n";