Download code

From LiteratePrograms

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 # Copyright (c) 2010 the authors listed at the following URL, and/or
 2 # the authors of referenced articles or incorporated external code:
 3 # http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python)?action=history&offset=20060430000847
 4 # 
 5 # Permission is hereby granted, free of charge, to any person obtaining
 6 # a copy of this software and associated documentation files (the
 7 # "Software"), to deal in the Software without restriction, including
 8 # without limitation the rights to use, copy, modify, merge, publish,
 9 # distribute, sublicense, and/or sell copies of the Software, and to
10 # permit persons to whom the Software is furnished to do so, subject to
11 # the following conditions:
12 # 
13 # The above copyright notice and this permission notice shall be
14 # included in all copies or substantial portions of the Software.
15 # 
16 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
19 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
20 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
21 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 # 
24 # Retrieved from: http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python)?oldid=4313
25 
26 from decimal import *
27 
28 def exp(x):
29     getcontext().prec += 3
30     if x < 0:
31         s = 1 / exp(abs(x))
32     else:
33         xpower = 1
34         ns = 0
35         s = 1
36         n = 0
37         factorial = 1
38         while s != ns:
39             s = ns
40             term = Decimal(xpower) / factorial
41             ns = s + term
42             xpower *= x
43             n += 1
44             factorial *= n
45     getcontext().prec -= 3
46     return +s
47 
48 def test_exp():
49     assert exp(4) == +Decimal("54.59815003314423907811026120286")
50     assert exp(0) == 1
51     assert exp(-8) == +Decimal("0.0003354626279025118388213891257809")
52     assert exp(Decimal("0.6931471805599453094172321214582")) == 2
53 
54 from math import log as _flog
55 
56 def log(x):
57     if x < 0:
58         return Decimal("NaN")
59     if x == 0:
60         return Decimal("-inf")
61     getcontext().prec += 3
62     eps = Decimal("10")**(-getcontext().prec+2)
63     # A good initial estimate is needed

64     r = Decimal(repr(_flog(float(x))))
65     while 1:
66         r2 = r - 1 + x/exp(r)
67         if abs(r2-r) < eps:
68             break
69         else:
70             r = r2
71     getcontext().prec -= 3
72     return +r
73 
74 def test_log():
75     A = [4, 1, Decimal("0.67753892")]
76     for p in [5, 40, 28]:
77         getcontext().prec = p
78         eps = Decimal("10")**(-getcontext().prec + 2)
79         for a in A:
80             w = exp(log(a))
81             assert abs(w - a) < eps
82 
83 
84 
85 if __name__ == "__main__":
86     test_exp()
87     test_log()
88 
89 


Views
Personal tools