Download code

Jump to: navigation, search

Back to Arbitrary-precision_elementary_mathematical_functions_(Python)

Download for Windows: single file, zip

Download for UNIX: single file, zip, tar.gz, tar.bz2

hpfunc.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python)?oldid=19016
14 
15 from decimal import *
16 
17 def exp(x):
18     getcontext().prec += 3
19     if x < 0:
20         s = 1 / exp(abs(x))
21     else:
22         xpower = 1
23         ns = 0
24         s = 1
25         n = 0
26         factorial = 1
27         while s != ns:
28             s = ns
29             term = Decimal(xpower) / factorial
30             ns = s + term
31             xpower *= x
32             n += 1
33             factorial *= n
34     getcontext().prec -= 3
35     return +s
36 
37 def test_exp():
38     assert exp(4) == +Decimal("54.59815003314423907811026120286")
39     assert exp(0) == 1
40     assert exp(-8) == +Decimal("0.0003354626279025118388213891257809")
41     assert exp(Decimal("0.6931471805599453094172321214582")) == 2
42 
43 from math import log as _flog
44 
45 def log(x):
46     if x < 0:
47         return Decimal("NaN")
48     if x == 0:
49         return Decimal("-inf")
50     getcontext().prec += 3
51     eps = Decimal("10")**(-getcontext().prec+2)
52     # A good initial estimate is needed
53     r = Decimal(repr(_flog(float(x))))
54     while 1:
55         r2 = r - 1 + x/exp(r)
56         if abs(r2-r) < eps:
57             break
58         else:
59             r = r2
60     getcontext().prec -= 3
61     return +r
62 def test_log():
63     A = [4, 1, Decimal("0.67753892")]
64     for p in [5, 40, 28]:
65         getcontext().prec = p
66         eps = Decimal("10")**(-getcontext().prec + 2)
67         for a in A:
68             w = exp(log(a))
69             assert abs(w - a) < eps
70 
71 if __name__ == "__main__":
72     test_exp()
73     test_log()
74 


hijacker
hijacker
hijacker
hijacker