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) * |
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) * |
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) * |
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) * |
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) * |
Date: 2010-11-21 02:47 |
Added a comment to the docstring.
See r86631.
|
|
Date |
User |
Action |
Args |
2022-04-11 14:56:56 | admin | set | github: 52018 |
2010-11-21 02:50:49 | rhettinger | set | status: open -> closed resolution: fixed |
2010-11-21 02:47:35 | rhettinger | set | messages:
+ msg121841 |
2010-01-29 20:33:22 | Steve Howell | set | nosy:
+ Steve Howell messages:
+ msg98536
|
2010-01-24 19:58:56 | ahojnnes | set | messages:
+ msg98249 |
2010-01-24 19:58:07 | ahojnnes | set | messages:
+ msg98248 |
2010-01-24 19:39:00 | rhettinger | set | priority: normal -> low assignee: georg.brandl -> rhettinger messages:
+ msg98246
|
2010-01-24 14:46:41 | mark.dickinson | set | files:
+ 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:30 | ahojnnes | set | messages:
+ msg98227 |
2010-01-24 13:48:48 | mark.dickinson | set | messages:
+ msg98226 |
2010-01-24 13:46:11 | ahojnnes | set | messages:
+ msg98225 |
2010-01-24 13:26:39 | mark.dickinson | set | nosy:
+ rhettinger messages:
+ msg98222
|
2010-01-24 13:16:25 | ezio.melotti | set | priority: normal nosy:
+ mark.dickinson
stage: test needed |
2010-01-24 13:09:49 | ahojnnes | create | |