Download code
From LiteratePrograms
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
