Program Listing for File PolynomialRoots-Utils.hh

↰ Return to documentation for file (PolynomialRoots-Utils.hh)

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2014                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/

#ifndef RPOLY_HH
#define RPOLY_HH

#include <utility>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <iostream>
#include <limits>

/*
..
.. N. FLOCKE
.. Algorithm 954: An Accurate and Efficient Cubic and Quartic Equation Solver
.. for Physical Applications
.. ACM TOMS, Vol. 41, No. 4, 2015.
.. DOI: http://dx.doi.org/10.1145/2699468
..
*/

#include <cstdint>

namespace PolynomialRoots {

  typedef double valueType;
  typedef int    indexType;
  typedef std::complex<valueType> complexType;

  #ifndef DOXYGEN_SHOULD_SKIP_THIS

  static int       const bitsValueType = std::numeric_limits<valueType>::digits;
  static valueType const splitFactor   = static_cast<valueType>((std::uint64_t(1)<<(bitsValueType-2))+1); // one extra digit is implicitly 1

  /*
  ||         _   _ _
  ||   _   _| |_(_) |___
  ||  | | | | __| | / __|
  ||  | |_| | |_| | \__ \
  ||   \__,_|\__|_|_|___/
  */
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // a + b = x + err
  static
  inline
  void
  TwoSum(
    valueType   a,
    valueType   b,
    valueType & x,
    valueType & err
  ) {
    x = a+b;
    valueType z = x-a;
    err = (a-(x-z))+(b-z);
  }

  static
  inline
  void
  TwoSum(
    complexType   a,
    complexType   b,
    complexType & x,
    complexType & err
  ) {
    valueType s1, e1, s2, e2;
    TwoSum( a.real(), b.real(), s1, e1 );
    TwoSum( a.imag(), b.imag(), s2, e2 );
    x   = complexType(s1,s2);
    err = complexType(e1,e2);
  }

  // a = x + y
  static
  inline
  void
  Split( valueType a, valueType & x, valueType & y ) {
    valueType c = splitFactor*a;
    x = c-(c-a);
    y = a-x;
  }

  // a * b = x + err
  static
  inline
  void
  TwoProduct(
    valueType   a,
    valueType   b,
    valueType & x,
    valueType & err
  ) {
    valueType a1, a2, b1, b2;
    Split( a, a1, a2 );
    Split( b, b1, b2 );
    x   = a*b;
    err = a2*b2-(((x-a1*b1)-a2*b1)-a1*b2);
  }

  static
  inline
  void
  TwoProduct(
    complexType   a,
    complexType   b,
    complexType & p,
    complexType & e,
    complexType & f,
    complexType & g
  ) {
    valueType z1, z2, z3, z4, z5, z6, h1, h2, h3, h4, h5, h6;
    TwoProduct(a.real(), b.real(), z1, h1 );
    TwoProduct(a.imag(), b.imag(), z2, h2 );
    TwoProduct(a.real(), b.imag(), z3, h3 );
    TwoProduct(a.imag(), b.real(), z4, h4 );
    TwoSum(z1, -z2, z5, h5);
    TwoSum(z3, z4, z6, h6);
    p = complexType(z5,z6);
    e = complexType(h1,h3);
    f = complexType(-h2,h4);
    g = complexType(h5,h6);
  }

  #endif

}

#endif