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 <math_BracketedRoot.hxx>
26 #include <Standard_RangeError.hxx>
27 #include <StdFail_NotDone.hxx>
28 #include <TColStd_Array1OfReal.hxx>
37 static Standard_Boolean myDebug = 0;
38 static Standard_Integer nbsolve = 0;
41 class DerivFunction: public math_Function
43 math_FunctionWithDerivative *myF;
46 DerivFunction(math_FunctionWithDerivative& theF)
50 virtual Standard_Boolean Value(const Standard_Real theX, Standard_Real& theFval)
52 return myF->Derivative(theX, theFval);
56 static void AppendRoot(TColStd_SequenceOfReal& Sol,
57 TColStd_SequenceOfInteger& NbStateSol,
58 const Standard_Real X,
59 math_FunctionWithDerivative& F,
60 // const Standard_Real K,
62 const Standard_Real dX) {
64 Standard_Integer n=Sol.Length();
68 cout << " Ajout de la solution numero : " << n+1 << endl;
69 cout << " Valeur de la racine :" << X << endl;
76 NbStateSol.Append(F.GetStateNumber());
80 Standard_Integer pl=n+1;
96 NbStateSol.Append(F.GetStateNumber());
99 Sol.InsertBefore(pl,X);
101 NbStateSol.InsertBefore(pl,F.GetStateNumber());
106 static void Solve(math_FunctionWithDerivative& F,
107 const Standard_Real K,
108 const Standard_Real x1,
109 const Standard_Real y1,
110 const Standard_Real x2,
111 const Standard_Real y2,
112 const Standard_Real tol,
113 const Standard_Real dX,
114 TColStd_SequenceOfReal& Sol,
115 TColStd_SequenceOfInteger& NbStateSol) {
118 cout <<"--> Resolution :" << ++nbsolve << endl;
119 cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
120 cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
124 Standard_Integer iter=0;
125 Standard_Real tols2 = 0.5*tol;
126 Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
127 a=x1;b=c=x2;fa=y1; fb=fc=y2;
128 for(iter=1;iter<=ITMAX;iter++) {
129 if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
132 if(Abs(fc)<Abs(fb)) {
133 a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
135 tol1 = EPSEPS*Abs(b)+tols2;
137 if(Abs(xm)<tol1 || fb==0) {
138 //-- On tente une iteration de newton
139 Standard_Real Xp,Yp,Dp;
140 Standard_Integer itern=5;
144 Ok = F.Values(Xp,Yp,Dp);
147 if(Dp>1e-10 || Dp<-1e-10) {
150 if(Xp<=x2 && Xp>=x1) {
151 F.Value(Xp,Yp); Yp-=K;
152 if(Abs(Yp)<Abs(fb)) {
160 while(Ok && --itern>=0);
162 AppendRoot(Sol,NbStateSol,b,F,K,dX);
165 if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
174 p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
175 q=(q-1.0)*(r-1.0)*(s-1.0);
181 min1=3.0*xm*q-Abs(tol1*q);
183 if((p+p)<( (min1<min2) ? min1 : min2)) {
202 if(xm>=0) b+=Abs(tol1);
209 cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
215 #define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement
216 //#define MATH_FUNCTIONROOTS_OLDCODE // Ancien
217 //#define MATH_FUNCTIONROOTS_CHECK // Check
219 math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
220 const Standard_Real A,
221 const Standard_Real B,
222 const Standard_Integer NbSample,
223 const Standard_Real _EpsX,
224 const Standard_Real EpsF,
225 const Standard_Real EpsNull,
226 const Standard_Real K )
230 cout << "---- Debut de math_FunctionRoots ----" << endl;
236 TColStd_SequenceOfReal StaticSol;
240 #ifdef MATH_FUNCTIONROOTS_NEWCODE
242 Done = Standard_True;
245 Standard_Integer N=NbSample;
246 //-- ------------------------------------------------------------
247 //-- Verifications de bas niveau
256 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
257 Standard_Real EpsX = _EpsX;
258 Standard_Real DeltaU = Abs(X0)+Abs(XN);
259 Standard_Real NEpsX = 0.0000000001 * DeltaU;
264 //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
265 //-- A .............................................................. B
266 //-- X0 X1 X2 ........................................ Xn-1 Xn
270 double dx = (XN-X0)/N;
271 TColStd_Array1OfReal ptrval(0, N);
272 Standard_Integer Nvalid = -1;
273 Standard_Real aux = 0;
274 for(i=0; i<=N ; i++,X+=dx) {
277 if(Ok) ptrval(++Nvalid) = aux - K;
280 //-- Toute la fonction est nulle ?
283 Done = Standard_False;
287 AllNull=Standard_True;
288 // for(i=0;AllNull && i<=N;i++) {
289 for(i=0;AllNull && i<=N;i++) {
290 if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
291 AllNull=Standard_False;
295 //-- tous les points echantillons sont dans la tolerance
299 //-- Il y a des points hors tolerance
300 //-- on detecte les changements de signes STRICTS
301 Standard_Integer ip1;
302 // Standard_Boolean chgtsign=Standard_False;
303 Standard_Real tol = EpsX;
305 for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
309 if(ptrval(ip1)>0.0) {
310 //-- --------------------------------------------------
311 //-- changement de signe dans Xi Xi+1
312 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
316 if(ptrval(ip1)<0.0) {
317 //-- --------------------------------------------------
318 //-- changement de signe dans Xi Xi+1
319 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
323 //-- On detecte les cas ou la fct s annule sur des Xi et est
324 //-- non nulle au voisinage de Xi
326 //-- On prend 2 points u0,u1 au voisinage de Xi
327 //-- Si (F(u0)-K)*(F(u1)-K) <0 on lance une recherche
328 //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
329 for(i=0; i<=N; i++) {
331 // Standard_Real Val,Deriv;
335 u0=dx*0.5; u1=X+u0; u0+=X;
342 F.Value(u0,y0); y0-=K;
343 F.Value(u1,y1); y1-=K;
345 Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
348 if(y0!=0.0 || y1!=0.0) {
349 AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
354 //-- --------------------------------------------------------------------------------
355 //-- Il faut traiter differement le cas des points en bout :
356 if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
357 AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
359 if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
360 AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
363 //-- --------------------------------------------------------------------------------
364 //-- --------------------------------------------------------------------------------
365 //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
366 //-- un maximum avec f(x)<0
367 //-- On reprend une discretisation plus fine au voisinage de ces extremums
369 //-- Recherche d un minima positif
370 Standard_Real xm,ym,dym,xm1,xp1;
371 Standard_Real majdx = 5.0*dx;
372 Standard_Boolean Rediscr;
373 // Standard_Real ptrvalbis[MAXBIS];
374 Standard_Integer im1=0;
376 for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
377 Rediscr = Standard_False;
380 if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
381 //-- Peut on traverser l axe Ox
382 //-- -------------- Estimation a partir de Xim1
385 F.Values(xm1,ym,dym); ym-=K;
386 if(dym<-1e-10 || dym>1e-10) { // normalement dym < 0
387 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
388 if(t<majdx && t > -majdx) {
389 Rediscr = Standard_True;
392 //-- -------------- Estimation a partir de Xip1
393 if(Rediscr==Standard_False) {
395 if(xp1 > XN ) xp1=XN;
396 F.Values(xp1,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;
406 else if(ptrval(i)<0.0) {
407 if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
408 //-- Peut on traverser l axe Ox
409 //-- -------------- Estimation a partir de Xim1
412 F.Values(xm1,ym,dym); ym-=K;
413 if(dym>1e-10 || dym<-1e-10) { // normalement dym > 0
414 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
415 if(t<majdx && t > -majdx) {
416 Rediscr = Standard_True;
419 //-- -------------- Estimation a partir de Xim1
420 if(Rediscr==Standard_False) {
423 F.Values(xm1,ym,dym); ym-=K;
424 if(dym>1e-10 || dym<-1e-10) { // normalement dym < 0
425 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
426 if(t<majdx && t > -majdx) {
427 Rediscr = Standard_True;
434 Standard_Real x0 = xm - dx;
435 Standard_Real x3 = xm + dx;
436 if (x0 < X0) x0 = X0;
437 if (x3 > XN) x3 = XN;
438 Standard_Real aSolX1 = 0., aSolX2 = 0.;
439 Standard_Real aVal1 = 0., aVal2 = 0.;
440 Standard_Real aDer1 = 0., aDer2 = 0.;
441 Standard_Boolean isSol1 = Standard_False;
442 Standard_Boolean isSol2 = Standard_False;
443 //-- ----------------------------------------------------
444 //-- Find minimum of the function |F| between x0 and x3
445 //-- by searching for the zero of the function derivative
446 DerivFunction aDerF(F);
447 math_BracketedRoot aBR(aDerF, x0, x3, EpsX);
451 F.Value(aSolX1, aVal1);
455 isSol1 = Standard_True;
460 //-- --------------------------------------------------
461 //-- On recherche un extrema entre x0 et x3
462 //-- x1 et x2 sont tels que x0<x1<x2<x3
463 //-- et |f(x0)| > |f(x1)| et |f(x3)| > |f(x2)|
465 //-- En entree : a=xm-dx b=xm c=xm+dx
466 Standard_Real x1, x2, f0, f3;
467 Standard_Real R = 0.61803399;
468 Standard_Real C = 1.0 - R;
469 Standard_Real tolCR = NEpsX*10.0;
472 Standard_Boolean recherche_minimum = (f0 > 0.0);
474 if (Abs(x3 - xm) > Abs(x0 - xm)) { x1 = xm; x2 = xm + C*(x3 - xm); }
475 else { x2 = xm; x1 = xm - C*(xm - x0); }
476 Standard_Real f1, f2;
477 F.Value(x1, f1); f1 -= K;
478 F.Value(x2, f2); f2 -= K;
479 //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
480 Standard_Real tolX = 0.001 * NEpsX;
481 while (Abs(x3 - x0) > tolCR*(Abs(x1) + Abs(x2)) && (Abs(x1 - x2) > tolX)) {
482 //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
483 //-- x0,f0,x1,f1,x2,f2,x3,f3);
484 if (recherche_minimum) {
486 x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
487 f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
490 x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
491 f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
496 x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
497 f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
500 x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
501 f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
504 //-- On ne fait pas que chercher des extremas. Il faut verifier
505 //-- si on ne tombe pas sur une racine
507 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
508 Solve(F, K, x0, f0, x1, f1, tol, NEpsX, Sol, NbStateSol);
511 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
512 Solve(F, K, x2, f2, x3, f3, tol, NEpsX, Sol, NbStateSol);
515 if ((recherche_minimum && f1<f2) || (!recherche_minimum && f1>f2)) {
516 //-- x1,f(x1) minimum
517 if (Abs(f1) < EpsF) {
518 isSol2 = Standard_True;
524 //-- x2.f(x2) minimum
525 if (Abs(f2) < EpsF) {
526 isSol2 = Standard_True;
531 // Choose the best solution between aSolX1, aSolX2
532 if (isSol1 && isSol2)
534 if (aVal2 - aVal1 > EpsF)
535 AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
536 else if (aVal1 - aVal2 > EpsF)
537 AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
541 F.Derivative(aSolX2, aDer2);
544 AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
546 AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
550 AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
552 AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
553 } //-- Recherche d un extrema
558 #ifdef MATH_FUNCTIONROOTS_CHECK
561 Standard_Integer n=Sol.Length();
562 for(Standard_Integer ii=1;ii<=n;ii++) {
563 StaticSol.Append(Sol.Value(ii));
572 #ifdef MATH_FUNCTIONROOTS_OLDCODE
574 //-- ********************************************************************************
575 //-- ANCIEN TRAITEMENT
576 //-- ********************************************************************************
579 // calculate all the real roots of a function within the range
580 // A..B. whitout condition on A and B
581 // a solution X is found when
582 // abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
583 // The function is considered as null between A and B if
584 // abs(F-K) <= EpsNull within this range.
585 Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
586 //-- Il ne faut pas EpsX = 0.000...001 car dans ce cas
587 //-- U + Nn*EpsX == U
588 Standard_Real Lowr,Upp;
589 Standard_Real Increment;
591 Standard_Real FLowr,FUpp,DFLowr,DFUpp;
593 Standard_Real Fxu,DFxu,FFxu,DFFxu;
594 Standard_Real Fyu,DFyu,FFyu,DFFyu;
595 Standard_Boolean Finish;
597 Standard_Integer Nbiter = 30;
598 Standard_Integer Iter;
599 Standard_Real Ambda,T;
600 Standard_Real AA,BB,CC;
602 Standard_Real Alfa1=0,Alfa2;
603 Standard_Real OldDF = RealLast();
604 Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
607 Done = Standard_False;
609 StdFail_NotDone_Raise_if(NbSample <= 0, " ");
622 Increment = (Upp-Lowr)/NbSample;
623 StdFail_NotDone_Raise_if(Increment < EpsX, " ");
624 Done = Standard_True;
625 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
626 Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
627 Standard_Real NEpsX = 0.0000000001 * DeltaU;
630 //-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
633 Null2 = EpsNull*EpsNull;
635 Ok = F.Values(Lowr,FLowr,DFLowr);
638 Done = Standard_False;
644 Ok = F.Values(Upp,FUpp,DFUpp);
647 Done = Standard_False;
656 Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
659 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
660 AllNull = ( FFyu <= Null2 );
672 Fyu = FLowr + (U-Lowr)*DFLowr;
676 Fyu = FUpp + (U-Upp)*DFUpp;
680 Ok = F.Values(U,Fyu,DFyu);
683 Done = Standard_False;
690 DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
692 if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
694 if (AllNull) { //rechercher les vraix zeros depuis le debut
696 AllNull = Standard_False;
698 Fxu = FLowr-EpsX*DFLowr;
701 DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
703 Ok = F.Values(U,Fyu,DFyu);
706 Done = Standard_False;
712 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
714 Standard_Real FxuFyu=Fxu*Fyu;
716 if ( (DFFyu > 0. && DFFxu <= 0.)
717 || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
718 || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
719 || (FxuFyu <= 0.) ) {
720 // recherche d 1 minimun possible
721 Finish = Standard_False;
728 // chercher si f peut s annuler pour eviter
729 // des iterations inutiles
730 if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
731 Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
733 Finish = Standard_True;
734 FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
736 else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
737 (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
739 Finish = Standard_True;
741 FFi = FFyu; // pour recalculer yu
743 else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
744 (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
746 Finish = Standard_True;
748 FFi = FFxu; // pour recalculer U
751 else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
753 Finish = Standard_True;
757 else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
759 Finish = Standard_True;
765 // calcul des 2 solutions annulant la derivee de l interpolation cubique
766 // Ambda*t=(U-Xu) F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
767 // df=aa*t*t+2*bb*t+cc
769 AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
770 BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
773 if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
778 math_DirectPolynomialRoots Solv(AA,BB,CC);
779 if ( !Solv.InfiniteRoots() ) {
780 Nn=Solv.NbSolutions();
782 if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
783 else Finish = Standard_True;
786 Alfa1 = Solv.Value(1);
788 Alfa2 = Solv.Value(2);
794 if (Alfa1 > 1. || Alfa2 < 0.){
795 // resolution par dichotomie
796 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
797 else Finish = Standard_True;
799 else if ( Alfa1 < 0. ||
800 ( DFFxu > 0. && DFFyu >= 0.) ) {
802 //(cas changement de signe de la distance signee sans
803 // changement de signe de la derivee:
804 //cas de 'presque'tangence avec 2
805 // solutions proches) ,on prend la plus grane racine
807 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
808 else Finish = Standard_True;
815 else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
816 //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
817 else Finish = Standard_True;
820 // petits tests pour diminuer le nombre d iterations
821 if (Alfa1 <= EpsX ) {
824 else if (Alfa1 >= (1.-EpsX) ) {
825 Alfa1 = Alfa1+Alfa1-1.;
827 Alfa1 = Ambda*(1. - Alfa1);
830 AA = FLowr + (U - Lowr)*DFLowr;
833 else if ( U >= Upp ) {
834 AA = FUpp + (U - Upp)*DFUpp;
838 Ok = F.Values(U,AA,BB);
841 Done = Standard_False;
849 if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
856 if (Alfa1 > Ambda*0.5) {
858 // determination d 1 autre borne pour diviser
859 //le nouvel intervalle par 2 au -
862 AA = FLowr + (Xu - Lowr)*DFLowr;
865 else if ( Xu >= Upp ) {
866 AA = FUpp + (Xu - Upp)*DFUpp;
870 Ok = F.Values(Xu,AA,BB);
873 Done = Standard_False;
881 if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
890 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
891 // changement de signe sur l intervalle u,U
896 FFi = Min(FFxu,FFyu);
901 else { Ambda = Alfa1;}
903 else { Ambda = Alfa1;}
910 if ( (Ambda-Alfa1) > Ambda*0.5 ) {
912 Xu = U - (Ambda-Alfa1)*0.5;
914 AA = FLowr + (Xu - Lowr)*DFLowr;
917 else if ( Xu >= Upp ) {
918 AA = FUpp + (Xu - Upp)*DFUpp;
922 Ok = F.Values(Xu,AA,BB);
925 Done = Standard_False;
933 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
938 Ambda = ( Ambda - Alfa1 )*0.5;
941 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
946 Ambda = ( Ambda - Alfa1 )*0.5;
948 FFi = Min(FFxu , FFyu);
953 Ambda = Ambda - Alfa1;
958 Ambda = Ambda - Alfa1;
962 if (Abs(FFxu) <= Standard_Underflow ||
963 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
965 Finish = Standard_True;
966 if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
969 else if (Abs(FFyu) <= Standard_Underflow ||
970 (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
972 Finish = Standard_True;
973 if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
978 Finish = Iter >= Nbiter || (Ambda <= EpsX &&
979 ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
982 } // fin interpolation cubique
984 // restitution du meilleur resultat
986 if ( FFxu < FFi ) { U = U + T -Ambda;}
987 else if ( FFyu < FFi ) { U = U + T;}
989 if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
990 U = Max( Lowr , Min( U , Upp));
993 if ( Abs(FFi) < EpsF ) {
995 if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
996 else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
997 else if (Fxu*Fyu > 0.) { AA = 0.;}
998 else { AA = Fyu-Fxu;}
999 if (!Sol.IsEmpty()) {
1000 if (Abs(Sol.Last() - U) > 5.*EpsX
1001 || (OldDF != RealLast() && OldDF*AA < 0.) ) {
1003 NbStateSol.Append(F.GetStateNumber());
1008 NbStateSol.Append(F.GetStateNumber());
1016 while ( Nn < 1000000 && DFFyu <= 0.) {
1017 // on repart d 1 pouyem plus loin
1020 AA = FLowr + (U - Lowr)*DFLowr;
1023 else if ( U >= Upp ) {
1024 AA = FUpp + (U - Upp)*DFUpp;
1028 Ok = F.Values(U,AA,BB);
1039 DFFyu = 2.*DFyu*Fyu;
1047 #ifdef MATH_FUNCTIONROOTS_CHECK
1049 Standard_Integer n1=StaticSol.Length();
1050 Standard_Integer n2=Sol.Length();
1052 printf("\n mathFunctionRoots : n1=%d n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
1053 for(Standard_Integer x1=1; x1<=n1;x1++) {
1054 Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
1055 printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
1058 for(Standard_Integer x2=1; x2<=n2; x2++) {
1059 Standard_Real v; F.Value(Sol(x2),v); v-=K;
1060 printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
1064 Standard_Integer n=n1;
1066 for(Standard_Integer i=1;i<=n;i++) {
1067 Standard_Real t = Sol(i)-StaticSol(i);
1069 printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
1080 void math_FunctionRoots::Dump(Standard_OStream& o) const
1083 o << "math_FunctionRoots ";
1085 o << " Status = Done \n";
1086 o << " Number of solutions = " << Sol.Length() << endl;
1087 for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
1088 o << " Solution Number " << i << "= " << Sol.Value(i) << endl;
1092 o<< " Status = not Done \n";