fix a major inefficiency in cos_sin for huge arguments. add manydigits problems demo. tweaks to manual and some housekeeping
This commit is contained in:
Fredrik Johansson
@@ -41,6 +41,10 @@ To install, unpack the mpmath archive and run
|
||||
|
||||
python setup.py install
|
||||
|
||||
Mpmath can also be installed using
|
||||
|
||||
python -m easy_install mpmath
|
||||
|
||||
The latest development code is available from
|
||||
http://code.google.com/p/mpmath/source/checkout
|
||||
|
||||
|
||||
@@ -0,0 +1,104 @@
|
||||
"""
|
||||
This script calculates solutions to some of the problems from the
|
||||
"Many Digits" competition:
|
||||
http://www.cs.ru.nl/~milad/manydigits/problems.php
|
||||
|
||||
Run with:
|
||||
|
||||
python manydigits.py
|
||||
|
||||
"""
|
||||
|
||||
from mpmath import *
|
||||
from mpmath.lib import *
|
||||
|
||||
dps = 100
|
||||
mp.dps = dps + 10
|
||||
|
||||
def pr(x):
|
||||
"""Return the first dps digits after the decimal point"""
|
||||
x = x._mpf_
|
||||
p = int(dps*3.33 + 10)
|
||||
t = to_fixed(x, p)
|
||||
d = bin_to_radix(t, p, 10, dps)
|
||||
s = str(d).zfill(dps)[-dps:]
|
||||
return s[:dps//2] + "\n" + s[dps//2:]
|
||||
|
||||
print """
|
||||
This script prints answers to a selection of the "Many Digits"
|
||||
competition problems: http://www.cs.ru.nl/~milad/manydigits/problems.php
|
||||
|
||||
The output for each problem is the first 100 first digits after the
|
||||
decimal point in the result.
|
||||
"""
|
||||
|
||||
print "C01: sin(tan(cos(1)))"
|
||||
print pr(sin(tan(cos(1))))
|
||||
print
|
||||
|
||||
print "C02: sqrt(e/pi)"
|
||||
print pr(sqrt(e/pi))
|
||||
print
|
||||
|
||||
print "C03: sin((e+1)^3)"
|
||||
print pr(sin((e+1)**3))
|
||||
print
|
||||
|
||||
print "C04: exp(pi*sqrt(2011))"
|
||||
mp.dps += 65
|
||||
print pr(exp(pi*sqrt(2011)))
|
||||
mp.dps -= 65
|
||||
print
|
||||
|
||||
print "C05: exp(exp(exp(1/2)))"
|
||||
print pr(exp(exp(exp(0.5))))
|
||||
print
|
||||
|
||||
print "C06: arctanh(1-arctanh(1-arctanh(1-arctanh(1/pi))))"
|
||||
print pr(atanh(1-atanh(1-atanh(1-atanh(1/pi)))))
|
||||
print
|
||||
|
||||
print "C07: pi^1000"
|
||||
mp.dps += 505
|
||||
print pr(pi**1000)
|
||||
mp.dps -= 505
|
||||
print
|
||||
|
||||
print "C08: sin(6^(6^6))"
|
||||
print pr(sin(6**(6**6)))
|
||||
print
|
||||
|
||||
print "C09: sin(10*arctan(tanh(pi*(2011^(1/2))/3)))"
|
||||
mp.dps += 150
|
||||
print pr(sin(10*atan(tanh(pi*sqrt(2011)/3))))
|
||||
mp.dps -= 150
|
||||
print
|
||||
|
||||
print "C10: (7+2^(1/5)-5*(8^(1/5)))^(1/3) + 4^(1/5)-2^(1/5)"
|
||||
a = mpf(1)/5
|
||||
print pr(((7 + 2**a - 5*(8**a))**(mpf(1)/3) + 4**a - 2**a))
|
||||
print
|
||||
|
||||
print "C11: tan(2^(1/2))+arctanh(sin(1))"
|
||||
print pr((tan(sqrt(2)) + atanh(sin(1))))
|
||||
print
|
||||
|
||||
print "C12: arcsin(1/e^2) + arcsinh(e^2)"
|
||||
print pr(asin(1/exp(2)) + asinh(exp(2)))
|
||||
print
|
||||
|
||||
print "C17: S= -4*Zeta(2) - 2*Zeta(3) + 4*Zeta(2)*Zeta(3) + 2*Zeta(5)"
|
||||
print pr(-4*zeta(2) - 2*zeta(3) + 4*zeta(2)*zeta(3) + 2*zeta(5))
|
||||
print
|
||||
|
||||
print "C18: Catalan G = Sum{i=0}{\infty}(-1)^i/(2i+1)^2"
|
||||
print pr(catalan)
|
||||
print
|
||||
|
||||
print "C21: Equation exp(cos(x)) = x"
|
||||
print pr(secant(lambda x: exp(cos(x))-x, 1))
|
||||
print
|
||||
|
||||
print "C22: J = integral(sin(sin(sin(x)))), x=0..1"
|
||||
print pr(quadts(lambda x: sin(sin(sin(x))), 0, 1))
|
||||
print
|
||||
+5
-5
@@ -387,7 +387,7 @@ pre.literal-block, pre.doctest-block {
|
||||
</div>
|
||||
<div class="section">
|
||||
<h1><a class="toc-backref" href="#id2" id="basics" name="basics">2 Basics</a></h1>
|
||||
<p>For download and installation instructions, please refer to the README or the mpmath website (in most cases, installation should be as simple as running <tt class="docutils literal"><span class="pre">python</span> <span class="pre">easy_install</span> <span class="pre">mpmath</span></tt>). After the setup has completed, you can fire up the interactive Python interpreter and try the following:</p>
|
||||
<p>For download and installation instructions, please refer to the README or the mpmath website. After the setup has completed, you can fire up the interactive Python interpreter and try the following:</p>
|
||||
<pre class="literal-block">
|
||||
>>> from mpmath import *
|
||||
>>> mp.dps = 50
|
||||
@@ -412,7 +412,7 @@ mpf('2.0')
|
||||
>>> mpf("inf")
|
||||
mpf('+inf')
|
||||
</pre>
|
||||
<p>An <tt class="docutils literal"><span class="pre">mpc</span></tt> represents a complex number in rectangular form as a pair of <tt class="docutils literal"><span class="pre">mpf</span></tt> instances. It can be constructed from a Python <tt class="docutils literal"><span class="pre">complex</span></tt>, a real number, or a pair of real numbers:</p>
|
||||
<p>The <tt class="docutils literal"><span class="pre">mpc</span></tt> type represents a complex number in rectangular form as a pair of <tt class="docutils literal"><span class="pre">mpf</span></tt> instances. It can be constructed from a Python <tt class="docutils literal"><span class="pre">complex</span></tt>, a real number, or a pair of real numbers:</p>
|
||||
<pre class="literal-block">
|
||||
>>> mpc(2,3)
|
||||
mpc(real='2.0', imag='3.0')
|
||||
@@ -698,7 +698,7 @@ mpf('0.2500000000000000000000000000000000000000057')
|
||||
<td>Natural logarithm (optionally base-b logarithm)</td>
|
||||
</tr>
|
||||
<tr><td><tt class="docutils literal"><span class="pre">power(x,y)</span></tt></td>
|
||||
<td>Power, x^y</td>
|
||||
<td>Power, <tt class="docutils literal"><span class="pre">x**y</span></tt></td>
|
||||
</tr>
|
||||
<tr><td><tt class="docutils literal"><span class="pre">cos(x)</span></tt></td>
|
||||
<td>Cosine</td>
|
||||
@@ -901,14 +901,14 @@ mpf('0.2500000000000000000000000000000000000000057')
|
||||
</div>
|
||||
<div class="section">
|
||||
<h3><a id="error-detection" name="error-detection">3.1.3 Error detection</a></h3>
|
||||
<p>The tanh-sinh algorithm is not suitable for adaptive quadrature, and does not perform well if there are singularities between the endpoints or if the integrand is very bumpy or oscillatory (such integrals should manually be split into smaller pieces). If the <tt class="docutils literal"><span class="pre">error</span></tt> option is set, <tt class="docutils literal"><span class="pre">quadts</span></tt> will return an error estimate along with the result; although this estimate is not always correct, it can be useful for debugging. You can also pass <tt class="docutils literal"><span class="pre">quadts</span></tt> the option <tt class="docutils literal"><span class="pre">verbose=True</span></tt> to show detailed progress.</p>
|
||||
<p>The tanh-sinh algorithm is not suitable for adaptive quadrature, and does not perform well if there are singularities between the endpoints or if the integrand is oscillatory (such integrals should manually be split into smaller pieces). If the <tt class="docutils literal"><span class="pre">error</span></tt> option is set, <tt class="docutils literal"><span class="pre">quadts</span></tt> will return an error estimate along with the result; although this estimate is not always correct, it can be useful for debugging. You can also pass <tt class="docutils literal"><span class="pre">quadts</span></tt> the option <tt class="docutils literal"><span class="pre">verbose=True</span></tt> to show detailed progress.</p>
|
||||
<p>A simple example where the algorithm fails is the function f(<em>x</em>) = abs(sin(<em>x</em>)), which is not smooth at <em>x</em> = pi. In this case, a close value is calculated, but the result is nowhere near the target accuracy; however, <tt class="docutils literal"><span class="pre">quadts</span></tt> gives a good estimate of the magnitude of the error:</p>
|
||||
<pre class="literal-block">
|
||||
>>> mp.dps = 15
|
||||
>>> quadts(lambda x: abs(sin(x)), 0, 2*pi, error=True)
|
||||
(mpf('3.9990089417677899'), mpf('0.001'))
|
||||
</pre>
|
||||
<p>Attempting to evaluate oscillatory integrals on large intervals by means of the tanh-sinh method is generally futile. This integral should be pi/2 = 1.57:</p>
|
||||
<p>This highly oscillatory integral should be pi/2 = 1.57:</p>
|
||||
<pre class="literal-block">
|
||||
>>> print quadts(lambda x: sin(x)/x, 0, inf, error=True)
|
||||
(mpf('2.3840907358976544'), mpf('1.0'))
|
||||
|
||||
Binary file not shown.
+5
-5
@@ -28,7 +28,7 @@ This manual gives an introduction to mpmath's major features. Some supplementary
|
||||
Basics
|
||||
======
|
||||
|
||||
For download and installation instructions, please refer to the README or the mpmath website (in most cases, installation should be as simple as running ``python easy_install mpmath``). After the setup has completed, you can fire up the interactive Python interpreter and try the following::
|
||||
For download and installation instructions, please refer to the README or the mpmath website. After the setup has completed, you can fire up the interactive Python interpreter and try the following::
|
||||
|
||||
>>> from mpmath import *
|
||||
>>> mp.dps = 50
|
||||
@@ -55,7 +55,7 @@ Mpmath provides two main numerical types: ``mpf`` and ``mpc``. The ``mpf`` type
|
||||
>>> mpf("inf")
|
||||
mpf('+inf')
|
||||
|
||||
An ``mpc`` represents a complex number in rectangular form as a pair of ``mpf`` instances. It can be constructed from a Python ``complex``, a real number, or a pair of real numbers::
|
||||
The ``mpc`` type represents a complex number in rectangular form as a pair of ``mpf`` instances. It can be constructed from a Python ``complex``, a real number, or a pair of real numbers::
|
||||
|
||||
>>> mpc(2,3)
|
||||
mpc(real='2.0', imag='3.0')
|
||||
@@ -313,7 +313,7 @@ Function Description
|
||||
``hypot(x,y)`` Euclidean norm
|
||||
``exp(x)`` Exponential function
|
||||
``log(x,b)`` Natural logarithm (optionally base-b logarithm)
|
||||
``power(x,y)`` Power, x^y
|
||||
``power(x,y)`` Power, ``x**y``
|
||||
``cos(x)`` Cosine
|
||||
``sin(x)`` Sine
|
||||
``tan(x)`` Tangent
|
||||
@@ -454,7 +454,7 @@ While double integrals are reasonably fast, even a simple triple integral at ver
|
||||
Error detection
|
||||
...............
|
||||
|
||||
The tanh-sinh algorithm is not suitable for adaptive quadrature, and does not perform well if there are singularities between the endpoints or if the integrand is very bumpy or oscillatory (such integrals should manually be split into smaller pieces). If the ``error`` option is set, ``quadts`` will return an error estimate along with the result; although this estimate is not always correct, it can be useful for debugging. You can also pass ``quadts`` the option ``verbose=True`` to show detailed progress.
|
||||
The tanh-sinh algorithm is not suitable for adaptive quadrature, and does not perform well if there are singularities between the endpoints or if the integrand is oscillatory (such integrals should manually be split into smaller pieces). If the ``error`` option is set, ``quadts`` will return an error estimate along with the result; although this estimate is not always correct, it can be useful for debugging. You can also pass ``quadts`` the option ``verbose=True`` to show detailed progress.
|
||||
|
||||
A simple example where the algorithm fails is the function f(*x*) = abs(sin(*x*)), which is not smooth at *x* = pi. In this case, a close value is calculated, but the result is nowhere near the target accuracy; however, ``quadts`` gives a good estimate of the magnitude of the error::
|
||||
|
||||
@@ -462,7 +462,7 @@ A simple example where the algorithm fails is the function f(*x*) = abs(sin(*x*)
|
||||
>>> quadts(lambda x: abs(sin(x)), 0, 2*pi, error=True)
|
||||
(mpf('3.9990089417677899'), mpf('0.001'))
|
||||
|
||||
Attempting to evaluate oscillatory integrals on large intervals by means of the tanh-sinh method is generally futile. This integral should be pi/2 = 1.57::
|
||||
This highly oscillatory integral should be pi/2 = 1.57::
|
||||
|
||||
>>> print quadts(lambda x: sin(x)/x, 0, inf, error=True)
|
||||
(mpf('2.3840907358976544'), mpf('1.0'))
|
||||
|
||||
+14
-15
@@ -1641,14 +1641,6 @@ def sin_taylor(x, prec):
|
||||
k += 2
|
||||
return s
|
||||
|
||||
def trig_reduce(x, prec):
|
||||
pi_ = pi_fixed(prec)
|
||||
pi4 = pi_ >> 2
|
||||
pi2 = pi_ >> 1
|
||||
n, rem = divmod(x + pi4, pi2)
|
||||
rem -= pi4
|
||||
return n, rem
|
||||
|
||||
def cos_sin(x, prec, rounding):
|
||||
"""Simultaneously compute (cos(x), sin(x)) for real x."""
|
||||
|
||||
@@ -1661,6 +1653,7 @@ def cos_sin(x, prec, rounding):
|
||||
return fone, fzero
|
||||
|
||||
magnitude = bc + exp
|
||||
abs_mag = abs(magnitude)
|
||||
|
||||
# Very close to 0
|
||||
if magnitude < -prec:
|
||||
@@ -1682,17 +1675,23 @@ def cos_sin(x, prec, rounding):
|
||||
s = fadd(x, (1, 1, magnitude-prec-4, 1), prec, rounding)
|
||||
return c, s
|
||||
|
||||
bits_from_unit = abs(magnitude)
|
||||
|
||||
prec1 = prec + bits_from_unit + 15
|
||||
wp = prec1
|
||||
wp = wp1 = prec + 15
|
||||
|
||||
while 1:
|
||||
n, rx = trig_reduce(to_fixed(x, wp), wp)
|
||||
# Reduce modulo pi/4
|
||||
rp = wp + abs_mag
|
||||
a = to_fixed(x, rp)
|
||||
pi_ = pi_fixed(rp)
|
||||
pi4 = pi_ >> 2
|
||||
pi2 = pi_ >> 1
|
||||
n, rx = divmod(a+pi4, pi2)
|
||||
rx -= pi4
|
||||
rx >>= abs_mag
|
||||
|
||||
# If we're close to a root, we have to increase the
|
||||
# fixed-point precision to obtain full relative accuracy
|
||||
if abs(rx >> (prec1-8)) < 10:
|
||||
wp += prec1 - bitcount(abs(rx))
|
||||
if abs(rx >> (wp1 - 8)) < 10:
|
||||
wp += wp1 - bitcount(abs(rx))
|
||||
else:
|
||||
break
|
||||
|
||||
|
||||
@@ -7,6 +7,6 @@ setup(name='mpmath',
|
||||
author='Fredrik Johansson',
|
||||
author_email='fredrik.johansson@gmail.com',
|
||||
license = 'BSD',
|
||||
packages=['mpmath', 'mpmath/lib', 'mpmath/apps', 'mpmath/tests'],
|
||||
packages=['mpmath', 'demo', 'mpmath/tests'],
|
||||
classifiers=['Topic :: Scientific/Engineering :: Mathematics']
|
||||
)
|
||||
|
||||
Reference in New Issue
Block a user