1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
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.
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.
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.
20 #define No_Standard_RangeError
21 #define No_Standard_OutOfRange
22 #define No_Standard_DimensionError
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>
41 static Standard_Boolean myDebug = 0;
42 static Standard_Integer nbsolve = 0;
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,
51 const Standard_Real dX) {
53 Standard_Integer n=Sol.Length();
57 cout << " Ajout de la solution numero : " << n+1 << endl;
58 cout << " Valeur de la racine :" << X << endl;
65 NbStateSol.Append(F.GetStateNumber());
69 Standard_Integer pl=n+1;
85 NbStateSol.Append(F.GetStateNumber());
88 Sol.InsertBefore(pl,X);
90 NbStateSol.InsertBefore(pl,F.GetStateNumber());
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) {
107 cout <<"--> Resolution :" << ++nbsolve << endl;
108 cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
109 cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
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)) {
121 if(Abs(fc)<Abs(fb)) {
122 a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
124 tol1 = EPSEPS*Abs(b)+tols2;
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;
133 Ok = F.Values(Xp,Yp,Dp);
136 if(Dp>1e-10 || Dp<-1e-10) {
139 if(Xp<=x2 && Xp>=x1) {
140 F.Value(Xp,Yp); Yp-=K;
141 if(Abs(Yp)<Abs(fb)) {
149 while(Ok && --itern>=0);
151 AppendRoot(Sol,NbStateSol,b,F,K,dX);
154 if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
163 p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
164 q=(q-1.0)*(r-1.0)*(s-1.0);
170 min1=3.0*xm*q-Abs(tol1*q);
172 if((p+p)<( (min1<min2) ? min1 : min2)) {
191 if(xm>=0) b+=Abs(tol1);
198 cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
205 TColStd_SequenceOfReal StaticSol;
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 )
220 cout << "---- Debut de math_FunctionRoots ----" << endl;
225 static Standard_Integer methode = 1; //-- 1:(Nv Traitement) 3:(Nv + Ancien +check) 2:(Ancien)
229 Done = Standard_True;
232 Standard_Integer N=NbSample;
233 //-- ------------------------------------------------------------
234 //-- Verifications de bas niveau
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;
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
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) {
264 if(Ok) ptrval(++Nvalid) = aux - K;
267 //-- Toute la fonction est nulle ?
270 Done = Standard_False;
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;
282 //-- tous les points echantillons sont dans la tolerance
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;
292 for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
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);
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);
310 //-- On detecte les cas ou la fct s annule sur des Xi et est
311 //-- non nulle au voisinage de Xi
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++) {
318 // Standard_Real Val,Deriv;
322 u0=dx*0.5; u1=X+u0; u0+=X;
329 F.Value(u0,y0); y0-=K;
330 F.Value(u1,y1); y1-=K;
332 Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
335 if(y0!=0.0 || y1!=0.0) {
336 AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
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);
346 if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
347 AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
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
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;
363 for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
364 Rediscr = Standard_False;
367 if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
368 //-- Peut on traverser l axe Ox
369 //-- -------------- Estimation a partir de Xim1
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;
379 //-- -------------- Estimation a partir de Xip1
380 if(Rediscr==Standard_False) {
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;
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
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;
406 //-- -------------- Estimation a partir de Xim1
407 if(Rediscr==Standard_False) {
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;
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)|
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;
437 Standard_Boolean recherche_minimum = (f0>0.0);
439 if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); }
440 else { x2=xm; x1=xm-C*(xm-x0); }
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) {
450 x0=x1; x1=x2; x2=R*x1+C*x3;
451 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
454 x3=x2; x2=x1; x1=R*x2+C*x0;
455 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
460 x0=x1; x1=x2; x2=R*x1+C*x3;
461 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
464 x3=x2; x2=x1; x1=R*x2+C*x0;
465 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
468 //-- On ne fait pas que chercher des extremas. Il faut verifier
469 //-- si on ne tombe pas sur une racine
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);
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);
480 //-- x1,f(x1) minimum
482 AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
486 //-- x2.f(x2) minimum
488 AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
491 } //-- Recherche d un extrema
498 Standard_Integer n=Sol.Length();
499 for(Standard_Integer ii=1;ii<=n;ii++) {
500 StaticSol.Append(Sol.Value(ii));
508 //-- ********************************************************************************
509 //-- ANCIEN TRAITEMENT
510 //-- ********************************************************************************
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;
525 Standard_Real FLowr,FUpp,DFLowr,DFUpp;
527 Standard_Real Fxu,DFxu,FFxu,DFFxu;
528 Standard_Real Fyu,DFyu,FFyu,DFFyu;
529 Standard_Boolean Finish;
531 Standard_Integer Nbiter = 30;
532 Standard_Integer Iter;
533 Standard_Real Ambda,T;
534 Standard_Real AA,BB,CC;
536 Standard_Real Alfa1=0,Alfa2;
537 Standard_Real OldDF = RealLast();
538 Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
541 Done = Standard_False;
543 StdFail_NotDone_Raise_if(NbSample <= 0, " ");
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;
564 //-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
567 Null2 = EpsNull*EpsNull;
569 Ok = F.Values(Lowr,FLowr,DFLowr);
572 Done = Standard_False;
578 Ok = F.Values(Upp,FUpp,DFUpp);
581 Done = Standard_False;
590 Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
593 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
594 AllNull = ( FFyu <= Null2 );
606 Fyu = FLowr + (U-Lowr)*DFLowr;
610 Fyu = FUpp + (U-Upp)*DFUpp;
614 Ok = F.Values(U,Fyu,DFyu);
617 Done = Standard_False;
624 DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
626 if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
628 if (AllNull) { //rechercher les vraix zeros depuis le debut
630 AllNull = Standard_False;
632 Fxu = FLowr-EpsX*DFLowr;
635 DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
637 Ok = F.Values(U,Fyu,DFyu);
640 Done = Standard_False;
646 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
648 Standard_Real FxuFyu=Fxu*Fyu;
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;
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. ) {
667 Finish = Standard_True;
668 FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
670 else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
671 (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
673 Finish = Standard_True;
675 FFi = FFyu; // pour recalculer yu
677 else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
678 (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
680 Finish = Standard_True;
682 FFi = FFxu; // pour recalculer U
685 else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
687 Finish = Standard_True;
691 else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
693 Finish = Standard_True;
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
703 AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
704 BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
707 if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
712 math_DirectPolynomialRoots Solv(AA,BB,CC);
713 if ( !Solv.InfiniteRoots() ) {
714 Nn=Solv.NbSolutions();
716 if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
717 else Finish = Standard_True;
720 Alfa1 = Solv.Value(1);
722 Alfa2 = Solv.Value(2);
728 if (Alfa1 > 1. || Alfa2 < 0.){
729 // resolution par dichotomie
730 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
731 else Finish = Standard_True;
733 else if ( Alfa1 < 0. ||
734 ( DFFxu > 0. && DFFyu >= 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
741 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
742 else Finish = Standard_True;
749 else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
750 //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
751 else Finish = Standard_True;
754 // petits tests pour diminuer le nombre d iterations
755 if (Alfa1 <= EpsX ) {
758 else if (Alfa1 >= (1.-EpsX) ) {
759 Alfa1 = Alfa1+Alfa1-1.;
761 Alfa1 = Ambda*(1. - Alfa1);
764 AA = FLowr + (U - Lowr)*DFLowr;
767 else if ( U >= Upp ) {
768 AA = FUpp + (U - Upp)*DFUpp;
772 Ok = F.Values(U,AA,BB);
775 Done = Standard_False;
783 if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
790 if (Alfa1 > Ambda*0.5) {
792 // determination d 1 autre borne pour diviser
793 //le nouvel intervalle par 2 au -
796 AA = FLowr + (Xu - Lowr)*DFLowr;
799 else if ( Xu >= Upp ) {
800 AA = FUpp + (Xu - Upp)*DFUpp;
804 Ok = F.Values(Xu,AA,BB);
807 Done = Standard_False;
815 if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
824 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
825 // changement de signe sur l intervalle u,U
830 FFi = Min(FFxu,FFyu);
835 else { Ambda = Alfa1;}
837 else { Ambda = Alfa1;}
844 if ( (Ambda-Alfa1) > Ambda*0.5 ) {
846 Xu = U - (Ambda-Alfa1)*0.5;
848 AA = FLowr + (Xu - Lowr)*DFLowr;
851 else if ( Xu >= Upp ) {
852 AA = FUpp + (Xu - Upp)*DFUpp;
856 Ok = F.Values(Xu,AA,BB);
859 Done = Standard_False;
867 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
872 Ambda = ( Ambda - Alfa1 )*0.5;
875 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
880 Ambda = ( Ambda - Alfa1 )*0.5;
882 FFi = Min(FFxu , FFyu);
887 Ambda = Ambda - Alfa1;
892 Ambda = Ambda - Alfa1;
896 if (Abs(FFxu) <= Standard_Underflow ||
897 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
899 Finish = Standard_True;
900 if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
903 else if (Abs(FFyu) <= Standard_Underflow ||
904 (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
906 Finish = Standard_True;
907 if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
912 Finish = Iter >= Nbiter || (Ambda <= EpsX &&
913 ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
916 } // fin interpolation cubique
918 // restitution du meilleur resultat
920 if ( FFxu < FFi ) { U = U + T -Ambda;}
921 else if ( FFyu < FFi ) { U = U + T;}
923 if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
924 U = Max( Lowr , Min( U , Upp));
927 if ( Abs(FFi) < EpsF ) {
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
937 NbStateSol.Append(F.GetStateNumber());
942 NbStateSol.Append(F.GetStateNumber());
950 while ( Nn < 1000000 && DFFyu <= 0.) {
951 // on repart d 1 pouyem plus loin
954 AA = FLowr + (U - Lowr)*DFLowr;
957 else if ( U >= Upp ) {
958 AA = FUpp + (U - Upp)*DFUpp;
962 Ok = F.Values(U,AA,BB);
982 Standard_Integer n1=StaticSol.Length();
983 Standard_Integer n2=Sol.Length();
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);
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);
997 Standard_Integer n=n1;
999 for(Standard_Integer i=1;i<=n;i++) {
1000 Standard_Real t = Sol(i)-StaticSol(i);
1002 printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
1011 void math_FunctionRoots::Dump(Standard_OStream& o) const
1014 o << "math_FunctionRoots ";
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;
1023 o<< " Status = not Done \n";