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
21 #include <StdFail_NotDone.hxx>
22 #include <Standard_RangeError.hxx>
23 #include <math_DirectPolynomialRoots.hxx>
24 #include <math_FunctionRoots.ixx>
25 #include <math_FunctionWithDerivative.hxx>
26 #include <TColStd_Array1OfReal.hxx>
37 static Standard_Boolean myDebug = 0;
38 static Standard_Integer nbsolve = 0;
41 static void AppendRoot(TColStd_SequenceOfReal& Sol,
42 TColStd_SequenceOfInteger& NbStateSol,
43 const Standard_Real X,
44 math_FunctionWithDerivative& F,
45 // const Standard_Real K,
47 const Standard_Real dX) {
49 Standard_Integer n=Sol.Length();
53 cout << " Ajout de la solution numero : " << n+1 << endl;
54 cout << " Valeur de la racine :" << X << endl;
61 NbStateSol.Append(F.GetStateNumber());
65 Standard_Integer pl=n+1;
81 NbStateSol.Append(F.GetStateNumber());
84 Sol.InsertBefore(pl,X);
86 NbStateSol.InsertBefore(pl,F.GetStateNumber());
91 static void Solve(math_FunctionWithDerivative& F,
92 const Standard_Real K,
93 const Standard_Real x1,
94 const Standard_Real y1,
95 const Standard_Real x2,
96 const Standard_Real y2,
97 const Standard_Real tol,
98 const Standard_Real dX,
99 TColStd_SequenceOfReal& Sol,
100 TColStd_SequenceOfInteger& NbStateSol) {
103 cout <<"--> Resolution :" << ++nbsolve << endl;
104 cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
105 cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
109 Standard_Integer iter=0;
110 Standard_Real tols2 = 0.5*tol;
111 Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
112 a=x1;b=c=x2;fa=y1; fb=fc=y2;
113 for(iter=1;iter<=ITMAX;iter++) {
114 if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
117 if(Abs(fc)<Abs(fb)) {
118 a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
120 tol1 = EPSEPS*Abs(b)+tols2;
122 if(Abs(xm)<tol1 || fb==0) {
123 //-- On tente une iteration de newton
124 Standard_Real Xp,Yp,Dp;
125 Standard_Integer itern=5;
129 Ok = F.Values(Xp,Yp,Dp);
132 if(Dp>1e-10 || Dp<-1e-10) {
135 if(Xp<=x2 && Xp>=x1) {
136 F.Value(Xp,Yp); Yp-=K;
137 if(Abs(Yp)<Abs(fb)) {
145 while(Ok && --itern>=0);
147 AppendRoot(Sol,NbStateSol,b,F,K,dX);
150 if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
159 p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
160 q=(q-1.0)*(r-1.0)*(s-1.0);
166 min1=3.0*xm*q-Abs(tol1*q);
168 if((p+p)<( (min1<min2) ? min1 : min2)) {
187 if(xm>=0) b+=Abs(tol1);
194 cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
200 #define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement
201 //#define MATH_FUNCTIONROOTS_OLDCODE // Ancien
202 //#define MATH_FUNCTIONROOTS_CHECK // Check
204 math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
205 const Standard_Real A,
206 const Standard_Real B,
207 const Standard_Integer NbSample,
208 const Standard_Real _EpsX,
209 const Standard_Real EpsF,
210 const Standard_Real EpsNull,
211 const Standard_Real K )
215 cout << "---- Debut de math_FunctionRoots ----" << endl;
221 TColStd_SequenceOfReal StaticSol;
225 #ifdef MATH_FUNCTIONROOTS_NEWCODE
227 Done = Standard_True;
230 Standard_Integer N=NbSample;
231 //-- ------------------------------------------------------------
232 //-- Verifications de bas niveau
241 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
242 Standard_Real EpsX = _EpsX;
243 Standard_Real DeltaU = Abs(X0)+Abs(XN);
244 Standard_Real NEpsX = 0.0000000001 * DeltaU;
249 //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
250 //-- A .............................................................. B
251 //-- X0 X1 X2 ........................................ Xn-1 Xn
255 double dx = (XN-X0)/N;
256 TColStd_Array1OfReal ptrval(0, N);
257 Standard_Integer Nvalid = -1;
258 Standard_Real aux = 0;
259 for(i=0; i<=N ; i++,X+=dx) {
262 if(Ok) ptrval(++Nvalid) = aux - K;
265 //-- Toute la fonction est nulle ?
268 Done = Standard_False;
272 AllNull=Standard_True;
273 // for(i=0;AllNull && i<=N;i++) {
274 for(i=0;AllNull && i<=N;i++) {
275 if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
276 AllNull=Standard_False;
280 //-- tous les points echantillons sont dans la tolerance
284 //-- Il y a des points hors tolerance
285 //-- on detecte les changements de signes STRICTS
286 Standard_Integer ip1;
287 // Standard_Boolean chgtsign=Standard_False;
288 Standard_Real tol = EpsX;
290 for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
294 if(ptrval(ip1)>0.0) {
295 //-- --------------------------------------------------
296 //-- changement de signe dans Xi Xi+1
297 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
301 if(ptrval(ip1)<0.0) {
302 //-- --------------------------------------------------
303 //-- changement de signe dans Xi Xi+1
304 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
308 //-- On detecte les cas ou la fct s annule sur des Xi et est
309 //-- non nulle au voisinage de Xi
311 //-- On prend 2 points u0,u1 au voisinage de Xi
312 //-- Si (F(u0)-K)*(F(u1)-K) <0 on lance une recherche
313 //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
314 for(i=0; i<=N; i++) {
316 // Standard_Real Val,Deriv;
320 u0=dx*0.5; u1=X+u0; u0+=X;
327 F.Value(u0,y0); y0-=K;
328 F.Value(u1,y1); y1-=K;
330 Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
333 if(y0!=0.0 || y1!=0.0) {
334 AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
339 //-- --------------------------------------------------------------------------------
340 //-- Il faut traiter differement le cas des points en bout :
341 if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
342 AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
344 if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
345 AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
348 //-- --------------------------------------------------------------------------------
349 //-- --------------------------------------------------------------------------------
350 //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
351 //-- un maximum avec f(x)<0
352 //-- On reprend une discretisation plus fine au voisinage de ces extremums
354 //-- Recherche d un minima positif
355 Standard_Real xm,ym,dym,xm1,xp1;
356 Standard_Real majdx = 5.0*dx;
357 Standard_Boolean Rediscr;
358 // Standard_Real ptrvalbis[MAXBIS];
359 Standard_Integer im1=0;
361 for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
362 Rediscr = Standard_False;
365 if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
366 //-- Peut on traverser l axe Ox
367 //-- -------------- Estimation a partir de Xim1
370 F.Values(xm1,ym,dym); ym-=K;
371 if(dym<-1e-10 || dym>1e-10) { // normalement dym < 0
372 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
373 if(t<majdx && t > -majdx) {
374 Rediscr = Standard_True;
377 //-- -------------- Estimation a partir de Xip1
378 if(Rediscr==Standard_False) {
380 if(xp1 > XN ) xp1=XN;
381 F.Values(xp1,ym,dym); ym-=K;
382 if(dym<-1e-10 || dym>1e-10) { // normalement dym > 0
383 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
384 if(t<majdx && t > -majdx) {
385 Rediscr = Standard_True;
391 else if(ptrval(i)<0.0) {
392 if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
393 //-- Peut on traverser l axe Ox
394 //-- -------------- Estimation a partir de Xim1
397 F.Values(xm1,ym,dym); ym-=K;
398 if(dym>1e-10 || dym<-1e-10) { // normalement dym > 0
399 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
400 if(t<majdx && t > -majdx) {
401 Rediscr = Standard_True;
404 //-- -------------- Estimation a partir de Xim1
405 if(Rediscr==Standard_False) {
408 F.Values(xm1,ym,dym); ym-=K;
409 if(dym>1e-10 || dym<-1e-10) { // normalement dym < 0
410 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
411 if(t<majdx && t > -majdx) {
412 Rediscr = Standard_True;
419 //-- --------------------------------------------------
420 //-- On recherche un extrema entre x0 et x3
421 //-- x1 et x2 sont tels que x0<x1<x2<x3
422 //-- et |f(x0)| > |f(x1)| et |f(x3)| > |f(x2)|
424 //-- En entree : a=xm-dx b=xm c=xm+dx
425 Standard_Real x0,x1,x2,x3,f0,f3;
426 Standard_Real R=0.61803399;
427 Standard_Real C=1.0-R;
428 Standard_Real tolCR=NEpsX*10.0;
435 Standard_Boolean recherche_minimum = (f0>0.0);
437 if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); }
438 else { x2=xm; x1=xm-C*(xm-x0); }
440 F.Value(x1,f1); f1-=K;
441 F.Value(x2,f2); f2-=K;
442 //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
443 Standard_Real tolX = 0.001 * NEpsX;
444 while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > tolX)) {
445 //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
446 //-- x0,f0,x1,f1,x2,f2,x3,f3);
447 if(recherche_minimum) {
449 x0=x1; x1=x2; x2=R*x1+C*x3;
450 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
453 x3=x2; x2=x1; x1=R*x2+C*x0;
454 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
459 x0=x1; x1=x2; x2=R*x1+C*x3;
460 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
463 x3=x2; x2=x1; x1=R*x2+C*x0;
464 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
467 //-- On ne fait pas que chercher des extremas. Il faut verifier
468 //-- si on ne tombe pas sur une racine
470 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
471 Solve(F,K,x0,f0,x1,f1,tol,NEpsX,Sol,NbStateSol);
474 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
475 Solve(F,K,x2,f2,x3,f3,tol,NEpsX,Sol,NbStateSol);
479 //-- x1,f(x1) minimum
481 AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
485 //-- x2.f(x2) minimum
487 AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
490 } //-- Recherche d un extrema
495 #ifdef MATH_FUNCTIONROOTS_CHECK
498 Standard_Integer n=Sol.Length();
499 for(Standard_Integer ii=1;ii<=n;ii++) {
500 StaticSol.Append(Sol.Value(ii));
509 #ifdef MATH_FUNCTIONROOTS_OLDCODE
511 //-- ********************************************************************************
512 //-- ANCIEN TRAITEMENT
513 //-- ********************************************************************************
516 // calculate all the real roots of a function within the range
517 // A..B. whitout condition on A and B
518 // a solution X is found when
519 // abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
520 // The function is considered as null between A and B if
521 // abs(F-K) <= EpsNull within this range.
522 Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
523 //-- Il ne faut pas EpsX = 0.000...001 car dans ce cas
524 //-- U + Nn*EpsX == U
525 Standard_Real Lowr,Upp;
526 Standard_Real Increment;
528 Standard_Real FLowr,FUpp,DFLowr,DFUpp;
530 Standard_Real Fxu,DFxu,FFxu,DFFxu;
531 Standard_Real Fyu,DFyu,FFyu,DFFyu;
532 Standard_Boolean Finish;
534 Standard_Integer Nbiter = 30;
535 Standard_Integer Iter;
536 Standard_Real Ambda,T;
537 Standard_Real AA,BB,CC;
539 Standard_Real Alfa1=0,Alfa2;
540 Standard_Real OldDF = RealLast();
541 Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
544 Done = Standard_False;
546 StdFail_NotDone_Raise_if(NbSample <= 0, " ");
559 Increment = (Upp-Lowr)/NbSample;
560 StdFail_NotDone_Raise_if(Increment < EpsX, " ");
561 Done = Standard_True;
562 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
563 Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
564 Standard_Real NEpsX = 0.0000000001 * DeltaU;
567 //-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
570 Null2 = EpsNull*EpsNull;
572 Ok = F.Values(Lowr,FLowr,DFLowr);
575 Done = Standard_False;
581 Ok = F.Values(Upp,FUpp,DFUpp);
584 Done = Standard_False;
593 Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
596 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
597 AllNull = ( FFyu <= Null2 );
609 Fyu = FLowr + (U-Lowr)*DFLowr;
613 Fyu = FUpp + (U-Upp)*DFUpp;
617 Ok = F.Values(U,Fyu,DFyu);
620 Done = Standard_False;
627 DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
629 if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
631 if (AllNull) { //rechercher les vraix zeros depuis le debut
633 AllNull = Standard_False;
635 Fxu = FLowr-EpsX*DFLowr;
638 DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
640 Ok = F.Values(U,Fyu,DFyu);
643 Done = Standard_False;
649 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
651 Standard_Real FxuFyu=Fxu*Fyu;
653 if ( (DFFyu > 0. && DFFxu <= 0.)
654 || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
655 || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
656 || (FxuFyu <= 0.) ) {
657 // recherche d 1 minimun possible
658 Finish = Standard_False;
665 // chercher si f peut s annuler pour eviter
666 // des iterations inutiles
667 if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
668 Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
670 Finish = Standard_True;
671 FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
673 else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
674 (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
676 Finish = Standard_True;
678 FFi = FFyu; // pour recalculer yu
680 else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
681 (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
683 Finish = Standard_True;
685 FFi = FFxu; // pour recalculer U
688 else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
690 Finish = Standard_True;
694 else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
696 Finish = Standard_True;
702 // calcul des 2 solutions annulant la derivee de l interpolation cubique
703 // Ambda*t=(U-Xu) F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
704 // df=aa*t*t+2*bb*t+cc
706 AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
707 BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
710 if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
715 math_DirectPolynomialRoots Solv(AA,BB,CC);
716 if ( !Solv.InfiniteRoots() ) {
717 Nn=Solv.NbSolutions();
719 if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
720 else Finish = Standard_True;
723 Alfa1 = Solv.Value(1);
725 Alfa2 = Solv.Value(2);
731 if (Alfa1 > 1. || Alfa2 < 0.){
732 // resolution par dichotomie
733 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
734 else Finish = Standard_True;
736 else if ( Alfa1 < 0. ||
737 ( DFFxu > 0. && DFFyu >= 0.) ) {
739 //(cas changement de signe de la distance signee sans
740 // changement de signe de la derivee:
741 //cas de 'presque'tangence avec 2
742 // solutions proches) ,on prend la plus grane racine
744 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
745 else Finish = Standard_True;
752 else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
753 //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
754 else Finish = Standard_True;
757 // petits tests pour diminuer le nombre d iterations
758 if (Alfa1 <= EpsX ) {
761 else if (Alfa1 >= (1.-EpsX) ) {
762 Alfa1 = Alfa1+Alfa1-1.;
764 Alfa1 = Ambda*(1. - Alfa1);
767 AA = FLowr + (U - Lowr)*DFLowr;
770 else if ( U >= Upp ) {
771 AA = FUpp + (U - Upp)*DFUpp;
775 Ok = F.Values(U,AA,BB);
778 Done = Standard_False;
786 if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
793 if (Alfa1 > Ambda*0.5) {
795 // determination d 1 autre borne pour diviser
796 //le nouvel intervalle par 2 au -
799 AA = FLowr + (Xu - Lowr)*DFLowr;
802 else if ( Xu >= Upp ) {
803 AA = FUpp + (Xu - Upp)*DFUpp;
807 Ok = F.Values(Xu,AA,BB);
810 Done = Standard_False;
818 if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
827 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
828 // changement de signe sur l intervalle u,U
833 FFi = Min(FFxu,FFyu);
838 else { Ambda = Alfa1;}
840 else { Ambda = Alfa1;}
847 if ( (Ambda-Alfa1) > Ambda*0.5 ) {
849 Xu = U - (Ambda-Alfa1)*0.5;
851 AA = FLowr + (Xu - Lowr)*DFLowr;
854 else if ( Xu >= Upp ) {
855 AA = FUpp + (Xu - Upp)*DFUpp;
859 Ok = F.Values(Xu,AA,BB);
862 Done = Standard_False;
870 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
875 Ambda = ( Ambda - Alfa1 )*0.5;
878 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
883 Ambda = ( Ambda - Alfa1 )*0.5;
885 FFi = Min(FFxu , FFyu);
890 Ambda = Ambda - Alfa1;
895 Ambda = Ambda - Alfa1;
899 if (Abs(FFxu) <= Standard_Underflow ||
900 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
902 Finish = Standard_True;
903 if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
906 else if (Abs(FFyu) <= Standard_Underflow ||
907 (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
909 Finish = Standard_True;
910 if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
915 Finish = Iter >= Nbiter || (Ambda <= EpsX &&
916 ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
919 } // fin interpolation cubique
921 // restitution du meilleur resultat
923 if ( FFxu < FFi ) { U = U + T -Ambda;}
924 else if ( FFyu < FFi ) { U = U + T;}
926 if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
927 U = Max( Lowr , Min( U , Upp));
930 if ( Abs(FFi) < EpsF ) {
932 if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
933 else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
934 else if (Fxu*Fyu > 0.) { AA = 0.;}
935 else { AA = Fyu-Fxu;}
936 if (!Sol.IsEmpty()) {
937 if (Abs(Sol.Last() - U) > 5.*EpsX
940 NbStateSol.Append(F.GetStateNumber());
945 NbStateSol.Append(F.GetStateNumber());
953 while ( Nn < 1000000 && DFFyu <= 0.) {
954 // on repart d 1 pouyem plus loin
957 AA = FLowr + (U - Lowr)*DFLowr;
960 else if ( U >= Upp ) {
961 AA = FUpp + (U - Upp)*DFUpp;
965 Ok = F.Values(U,AA,BB);
984 #ifdef MATH_FUNCTIONROOTS_CHECK
986 Standard_Integer n1=StaticSol.Length();
987 Standard_Integer n2=Sol.Length();
989 printf("\n mathFunctionRoots : n1=%d n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
990 for(Standard_Integer x1=1; x1<=n1;x1++) {
991 Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
992 printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
995 for(Standard_Integer x2=1; x2<=n2; x2++) {
996 Standard_Real v; F.Value(Sol(x2),v); v-=K;
997 printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
1001 Standard_Integer n=n1;
1003 for(Standard_Integer i=1;i<=n;i++) {
1004 Standard_Real t = Sol(i)-StaticSol(i);
1006 printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
1017 void math_FunctionRoots::Dump(Standard_OStream& o) const
1020 o << "math_FunctionRoots ";
1022 o << " Status = Done \n";
1023 o << " Number of solutions = " << Sol.Length() << endl;
1024 for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
1025 o << " Solution Number " << i << "= " << Sol.Value(i) << endl;
1029 o<< " Status = not Done \n";