opensine.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file opensine.h
///
/// An open-interval version of the sine function.
/// This file contains the definition of the sine function in the interval
/// arithmetic based on open intervals, as opposed to closed intervals.
/// Based on this definition, also the other trigonometric functions
/// can be implemented.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2007-2010 by Pawel Pilarczyk.
//
// This file is part of my research program package.  This is free software;
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation;
// either version 2 of the License, or (at your option) any later version.
//
// This software 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on December 28, 2008. Last revision: December 30, 2008.


#ifndef _FINRESDYN_OPENSINE_H_
#define _FINRESDYN_OPENSINE_H_


// include local header files
#include "rounding.h"
#include "numberpi.h"


// --------------------------------------------------
// ---- a class that computes the sine function -----
// --------------------------------------------------

/// Computes the sine function using the expansion
/// into the Taylor series at pi/2.
template <class NumType, int nTerms>
class tOpenSine
{
public:
        /// Computes the sine function in the open interval arithmetic.
        static int compute (const NumType &xLeft, const NumType &xRight,
                NumType &yLeft, NumType &yRight);

private:
        /// Computes the provided number of terms in the expansion
        /// of sin(pi+x) in the Taylor series at pi.
        static void computeTaylorTerms (const NumType &x,
                NumType *lower_x2k, NumType *upper_x2k, int n);

        /// Returns a lower bound for the sine function.
        /// The argument must be in the interval (0,pi).
        static NumType lowerBoundTaylor (const NumType &x);

        /// Returns an upper bound for the sine function.
        /// The argument must be in the interval (0,pi).
        static NumType upperBoundTaylor (const NumType &x);

        /// Computes the sine function in the open interval arithmetic
        /// with the arguments shifted so that xLeft is in [0,pi].
        static int shifted1pi (const NumType &xLeft, const NumType &xRight,
                NumType &yLeft, NumType &yRight);

        /// Computes the sine function in the open interval arithmetic
        /// with the arguments shifted so that xLeft is in [0,2pi].
        static int shifted2pi (const NumType &xLeft, const NumType &xRight,
                NumType &yLeft, NumType &yRight);

}; /* class tOpenSine */

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

template <class NumType, int nTerms>
inline void tOpenSine<NumType,nTerms>::computeTaylorTerms
        (const NumType &x, NumType *lower_x2k, NumType *upper_x2k, int n)
{
        // prepare shortcuts for the rounding mode and the number pi
        typedef tRounding<NumType> rnd;
        typedef tNumberPi<NumType> pi;

        // prepare a lower bound and an upper bound for x^2
        NumType lowerDiff = (x <= pi::lowerBound () / 2) ?
                rnd::sub_down (pi::lowerBound () / 2, x) :
                rnd::sub_down (x, pi::upperBound () / 2);
        if (lowerDiff < 0)
                lowerDiff = 0;
        NumType lower_x2 = rnd::mul_down (lowerDiff, lowerDiff);
        NumType upperDiff =  (x <= pi::upperBound () / 2) ?
                rnd::sub_up (pi::upperBound () / 2, x) :
                rnd::sub_up (x, pi::lowerBound () / 2);
        NumType upper_x2 = rnd::mul_up (upperDiff, upperDiff);

        // compute the terms, including the factorial, but without the sign,
        // for the Taylor expansion as in the following formula:
        // \sin (\pi/2 + x) = \sum_{k = 0}^{\infty} (-1)^k x^{2k} / (2k)!
        for (int k = 0; k < n; ++ k)
        {
                if (k == 0)
                {
                        lower_x2k [0] = 1;
                        upper_x2k [0] = 1;
                }
                else if (k == 1)
                {
                        lower_x2k [1] = lower_x2 / 2;
                        upper_x2k [1] = upper_x2 / 2;
                }
                else
                {
                        long int factor = ((k << 1) - 1) * (k << 1);
                        lower_x2k [k] = rnd::mul_down (lower_x2k [k - 1],
                                lower_x2);
                        lower_x2k [k] = rnd::div_down (lower_x2k [k],
                                factor);
                        upper_x2k [k] = rnd::mul_up (upper_x2k [k - 1],
                                upper_x2);
                        upper_x2k [k] = rnd::div_up (upper_x2k [k],
                                factor);
                }
        }

        return;
} /* computeTaylorTerms */

template <class NumType, int nTerms>
inline NumType tOpenSine<NumType,nTerms>::lowerBoundTaylor (const NumType &x)
{
        // prepare a shortcut for the rounding mode
        typedef tRounding<NumType> rnd;

        // determine an even strict upper bound for the number of terms
        int n = (nTerms & 1) ? (nTerms + 1) : (nTerms);

        // prepare all the terms to add with factorials included
        NumType lower_x2k [nTerms + 1];
        NumType upper_x2k [nTerms + 1];
        computeTaylorTerms (x, lower_x2k, upper_x2k, n);

        // compute the sum for the given number of terms
        NumType result = 0;
        for (int k = n - 1; k >= 0; -- k)
        {
                // subtract the term if its number is odd
                if (k & 1)
                        result = rnd::sub_down (result, upper_x2k [k]);

                // add the term if its number is even
                else
                        result = rnd::add_down (result, lower_x2k [k]);
        }

        return result;
} /* lowerBoundTaylor */

