mirror of
https://github.com/Open-Cascade-SAS/OCCT.git
synced 2026-06-02 18:06:43 +08:00
3ed30348aa
Resolved some C4702 (unreachable code) and C4127 (conditional expression is constant). small corrections in NoSuchObject invoking. Macros names changing, deadcode deleting, re-writing "for" loops into equivalent "if" structures. changed condition in "if" block, deadcode deleted. Small changes in else statement.
1035 lines
26 KiB
C++
Executable File
1035 lines
26 KiB
C++
Executable File
// Copyright (c) 1997-1999 Matra Datavision
|
|
// Copyright (c) 1999-2012 OPEN CASCADE SAS
|
|
//
|
|
// The content of this file is subject to the Open CASCADE Technology Public
|
|
// License Version 6.5 (the "License"). You may not use the content of this file
|
|
// except in compliance with the License. Please obtain a copy of the License
|
|
// at http://www.opencascade.org and read it completely before using this file.
|
|
//
|
|
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
|
|
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
|
|
//
|
|
// The Original Code and all software distributed under the License is
|
|
// distributed on an "AS IS" basis, without warranty of any kind, and the
|
|
// Initial Developer hereby disclaims all such warranties, including without
|
|
// limitation, any warranties of merchantability, fitness for a particular
|
|
// purpose or non-infringement. Please see the License for the specific terms
|
|
// and conditions governing the rights and limitations under the License.
|
|
|
|
//#ifndef DEB
|
|
#define No_Standard_RangeError
|
|
#define No_Standard_OutOfRange
|
|
#define No_Standard_DimensionError
|
|
//#endif
|
|
|
|
#include <StdFail_NotDone.hxx>
|
|
#include <Standard_RangeError.hxx>
|
|
#include <math_DirectPolynomialRoots.hxx>
|
|
#include <math_FunctionRoots.ixx>
|
|
#include <math_FunctionWithDerivative.hxx>
|
|
#include <TColStd_Array1OfReal.hxx>
|
|
//#ifdef WNT
|
|
#include <stdio.h>
|
|
//#endif
|
|
|
|
#define ITMAX 100
|
|
#define EPS 1e-14
|
|
#define EPSEPS 2e-14
|
|
#define MAXBIS 100
|
|
|
|
# ifdef DEB
|
|
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 DEB
|
|
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 DEB
|
|
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 DEB
|
|
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 DEB
|
|
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");
|
|
while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > 0)) {
|
|
//-- 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*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";
|
|
}
|
|
}
|