[occt.git] / src / math / math_FunctionAllRoots.cxx

// 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 <math_FunctionAllRoots.ixx>

#include <Standard_NumericError.hxx>
#include <Standard_OutOfRange.hxx>
#include <math_FunctionRoots.hxx>
#include <math_FunctionWithDerivative.hxx>
#include <math_FunctionSample.hxx>

math_FunctionAllRoots::math_FunctionAllRoots (
                           math_FunctionWithDerivative& F,
			   const math_FunctionSample& S,
			   const Standard_Real EpsX, const Standard_Real EpsF,
			   const Standard_Real EpsNul) {

    done=Standard_False;


    Standard_Boolean Nul, PNul, InterNul, Nuld, Nulf;
    Standard_Real DebNul,FinNul;
    Standard_Integer Indd,Indf;
    Standard_Real cst,val,valsav=0,valbid;
    Standard_Boolean bid,fini;
    Standard_Integer Nbp,i;

    Nbp=S.NbPoints();
    bid=F.Value(S.GetParameter(1),val);
    PNul=Abs(val)<=EpsNul;
    if (!PNul) {valsav=val;}
    InterNul=Standard_False;
    Nuld=Standard_False;
    Nulf=Standard_False;

    i=2;
    fini=(i>Nbp);

    while (!fini) {

      bid=F.Value(S.GetParameter(i),val);
      Nul=Abs(val)<=EpsNul;
      if (!Nul) {
	valsav=val;
      }
      if (InterNul && (!Nul)) {
	InterNul=Standard_False;
	pdeb.Append(DebNul);
	ideb.Append(Indd);
	if (val>0.0) {
	  cst=EpsNul;
	}
	else {
	  cst=-EpsNul;
	}
	math_FunctionRoots Res1(F,S.GetParameter(i-1),S.GetParameter(i),10,
				EpsX,EpsF,0.0,cst);
	Standard_NumericError_Raise_if((!Res1.IsDone())   ||
				  (Res1.IsAllNull()) ||
				  (Res1.NbSolutions()==0), " ");

	FinNul=Res1.Value(1);
	Indf=Res1.StateNumber(1);

	cst=-cst;
	math_FunctionRoots Res2(F,S.GetParameter(i-1),S.GetParameter(i),10,
				EpsX,EpsF,0.0,cst);
	Standard_NumericError_Raise_if((!Res2.IsDone())   ||
				       (Res2.IsAllNull())," ");
	  
				   //-- || (Res2.NbSolutions()!=0), " ");  lbr le 13 mai 87 (!=0 -> ==0)
	if (Res2.NbSolutions()!=0) {
	  if (Res2.Value(1) < FinNul) {
	    FinNul = Res2.Value(1);
	    Indf = Res2.StateNumber(1);
	  }
	}
	pfin.Append(FinNul);
	ifin.Append(Indf);
      }
      else if ((!InterNul) && PNul && Nul) {
	InterNul=Standard_True;
	if (i==2) {
	  DebNul=S.GetParameter(1);
	  bid = F.Value(DebNul,valbid);
	  Indd = F.GetStateNumber();
	  Nuld=Standard_True;
	}
	else {
	  if (valsav>0.0) {
	    cst=EpsNul;
	  }
	  else {
	    cst=-EpsNul;
	  }
	  math_FunctionRoots Res1(F,S.GetParameter(i-2),S.GetParameter(i-1),10,
				  EpsX,EpsF,0.0,cst);
	  Standard_NumericError_Raise_if((!Res1.IsDone())   ||
				  (Res1.IsAllNull()) ||
				  (Res1.NbSolutions()==0), " ");
	  DebNul=Res1.Value(Res1.NbSolutions());
	  Indd = Res1.StateNumber(Res1.NbSolutions());

	  cst=-cst;
	  math_FunctionRoots Res3(F,S.GetParameter(i-2),S.GetParameter(i-1),10,
				  EpsX,EpsF,0.0,cst);
	  Standard_NumericError_Raise_if((!Res3.IsDone()) ||
				    (Res3.IsAllNull()), " ");

	  if (Res3.NbSolutions()!=0) {
	    if (Res3.Value(Res3.NbSolutions()) > DebNul) {
	      DebNul = Res3.Value(Res3.NbSolutions());
	      Indd = Res3.StateNumber(Res3.NbSolutions());
	    }
	  }
	}
      }
      i=i+1;
      PNul=Nul;
      fini=(i>Nbp);
    }

    if (InterNul) {            // rajouter l intervalle finissant au dernier pt
      pdeb.Append(DebNul);
      ideb.Append(Indd);
      FinNul = S.GetParameter(Nbp);
      bid = F.Value(FinNul,valbid);
      Indf = F.GetStateNumber();
      pfin.Append(FinNul);
      ifin.Append(Indf);
      Nulf=Standard_True;
    }

    if (pdeb.Length()==0) {  // Pas d intervalle nul

      math_FunctionRoots Res(F,S.GetParameter(1),S.GetParameter(Nbp),Nbp,
			     EpsX,EpsF,0.0);
      Standard_NumericError_Raise_if((!Res.IsDone()) ||
				(Res.IsAllNull()), " ");

      for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
	piso.Append(Res.Value(j));
	iiso.Append(Res.StateNumber(j));
      }
    }
    else {
      Standard_Integer NbpMin = 3;
      Standard_Integer Nbrpt;
      if (!Nuld) {          // Recherche des solutions entre S.GetParameter(1) 
                            // et le debut du 1er intervalle nul

	Nbrpt=(Standard_Integer ) IntegerPart(
                              Abs((pdeb.Value(1)-S.GetParameter(1))/
			      (S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
        math_FunctionRoots Res(F,S.GetParameter(1),pdeb.Value(1),
			       Max(Nbrpt, NbpMin), EpsX,EpsF,0.0);
	Standard_NumericError_Raise_if((!Res.IsDone()) ||
				  (Res.IsAllNull()), " ");

	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
	  piso.Append(Res.Value(j));
	  iiso.Append(Res.StateNumber(j));
	}
      }
      for (Standard_Integer k=2; k<=pdeb.Length(); k++) {
	Nbrpt=(Standard_Integer )
                              IntegerPart(Abs((pdeb.Value(k)-pfin.Value(k-1))/
			      (S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
        math_FunctionRoots Res(F,pfin.Value(k-1),pdeb.Value(k),
			       Max(Nbrpt, NbpMin), EpsX,EpsF,0.0);
	Standard_NumericError_Raise_if((!Res.IsDone()) ||
				  (Res.IsAllNull()), " ");

	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
	  piso.Append(Res.Value(j));
	  iiso.Append(Res.StateNumber(j));
	}
      }
      if (!Nulf) {            // Recherche des solutions entre la fin du
                              // dernier intervalle nul et Value(Nbp).

	Nbrpt=(Standard_Integer )
                              IntegerPart(Abs((S.GetParameter(Nbp)-
			       pfin.Value(pdeb.Length()))/
			      (S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
        math_FunctionRoots Res(F,pfin.Value(pdeb.Length()),
			       S.GetParameter(Nbp), Max(Nbrpt, NbpMin),
			       EpsX,EpsF,0.0);
	Standard_NumericError_Raise_if((!Res.IsDone()) ||
				  (Res.IsAllNull()), " ");

	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
	  piso.Append(Res.Value(j));
	  iiso.Append(Res.StateNumber(j));
	}
      }
    }
    done=Standard_True;
  }


void math_FunctionAllRoots::Dump(Standard_OStream& o) const {
  
  o<< "math_FunctionAllRoots ";
  if(done) {
    o<< " Status = Done \n";
    o << " Number of null intervals = " << pdeb.Length() << endl;
    o << " Number of points where the function is null: " << piso.Length() << endl;
  }
  else {
    o<< " Status = not Done \n";
  }
}
Commit	Line	Data
b311480e	1	// Copyright (c) 1997-1999 Matra Datavision
	2	// Copyright (c) 1999-2012 OPEN CASCADE SAS
	3	//
	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.
	8	//
	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.
	11	//
	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.
	18
7fd59977	19	//#ifndef DEB
	20	#define No_Standard_RangeError
	21	#define No_Standard_OutOfRange
	22	#define No_Standard_DimensionError
	23	//#endif
	24
	25	#include <math_FunctionAllRoots.ixx>
	26
	27	#include <Standard_NumericError.hxx>
	28	#include <Standard_OutOfRange.hxx>
	29	#include <math_FunctionRoots.hxx>
	30	#include <math_FunctionWithDerivative.hxx>
	31	#include <math_FunctionSample.hxx>
	32
b311480e	33	math_FunctionAllRoots::math_FunctionAllRoots (
7fd59977	34	math_FunctionWithDerivative& F,
	35	const math_FunctionSample& S,
	36	const Standard_Real EpsX, const Standard_Real EpsF,
	37	const Standard_Real EpsNul) {
	38
	39	done=Standard_False;
	40
	41
	42
	43	Standard_Boolean Nul, PNul, InterNul, Nuld, Nulf;
	44	Standard_Real DebNul,FinNul;
	45	Standard_Integer Indd,Indf;
	46	Standard_Real cst,val,valsav=0,valbid;
	47	Standard_Boolean bid,fini;
	48	Standard_Integer Nbp,i;
	49
	50	Nbp=S.NbPoints();
	51	bid=F.Value(S.GetParameter(1),val);
	52	PNul=Abs(val)<=EpsNul;
	53	if (!PNul) {valsav=val;}
	54	InterNul=Standard_False;
	55	Nuld=Standard_False;
	56	Nulf=Standard_False;
	57
	58	i=2;
	59	fini=(i>Nbp);
	60
	61	while (!fini) {
	62
	63	bid=F.Value(S.GetParameter(i),val);
	64	Nul=Abs(val)<=EpsNul;
	65	if (!Nul) {
	66	valsav=val;
	67	}
	68	if (InterNul && (!Nul)) {
	69	InterNul=Standard_False;
	70	pdeb.Append(DebNul);
	71	ideb.Append(Indd);
	72	if (val>0.0) {
	73	cst=EpsNul;
	74	}
	75	else {
	76	cst=-EpsNul;
	77	}
	78	math_FunctionRoots Res1(F,S.GetParameter(i-1),S.GetParameter(i),10,
	79	EpsX,EpsF,0.0,cst);
	80	Standard_NumericError_Raise_if((!Res1.IsDone()) \|\|
	81	(Res1.IsAllNull()) \|\|
	82	(Res1.NbSolutions()==0), " ");
	83
	84	FinNul=Res1.Value(1);
	85	Indf=Res1.StateNumber(1);
	86
	87	cst=-cst;
	88	math_FunctionRoots Res2(F,S.GetParameter(i-1),S.GetParameter(i),10,
	89	EpsX,EpsF,0.0,cst);
	90	Standard_NumericError_Raise_if((!Res2.IsDone()) \|\|
	91	(Res2.IsAllNull())," ");
	92
	93	//-- \|\| (Res2.NbSolutions()!=0), " "); lbr le 13 mai 87 (!=0 -> ==0)
	94	if (Res2.NbSolutions()!=0) {
	95	if (Res2.Value(1) < FinNul) {
	96	FinNul = Res2.Value(1);
	97	Indf = Res2.StateNumber(1);
98	}
99	}
100	pfin.Append(FinNul);
101	ifin.Append(Indf);
102	}
103	else if ((!InterNul) && PNul && Nul) {
104	InterNul=Standard_True;
105	if (i==2) {
106	DebNul=S.GetParameter(1);
107	bid = F.Value(DebNul,valbid);
108	Indd = F.GetStateNumber();
109	Nuld=Standard_True;
110	}
111	else {
112	if (valsav>0.0) {
113	cst=EpsNul;
114	}
115	else {
116	cst=-EpsNul;
117	}
118	math_FunctionRoots Res1(F,S.GetParameter(i-2),S.GetParameter(i-1),10,
119	EpsX,EpsF,0.0,cst);
120	Standard_NumericError_Raise_if((!Res1.IsDone()) \|\|
121	(Res1.IsAllNull()) \|\|
122	(Res1.NbSolutions()==0), " ");
123	DebNul=Res1.Value(Res1.NbSolutions());
124	Indd = Res1.StateNumber(Res1.NbSolutions());
125
126	cst=-cst;
127	math_FunctionRoots Res3(F,S.GetParameter(i-2),S.GetParameter(i-1),10,
128	EpsX,EpsF,0.0,cst);
129	Standard_NumericError_Raise_if((!Res3.IsDone()) \|\|
130	(Res3.IsAllNull()), " ");
131
132	if (Res3.NbSolutions()!=0) {
133	if (Res3.Value(Res3.NbSolutions()) > DebNul) {
134	DebNul = Res3.Value(Res3.NbSolutions());
135	Indd = Res3.StateNumber(Res3.NbSolutions());
136	}
137	}
138	}
139	}
140	i=i+1;
141	PNul=Nul;
142	fini=(i>Nbp);
143	}
144
145	if (InterNul) { // rajouter l intervalle finissant au dernier pt
146	pdeb.Append(DebNul);
147	ideb.Append(Indd);
148	FinNul = S.GetParameter(Nbp);
149	bid = F.Value(FinNul,valbid);
150	Indf = F.GetStateNumber();
151	pfin.Append(FinNul);
152	ifin.Append(Indf);
153	Nulf=Standard_True;
154	}
155
156	if (pdeb.Length()==0) { // Pas d intervalle nul
157
158	math_FunctionRoots Res(F,S.GetParameter(1),S.GetParameter(Nbp),Nbp,
159	EpsX,EpsF,0.0);
160	Standard_NumericError_Raise_if((!Res.IsDone()) \|\|
161	(Res.IsAllNull()), " ");
162
163	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
164	piso.Append(Res.Value(j));
165	iiso.Append(Res.StateNumber(j));
166	}
167	}
168	else {
169	Standard_Integer NbpMin = 3;
170	Standard_Integer Nbrpt;
171	if (!Nuld) { // Recherche des solutions entre S.GetParameter(1)
172	// et le debut du 1er intervalle nul
173
174	Nbrpt=(Standard_Integer ) IntegerPart(
175	Abs((pdeb.Value(1)-S.GetParameter(1))/
176	(S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
177	math_FunctionRoots Res(F,S.GetParameter(1),pdeb.Value(1),
178	Max(Nbrpt, NbpMin), EpsX,EpsF,0.0);
179	Standard_NumericError_Raise_if((!Res.IsDone()) \|\|
180	(Res.IsAllNull()), " ");
181
182	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
183	piso.Append(Res.Value(j));
184	iiso.Append(Res.StateNumber(j));
185	}
186	}
187	for (Standard_Integer k=2; k<=pdeb.Length(); k++) {
188	Nbrpt=(Standard_Integer )
189	IntegerPart(Abs((pdeb.Value(k)-pfin.Value(k-1))/
190	(S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
191	math_FunctionRoots Res(F,pfin.Value(k-1),pdeb.Value(k),
192	Max(Nbrpt, NbpMin), EpsX,EpsF,0.0);
193	Standard_NumericError_Raise_if((!Res.IsDone()) \|\|
194	(Res.IsAllNull()), " ");
195
196	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
197	piso.Append(Res.Value(j));
198	iiso.Append(Res.StateNumber(j));
199	}
200	}
201	if (!Nulf) { // Recherche des solutions entre la fin du
202	// dernier intervalle nul et Value(Nbp).
203
204	Nbrpt=(Standard_Integer )
205	IntegerPart(Abs((S.GetParameter(Nbp)-
206	pfin.Value(pdeb.Length()))/
207	(S.GetParameter(Nbp)-S.GetParameter(1)))*Nbp);
208	math_FunctionRoots Res(F,pfin.Value(pdeb.Length()),
209	S.GetParameter(Nbp), Max(Nbrpt, NbpMin),
210	EpsX,EpsF,0.0);
211	Standard_NumericError_Raise_if((!Res.IsDone()) \|\|
212	(Res.IsAllNull()), " ");
213
214	for (Standard_Integer j=1; j<=Res.NbSolutions(); j++) {
215	piso.Append(Res.Value(j));
216	iiso.Append(Res.StateNumber(j));
217	}
218	}
219	}
220	done=Standard_True;
221	}
222
223
224	void math_FunctionAllRoots::Dump(Standard_OStream& o) const {
225
226	o<< "math_FunctionAllRoots ";
227	if(done) {
228	o<< " Status = Done \n";
229	o << " Number of null intervals = " << pdeb.Length() << endl;
230	o << " Number of points where the function is null: " << piso.Length() << endl;
231	}
232	else {
233	o<< " Status = not Done \n";
234	}
235	}