classification
Title: sin/cos function in decimal-docs
Type: behavior Stage: patch review
Components: Documentation Versions: Python 3.1, Python 3.2, Python 2.7, Python 2.6
process
Status: closed Resolution: fixed
Dependencies: Superseder:
Assigned To: rhettinger Nosy List: Steve Howell, ahojnnes, georg.brandl, mark.dickinson, rhettinger
Priority: low Keywords: patch

Created on 2010-01-24 13:09 by ahojnnes, last changed 2010-11-21 02:50 by rhettinger. This issue is now closed.

Files
File name Uploaded Description Edit
issue7770.patch mark.dickinson, 2010-01-24 14:46
Messages (11)
msg98221 - (view) Author: Johannes Schönberger (ahojnnes) Date: 2010-01-24 13:09
Unfortunately the sin/cos function in the decimal docs seems to return false results.

In [33]: sin(decimal.Decimal('75'))
Out[33]: Decimal('0.377879483645203210442266845614')

In [34]: sin(decimal.Decimal('76'))
Out[34]: Decimal('-2.08828460009724889326220807212')

In [42]: sin(decimal.Decimal('100'))
Out[42]: Decimal('-58433045378.5877442230422602000')

#####

In [37]: cos(decimal.Decimal('79'))
Out[37]: Decimal('-14.3603762068086273628189466621')

In [38]: cos(decimal.Decimal('77'))
Out[38]: Decimal('-13.6219693138941571794243404126')

In [39]: cos(decimal.Decimal('75'))
Out[39]: Decimal('1.25761570869417008177315814688')

In [40]: cos(decimal.Decimal('72'))
Out[40]: Decimal('-1.02323683857735456186757099584')
msg98222 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-01-24 13:26
I think that recipe is meant more as a simple example of what's possible than as a bullet-proof library function.  Can you suggest changes that would improve accuracy while also keeping the recipe clear and simple?

