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