template <class NumType, int nTerms>
inline NumType tOpenSine<NumType,nTerms>::upperBoundTaylor (const NumType &x)
{
        // prepare a shortcut for the rounding mode
        typedef tRounding<NumType> rnd;

        // determine an odd strict upper bound for the number of terms
        int n = (nTerms & 1) ? (nTerms) : (nTerms + 1);

        // prepare all the terms to add with factorials included
        NumType lower_x2k [nTerms + 1];
        NumType upper_x2k [nTerms + 1];
        computeTaylorTerms (x, lower_x2k, upper_x2k, n);

        // compute the sum for the given number of terms
        NumType result = 0;
        for (int k = n - 1; k >= 0; -- k)
        {
                // subtract the term if its number is odd
                if (k & 1)
                        result = rnd::sub_up (result, lower_x2k [k]);

                // add the term if its number is even
                else
                        result = rnd::add_up (result, upper_x2k [k]);
        }

        return result;
} /* upperBoundTaylor */

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

template <class NumType, int nTerms>
inline int tOpenSine<NumType,nTerms>::compute
        (const NumType &xLeft, const NumType &xRight,
        NumType &yLeft, NumType &yRight)
{
        // prepare shortcuts for the rounding mode and the number pi
        typedef tRounding<NumType> rnd;
        typedef tNumberPi<NumType> pi;

        // prepare bounds for the number 2pi
        NumType lower2pi = pi::lowerBound () * 2;
        NumType upper2pi = pi::upperBound () * 2;

        // shift the interval so that the left bound is in [0,2pi]
        if (xLeft < 0)
        {
                int periods = static_cast<int> ((-xLeft) / upper2pi) + 1;
                NumType shiftMin = rnd::mul_down (periods, lower2pi);
                NumType shiftMax = rnd::mul_up (periods, upper2pi);
                NumType newLeft = rnd::add_down (xLeft, shiftMin);
                NumType newRight = rnd::add_up (xRight, shiftMax);
                while (newLeft < 0)
                {
                //      std::cout << "Debug: sine add (" << newLeft << ").\n";
                        newLeft = rnd::add_down (newLeft, lower2pi);
                        newRight = rnd::add_up (newRight, upper2pi);
                }
                while (newLeft > lower2pi)
                {
                //      std::cout << "Debug: sine sub (" << newLeft << ").\n";
                        newLeft = rnd::sub_down (newLeft, lower2pi);
                        newRight = rnd::sub_up (newRight, upper2pi);
                }
                return shifted2pi (newLeft, newRight, yLeft, yRight);
        }
        else if (xLeft > upper2pi)
        {
                int periods = static_cast<int> (xLeft / upper2pi);
                NumType shiftMin = rnd::mul_down (periods, lower2pi);
                NumType shiftMax = rnd::mul_up (periods, upper2pi);
                NumType newLeft = rnd::sub_down (xLeft, shiftMax);
                NumType newRight = rnd::sub_up (xRight, shiftMin);
                while (newLeft > upper2pi)
                {
                //      std::cout << "Debug: sine+sub (" << newLeft << ").\n";
                        newLeft = rnd::sub_down (newLeft, upper2pi);
                        newRight = rnd::sub_up (newRight, lower2pi);
                }
                return shifted2pi (newLeft, newRight, yLeft, yRight);
        }
        else
        {
                return shifted2pi (xLeft, xRight, yLeft, yRight);
        }
} /* compute */