The usual method for sin/cos is to do an initial reduction by an integral multiple of pi/2 so that the argument lies in [-pi/4, pi/4].  (Though that's still problematic for very large arguments, where pi needs to be computed to many more places than the desired output accuracy.)

Perhaps it would be best to just add a note to the docstring indicating that this is an implementation that's suitable for small arguments; for large arguments, more sophisticated algorithms are needed.
msg98225 - (view) Author: Johannes Schönberger (ahojnnes) Date: 2010-01-24 13:46
This is the version I would suggest and which is quite accurate unless the numbers get extremely large.

def sin(x):
    """Return the sine of x as measured in radians.

    >>> print sin(Decimal('0.5'))
    0.4794255386042030002732879352
    >>> print sin(0.5)
    0.479425538604
    >>> print sin(0.5+0j)
    (0.479425538604+0j)

    """
    x = x - pi()*int(x / pi())
    getcontext().prec += 2
    i, lasts, s, fact, num, sign = 1, 0, x, 1, x, 1
    while s != lasts:
        lasts = s
        i += 2
        fact *= i * (i-1)
        num *= x * x
        sign *= -1
        s += num / fact * sign
    getcontext().prec -= 2
    return +s
msg98226 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-01-24 13:48
Hmm.  Isn't sine periodic with period 2*pi, not pi?

I'd also suggest making use of remainder_near.
msg98227 - (view) Author: Johannes Schönberger (ahojnnes) Date: 2010-01-24 13:56
stupid, it is much better to just use the modulo operator. The same works with cos...

def sin(x):
    """Return the sine of x as measured in radians.

    >>> print sin(Decimal('0.5'))
    0.4794255386042030002732879352
    >>> print sin(0.5)
    0.479425538604
    >>> print sin(0.5+0j)
    (0.479425538604+0j)

    """
    x %= pi()
    getcontext().prec += 2
    i, lasts, s, fact, num, sign = 1, 0, x, 1, x, 1
    while s != lasts:
        lasts = s
        i += 2
        fact *= i * (i-1)
        num *= x * x
        sign *= -1
        s += num / fact * sign
    getcontext().prec -= 2
    return +s

def cos(x):
    """Return the cosine of x as measured in radians.

    >>> print cos(Decimal('0.5'))
    0.8775825618903727161162815826
    >>> print cos(0.5)
    0.87758256189
    >>> print cos(0.5+0j)
    (0.87758256189+0j)

    """
    x %= pi()
    getcontext().prec += 2
    i, lasts, s, fact, num, sign = 0, 0, 1, 1, 1, 1
    while s != lasts:
        lasts = s
        i += 2
        fact *= i * (i-1)
        num *= x * x
        sign *= -1
        s += num / fact * sign
    getcontext().prec -= 2
    return +s
msg98228 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-01-24 14:46
Johannes: note that the reduction needs to be by 2*pi, not pi.  The remainder_near method is slightly better than the modulo operator here:
remainder_near reduces to the range [-pi, pi], while the modulo operator reduces to the range (-2*pi, 2*pi), so ends up passing larger arguments to the Taylor series.  (And it's a nice excuse to show what remainder_near is good for!)  Also, the reduction should happen *after* the temporary increase in context precision.

Here's a doc patch, that just adds a single line:

x = x.remainder_near(2*pi())

after getcontext().prec += 2, for both the sin and cos recipes.  It also removes the float and complex examples in the docstrings, since remainder_near doesn't exist for those types.

Raymond, is this change okay with you?
msg98246 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2010-01-24 19:38
I'll look at it when I get a chance.  Want to continue to show how the same recipe can work for different input types.  I may just document the function's range of inputs.
msg98248 - (view) Author: Johannes Schönberger (ahojnnes) Date: 2010-01-24 19:58
This adds two further lines, so you can pass int's, float's or Decimal types to the function.

def sin(x):
    """Return the sine of x as measured in radians.

    >>> print sin(Decimal('0.5'))
    0.4794255386042030002732879352
    >>> print sin(0.5)
    0.479425538604
    >>> print sin(0.5+0j)
    (0.479425538604+0j)

    """
    getcontext().prec += 2
    if not isinstance(x, Decimal):
        x = Decimal(x)
    x = x.remainder_near(2*pi())
    i, lasts, s, fact, num, sign = 1, 0, x, 1, x, 1
    while s != lasts:
        lasts = s
        i += 2
        fact *= i * (i-1)
        num *= x * x
        sign *= -1
        s += num / fact * sign
    getcontext().prec -= 2
    return +s
msg98249 - (view) Author: Johannes Schönberger (ahojnnes) Date: 2010-01-24 19:58
sorry, forgot the str:

def rsin(x):
    """Return the sine of x as measured in radians.

    >>> print sin(Decimal('0.5'))
    0.4794255386042030002732879352
    >>> print sin(0.5)
    0.479425538604
    >>> print sin(0.5+0j)
    (0.479425538604+0j)

    """
    getcontext().prec += 2
    if not isinstance(x, Decimal):
        x = Decimal(str(x))
    x = x.remainder_near(2*pi())
    i, lasts, s, fact, num, sign = 1, 0, x, 1, x, 1
    while s != lasts:
        lasts = s
        i += 2
        fact *= i * (i-1)
        num *= x * x
        sign *= -1
        s += num / fact * sign
    getcontext().prec -= 2
    return +s
msg98536 - (view) Author: showell (Steve Howell) Date: 2010-01-29 20:33
+1 on showing off remainder_near.  I recently wrote a program where I reinvented the logic (on a unit circle too), not knowing it was already implemented in Decimal.
msg121841 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2010-11-21 02:47
Added a comment to the docstring.
See r86631.
History
Date User Action Args
2010-11-21 02:50:49rhettingersetstatus: open -> closed
resolution: fixed
2010-11-21 02:47:35rhettingersetmessages: + msg121841
2010-01-29 20:33:22Steve Howellsetnosy: + Steve Howell
messages: + msg98536
2010-01-24 19:58:56ahojnnessetmessages: + msg98249
2010-01-24 19:58:07ahojnnessetmessages: + msg98248
2010-01-24 19:39:00rhettingersetpriority: normal -> low
assignee: georg.brandl -> rhettinger
messages: + msg98246
2010-01-24 14:46:41mark.dickinsonsetfiles: + issue7770.patch
versions: + Python 3.1, Python 2.7, Python 3.2
messages: + msg98228

keywords: + patch
stage: test needed -> patch review
2010-01-24 13:56:30ahojnnessetmessages: + msg98227
2010-01-24 13:48:48mark.dickinsonsetmessages: + msg98226
2010-01-24 13:46:11ahojnnessetmessages: + msg98225
2010-01-24 13:26:39mark.dickinsonsetnosy: + rhettinger
messages: + msg98222
2010-01-24 13:16:25ezio.melottisetpriority: normal
nosy: + mark.dickinson

stage: test needed
2010-01-24 13:09:49ahojnnescreate