p.lang.pascal.borland:412   
   From: WE@completely.invalid (Wolfgang Ehrhardt)   
      
   On Sun, 16 Jan 2011 20:13:26 +0000, Robert AH Prins   
    wrote:   
   >Two more replies from Norbert Juffa:   
   ...   
   >I got repro. The failing input is 4003_A65AF678_54B28211. With rounding    
   >to positive infinity, multiplying the input by L2E gives 30.0 + 1 ulp.    
   >FRNDINT then rounds this up further to 31.0. Things are already    
   >irreparably screwed up at this point. In any event, the input to F2XM1    
   >becomes -(1.0 + 1 ulp), which is clearly out of bounds. On my CPU it    
   >returns the unmodified input operand in this case.   
      
   Meanwhile I found that the unmodified code has a lot more errors for   
   rmTruncate with negative x arguments: here are the smallest values   
   found (also with hex dumps for better reference), the positive   
   arguments give wrong results with rmUp the negative with rmTruncate.   
      
           x               exp(x)                 x as LSB Hex   
   -6.238324625040 -1.27054942088145051E-0021  C001C7A05AF6CC0968E2   
   -7.624618986159 -3.17637355220362627E-0022  C001F3FCE0F4C07D474D   
   -8.317766166719 -5.29395592033937712E-0023  C002851591F9DD5B9B41   
   -9.010913347279 -7.94093388050906568E-0023  C002902CB3795A7892DC   
      
    6.238324625040 -1.11022302462515654E-0016  4001C7A05AF6CC0968E1   
    7.624618986159 -4.44089209850062616E-0016  4001F3FCE0F4C07D474C   
   12.476649250079 -1.13686837721616030E-0013  4002C7A05AF6CC0968E1   
   15.249237972319 -1.81898940354585648E-0012  4002F3FCE0F4C07D474C   
      
   >The only reasonable workaround for this that I can see is to reduce the    
   >input to F2XM1 further, to make sure we are definitely inside the    
   >interval supported by the hardware. In other words, replace F2MX1 by the    
   >following instruction sequence:   
   >   
   >  fmul  dword ptr [half]   ; 0.5*frac(z) int(z)   
   >  f2mx1                    ; 2^(0.5*frac(z))-1 int(z)   
   >  fld   st(0)              ; 2^(0.5*frac(z))-1 2^(0.5*frac(z))-1 int(z)   
   >  fadd  dword ptr [two]    ; 2^(0.5*frac(z))+1 2^(0.5*frac(z))-1 int(z)   
   >  fmulp st(1),st           ; 2^frac(z)-1 int(z)   
   >   
   This code does not show the errors and seems to work ok (Note that   
   there is a small type f2mx1 must be changed to f2xm1; and the   
   constants should be single). The priced to pay is a small decrease in   
   accuracy, here the result for 20000 random arguments with standard   
   round to nearest (the error is calculated with my MPArith routines):   
      
    ------ Old version ------   
   Test of amath.exp   
   at 10000 random values in [-100.0000000000 .. 100.0000000000]   
   RMS = 0.24, max rel = 0.71 eps at    
    x(ext) =  7.24700653776128467E+0001 = $400590F0AC68BFA87B80   
    y(ext) =  2.97405843003663398E+0031 = $4067BBB0815EA0230490   
    y(mpf) = +2.97405843003663397839071382253654488039888365036E+31   
      
   Test of amath.exp   
   at 10000 random values in [100.0000000000 .. 10000.0000000000]   
   RMS = 0.24, max rel = 0.72 eps at    
    x(ext) =  4.80456952900561126E+0003 = $400B96248E653929CA47   
    y(ext) =  3.96309394729273526E+2086 = $5B12B8A5DD42BAF8E736   
    y(mpf) = +3.96309394729273526538566268155478642708748950161E+2086   
      
    ------ New version ------   
   Test of amath.exp   
   at 10000 random values in [-100.0000000000 .. 100.0000000000]   
   RMS = 0.25, max rel = 0.83 eps at   
    x(ext) =  1.83770128884148676E+0001 = $400393041F554F4C8000   
    y(ext) =  9.57271857222288938E+0007 = $4019B695CA371C7FC4DE   
    y(mpf) = +9.57271857222288937714622943945524368866640537394E+7   
      
   Test of amath.exp   
   at 10000 random values in [100.0000000000 .. 10000.0000000000]   
   RMS = 0.24, max rel = 0.83 eps at   
    x(ext) =  5.98844546289204336E+0003 = $400BBB23904ED947432E   
    y(ext) =  5.60815118133691577E+2600 = $61BEB51753447EF41B4A   
    y(mpf) = +5.60815118133691576592733747663153598713873826468E+2600   
      
   >If there is consensus that Pascal math functions shall work properly    
   >when the dynamic rounding mode is not round-to-nearest-even, I believe    
   >this is the way to go (least impact on execution time).   
   >   
   >I dug out my old Pascal manuals and can't find anything about the    
   >prescribed behavior of math functions when a user dials in some rounding    
   >mode different from the default. Given the naive approach to exp()    
   >computation in Borland's original library, it is quite possible that it    
   >happened to work as desired in round-up mode. My question is: is that to    
   >be considered an implementation artifact or a feature?   
      
   I guess it is a pure accident, especially as Borland Pascal has no   
   library functions for getting/setting rounding modes. The situation is   
   a bit different for Virtual Pascal and Delphi. But IMO at least the   
   Delphi math unit designer/implementers did not understand there task.   
   My AMath unit is an attempt to overcome these Delphi shortcomings.    
      
   Robert: Please forward a big "Thank you" to Norbert.   
      
   Wolfgang   
   --    
   Do not use the e-mail address from the header! You can get a   
   valid address from http://home.netsurf.de/wolfgang.ehrhardt   
   (Free open source Crypto, CRC/Hash, MPArith for Pascal/Delphi)   
      
   --- Internet Rex 2.31   
    * Origin: The gateway at Omicron Theta (1:261/20.999)