asin.h

Go to the documentation of this file.

/*
 * aasin.h
 * The basic idea is to exploit Pade' polynomials.
 * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
 * moshier@na-net.ornl.gov) as well as actual code.
 * The Cephes library can be found here:  http://www.netlib.org/cephes/
 *
 *  Created on: Jun 23, 2012
 *      Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
 */

/*
 * VDT is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser Public License for more details.
 *
 * You should have received a copy of the GNU Lesser Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef ASIN_H_
#define ASIN_H_

#include "vdtcore_common.h"

namespace vdt{

namespace details{

const double RX1asin = 2.967721961301243206100E-3;
const double RX2asin = -5.634242780008963776856E-1;
const double RX3asin = 6.968710824104713396794E0;
const double RX4asin = -2.556901049652824852289E1;
const double RX5asin = 2.853665548261061424989E1;

const double SX1asin = -2.194779531642920639778E1;
const double SX2asin =  1.470656354026814941758E2;
const double SX3asin = -3.838770957603691357202E2;
const double SX4asin = 3.424398657913078477438E2;

const double PX1asin = 4.253011369004428248960E-3;
const double PX2asin = -6.019598008014123785661E-1;
const double PX3asin = 5.444622390564711410273E0;
const double PX4asin = -1.626247967210700244449E1;
const double PX5asin = 1.956261983317594739197E1;
const double PX6asin = -8.198089802484824371615E0;

const double QX1asin = -1.474091372988853791896E1;
const double QX2asin =  7.049610280856842141659E1;
const double QX3asin = -1.471791292232726029859E2;
const double QX4asin = 1.395105614657485689735E2;
const double QX5asin = -4.918853881490881290097E1;

inline double getRX(const double x){
   double rx = RX1asin;
   rx*= x;
   rx+= RX2asin;
   rx*= x;
   rx+= RX3asin;
   rx*= x;
   rx+= RX4asin;
   rx*= x;
   rx+= RX5asin;
   return rx;
}
inline double getSX(const double x){
   double sx = x;
   sx+= SX1asin;
   sx*= x;
   sx+= SX2asin;
   sx*= x;
   sx+= SX3asin;
   sx*= x;
   sx+= SX4asin;
   return sx;
}

inline double getPX(const double x){
   double px = PX1asin;
   px*= x;
   px+= PX2asin;
   px*= x;
   px+= PX3asin;
   px*= x;
   px+= PX4asin;
   px*= x;
   px+= PX5asin;
   px*= x;
   px+= PX6asin;
   return px;
}

inline double getQX(const double x){
   double qx = x;
   qx+= QX1asin;
   qx*= x;
   qx+= QX2asin;
   qx*= x;
   qx+= QX3asin;
   qx*= x;
   qx+= QX4asin;
   qx*= x;
   qx+= QX5asin;
   return qx;
   }
}

}

namespace vdt{

// asin double precision --------------------------------------------------------
/// Double Precision asin
inline double fast_asin(double x){

   const uint64_t sign_mask = details::getSignMask(x);
   x = std::fabs(x);
   const double a = x;


   double zz = 1.0 - a;
   double px = details::getRX(zz);
   double qx = details::getSX(zz);

   const double p = zz * px/qx;

   zz = std::sqrt(zz+zz);
   double z = details::PIO4 - zz;
   zz = zz * p - details::MOREBITS;
   z -= zz;
   z += details::PIO4;

   if( a < 0.625 ){
      zz = a * a;
      px = details::getPX(zz);
      qx = details::getQX(zz);
      z = zz*px/qx;
      z = a * z + a;
   }


   // Linear approx, not sooo needed but seable. Price is cheap though
   double res = a < 1e-8? a : z ;
        // Restore Sign
   return details::dpORuint64(res,sign_mask);

}

//------------------------------------------------------------------------------
/// Single Precision asin
inline float fast_asinf(float x){


    uint32_t flag=0;

    const uint32_t sign_mask = details::getSignMask(x);
    const float a = std::fabs(x);

    float z;
    if( a > 0.5f )
        {
        z = 0.5f * (1.0f - a);
        x = sqrtf( z );
        flag = 1;
        }
    else
        {
        x = a;
        z = x * x;
        }

    z = (((( 4.2163199048E-2f * z
            + 2.4181311049E-2f) * z
            + 4.5470025998E-2f) * z
            + 7.4953002686E-2f) * z
            + 1.6666752422E-1f) * z * x
            + x;

//     if( flag != 0 )
//       {
//       z = z + z;
//       z = PIO2F - z;
//       }

    // No branch with the two coefficients

    float tmp = z + z;
    tmp = details::PIO2F - tmp;

    // Linear approx, not sooo needed but seable. Price is cheap though
    float res = a < 1e-4f? a : tmp * flag + (1-flag) * z ;

    // Restore Sign
    return details::spORuint32(res,sign_mask);

}

//------------------------------------------------------------------------------
// The cos is in this file as well

inline double fast_acos( double x ){return details::PIO2  - fast_asin(x);}

//------------------------------------------------------------------------------

inline float fast_acosf( float x ){return details::PIO2F  - fast_asinf(x);}

//------------------------------------------------------------------------------

// // Vector signatures
//
// void asinv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void fast_asinv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void asinfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
// void fast_asinfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
//
// void acosv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void fast_acosv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void acosfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
// void fast_acosfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);

} //vdt namespace

#endif /* ASIN_H_ */