aimstil  5.0.5
poly_solver.tpp
Go to the documentation of this file.
1 
2 namespace til { namespace math
3 {
4 
5  //---------------------------------------------------------------------------
6 
7  template < typename TPrec, typename TInfinitySolutionsPolicy >
8  void
9  PolySolver_real< TPrec, TInfinitySolutionsPolicy >::
10  solve(TPrec a, TPrec b)
11  {
12  const TPrec EPSILON = 128*std::numeric_limits<TPrec>::epsilon();
13 
14  if (std::abs(a) < EPSILON)
15  {
16  if (std::abs(b) < EPSILON)
17  {
18  m_infSolPolicy(*this);
19  }
20  else
21  {
22  m_nSols = 0;
23  }
24  }
25  else
26  {
27  m_nSols = 1;
28  m_sols[0] = -b / a;
29  }
30  }
31 
32  //---------------------------------------------------------------------------
33 
34  template < typename TPrec, typename TInfinitySolutionsPolicy >
35  void
36  PolySolver_real< TPrec, TInfinitySolutionsPolicy >::
37  solve(TPrec a, TPrec b, TPrec c)
38  {
39  const TPrec EPSILON = 128*std::numeric_limits<TPrec>::epsilon();
40 
41  // Check that a is not too small -- otherwise we'll have numerical problems, that
42  // are solved here by considering that a small a makes the polynomial first degree
43  if (std::abs(a) < EPSILON)
44  {
45  this->solve(b,c);
46  return;
47  }
48 
49  TPrec delta = b*b-4*a*c;
50 
51  // Start with this case, as it is probably the most frequent case
52  if (delta > EPSILON)
53  {
54  m_nSols = 2;
55  delta = std::sqrt(delta);
56  a *= 2;
57  m_sols[0] = (-b-delta)/a;
58  m_sols[1] = (-b+delta)/a;
59  }
60  // TODO: I am not sure whether this should not be delta < 0 instead...
61  else if (delta < -EPSILON)
62  {
63  m_nSols = 0;
64  }
65  else
66  {
67  m_nSols = 1;
68  m_sols[0] = -b / (2*a);
69  }
70  }
71 
72  //---------------------------------------------------------------------------
73 
74 }} // namespace til::math
75