src/math/math_FunctionAllRoots.cxx

   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
  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
  33 math_FunctionAllRoots::math_FunctionAllRoots (
  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 }