src/math/math_FunctionAllRoots.cxx

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