precisionarbitrary-precisionmpmath

mpmath and mathmatica discrepancy in output

I need to eventually do very precise calculations with complex numbers, which is why I am trying to use the mp module from mp.math. (Admittedly I have read the mp math documentation, but am new to it) After running into incorrect answers, I traced my issue back to the fact that when I do a calculation on mpmath (see bellow) I get a different result than when I do the same with mathematica. Could someone help me figure out what I am doing wrong?

mp.math: from mpmath import mp mp.dps = 300;

mp.exp(mp.fmul(mp.mpf(str(-64)),mp.mpf(str(.2))))

Output: mpf('0.00000276077257203720079298503881941442964716356235195837374540868148586730700359076564381986686089035467149986363953794158687758910992470828322377833660309685941797430472579604476223424097045454723393291529826485166261076775808555479838036656035924269261878257522607849410841260713806720383331539818973042529')

mathmatica SetPrecision[Exp[-64*.2], 300]

Output: 2.76077257203719864675777080631480231431851279921829700469970703125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^-6

I would expect the results to be the same, but by around the 20th decimal place they start to differ.

Solution

  • In this case I think you appear to have found a bug in Mathematica or not asked the question of it in quite the right way. Its answer for exp(-12.8) doesn't look at all plausible to me with all those trailing zeroes.

    When I throw the calculation at Julia and increase the precision sufficiently it agrees with the result that you obtained all the way when I increased precision to 1000 bits. I'm not sure of the implementation details of Julia's BigFloats but they have never let me down before.

    julia> exp(-BigFloat(64)/BigFloat(5))

    2.760772572037200792985038819414429647163562351958373745 408681485867307003590676e-06

    julia> setprecision(BigFloat, 500)
500

julia> exp(-BigFloat(64)/BigFloat(5))

    2.7607725720372007929850388194144296471635623519583737454086814858673070035907656438198668608903546714998636395379415868775891099247082832237783366030938e-06

    julia> setprecision(BigFloat, 1000)
1000

julia> exp(-BigFloat(64)/BigFloat(5))

    2.76077257203720079298503881941442964716356235195837374540868148586730700359076564381986686089035467149986363953794158687758910992470828322377833660309685941797430472579604476223424097045454723393291529826485166261076775808555479838036656035924269261878257522607849410841260713806720383331539818973042529e-06

    log(BigFloat("2.76077257203720079298503881941442964716356235195837374540868148586730700359076564381986686089035467149986363953794158687758910992470828322377833660309685941797430472579604476223424097045454723393291529826485166261076775808555479838036656035924269261878257522607849410841260713806720383331539818973042529e-06"))

    -12.8000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003

    Taking the log of the Mathematica answer in Julia it is clear that Mathematica has computed exp using standard 64bit double arithmetic. Oops!

    julia> log(BigFloat("2.76077257203719864675777080631480231431851279921829700469970703125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-6"))

    -12.8000000000000007774009673058214261052537582953257100770542758344316669794447948750058707481330343868028670327576657852816894936980836750044670968536044806438802211617596023783929373187330881370009468791968179640611989465169433207621738713056795607656520114787722487328419418068321396252206120875122593

    I think your answers from mpmath are in fact entirely correct and Mathematica is giving you misleading and incorrect answers (looking suspiciously like 64bit double 1/2^N truncated). PS Sorry about the formatting above. I don't seem to be able to indent the very long output lines (and breaking them up is too error prone)