292 lines
10 KiB
Python
292 lines
10 KiB
Python
|
"""
|
||
|
|
Check that the output from irrational functions is accurate for
|
||
|
|
high-precision input, from 5 to 200 digits. The reference values were
|
||
|
|
verified with Mathematica.
|
||
|
|
"""
|
||
|
|
|
||
|
|
import time
|
||
|
|
from mpmath import *
|
||
|
|
|
||
|
|
precs = [5, 15, 28, 35, 57, 80, 100, 150, 200]
|
||
|
|
|
||
|
|
# sqrt(3) + pi/2
|
||
|
|
a = \
|
||
|
|
"3.302847134363773912758768033145623809041389953497933538543279275605"\
|
||
|
|
"841220051904536395163599428307109666700184672047856353516867399774243594"\
|
||
|
|
"67433521615861420725323528325327484262075464241255915238845599752675"
|
||
|
|
|
||
|
|
# e + 1/euler**2
|
||
|
|
b = \
|
||
|
|
"5.719681166601007617111261398629939965860873957353320734275716220045750"\
|
||
|
|
"31474116300529519620938123730851145473473708966080207482581266469342214"\
|
||
|
|
"824842256999042984813905047895479210702109260221361437411947323431"
|
||
|
|
|
||
|
|
# sqrt(a)
|
||
|
|
sqrt_a = \
|
||
|
|
"1.817373691447021556327498239690365674922395036495564333152483422755"\
|
||
|
|
"144321726165582817927383239308173567921345318453306994746434073691275094"\
|
||
|
|
"484777905906961689902608644112196725896908619756404253109722911487"
|
||
|
|
|
||
|
|
# sqrt(a+b*i).real
|
||
|
|
sqrt_abi_real = \
|
||
|
|
"2.225720098415113027729407777066107959851146508557282707197601407276"\
|
||
|
|
"89160998185797504198062911768240808839104987021515555650875977724230130"\
|
||
|
|
"3584116233925658621288393930286871862273400475179312570274423840384"
|
||
|
|
|
||
|
|
# sqrt(a+b*i).imag
|
||
|
|
sqrt_abi_imag = \
|
||
|
|
"1.2849057639084690902371581529110949983261182430040898147672052833653668"\
|
||
|
|
"0629534491275114877090834296831373498336559849050755848611854282001250"\
|
||
|
|
"1924311019152914021365263161630765255610885489295778894976075186"
|
||
|
|
|
||
|
|
# log(a)
|
||
|
|
log_a = \
|
||
|
|
"1.194784864491089550288313512105715261520511949410072046160598707069"\
|
||
|
|
"4336653155025770546309137440687056366757650909754708302115204338077595203"\
|
||
|
|
"83005773986664564927027147084436553262269459110211221152925732612"
|
||
|
|
|
||
|
|
# log(a+b*i).real
|
||
|
|
log_abi_real = \
|
||
|
|
"1.8877985921697018111624077550443297276844736840853590212962006811663"\
|
||
|
|
"04949387789489704203167470111267581371396245317618589339274243008242708"\
|
||
|
|
"014251531496104028712866224020066439049377679709216784954509456421"
|
||
|
|
|
||
|
|
# log(a+b*i).imag
|
||
|
|
log_abi_imag = \
|
||
|
|
"1.0471204952840802663567714297078763189256357109769672185219334169734948"\
|
||
|
|
"4265809854092437285294686651806426649541504240470168212723133326542181"\
|
||
|
|
"8300136462287639956713914482701017346851009323172531601894918640"
|
||
|
|
|
||
|
|
# exp(a)
|
||
|
|
exp_a = \
|
||
|
|
"27.18994224087168661137253262213293847994194869430518354305430976149"\
|
||
|
|
"382792035050358791398632888885200049857986258414049540376323785711941636"\
|
||
|
|
"100358982497583832083513086941635049329804685212200507288797531143"
|
||
|
|
|
||
|
|
# exp(a+b*i).real
|
||
|
|
exp_abi_real = \
|
||
|
|
"22.98606617170543596386921087657586890620262522816912505151109385026"\
|
||
|
|
"40160179326569526152851983847133513990281518417211964710397233157168852"\
|
||
|
|
"4963130831190142571659948419307628119985383887599493378056639916701"
|
||
|
|
|
||
|
|
# exp(a+b*i).imag
|
||
|
|
exp_abi_imag = \
|
||
|
|
"-14.523557450291489727214750571590272774669907424478129280902375851196283"\
|
||
|
|
"3377162379031724734050088565710975758824441845278120105728824497308303"\
|
||
|
|
"6065619788140201636218705414429933685889542661364184694108251449"
|
||
|
|
|
||
|
|
# a**b
|
||
|
|
pow_a_b = \
|
||
|
|
"928.7025342285568142947391505837660251004990092821305668257284426997"\
|
||
|
|
"361966028275685583421197860603126498884545336686124793155581311527995550"\
|
||
|
|
"580229264427202446131740932666832138634013168125809402143796691154"
|
||
|
|
|
||
|
|
# (a**(a+b*i)).real
|
||
|
|
pow_a_abi_real = \
|
||
|
|
"44.09156071394489511956058111704382592976814280267142206420038656267"\
|
||
|
|
"67707916510652790502399193109819563864568986234654864462095231138500505"\
|
||
|
|
"8197456514795059492120303477512711977915544927440682508821426093455"
|
||
|
|
|
||
|
|
# (a**(a+b*i)).imag
|
||
|
|
pow_a_abi_imag = \
|
||
|
|
"27.069371511573224750478105146737852141664955461266218367212527612279886"\
|
||
|
|
"9322304536553254659049205414427707675802193810711302947536332040474573"\
|
||
|
|
"8166261217563960235014674118610092944307893857862518964990092301"
|
||
|
|
|
||
|
|
# ((a+b*i)**(a+b*i)).real
|
||
|
|
pow_abi_abi_real = \
|
||
|
|
"-0.15171310677859590091001057734676423076527145052787388589334350524"\
|
||
|
|
"8084195882019497779202452975350579073716811284169068082670778986235179"\
|
||
|
|
"0813026562962084477640470612184016755250592698408112493759742219150452"\
|
||
|
|
|
||
|
|
# ((a+b*i)**(a+b*i)).imag
|
||
|
|
pow_abi_abi_imag = \
|
||
|
|
"1.2697592504953448936553147870155987153192995316950583150964099070426"\
|
||
|
|
"4736837932577176947632535475040521749162383347758827307504526525647759"\
|
||
|
|
"97547638617201824468382194146854367480471892602963428122896045019902"
|
||
|
|
|
||
|
|
# sin(a)
|
||
|
|
sin_a = \
|
||
|
|
"-0.16055653857469062740274792907968048154164433772938156243509084009"\
|
||
|
|
"38437090841460493108570147191289893388608611542655654723437248152535114"\
|
||
|
|
"528368009465836614227575701220612124204622383149391870684288862269631"
|
||
|
|
|
||
|
|
# sin(1000*a)
|
||
|
|
sin_1000a = \
|
||
|
|
"-0.85897040577443833776358106803777589664322997794126153477060795801"\
|
||
|
|
"09151695416961724733492511852267067419573754315098042850381158563024337"\
|
||
|
|
"216458577140500488715469780315833217177634490142748614625281171216863"
|
||
|
|
|
||
|
|
# sin(a+b*i)
|
||
|
|
sin_abi_real = \
|
||
|
|
"-24.4696999681556977743346798696005278716053366404081910969773939630"\
|
||
|
|
"7149215135459794473448465734589287491880563183624997435193637389884206"\
|
||
|
|
"02151395451271809790360963144464736839412254746645151672423256977064"
|
||
|
|
|
||
|
|
sin_abi_imag = \
|
||
|
|
"-150.42505378241784671801405965872972765595073690984080160750785565810981"\
|
||
|
|
"8314482499135443827055399655645954830931316357243750839088113122816583"\
|
||
|
|
"7169201254329464271121058839499197583056427233866320456505060735"
|
||
|
|
|
||
|
|
# cos
|
||
|
|
cos_a = \
|
||
|
|
"-0.98702664499035378399332439243967038895709261414476495730788864004"\
|
||
|
|
"05406821549361039745258003422386169330787395654908532996287293003581554"\
|
||
|
|
"257037193284199198069707141161341820684198547572456183525659969145501"
|
||
|
|
|
||
|
|
cos_1000a = \
|
||
|
|
"-0.51202523570982001856195696460663971099692261342827540426136215533"\
|
||
|
|
"52686662667660613179619804463250686852463876088694806607652218586060613"\
|
||
|
|
"951310588158830695735537073667299449753951774916401887657320950496820"
|
||
|
|
|
||
|
|
# tan
|
||
|
|
tan_a = \
|
||
|
|
"0.162666873675188117341401059858835168007137819495998960250142156848"\
|
||
|
|
"639654718809412181543343168174807985559916643549174530459883826451064966"\
|
||
|
|
"7996119428949951351938178809444268785629011625179962457123195557310"
|
||
|
|
|
||
|
|
tan_abi_real = \
|
||
|
|
"6.822696615947538488826586186310162599974827139564433912601918442911"\
|
||
|
|
"1026830824380070400102213741875804368044342309515353631134074491271890"\
|
||
|
|
"467615882710035471686578162073677173148647065131872116479947620E-6"
|
||
|
|
|
||
|
|
tan_abi_imag = \
|
||
|
|
"0.9999795833048243692245661011298447587046967777739649018690797625964167"\
|
||
|
|
"1446419978852235960862841608081413169601038230073129482874832053357571"\
|
||
|
|
"62702259309150715669026865777947502665936317953101462202542168429"
|
||
|
|
|
||
|
|
|
||
|
|
def test_hp():
|
||
|
|
for dps in precs:
|
||
|
|
mp.dps = dps + 8
|
||
|
|
aa = mpf(a)
|
||
|
|
bb = mpf(b)
|
||
|
|
a1000 = 1000*mpf(a)
|
||
|
|
abi = mpc(aa, bb)
|
||
|
|
mp.dps = dps
|
||
|
|
assert (sqrt(3) + pi/2).ae(aa)
|
||
|
|
assert (e + 1/euler**2).ae(bb)
|
||
|
|
|
||
|
|
assert sqrt(aa).ae(mpf(sqrt_a))
|
||
|
|
assert sqrt(abi).ae(mpc(sqrt_abi_real, sqrt_abi_imag))
|
||
|
|
|
||
|
|
assert log(aa).ae(mpf(log_a))
|
||
|
|
assert log(abi).ae(mpc(log_abi_real, log_abi_imag))
|
||
|
|
|
||
|
|
assert exp(aa).ae(mpf(exp_a))
|
||
|
|
assert exp(abi).ae(mpc(exp_abi_real, exp_abi_imag))
|
||
|
|
|
||
|
|
assert (aa**bb).ae(mpf(pow_a_b))
|
||
|
|
assert (aa**abi).ae(mpc(pow_a_abi_real, pow_a_abi_imag))
|
||
|
|
assert (abi**abi).ae(mpc(pow_abi_abi_real, pow_abi_abi_imag))
|
||
|
|
|
||
|
|
assert sin(a).ae(mpf(sin_a))
|
||
|
|
assert sin(a1000).ae(mpf(sin_1000a))
|
||
|
|
assert sin(abi).ae(mpc(sin_abi_real, sin_abi_imag))
|
||
|
|
|
||
|
|
assert cos(a).ae(mpf(cos_a))
|
||
|
|
assert cos(a1000).ae(mpf(cos_1000a))
|
||
|
|
|
||
|
|
assert tan(a).ae(mpf(tan_a))
|
||
|
|
assert tan(abi).ae(mpc(tan_abi_real, tan_abi_imag))
|
||
|
|
|
||
|
|
# check that complex cancellation is avoided so that both
|
||
|
|
# real and imaginary parts have high relative accuracy.
|
||
|
|
# abs_eps should be 0, but has to be set to 1e-205 to pass the
|
||
|
|
# 200-digit case, probably due to slight inaccuracy in the
|
||
|
|
# precomputed input
|
||
|
|
assert (tan(abi).real).ae(mpf(tan_abi_real), abs_eps=1e-205)
|
||
|
|
assert (tan(abi).imag).ae(mpf(tan_abi_imag), abs_eps=1e-205)
|
||
|
|
mp.dps = 460
|
||
|
|
assert str(log(3))[-20:] == '02166121184001409826'
|
||
|
|
mp.dps = 15
|
||
|
|
|
||
|
|
# Since str(a) can differ in the last digit from rounded a, and I want
|
||
|
|
# to compare the last digits of big numbers with the results in Mathematica,
|
||
|
|
# I made this hack to get the last 20 digits of rounded a
|
||
|
|
|
||
|
|
def last_digits(a):
|
||
|
|
r = repr(a)
|
||
|
|
s = str(a)
|
||
|
|
#dps = mp.dps
|
||
|
|
#mp.dps += 3
|
||
|
|
m = 10
|
||
|
|
r = r.replace(s[:-m],'')
|
||
|
|
r = r.replace("mpf('",'').replace("')",'')
|
||
|
|
num0 = 0
|
||
|
|
for c in r:
|
||
|
|
if c == '0':
|
||
|
|
num0 += 1
|
||
|
|
else:
|
||
|
|
break
|
||
|
|
b = float(int(r))/10**(len(r) - m)
|
||
|
|
if b >= 10**m - 0.5: # pragma: no cover
|
||
|
|
raise NotImplementedError
|
||
|
|
n = int(round(b))
|
||
|
|
sn = str(n)
|
||
|
|
s = s[:-m] + '0'*num0 + sn
|
||
|
|
return s[-20:]
|
||
|
|
|
||
|
|
# values checked with Mathematica
|
||
|
|
def test_log_hp():
|
||
|
|
mp.dps = 2000
|
||
|
|
a = mpf(10)**15000/3
|
||
|
|
r = log(a)
|
||
|
|
res = last_digits(r)
|
||
|
|
# Mathematica N[Log[10^15000/3], 2000]
|
||
|
|
# ...7443804441768333470331
|
||
|
|
assert res == '43804441768333470331'
|
||
|
|
|
||
|
|
# see issue 145
|
||
|
|
r = log(mpf(3)/2)
|
||
|
|
# Mathematica N[Log[3/2], 2000]
|
||
|
|
# ...69653749808140753263288
|
||
|
|
res = last_digits(r)
|
||
|
|
assert res == '53749808140753263288'
|
||
|
|
|
||
|
|
mp.dps = 10000
|
||
|
|
r = log(2)
|
||
|
|
res = last_digits(r)
|
||
|
|
# Mathematica N[Log[2], 10000]
|
||
|
|
# ...695615913401856601359655561
|
||
|
|
assert res == '13401856601359655561'
|
||
|
|
r = log(mpf(10)**10/3)
|
||
|
|
res = last_digits(r)
|
||
|
|
# Mathematica N[Log[10^10/3], 10000]
|
||
|
|
# ...587087654020631943060007154
|
||
|
|
assert res == '54020631943060007154', res
|
||
|
|
r = log(mpf(10)**100/3)
|
||
|
|
res = last_digits(r)
|
||
|
|
# Mathematica N[Log[10^100/3], 10000]
|
||
|
|
# ,,,59246336539088351652334666
|
||
|
|
assert res == '36539088351652334666', res
|
||
|
|
mp.dps += 10
|
||
|
|
a = 1 - mpf(1)/10**10
|
||
|
|
mp.dps -= 10
|
||
|
|
r = log(a)
|
||
|
|
res = last_digits(r)
|
||
|
|
# ...3310334360482956137216724048322957404
|
||
|
|
# 372167240483229574038733026370
|
||
|
|
# Mathematica N[Log[1 - 10^-10]*10^10, 10000]
|
||
|
|
# ...60482956137216724048322957404
|
||
|
|
assert res == '37216724048322957404', res
|
||
|
|
mp.dps = 10000
|
||
|
|
mp.dps += 100
|
||
|
|
a = 1 + mpf(1)/10**100
|
||
|
|
mp.dps -= 100
|
||
|
|
|
||
|
|
r = log(a)
|
||
|
|
res = last_digits(+r)
|
||
|
|
# Mathematica N[Log[1 + 10^-100]*10^10, 10030]
|
||
|
|
# ...3994733877377412241546890854692521568292338268273 10^-91
|
||
|
|
assert res == '39947338773774122415', res
|
||
|
|
|
||
|
|
mp.dps = 15
|
||
|
|
|
||
|
|
def test_exp_hp():
|
||
|
|
mp.dps = 4000
|
||
|
|
r = exp(mpf(1)/10)
|
||
|
|
# IntegerPart[N[Exp[1/10] * 10^4000, 4000]]
|
||
|
|
# ...92167105162069688129
|
||
|
|
assert int(r * 10**mp.dps) % 10**20 == 92167105162069688129
|