template <class NumType, int nTerms>
inline int tOpenSine<NumType,nTerms>::shifted1pi
        (const NumType &xLeft, const NumType &xRight,
        NumType &yLeft, NumType &yRight)
{
        // prepare shortcuts for the rounding mode and the number pi
        typedef tRounding<NumType> rnd;
        typedef tNumberPi<NumType> pi;

        // prepare bounds for the number pi
        const NumType lowerPi = pi::lowerBound ();
        const NumType upperPi = pi::upperBound ();

        // prepare bounds for the number 2pi
        const NumType lower2pi = lowerPi * 2;
        const NumType upper2pi = upperPi * 2;

        // check if the width of the interval exceeds the period of sine
//      if (rnd::sub_up (xRight, xLeft) >= lower2pi)
//      {
//              yLeft = rnd::sub_down (-1, rnd::min_number ());
//              yRight = rnd::add_up (1, rnd::min_number ());
//              return 0;
//      }

        // prepare bounds for the number pi/2
        const NumType lowerPi2 = lowerPi / 2;
        const NumType upperPi2 = upperPi / 2;

        // if both bounds are to the left from pi/2
        if (xRight <= upperPi2)
        {
                yLeft = lowerBoundTaylor (xLeft);
                yRight = upperBoundTaylor (xRight);
                return 0;
        }

        // prepare bounds for the number 3pi/2
        NumType upper3pi2 = rnd::add_up (upperPi, upperPi2);
//      NumType lower3pi2 = rnd::add_down (lowerPi, lowerPi2);

        // if the left bound is below pi/2, but the right bound is larger
        if (xLeft <= upperPi2)
        {
                // if the right bound is in [pi/2, pi]
                if (xRight <= upperPi)
                {
                        if (rnd::sub_up (upperPi2, xLeft) <=
                                rnd::sub_down (xRight, upperPi2))
                        {
                                yLeft = lowerBoundTaylor (xRight);
                        }
                        else if (rnd::sub_down (lowerPi2, xLeft) >=
                                rnd::sub_up (xRight, lowerPi2))
                        {
                                yLeft = lowerBoundTaylor (xLeft);
                        }
                        else
                        {
                                NumType leftVal = lowerBoundTaylor (xLeft);
                                NumType rightVal = lowerBoundTaylor (xRight);
                                yLeft = (leftVal < rightVal) ?
                                        leftVal : rightVal;
                        }
                        yRight = rnd::add_up (1, rnd::min_number ());
                        return 0;
                }

                // if the right bound is in [pi, 3/2 pi]
                if (xRight <= upper3pi2)
                {
                        yLeft = -upperBoundTaylor
                                (rnd::sub_up (xRight, lowerPi));
                        yRight = rnd::add_up (1, rnd::min_number ());
                        return 0;
                }

                // if the right bound is beyond 3/2 pi
                yLeft = rnd::sub_down (-1, rnd::min_number ());
                yRight = rnd::add_up (1, rnd::min_number ());
                return 0;
        }

        // if the left bound is between pi/2 and pi
//      if (xLeft <= upperPi)
        {
                // if the right bound is also between pi/2 and pi
                if (xRight <= upperPi)
                {
                        yLeft = lowerBoundTaylor (xRight);
                        yRight = upperBoundTaylor (xLeft);
                        return 0;
                }

                // if the right bound is in [pi, 3/2 pi]
                if (xRight <= upper3pi2)
                {
                        yLeft = -upperBoundTaylor
                                (rnd::sub_up (xRight, lowerPi));
                        yRight = upperBoundTaylor (xLeft);
                        return 0;
                }

                // if the right bound is in [3/2 pi, 2pi]
                if (xRight <= upper2pi)
                {
                        yLeft = rnd::sub_down (-1, rnd::min_number ());
                        yRight = upperBoundTaylor (xLeft);
                        return 0;
                }

                // if the right bound is in [2pi, 2pi+pi/2]
                if (xRight <= rnd::add_up (upper2pi, upperPi2))
                {
                        if (rnd::sub_up (upperPi, xLeft) <=
                                rnd::sub_down (xRight, upper2pi))
                        {
                                yRight = upperBoundTaylor
                                        (rnd::sub_up (xRight, lower2pi));
                        }
                        else if (rnd::sub_down (lowerPi, xLeft) >=
                                rnd::sub_up (xRight, lower2pi))
                        {
                                yRight = upperBoundTaylor (xLeft);
                        }
                        else
                        {
                                NumType leftVal = upperBoundTaylor (xLeft);
                                NumType rightVal = upperBoundTaylor
                                        (rnd::sub_up (xRight, lower2pi));
                                yRight = (leftVal < rightVal) ?
                                        rightVal : leftVal;
                        }
                        yLeft = rnd::sub_down (-1, rnd::min_number ());
                        return 0;
                }

                // otherwise the full range of values is covered
                yLeft = rnd::sub_down (-1, rnd::min_number ());
                yRight = rnd::add_up (1, rnd::min_number ());
                return 0;
        }
} /* shifted1pi */

template <class NumType, int nTerms>
inline int tOpenSine<NumType,nTerms>::shifted2pi
        (const NumType &xLeft, const NumType &xRight,
        NumType &yLeft, NumType &yRight)
{
        // prepare shortcuts for the rounding mode and the number pi
        typedef tRounding<NumType> rnd;
        typedef tNumberPi<NumType> pi;

        // if the left endpoint is below pi then fine
        if (xLeft <= pi::upperBound ())
                return shifted1pi (xLeft, xRight, yLeft, yRight);

        // compute the function with the arguments shifted by pi
        // (note the reversed order of results as they will be negated)
        shifted1pi (rnd::sub_down (xLeft, pi::upperBound ()),
                rnd::sub_up (xRight, pi::upperBound ()),
                yRight, yLeft);

        // flip the resulting interval and finish
        yLeft = -yLeft;
        yRight = -yRight;
        return 0;
} /* shifted2pi */


#endif // _FINRESDYN_OPENSINE_H_