mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-05-01 10:26:12 +03:00
f4dee9bb20
Added DRAW command dsetsignal, resetting OSD signal handler with either armed or disabled FPE handler, according to an option. If called without arguments, it sets FPE handler only if environment variable OSD_FPE is defined (with value different from 0). On start, DRAW calls dsetsignal to set FPE signal if CSF_FPE is defined. Test bugs fclasses bug6143 uses dsetsignal to set FPE handler unconditionally before the test command, and resets it to default at the end. A number of changes in the code have been done in order to fix floating point exceptions that became generated after enabling signals: - Global functions Sinh() and Cosh() defined in Standard_Real.hxx are improved to raise Standard_NumericError exception if argument is too big (greater than 710.47586), instead of relying on system treatment of floating point overflow. These functions are used instead of sinh and cosh in ElCLib.cxx. - Maximal value of parameter on hyperbola is restricted by 23 (corresponding to ~1e10 in 3d) in order to avoid FP overflow in Extrema_GenExtCS.cxx, ShapeFix_EdgeProjAux.cxx. - Interface of the root curve adaptor class Adaptor3d_Curve has been updated to add new virtual methods BasisCurve and OffsetValue. They complement the adaptor for the case of offset curves. These methods are used in Extrema_GenExtCS.cxx to restrict domain search in the case of offset of hyperbola, in order to get rid of floating point overflow. All classes inheriting Adaptor3d_Curve have been changed to implement the new virtual methods. - Protection against division by zero has been implemented in ApproxInt_KnotTools.cxx, BRepClass3d_SClassifier.cxx, BRepGProp_Face.cxx, BRepMesh_FastDiscretFace.cxx, Geom2dGcc_Circ2d2TanOnIter.cxx, Geom2dInt_Geom2dCurveTool.cxx, IntPolyh_MaillageAffinage.cxx. - Protection against calling of math functions of infinite arguments has been added in BRepCheck_Edge.cxx, BRepLib.cxx, CSLib_NormalPolyDef.cxx, Extrema_FuncExtPC.gxx, Extrema_GExtPC.gxx, Extrema_GLocateExtPC.gxx, Intf_InterferencePolygonPolyhedron.gxx, ShapeAnalysis_Surface.cxx, ShapeAnalysis_TransferParametersProj.cxx, ShapeAnalysis_Wire.cxx, math_FunctionSetRoot.cxx. - Proper initialization of local variables is done in BOPAlgo_PaveFiller_6.cxx, XSDRAWSTLVRML.cxx. - Inconsistent usage of Standard_Boolean* to access integer data in HLR (caused by #27772) is corrected Some test cases have been updated to actual state.
1031 lines
26 KiB
C++
1031 lines
26 KiB
C++
// Copyright (c) 1997-1999 Matra Datavision
|
|
// Copyright (c) 1999-2014 OPEN CASCADE SAS
|
|
//
|
|
// This file is part of Open CASCADE Technology software library.
|
|
//
|
|
// This library is free software; you can redistribute it and/or modify it under
|
|
// the terms of the GNU Lesser General Public License version 2.1 as published
|
|
// by the Free Software Foundation, with special exception defined in the file
|
|
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
|
|
// distribution for complete text of the license and disclaimer of any warranty.
|
|
//
|
|
// Alternatively, this file may be used under the terms of Open CASCADE
|
|
// commercial license or contractual agreement.
|
|
|
|
//#ifndef OCCT_DEBUG
|
|
#define No_Standard_RangeError
|
|
#define No_Standard_OutOfRange
|
|
#define No_Standard_DimensionError
|
|
|
|
//#endif
|
|
|
|
#include <math_DirectPolynomialRoots.hxx>
|
|
#include <math_FunctionRoots.hxx>
|
|
#include <math_FunctionWithDerivative.hxx>
|
|
#include <Standard_RangeError.hxx>
|
|
#include <StdFail_NotDone.hxx>
|
|
#include <TColStd_Array1OfReal.hxx>
|
|
|
|
#include <stdio.h>
|
|
#define ITMAX 100
|
|
#define EPS 1e-14
|
|
#define EPSEPS 2e-14
|
|
#define MAXBIS 100
|
|
|
|
#ifdef OCCT_DEBUG
|
|
static Standard_Boolean myDebug = 0;
|
|
static Standard_Integer nbsolve = 0;
|
|
#endif
|
|
|
|
static void AppendRoot(TColStd_SequenceOfReal& Sol,
|
|
TColStd_SequenceOfInteger& NbStateSol,
|
|
const Standard_Real X,
|
|
math_FunctionWithDerivative& F,
|
|
// const Standard_Real K,
|
|
const Standard_Real ,
|
|
const Standard_Real dX) {
|
|
|
|
Standard_Integer n=Sol.Length();
|
|
Standard_Real t;
|
|
#ifdef OCCT_DEBUG
|
|
if (myDebug) {
|
|
cout << " Ajout de la solution numero : " << n+1 << endl;
|
|
cout << " Valeur de la racine :" << X << endl;
|
|
}
|
|
#endif
|
|
|
|
if(n==0) {
|
|
Sol.Append(X);
|
|
F.Value(X,t);
|
|
NbStateSol.Append(F.GetStateNumber());
|
|
}
|
|
else {
|
|
Standard_Integer i=1;
|
|
Standard_Integer pl=n+1;
|
|
while(i<=n) {
|
|
t=Sol.Value(i);
|
|
if(t >= X) {
|
|
pl=i;
|
|
i=n;
|
|
}
|
|
if(Abs(X-t) <= dX) {
|
|
pl=0;
|
|
i=n;
|
|
}
|
|
i++;
|
|
} //-- while
|
|
if(pl>n) {
|
|
Sol.Append(X);
|
|
F.Value(X,t);
|
|
NbStateSol.Append(F.GetStateNumber());
|
|
}
|
|
else if(pl>0) {
|
|
Sol.InsertBefore(pl,X);
|
|
F.Value(X,t);
|
|
NbStateSol.InsertBefore(pl,F.GetStateNumber());
|
|
}
|
|
}
|
|
}
|
|
|
|
static void Solve(math_FunctionWithDerivative& F,
|
|
const Standard_Real K,
|
|
const Standard_Real x1,
|
|
const Standard_Real y1,
|
|
const Standard_Real x2,
|
|
const Standard_Real y2,
|
|
const Standard_Real tol,
|
|
const Standard_Real dX,
|
|
TColStd_SequenceOfReal& Sol,
|
|
TColStd_SequenceOfInteger& NbStateSol) {
|
|
#ifdef OCCT_DEBUG
|
|
if (myDebug) {
|
|
cout <<"--> Resolution :" << ++nbsolve << endl;
|
|
cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
|
|
cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
|
|
}
|
|
#endif
|
|
|
|
Standard_Integer iter=0;
|
|
Standard_Real tols2 = 0.5*tol;
|
|
Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
|
|
a=x1;b=c=x2;fa=y1; fb=fc=y2;
|
|
for(iter=1;iter<=ITMAX;iter++) {
|
|
if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
|
|
c=a; fc=fa; e=d=b-a;
|
|
}
|
|
if(Abs(fc)<Abs(fb)) {
|
|
a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
|
|
}
|
|
tol1 = EPSEPS*Abs(b)+tols2;
|
|
xm=0.5*(c-b);
|
|
if(Abs(xm)<tol1 || fb==0) {
|
|
//-- On tente une iteration de newton
|
|
Standard_Real Xp,Yp,Dp;
|
|
Standard_Integer itern=5;
|
|
Standard_Boolean Ok;
|
|
Xp=b;
|
|
do {
|
|
Ok = F.Values(Xp,Yp,Dp);
|
|
if(Ok) {
|
|
Ok=Standard_False;
|
|
if(Dp>1e-10 || Dp<-1e-10) {
|
|
Xp = Xp-(Yp-K)/Dp;
|
|
}
|
|
if(Xp<=x2 && Xp>=x1) {
|
|
F.Value(Xp,Yp); Yp-=K;
|
|
if(Abs(Yp)<Abs(fb)) {
|
|
b=Xp;
|
|
fb=Yp;
|
|
Ok=Standard_True;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
while(Ok && --itern>=0);
|
|
|
|
AppendRoot(Sol,NbStateSol,b,F,K,dX);
|
|
return;
|
|
}
|
|
if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
|
|
s=fb/fa;
|
|
if(a==c) {
|
|
p=xm*s; p+=p;
|
|
q=1.0-s;
|
|
}
|
|
else {
|
|
q=fa/fc;
|
|
r=fb/fc;
|
|
p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
|
|
q=(q-1.0)*(r-1.0)*(s-1.0);
|
|
}
|
|
if(p>0.0) {
|
|
q=-q;
|
|
}
|
|
p=Abs(p);
|
|
min1=3.0*xm*q-Abs(tol1*q);
|
|
min2=Abs(e*q);
|
|
if((p+p)<( (min1<min2) ? min1 : min2)) {
|
|
e=d;
|
|
d=p/q;
|
|
}
|
|
else {
|
|
d=xm;
|
|
e=d;
|
|
}
|
|
}
|
|
else {
|
|
d=xm;
|
|
e=d;
|
|
}
|
|
a=b;
|
|
fa=fb;
|
|
if(Abs(d)>tol1) {
|
|
b+=d;
|
|
}
|
|
else {
|
|
if(xm>=0) b+=Abs(tol1);
|
|
else b+=-Abs(tol1);
|
|
}
|
|
F.Value(b,fb);
|
|
fb-=K;
|
|
}
|
|
#ifdef OCCT_DEBUG
|
|
cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
|
|
#endif
|
|
}
|
|
|
|
#define NEWSEQ 1
|
|
|
|
#define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement
|
|
//#define MATH_FUNCTIONROOTS_OLDCODE // Ancien
|
|
//#define MATH_FUNCTIONROOTS_CHECK // Check
|
|
|
|
math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
|
|
const Standard_Real A,
|
|
const Standard_Real B,
|
|
const Standard_Integer NbSample,
|
|
const Standard_Real _EpsX,
|
|
const Standard_Real EpsF,
|
|
const Standard_Real EpsNull,
|
|
const Standard_Real K )
|
|
{
|
|
#ifdef OCCT_DEBUG
|
|
if (myDebug) {
|
|
cout << "---- Debut de math_FunctionRoots ----" << endl;
|
|
nbsolve = 0;
|
|
}
|
|
#endif
|
|
|
|
#if NEWSEQ
|
|
TColStd_SequenceOfReal StaticSol;
|
|
#endif
|
|
Sol.Clear();
|
|
NbStateSol.Clear();
|
|
#ifdef MATH_FUNCTIONROOTS_NEWCODE
|
|
{
|
|
Done = Standard_True;
|
|
Standard_Real X0=A;
|
|
Standard_Real XN=B;
|
|
Standard_Integer N=NbSample;
|
|
//-- ------------------------------------------------------------
|
|
//-- Verifications de bas niveau
|
|
if(B<A) {
|
|
X0=B;
|
|
XN=A;
|
|
}
|
|
N*=2;
|
|
if(N < 20) {
|
|
N=20;
|
|
}
|
|
//-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
|
|
Standard_Real EpsX = _EpsX;
|
|
Standard_Real DeltaU = Abs(X0)+Abs(XN);
|
|
Standard_Real NEpsX = 0.0000000001 * DeltaU;
|
|
if(EpsX < NEpsX) {
|
|
EpsX = NEpsX;
|
|
}
|
|
|
|
//-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
|
|
//-- A .............................................................. B
|
|
//-- X0 X1 X2 ........................................ Xn-1 Xn
|
|
Standard_Integer i;
|
|
Standard_Real X=X0;
|
|
Standard_Boolean Ok;
|
|
double dx = (XN-X0)/N;
|
|
TColStd_Array1OfReal ptrval(0, N);
|
|
Standard_Integer Nvalid = -1;
|
|
Standard_Real aux = 0;
|
|
for(i=0; i<=N ; i++,X+=dx) {
|
|
if( X > XN) X=XN;
|
|
Ok=F.Value(X,aux);
|
|
if(Ok) ptrval(++Nvalid) = aux - K;
|
|
// ptrval(i)-=K;
|
|
}
|
|
//-- Toute la fonction est nulle ?
|
|
|
|
if( Nvalid < N ) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
AllNull=Standard_True;
|
|
// for(i=0;AllNull && i<=N;i++) {
|
|
for(i=0;AllNull && i<=N;i++) {
|
|
if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
|
|
AllNull=Standard_False;
|
|
}
|
|
}
|
|
if(AllNull) {
|
|
//-- tous les points echantillons sont dans la tolerance
|
|
|
|
}
|
|
else {
|
|
//-- Il y a des points hors tolerance
|
|
//-- on detecte les changements de signes STRICTS
|
|
Standard_Integer ip1;
|
|
// Standard_Boolean chgtsign=Standard_False;
|
|
Standard_Real tol = EpsX;
|
|
Standard_Real X2;
|
|
for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
|
|
X2=X+dx;
|
|
if(X2 > XN) X2 = XN;
|
|
if(ptrval(i)<0.0) {
|
|
if(ptrval(ip1)>0.0) {
|
|
//-- --------------------------------------------------
|
|
//-- changement de signe dans Xi Xi+1
|
|
Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
|
|
}
|
|
}
|
|
else {
|
|
if(ptrval(ip1)<0.0) {
|
|
//-- --------------------------------------------------
|
|
//-- changement de signe dans Xi Xi+1
|
|
Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
|
|
}
|
|
}
|
|
}
|
|
//-- On detecte les cas ou la fct s annule sur des Xi et est
|
|
//-- non nulle au voisinage de Xi
|
|
//--
|
|
//-- On prend 2 points u0,u1 au voisinage de Xi
|
|
//-- Si (F(u0)-K)*(F(u1)-K) <0 on lance une recherche
|
|
//-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
|
|
for(i=0; i<=N; i++) {
|
|
if(ptrval(i)==0) {
|
|
// Standard_Real Val,Deriv;
|
|
X=X0+i*dx;
|
|
if (X>XN) X=XN;
|
|
Standard_Real u0,u1;
|
|
u0=dx*0.5; u1=X+u0; u0+=X;
|
|
if(u0<X0) u0=X0;
|
|
if(u0>XN) u0=XN;
|
|
if(u1<X0) u1=X0;
|
|
if(u1>XN) u1=XN;
|
|
|
|
Standard_Real y0,y1;
|
|
F.Value(u0,y0); y0-=K;
|
|
F.Value(u1,y1); y1-=K;
|
|
if(y0*y1 < 0.0) {
|
|
Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
|
|
}
|
|
else {
|
|
if(y0!=0.0 || y1!=0.0) {
|
|
AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//-- --------------------------------------------------------------------------------
|
|
//-- Il faut traiter differement le cas des points en bout :
|
|
if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
|
|
AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
|
|
}
|
|
if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
|
|
AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
|
|
}
|
|
|
|
//-- --------------------------------------------------------------------------------
|
|
//-- --------------------------------------------------------------------------------
|
|
//-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
|
|
//-- un maximum avec f(x)<0
|
|
//-- On reprend une discretisation plus fine au voisinage de ces extremums
|
|
//--
|
|
//-- Recherche d un minima positif
|
|
Standard_Real xm,ym,dym,xm1,xp1;
|
|
Standard_Real majdx = 5.0*dx;
|
|
Standard_Boolean Rediscr;
|
|
// Standard_Real ptrvalbis[MAXBIS];
|
|
Standard_Integer im1=0;
|
|
ip1=2;
|
|
for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
|
|
Rediscr = Standard_False;
|
|
if (xm > XN) xm=XN;
|
|
if(ptrval(i)>0.0) {
|
|
if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
|
|
//-- Peut on traverser l axe Ox
|
|
//-- -------------- Estimation a partir de Xim1
|
|
xm1=xm-dx;
|
|
if(xm1 < X0) xm1=X0;
|
|
F.Values(xm1,ym,dym); ym-=K;
|
|
if(dym<-1e-10 || dym>1e-10) { // normalement dym < 0
|
|
Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
|
|
if(t<majdx && t > -majdx) {
|
|
Rediscr = Standard_True;
|
|
}
|
|
}
|
|
//-- -------------- Estimation a partir de Xip1
|
|
if(Rediscr==Standard_False) {
|
|
xp1=xm+dx;
|
|
if(xp1 > XN ) xp1=XN;
|
|
F.Values(xp1,ym,dym); ym-=K;
|
|
if(dym<-1e-10 || dym>1e-10) { // normalement dym > 0
|
|
Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
|
|
if(t<majdx && t > -majdx) {
|
|
Rediscr = Standard_True;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else if(ptrval(i)<0.0) {
|
|
if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
|
|
//-- Peut on traverser l axe Ox
|
|
//-- -------------- Estimation a partir de Xim1
|
|
xm1=xm-dx;
|
|
if(xm1 < X0) xm1=X0;
|
|
F.Values(xm1,ym,dym); ym-=K;
|
|
if(dym>1e-10 || dym<-1e-10) { // normalement dym > 0
|
|
Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
|
|
if(t<majdx && t > -majdx) {
|
|
Rediscr = Standard_True;
|
|
}
|
|
}
|
|
//-- -------------- Estimation a partir de Xim1
|
|
if(Rediscr==Standard_False) {
|
|
xm1=xm-dx;
|
|
if(xm1 < X0) xm1=X0;
|
|
F.Values(xm1,ym,dym); ym-=K;
|
|
if(dym>1e-10 || dym<-1e-10) { // normalement dym < 0
|
|
Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
|
|
if(t<majdx && t > -majdx) {
|
|
Rediscr = Standard_True;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if(Rediscr) {
|
|
//-- --------------------------------------------------
|
|
//-- On recherche un extrema entre x0 et x3
|
|
//-- x1 et x2 sont tels que x0<x1<x2<x3
|
|
//-- et |f(x0)| > |f(x1)| et |f(x3)| > |f(x2)|
|
|
//--
|
|
//-- En entree : a=xm-dx b=xm c=xm+dx
|
|
Standard_Real x0,x1,x2,x3,f0,f3;
|
|
Standard_Real R=0.61803399;
|
|
Standard_Real C=1.0-R;
|
|
Standard_Real tolCR=NEpsX*10.0;
|
|
f0=ptrval(im1);
|
|
f3=ptrval(ip1);
|
|
x0=xm-dx;
|
|
x3=xm+dx;
|
|
if(x0 < X0) x0=X0;
|
|
if(x3 > XN) x3=XN;
|
|
Standard_Boolean recherche_minimum = (f0>0.0);
|
|
|
|
if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); }
|
|
else { x2=xm; x1=xm-C*(xm-x0); }
|
|
Standard_Real f1,f2;
|
|
F.Value(x1,f1); f1-=K;
|
|
F.Value(x2,f2); f2-=K;
|
|
//-- printf("\n *************** RECHERCHE MINIMUM **********\n");
|
|
Standard_Real tolX = 0.001 * NEpsX;
|
|
while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > tolX)) {
|
|
//-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
|
|
//-- x0,f0,x1,f1,x2,f2,x3,f3);
|
|
if(recherche_minimum) {
|
|
if(f2<f1) {
|
|
x0=x1; x1=x2; x2=R*x1+C*x3;
|
|
f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
|
|
}
|
|
else {
|
|
x3=x2; x2=x1; x1=R*x2+C*x0;
|
|
f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
|
|
}
|
|
}
|
|
else {
|
|
if(f2>f1) {
|
|
x0=x1; x1=x2; x2=R*x1+C*x3;
|
|
f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
|
|
}
|
|
else {
|
|
x3=x2; x2=x1; x1=R*x2+C*x0;
|
|
f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
|
|
}
|
|
}
|
|
//-- On ne fait pas que chercher des extremas. Il faut verifier
|
|
//-- si on ne tombe pas sur une racine
|
|
if(f1*f0 <0.0) {
|
|
//-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
|
|
Solve(F,K,x0,f0,x1,f1,tol,NEpsX,Sol,NbStateSol);
|
|
}
|
|
if(f2*f3 <0.0) {
|
|
//-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
|
|
Solve(F,K,x2,f2,x3,f3,tol,NEpsX,Sol,NbStateSol);
|
|
}
|
|
}
|
|
if(f1<f2) {
|
|
//-- x1,f(x1) minimum
|
|
if(Abs(f1) < EpsF) {
|
|
AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
|
|
}
|
|
}
|
|
else {
|
|
//-- x2.f(x2) minimum
|
|
if(Abs(f2) < EpsF) {
|
|
AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
|
|
}
|
|
}
|
|
} //-- Recherche d un extrema
|
|
} //-- for
|
|
}
|
|
|
|
#if NEWSEQ
|
|
#ifdef MATH_FUNCTIONROOTS_CHECK
|
|
{
|
|
StaticSol.Clear();
|
|
Standard_Integer n=Sol.Length();
|
|
for(Standard_Integer ii=1;ii<=n;ii++) {
|
|
StaticSol.Append(Sol.Value(ii));
|
|
}
|
|
Sol.Clear();
|
|
NbStateSol.Clear();
|
|
}
|
|
#endif
|
|
#endif
|
|
#endif
|
|
}
|
|
#ifdef MATH_FUNCTIONROOTS_OLDCODE
|
|
{
|
|
//-- ********************************************************************************
|
|
//-- ANCIEN TRAITEMENT
|
|
//-- ********************************************************************************
|
|
|
|
|
|
// calculate all the real roots of a function within the range
|
|
// A..B. whitout condition on A and B
|
|
// a solution X is found when
|
|
// abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
|
|
// The function is considered as null between A and B if
|
|
// abs(F-K) <= EpsNull within this range.
|
|
Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
|
|
//-- Il ne faut pas EpsX = 0.000...001 car dans ce cas
|
|
//-- U + Nn*EpsX == U
|
|
Standard_Real Lowr,Upp;
|
|
Standard_Real Increment;
|
|
Standard_Real Null2;
|
|
Standard_Real FLowr,FUpp,DFLowr,DFUpp;
|
|
Standard_Real U,Xu;
|
|
Standard_Real Fxu,DFxu,FFxu,DFFxu;
|
|
Standard_Real Fyu,DFyu,FFyu,DFFyu;
|
|
Standard_Boolean Finish;
|
|
Standard_Real FFi;
|
|
Standard_Integer Nbiter = 30;
|
|
Standard_Integer Iter;
|
|
Standard_Real Ambda,T;
|
|
Standard_Real AA,BB,CC;
|
|
Standard_Integer Nn;
|
|
Standard_Real Alfa1=0,Alfa2;
|
|
Standard_Real OldDF = RealLast();
|
|
Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
|
|
Standard_Boolean Ok;
|
|
|
|
Done = Standard_False;
|
|
|
|
StdFail_NotDone_Raise_if(NbSample <= 0, " ");
|
|
|
|
// initialisation
|
|
|
|
if (A > B) {
|
|
Lowr=B;
|
|
Upp=A;
|
|
}
|
|
else {
|
|
Lowr=A;
|
|
Upp=B;
|
|
}
|
|
|
|
Increment = (Upp-Lowr)/NbSample;
|
|
StdFail_NotDone_Raise_if(Increment < EpsX, " ");
|
|
Done = Standard_True;
|
|
//-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
|
|
Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
|
|
Standard_Real NEpsX = 0.0000000001 * DeltaU;
|
|
if(EpsX < NEpsX) {
|
|
EpsX = NEpsX;
|
|
//-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
|
|
}
|
|
//--
|
|
Null2 = EpsNull*EpsNull;
|
|
|
|
Ok = F.Values(Lowr,FLowr,DFLowr);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
FLowr = FLowr - K;
|
|
|
|
Ok = F.Values(Upp,FUpp,DFUpp);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
FUpp = FUpp - K;
|
|
|
|
// Calcul sur U
|
|
|
|
U = Lowr-EpsX;
|
|
Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
|
|
DFyu = DFLowr;
|
|
FFyu = Fyu*Fyu;
|
|
DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
|
|
AllNull = ( FFyu <= Null2 );
|
|
|
|
while ( U < Upp) {
|
|
|
|
Xu = U;
|
|
Fxu = Fyu;
|
|
DFxu = DFyu;
|
|
FFxu = FFyu;
|
|
DFFxu = DFFyu;
|
|
|
|
U = Xu + Increment;
|
|
if (U <= Lowr) {
|
|
Fyu = FLowr + (U-Lowr)*DFLowr;
|
|
DFyu = DFLowr;
|
|
}
|
|
else if (U >= Upp) {
|
|
Fyu = FUpp + (U-Upp)*DFUpp;
|
|
DFyu = DFUpp;
|
|
}
|
|
else {
|
|
Ok = F.Values(U,Fyu,DFyu);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
Fyu = Fyu - K;
|
|
}
|
|
FFyu = Fyu*Fyu;
|
|
DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
|
|
|
|
if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
|
|
|
|
if (AllNull) { //rechercher les vraix zeros depuis le debut
|
|
|
|
AllNull = Standard_False;
|
|
Xu = Lowr-EpsX;
|
|
Fxu = FLowr-EpsX*DFLowr;
|
|
DFxu = DFLowr;
|
|
FFxu = Fxu*Fxu;
|
|
DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
|
|
U = Xu + Increment;
|
|
Ok = F.Values(U,Fyu,DFyu);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
Fyu = Fyu - K;
|
|
FFyu = Fyu*Fyu;
|
|
DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
|
|
}
|
|
Standard_Real FxuFyu=Fxu*Fyu;
|
|
|
|
if ( (DFFyu > 0. && DFFxu <= 0.)
|
|
|| (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
|
|
|| (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
|
|
|| (FxuFyu <= 0.) ) {
|
|
// recherche d 1 minimun possible
|
|
Finish = Standard_False;
|
|
Ambda = Increment;
|
|
T = 0.;
|
|
Iter=0;
|
|
FFi=0.;
|
|
|
|
if (FxuFyu >0.) {
|
|
// chercher si f peut s annuler pour eviter
|
|
// des iterations inutiles
|
|
if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
|
|
Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
|
|
|
|
Finish = Standard_True;
|
|
FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
|
|
}
|
|
else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
|
|
(FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
|
|
|
|
Finish = Standard_True;
|
|
FFxu = 0.0;
|
|
FFi = FFyu; // pour recalculer yu
|
|
}
|
|
else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
|
|
(FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
|
|
|
|
Finish = Standard_True;
|
|
FFyu =0.0;
|
|
FFi = FFxu; // pour recalculer U
|
|
}
|
|
}
|
|
else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
|
|
|
|
Finish = Standard_True;
|
|
FFxu = 0.0;
|
|
FFi = FFyu;
|
|
}
|
|
else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
|
|
|
|
Finish = Standard_True;
|
|
FFyu =0.0;
|
|
FFi = FFxu;
|
|
}
|
|
while (!Finish) {
|
|
|
|
// calcul des 2 solutions annulant la derivee de l interpolation cubique
|
|
// Ambda*t=(U-Xu) F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
|
|
// df=aa*t*t+2*bb*t+cc
|
|
|
|
AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
|
|
BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
|
|
CC = Ambda*DFFxu;
|
|
|
|
if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
|
|
AA=BB=CC=0;
|
|
}
|
|
|
|
|
|
math_DirectPolynomialRoots Solv(AA,BB,CC);
|
|
if ( !Solv.InfiniteRoots() ) {
|
|
Nn=Solv.NbSolutions();
|
|
if (Nn <= 0) {
|
|
if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
|
|
else Finish = Standard_True;
|
|
}
|
|
else {
|
|
Alfa1 = Solv.Value(1);
|
|
if (Nn == 2) {
|
|
Alfa2 = Solv.Value(2);
|
|
if (Alfa1 > Alfa2){
|
|
AA = Alfa1;
|
|
Alfa1 = Alfa2;
|
|
Alfa2 = AA;
|
|
}
|
|
if (Alfa1 > 1. || Alfa2 < 0.){
|
|
// resolution par dichotomie
|
|
if (Fxu*Fyu < 0.) Alfa1 = 0.5;
|
|
else Finish = Standard_True;
|
|
}
|
|
else if ( Alfa1 < 0. ||
|
|
( DFFxu > 0. && DFFyu >= 0.) ) {
|
|
// si 2 derivees >0
|
|
//(cas changement de signe de la distance signee sans
|
|
// changement de signe de la derivee:
|
|
//cas de 'presque'tangence avec 2
|
|
// solutions proches) ,on prend la plus grane racine
|
|
if (Alfa2 > 1.) {
|
|
if (Fxu*Fyu < 0.) Alfa1 = 0.5;
|
|
else Finish = Standard_True;
|
|
}
|
|
else Alfa1 = Alfa2;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
|
|
//-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
|
|
else Finish = Standard_True;
|
|
|
|
if (!Finish) {
|
|
// petits tests pour diminuer le nombre d iterations
|
|
if (Alfa1 <= EpsX ) {
|
|
Alfa1+=Alfa1;
|
|
}
|
|
else if (Alfa1 >= (1.-EpsX) ) {
|
|
Alfa1 = Alfa1+Alfa1-1.;
|
|
}
|
|
Alfa1 = Ambda*(1. - Alfa1);
|
|
U = U + T - Alfa1;
|
|
if ( U <= Lowr ) {
|
|
AA = FLowr + (U - Lowr)*DFLowr;
|
|
BB = DFLowr;
|
|
}
|
|
else if ( U >= Upp ) {
|
|
AA = FUpp + (U - Upp)*DFUpp;
|
|
BB = DFUpp;
|
|
}
|
|
else {
|
|
Ok = F.Values(U,AA,BB);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
AA = AA - K;
|
|
}
|
|
FFi = AA*AA;
|
|
CC = AA*BB; CC+=CC;
|
|
if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
|
|
&& AA*Fxu > 0. ) {
|
|
FFxu = FFi;
|
|
DFFxu = CC;
|
|
Fxu = AA;
|
|
DFxu = BB;
|
|
T = Alfa1;
|
|
if (Alfa1 > Ambda*0.5) {
|
|
// remarque (1)
|
|
// determination d 1 autre borne pour diviser
|
|
//le nouvel intervalle par 2 au -
|
|
Xu = U + Alfa1*0.5;
|
|
if ( Xu <= Lowr ) {
|
|
AA = FLowr + (Xu - Lowr)*DFLowr;
|
|
BB = DFLowr;
|
|
}
|
|
else if ( Xu >= Upp ) {
|
|
AA = FUpp + (Xu - Upp)*DFUpp;
|
|
BB = DFUpp;
|
|
}
|
|
else {
|
|
Ok = F.Values(Xu,AA,BB);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
AA = AA -K;
|
|
}
|
|
FFi = AA*AA;
|
|
CC = AA*BB; CC+=CC;
|
|
if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
|
|
|| Fxu * AA < 0. ) {
|
|
Fyu = AA;
|
|
DFyu = BB;
|
|
FFyu = FFi;
|
|
DFFyu = CC;
|
|
T = Alfa1*0.5;
|
|
Ambda = Alfa1*0.5;
|
|
}
|
|
else if ( AA*Fyu < 0. && AA*Fxu >0.) {
|
|
// changement de signe sur l intervalle u,U
|
|
Fxu = AA;
|
|
DFxu = BB;
|
|
FFxu = FFi;
|
|
DFFxu = CC;
|
|
FFi = Min(FFxu,FFyu);
|
|
T = Alfa1*0.5;
|
|
Ambda = Alfa1*0.5;
|
|
U = Xu;
|
|
}
|
|
else { Ambda = Alfa1;}
|
|
}
|
|
else { Ambda = Alfa1;}
|
|
}
|
|
else {
|
|
Fyu = AA;
|
|
DFyu = BB;
|
|
FFyu = FFi;
|
|
DFFyu = CC;
|
|
if ( (Ambda-Alfa1) > Ambda*0.5 ) {
|
|
// meme remarque (1)
|
|
Xu = U - (Ambda-Alfa1)*0.5;
|
|
if ( Xu <= Lowr ) {
|
|
AA = FLowr + (Xu - Lowr)*DFLowr;
|
|
BB = DFLowr;
|
|
}
|
|
else if ( Xu >= Upp ) {
|
|
AA = FUpp + (Xu - Upp)*DFUpp;
|
|
BB = DFUpp;
|
|
}
|
|
else {
|
|
Ok = F.Values(Xu,AA,BB);
|
|
|
|
if(!Ok) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
|
|
AA = AA - K;
|
|
}
|
|
FFi = AA*AA;
|
|
CC = AA*BB; CC+=CC;
|
|
if ( AA*Fyu <=0. && AA*Fxu > 0.) {
|
|
FFxu = FFi;
|
|
DFFxu = CC;
|
|
Fxu = AA;
|
|
DFxu = BB;
|
|
Ambda = ( Ambda - Alfa1 )*0.5;
|
|
T = 0.;
|
|
}
|
|
else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
|
|
FFyu = FFi;
|
|
DFFyu = CC;
|
|
Fyu = AA;
|
|
DFyu = BB;
|
|
Ambda = ( Ambda - Alfa1 )*0.5;
|
|
T = 0.;
|
|
FFi = Min(FFxu , FFyu);
|
|
U = Xu;
|
|
}
|
|
else {
|
|
T =0.;
|
|
Ambda = Ambda - Alfa1;
|
|
}
|
|
}
|
|
else {
|
|
T =0.;
|
|
Ambda = Ambda - Alfa1;
|
|
}
|
|
}
|
|
// tests d arrets
|
|
if (Abs(FFxu) <= Standard_Underflow ||
|
|
(Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
|
|
){
|
|
Finish = Standard_True;
|
|
if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
|
|
FFi = FFyu;
|
|
}
|
|
else if (Abs(FFyu) <= Standard_Underflow ||
|
|
(Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
|
|
){
|
|
Finish = Standard_True;
|
|
if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
|
|
FFi = FFxu;
|
|
}
|
|
else {
|
|
Iter = Iter + 1;
|
|
Finish = Iter >= Nbiter || (Ambda <= EpsX &&
|
|
( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
|
|
}
|
|
}
|
|
} // fin interpolation cubique
|
|
|
|
// restitution du meilleur resultat
|
|
|
|
if ( FFxu < FFi ) { U = U + T -Ambda;}
|
|
else if ( FFyu < FFi ) { U = U + T;}
|
|
|
|
if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
|
|
U = Max( Lowr , Min( U , Upp));
|
|
Ok = F.Value(U,FFi);
|
|
FFi = FFi - K;
|
|
if ( Abs(FFi) < EpsF ) {
|
|
//coherence
|
|
if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
|
|
else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
|
|
else if (Fxu*Fyu > 0.) { AA = 0.;}
|
|
else { AA = Fyu-Fxu;}
|
|
if (!Sol.IsEmpty()) {
|
|
if (Abs(Sol.Last() - U) > 5.*EpsX
|
|
|| (OldDF != RealLast() && OldDF*AA < 0.) ) {
|
|
Sol.Append(U);
|
|
NbStateSol.Append(F.GetStateNumber());
|
|
}
|
|
}
|
|
else {
|
|
Sol.Append(U);
|
|
NbStateSol.Append(F.GetStateNumber());
|
|
}
|
|
OldDF = AA ;
|
|
}
|
|
}
|
|
DFFyu = 0.;
|
|
Fyu = 0.;
|
|
Nn = 1;
|
|
while ( Nn < 1000000 && DFFyu <= 0.) {
|
|
// on repart d 1 pouyem plus loin
|
|
U = U + Nn*EpsX;
|
|
if ( U <= Lowr ) {
|
|
AA = FLowr + (U - Lowr)*DFLowr;
|
|
BB = DFLowr;
|
|
}
|
|
else if ( U >= Upp ) {
|
|
AA = FUpp + (U - Upp)*DFUpp;
|
|
BB = DFUpp;
|
|
}
|
|
else {
|
|
Ok = F.Values(U,AA,BB);
|
|
AA = AA - K;
|
|
}
|
|
if ( AA*Fyu < 0.) {
|
|
U = U - Nn*EpsX;
|
|
Nn = 1000001;
|
|
}
|
|
else {
|
|
Fyu = AA;
|
|
DFyu = BB;
|
|
FFyu = AA*AA;
|
|
DFFyu = 2.*DFyu*Fyu;
|
|
Nn = 2*Nn;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#if NEWSEQ
|
|
#ifdef MATH_FUNCTIONROOTS_CHECK
|
|
{
|
|
Standard_Integer n1=StaticSol.Length();
|
|
Standard_Integer n2=Sol.Length();
|
|
if(n1!=n2) {
|
|
printf("\n mathFunctionRoots : n1=%d n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
|
|
for(Standard_Integer x1=1; x1<=n1;x1++) {
|
|
Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
|
|
printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
|
|
}
|
|
printf("\n");
|
|
for(Standard_Integer x2=1; x2<=n2; x2++) {
|
|
Standard_Real v; F.Value(Sol(x2),v); v-=K;
|
|
printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
|
|
}
|
|
printf("\n");
|
|
}
|
|
Standard_Integer n=n1;
|
|
if(n1>n2) n=n2;
|
|
for(Standard_Integer i=1;i<=n;i++) {
|
|
Standard_Real t = Sol(i)-StaticSol(i);
|
|
if(Abs(t)>NEpsX) {
|
|
printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
#endif
|
|
}
|
|
#endif
|
|
}
|
|
|
|
|
|
void math_FunctionRoots::Dump(Standard_OStream& o) const
|
|
{
|
|
|
|
o << "math_FunctionRoots ";
|
|
if(Done) {
|
|
o << " Status = Done \n";
|
|
o << " Number of solutions = " << Sol.Length() << endl;
|
|
for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
|
|
o << " Solution Number " << i << "= " << Sol.Value(i) << endl;
|
|
}
|
|
}
|
|
else {
|
|
o<< " Status = not Done \n";
|
|
}
|
|
}
|