quadratic equation solver

In this exercises you will experiment catastrofic cancellations and verify the accuracy of the solution of a quadratic equation

1. here are two of the possible implementations of the solution of the quadratic equation

    #include<cmath>
    #include<tuple>
    #include<limits>
      
    template<typename T>
    inline T det(T a, T b, T c) {
      // compute determinant for equation ax^2 + 2bx + c= 0
      return std::sqrt(b*b-a*c);
    }
   
    template<typename T>
    inline std::tuple<T,T> quadSolverNaive(T a, T b, T c) {
      // solve equation ax^2 + 2bx + c= 0
      // using naive solution (as at college)
      auto d = -T(1)/a;
      return std::make_tuple(d*(b-det(a,b,c)),d*(b+det(a,b,c)));  
    }
   
   
    template<typename T>
    inline std::tuple<T,T> quadSolverOpt(T a, T b, T c) {
      // solve equation ax^2 + 2bx + c= 0
      // using stable algorithm
      auto q = -(std::copysign(det(a,b,c),b)+b);
      return std::make_tuple(q/a,c/q);
    }
   

2. compute radius and center for a circle given chord and sagita (simplify orientation as in the slides so y0=0)
3. use standard circle equation to recompute the value of xt. See for which value of the sagita the equation breaks-down

4. use above solutions to recompute the value of xt using the “local” circle equation. See for which value of the sagita the equation breaks-down

    #include<iostream>
    #include<iomanip>
   
    template<typename T> 
    void print(T x) {
      std::cout<< std::hexfloat << x <<' ' <<  std::scientific
      << std::setprecision(std::numeric_limits<T>::digits10+3) << x << std::endl;
    }
   
    template<typename T>
    void circle() {
      std::cout <<' '<< std::endl;
   
      constexpr T micron = 1.e-3;
      constexpr T meter = 1000.;
      constexpr T halfChord = meter/2;
    
      
      T xt = std::sqrt(T(77.));
      std::cout << "xt "; print(xt);
    
      auto sagita = T(10.);  // make it smaller and smaller...
      auto yt=halfChord;
      // compute xm, radius and x0
     
      // recompute xt
      // y0 = 0
         
      // standard circle eq. (x-x0)^2 +(yt-y0)^2 = r^2 ->  x = x0 - sqrt(r^2-yt^2)
      // print(xt-x);  // compare solution with actual value to get the accuracy
   
      // circle w/r/t local point (x-xm)^2 + (y-ym)^2  -2ar(x-xm) -2br(y-ym) = 0 ->
      // (x-xm)^2  - 2(x-xm)r + yt^2  = 0  -> x = xm + ....
         
      // print(xt-(xm+std::get<1>(solution)));  // compare solution with actual value to get the accuracy
    }
   
    int main() {
         
      std::cout << "\nfloat\n" << std::endl;
      circle<float>();
      std::cout << "\ndouble\n" << std::endl;
      circle<double>();
         
      return 0;
    }