//main.cpp
#pragma hdrstop
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include "UnitSine.h"
//From htpp://www.richelbilderbeek.nl/CppGetRandomUniform.htm
const double GetRandomUniform()
{
return static_cast<double>(std::rand())/static_cast<double>(RAND_MAX);
}
//Range from [ - 3.0 * M_PI , 3.0 * M_PI ]
const double GetRandomNumber()
{
return (-3.0 * M_PI) + (GetRandomUniform() * 6.0 * M_PI);
}
int main()
{
Sine s;
//Check special values
s( 0.0 * M_PI);
s( 0.5 * M_PI);
s( 1.0 * M_PI);
s( 1.5 * M_PI);
s( 2.0 * M_PI);
s(-0.0 * M_PI);
s(-0.5 * M_PI);
s(-1.0 * M_PI);
s(-1.5 * M_PI);
s(-2.0 * M_PI);
for (int i=0; i!=1000000; ++i)
{
if (i % 100000 == 0) std::cout << ".";
const double x = GetRandomNumber();
const double y = s(x);
const double yToo = std::sin(x);
const double error = std::fabs(y - yToo);
assert( error < s.GetMaxError());
}
std::cout << "Sine tested as safe" << std::endl;
}
UnitSine.h
#ifndef UnitSineH
#define UnitSineH
#include <vector>
struct Sine
{
Sine(const int size = 1000, const double notInitializedValue = 123.0);
double operator()(const double x);
const double GetMaxError() const;
private:
std::vector<double> mLut; //Look-up table
const double mNotInitializedValue;
const double Transform(const double f) const;
const int CalculateIndex(const double f) const;
};
#endif
UnitSine.cpp
#pragma hdrstop
#include <cassert>
#include <cmath>
#include "UnitSine.h"
Sine::Sine(const int size, const double notInitializedValue)
: mLut(size,notInitializedValue),
mNotInitializedValue(notInitializedValue)
{
//Calculate the value of mLut[1]
//for the GetErrorMax method
this->operator()( (M_PI * 0.5) / static_cast<double>(mLut.size()));
assert(mLut[1] != mNotInitializedValue);
}
double Sine::operator()(const double x)
{
const double di = 2.0 * M_PI;
//preY e [ -di, di >
const double preY = std::fmod(x,di);
//y e [ 0.0 , di >
const double y = (preY < 0.0 ? preY + di : preY);
assert( y >= 0.0 );
assert( y < di );
//f = fractional part e [ 0.0, 4.0 >
const double f = 4.0 * (y / di);
assert( f >= 0.0 );
assert( f < 4.0 );
const bool addMinus = ( f < 2.0 ? false : true);
const int index = CalculateIndex(f);
assert(index >= 0);
if (index == static_cast<int>(mLut.size()) )
return ( addMinus == true ? -1.0 : 1.0);
assert(index < static_cast<int>(mLut.size()) );
//Check if already calculated
if (mLut[index] == mNotInitializedValue)
{
mLut[index] = std::sin( std::fmod(y,M_PI) );
assert(mLut[index] >= 0.0);
assert(mLut[index] <= 1.0);
}
if (addMinus == true)
{
return -mLut[index];
}
else
{
return mLut[index];
}
}
const double Sine::GetMaxError() const
{
//The maximal error is twice the distance between
//mLut[0] and mLut[1], because the sine is steepest there
//Because mLut[0] == 0.0, it is twice the value of mLut[1]
assert(mLut[1] != mNotInitializedValue);
return mLut[1] * 2.0;
}
const int Sine::CalculateIndex(const double f) const
{
const int size = mLut.size();
assert( f >= 0.0 );
assert( f < 4.0 );
const double fTransformed = Transform(f);
assert(fTransformed >= 0.0);
assert(fTransformed <= 1.0);
const int index = static_cast<int>(fTransformed * static_cast<double>(size));
assert(index >= 0);
assert(index <= size); // index == size is handled by caller (i.e. operator() )
return index;
}
const double Sine::Transform(const double f) const
{
assert( f >= 0.0 );
assert( f < 4.0 );
if (f < 1.0) return f;
if (f < 2.0) return 1.0 - (f - 1.0);
return Transform(f - 2.0);
}
#pragma package(smart_init)
|
