flat assembler
Message board for the users of flat assembler.

Index > Windows > stumped

Author
Thread Post new topic Reply to topic
jack2



Joined: 06 Jul 2008
Posts: 33
jack2 22 Sep 2010, 18:08
Ok just for fun I was translating Alan Miller's Quad Module ( http://jblevins.org/mirror/amiller/ ) to FreeBasic and came across a problem,
in the function exactmul2 there is this line of code
Code:
a1 = (a - t) + t
    

which FreeBasic trnslates to
Code:
  mov eax, dword ptr [a]
      fld qword ptr [eax]
 fsub qword ptr [t]
  fadd qword ptr [t]
  fstp qword ptr [a1]
    

which look OK to me, however it gives the wrong result, the following code works
Code:
    a1 = a - t
    a1 = a1 + t
    

which is translated into
Code:
;##    a1 = a - t
 mov eax, dword ptr [a]
      fld qword ptr [eax]
 fsub qword ptr [t]
  fstp qword ptr [a1]

;##    a1 = a1 + t
       fld qword ptr [t]
   fadd qword ptr [a1]
 fstp qword ptr [a1]
    

my queston is:
why do they give different results?
Post 22 Sep 2010, 18:08
View user's profile Send private message Reply with quote
LocoDelAssembly
Your code has a bug


Joined: 06 May 2005
Posts: 4624
Location: Argentina
LocoDelAssembly 22 Sep 2010, 18:18
Perhaps the internal FPU precision is set to 80-bit? In that case storing and loading from ram would have the effect of rounding and clearing of the least significant bits, causing the results to differ from that of your first code.

PS: BTW, by wrong what do you mean exactly? Do you have an example?
Post 22 Sep 2010, 18:18
View user's profile Send private message Reply with quote
jack2



Joined: 06 Jul 2008
Posts: 33
jack2 22 Sep 2010, 18:38
LocoDelAssembly I can give you an example but it's in FreeBasic, nonetheless here's the code
Code:

'MODULE quadruple_precision
'
' This version is for Compaq Fortran, Lahey/Fujitsu LF95 and
' Absoft ProFortran 6.0.
' N.B. The -o0 option (no optimization) must be used with LF95.
'
' N.B. This is quadruple precision implemented in SOFTWARE, hence it is SLOW.
' This package has NOT been tested with other Fortran 90 compilers.
' The operations +-*/ & sqrt will probably work correctly with other compilers
' for PCs.   Functions and routines which initialize quadruple-precision
' constants, particularly the trigonometric and exponential functions will
' only work if the compiler initializes double-precision constants from
' decimal code in exactly the same way as the Lahey compilers.

' The basic algorithms come from:
' Linnainmma, S.  Software for doubled-precision floating-point computations,
'    ACM Trans, Math. Software, vol. 7, pp.272-283, 1981.

' Programmer : Alan Miller
' e-mail: amiller @ bigpond.net.au   URL:  http://users.bigpond.net.au/amiller

' Latest revision - 18 September 2002

' 4 Aug 97.  Added new algorithm for exponential.  Takes half the time of the
'            previous Taylor series, but errors 2.5 times larger.   The old
'            exponential is included here under the name exp_taylor for anyone
'            who needs the extra accuracy.  To use it instead of the new
'            algorithm, change the module procedure name in INTERFACE exp from
'            longexp to exp_taylor.
' 5 Aug 97.  Found way to reduce cancellation errors in yesterday's algorithm
'            for the exponential.   Removed exp_taylor.
' 18 Aug 97. Added table of quadruple-precision constants.
' 8 Sept 97. Added str_quad which reads a character string and converts it to
'            quadruple-precision form.
' 15 Oct 97. Added quad_str which takes a quadruple-precision value and outputs
'            a character string containing its decimal value to 30 significant
'            digits.
'            Added overlays to the ** operator for quadruple-precision values.
' 15 Jan 98. Added ATAN2.
' 19 Jan 98. Added <, <=, ==, >= and >.
' 27 Dec 98. Added quad/real, quad/integer, integer/quad, real/quad, dp/quad,
'            int+quad, real+quad, dp+quad, int-quad, real-quad, dp-quad,
'            SUM, DOT_PRODUCT & MATMUL.
' 29 Dec 98. Correction to routine string_quad for strings of < 5 characters.
' 10 May 99. Added EPSILON for quad type.
' 5 Oct 99.  log10e corrected.
' 15 Oct 99. Corrected function quad_pow_int.
' 18 Oct 99. Rewrote quad_pow_int to use the binary power method.
' 11 Nov 99. Added overlaying of assignments, e.g. quad = int, etc.
' 21 Jan 00. Added inreface for EPSILON.
' 17 Feb 02. Corrected ACOS & ASIN for arguments close to +1 or -1.
' 18 Feb 02. Further improvement to ACOS.
' 18 Sep 02. Added reference to Linnainmaa in comments at beginning.


#include "string.bi"

'dp = 8 = double
Dim Shared As Integer nbits = 53
Dim Shared As Double constant = 134217729.0 '=2^(nbits - nbits\2) + 1
Dim Shared As Double zero = 0.0
Dim Shared As Double one  = 1.0

Type quad
  As Double hi, lo
End Type

Declare Function longadd(Byref a As quad, Byref c As quad) As quad

Function cquad(Byref a As Double, Byref b As Double) As quad
    Dim As quad c
    c.hi=a
    c.lo=b
    Return c
End Function

Dim Shared As quad _
             pi =  ( 0.3141592653589793D+01, 0.1224646799147353D-15 ),  _
          piby2 =  ( 0.1570796326794897D+01, -.3828568698926950D-15 ),  _
          piby3 =  ( 0.1047197551196598D+01, -.3292527815701405D-15 ),  _
          piby4 =  ( 0.7853981633974484D+00, -.8040613248383182D-16 ),  _
          piby6 =  ( 0.5235987755982990D+00, -.1646263907850702D-15 ),  _
          twopi =  ( 0.6283185307179586D+01, 0.2449293598294707D-15 ),  _
          ln_pi =  ( 0.1144729885849400D+01, 0.2323105560877391D-15 ),  _
         sqrtpi =  ( 0.1772453850905516D+01, -.7666586499825800D-16 ),  _
       fact_pt5 =  ( 0.8862269254527582D+00, -.1493552349616447D-15 ),  _
        sqrt2pi =  ( 0.2506628274631001D+01, -.6273750096546544D-15 ),  _
      lnsqrt2pi =  ( 0.9189385332046728D+00, -.3878294158067242D-16 ),  _
      one_on2pi =  ( 0.1591549430918953D+00, 0.4567181289366658D-16 ),  _
    two_on_rtpi =  ( 0.1128379167095513D+01, -.4287537502368968D-15 ),  _
        deg2rad =  ( 0.1745329251994330D-01, -.3174581724866598D-17 ),  _
        rad2deg =  ( 0.5729577951308232D+02, -.1987849567057628D-14 ),  _
            ln2 =  ( 0.6931471805599454D+00, -.8783183432405266D-16 ),  _
           ln10 =  ( 0.2302585092994046D+01, -.2170756223382249D-15 ),  _
          log2e =  ( 0.1442695040888964D+01, -.6457785410341630D-15 ),  _
         log10e =  ( 0.4342944819032519D+00, -.5552037773430574D-16 ),  _
        log2_10 =  ( 0.3321928094887362D+01, 0.1661617516973592D-15 ),  _
        log10_2 =  ( 0.3010299956639812D+00, -.2803728127785171D-17 ),  _
          euler =  ( 0.5772156649015330D+00, -.1159652176149463D-15 ),  _
              e =  ( 0.2718281828459045D+01, 0.1445646891729250D-15 ),  _
          sqrt2 =  ( 0.1414213562373095D+01, 0.1253716717905022D-15 ),  _
          sqrt3 =  ( 0.1732050807568877D+01, 0.3223954471431004D-15 ),  _
         sqrt10 =  ( 0.3162277660168380D+01, -.6348773795572286D-15 )


Function exactmul2(Byref a As Double, Byref c As Double) As quad
    '  Procedure exactmul2, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As quad ac
    Dim As Double a1, a2, c1, c2, t

    t = constant * a

    ' a1 = (a - t) + t ' FreeBasic's optimization removes the brackets
    ' and sets a1 = a which defeats the whole point.
    ' splitting the expression makes it work
    a1 = a - t
    a1 = a1 + t
'asm
''##    a1 = a - t
'    mov eax, dword ptr [a]
'        fld qword ptr [eax]
'       fsub qword ptr [t]
'        fstp qword ptr [a1]
'
''##    a1 = a1 + t
'        fld qword ptr [t]
' fadd qword ptr [a1]
'       fstp qword ptr [a1]
'end asm
'asm
''##     a1 = (a - t) + t
'    mov eax, dword ptr [a]
'  fld qword ptr [eax]
'       fsub qword ptr [t]
'        fadd qword ptr [t]
'        fstp qword ptr [a1]
'end asm
'    asm
'        mov       eax, DWORD PTR [a]               
'        movsd     xmm0, QWORD PTR [eax]            
'        movsd     xmm1, QWORD PTR [t]              
'        subsd     xmm0, xmm1                       
'        movsd     xmm1, QWORD PTR [t]              
'        addsd     xmm0, xmm1                       
'        movsd     QWORD PTR [a1], xmm0             
'    end asm
    
    a2 = a - a1
    t = constant * c
    c1 = (c - t) + t
    c2 = c - c1
    ac.hi = a * c
    ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2

    Return ac
End Function 'exactmul2


Function longmul(Byref a As quad, Byref c As quad) As quad
    '  Procedure longmul, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As Double zz
    Dim As quad z, ac

    z = exactmul2(a.hi, c.hi)
    zz = ((a.hi + a.lo) * c.lo + a.lo * c.hi) + z.lo
    ac.hi = z.hi + zz
    ac.lo = (z.hi - ac.hi) + zz

    Return ac
End Function 'longmul



Function mult_quad_int(Byref a As quad, Byval i As Integer) As quad
    '     Multiply quadruple-precision number (a) by an integer (i).

    Dim As quad c

    If (i = 0) Then
      c = cquad(zero, zero)
    Elseif (i = 1) Then
      c = a
    Elseif (i = -1) Then
      c.hi = -a.hi
      c.lo = -a.lo
    Else
      c = longadd(exactmul2(a.hi, Cdbl(i)) , exactmul2(a.lo, Cdbl(i)))
    End If

    Return c
End Function 'mult_quad_int



Function mult_int_quad(Byval i As Integer, Byref a As quad) As quad
    '     Multiply quadruple-precision number (a) by an integer (i).

    Dim As quad c

    If (i = 0) Then
      c = cquad(zero, zero)
    Elseif (i = 1) Then
      c = a
    Elseif (i = -1) Then
      c.hi = -a.hi
      c.lo = -a.lo
    Else
      c = longadd(exactmul2(a.hi, Cdbl(i)) , exactmul2(a.lo, Cdbl(i)))
    End If

    Return c
End Function 'mult_int_quad



Function mult_quad_dp(Byref a As quad, Byref b As Double) As quad
    '  Multiply a quadruple-precision number (a) by a double-precision number (b).

    '     Local variables
    Dim As quad z, c
    Dim As Double zz

    z = exactmul2(a.hi, b)
    zz = a.lo * b + z.lo
    c.hi = z.hi + zz
    c.lo = (z.hi - c.hi) + zz

    Return c
End Function 'mult_quad_dp



Function mult_quad_real(Byref a As quad, Byref b As Double) As quad
    '  Multiply a quadruple-precision number (a) by a real number (b).

    '     Local variables
    Dim As quad z, c
    Dim As Double zz

    z = exactmul2(a.hi, b)
    zz = a.lo * b + z.lo
    c.hi = z.hi + zz
    c.lo = (z.hi - c.hi) + zz

    Return c
End Function 'mult_quad_real



Function mult_dp_quad(Byref b As Double, Byref a As quad) As quad
    '  Multiply a quadruple-precision number (a) by a double-precision number (b).

    '     Local variables
    Dim As quad z, c
    Dim As Double zz

    z = exactmul2(a.hi, b)
    zz = a.lo * b + z.lo
    c.hi = z.hi + zz
    c.lo = (z.hi - c.hi) + zz

    Return c
End Function 'mult_dp_quad



Function mult_real_quad(Byref a As Double, Byref b As quad) As quad
    '  Multiply a quadruple-precision number (a) by a single-precision number (b).

    '     Local variables
    Dim As quad z, c
    Dim As Double zz

    z = exactmul2(b.hi, a)
    zz = b.lo * a + z.lo
    c.hi = z.hi + zz
    c.lo = (z.hi - c.hi) + zz

    Return c
End Function 'mult_real_quad



Function longdiv(Byref a As quad, Byref c As quad) As quad
    '  Procedure longdiv, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, ac

    z = a.hi / c.hi
    q = exactmul2(c.hi, z)
    zz = ((((a.hi - q.hi) - q.lo) + a.lo) - z*c.lo) / (c.hi + c.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'longdiv



Function div_quad_dp(Byref a As quad, Byref b As Double) As quad
    '     Divide a quadruple-precision number (a) by a double-precision number (b)

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, c

    z = a.hi / b
    q = exactmul2(b, z)
    zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
    c.hi = z + zz
    c.lo = (z - c.hi) + zz

    Return c
End Function 'div_quad_dp



Function div_quad_real(Byref a As quad, Byref b As Double) As quad
    'Divide a quadruple-precision number (a) by a single-precision number (b)

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, c

    z = a.hi / b
    q = exactmul2(b, z)
    zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
    c.hi = z + zz
    c.lo = (z - c.hi) + zz

    Return c
End Function 'div_quad_real



Function div_quad_int(Byref a As quad, Byval b As Integer) As quad
    'Divide a quadruple-precision number (a) by an integer (b)

    'Local variables
    Dim As Double z, zz
    Dim As quad q, c

    z = a.hi / b
    q = exactmul2(Cdbl(b), z)
    zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
    c.hi = z + zz
    c.lo = (z - c.hi) + zz

    Return c
End Function 'div_quad_int



Function div_int_quad(Byval a As Integer, Byref c As quad) As quad
    '  Procedure longdiv, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, ac

    z = Cdbl(a) / c.hi
    q = exactmul2(c.hi, z)
    zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'div_int_quad



Function div_real_quad(Byref a As Double, Byref c As quad) As quad
    '  Procedure longdiv, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, ac

    z = a / c.hi
    q = exactmul2(c.hi, z)
    zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'div_real_quad



Function div_dp_quad(Byref a As Double, Byref c As quad) As quad
    '  Procedure longdiv, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As Double z, zz
    Dim As quad q, ac

    z = a / c.hi
    q = exactmul2(c.hi, z)
    zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'div_dp_quad



Function longadd(Byref a As quad, Byref c As quad) As quad
    '  Procedure longadd, translated from Pascal, from:
    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    '     Local variables
    Dim As quad ac
    Dim As Double z, q, zz

    z = a.hi + c.hi
    q = a.hi - z
    zz = (((q + c.hi) + (a.hi - (q + z))) + a.lo) + c.lo
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'longadd



Function quad_add_int(Byref a As quad, Byval c As Integer) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac
    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_add_int



Function quad_add_real(Byref a As quad, Byref c As Double) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_add_real



Function quad_add_dp(Byref a As quad, Byref c As Double) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_add_dp



Function int_add_quad(Byval c As Integer, Byref a As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'int_add_quad



Function real_add_quad(Byref c As Double, Byref a As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'real_add_quad



Function dp_add_quad(Byref c As Double, Byref a As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi + c
    q = a.hi - z
    zz = (((q + c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'dp_add_quad



'FUNCTION longsub(byref a as quad, byref c as quad) as quad
'    '  Adapted from longadd by changing signs of c.
'    '  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
'    '  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.
'
'    '     Local variables
'    dim as double z, q, zz
'    dim as quad ac
'
'    z = a.hi - c.hi
'    q = a.hi - z
'    zz = (((q - c.hi) + (a.hi - (q + z))) + a.lo) - c.lo
'    ac.hi = z + zz
'    ac.lo = (z - ac.hi) + zz
'
'    RETURN ac
'END FUNCTION 'longsub

Function longsub(Byref a As quad, Byref c As quad) As quad
'  Adapted from longadd by changing signs of c.
'  Linnainmaa, Seppo (1981).   Software for doubled-precision floating-point
'  computations.   ACM Trans. on Math. Software (TOMS), 7, 272-283.

    Dim As quad ac

    '     Local variables
    Dim As Double z, q, zz, temp

'    asm
 z = a.hi - c.hi

'        mov       eax, DWORD PTR [a]                    
'        mov       edx, DWORD PTR [c]                    
'        movsd     xmm0, QWORD PTR [eax]         
'        movsd     xmm1, QWORD PTR [edx]         
'        subsd     xmm0, xmm1                                 
'        movsd     QWORD PTR [z], xmm0

 q = a.hi - z

'        mov       eax, DWORD PTR [a]                    
'        movsd     xmm0, QWORD PTR [eax]         
'        movsd     xmm1, QWORD PTR [z]              
'        subsd     xmm0, xmm1                                 
'        movsd     QWORD PTR [q], xmm0             

' zz = (((q - c.hi) + (a.hi - (q + z))) + a.lo) - c.lo
'line above rewritten in simple sub-expressions to make it work in FB
temp = q + z
temp = a.hi - temp
zz = q - c.hi
zz = zz + temp
zz = zz + a.lo
zz = zz - c.lo
 
'asm
'        mov       eax, DWORD Ptr [c]                    
'        movsd     xmm0, QWORD Ptr [q]             
'        movsd     xmm1, QWORD Ptr [eax]         
'        subsd     xmm0, xmm1                                 
'        mov       eax, DWORD Ptr [a]                    
'        movsd     xmm1, QWORD Ptr [q]             
'        movsd     xmm2, QWORD Ptr [z]              
'        addsd     xmm1, xmm2                                 
'        movsd     xmm2, QWORD Ptr [eax]         
'        subsd     xmm2, xmm1                                 
'        addsd     xmm0, xmm2                                 
'        mov       eax, DWORD Ptr [a]                    
'        movsd     xmm1, QWORD Ptr [8+eax]    
'        addsd     xmm0, xmm1                                 
'        mov       eax, DWORD Ptr [c]                    
'        movsd     xmm1, QWORD Ptr [8+eax]    
'        subsd     xmm0, xmm1                                 
'        movsd     QWORD Ptr [zz], xmm0            
'End asm
 ac.hi = z + zz

'        movsd     xmm0, QWORD PTR [z]              
'        movsd     xmm1, QWORD PTR [zz]            
'        addsd     xmm0, xmm1                                 
'        lea       eax, DWORD PTR [ac]                    
'        movsd     QWORD PTR [eax], xmm0         

 ac.lo = (z - ac.hi) + zz

'        lea       eax, DWORD PTR [ac]                    
'        movsd     xmm0, QWORD PTR [z]              
'        movsd     xmm1, QWORD PTR [eax]         
'        subsd     xmm0, xmm1                                 
'        movsd     xmm1, QWORD PTR [zz]            
'        addsd     xmm0, xmm1                                 
'        lea       eax, DWORD PTR [ac]                    
'        movsd     QWORD PTR [8+eax], xmm0    
'
'    end asm

    Return ac
End Function



Function quad_sub_int(Byref a As quad, Byval c As Integer) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi - c
    q = a.hi - z
    zz = (((q - c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_sub_int



Function quad_sub_real(Byref a As quad, Byval c As Double) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi - c
    q = a.hi - z
    zz = (((q - c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_sub_real



Function quad_sub_dp(Byref a As quad, Byref c As Double) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a.hi - c
    q = a.hi - z
    zz = (((q - c) + (a.hi - (q + z))) + a.lo)
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'quad_sub_dp



Function int_sub_quad(Byval a As Integer, Byref c As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = Cdbl(a) - c.hi
    q = Cdbl(a) - z
    zz = ((q - c.hi) + (Cdbl(a) - (q + z))) - c.lo
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'int_sub_quad



Function real_sub_quad(Byref a As Double, Byref c As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a - c.hi
    q = a - z
    zz = ((q - c.hi) + (a - (q + z))) - c.lo
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'real_sub_quad



Function dp_sub_quad(Byref a As Double, Byref c As quad) As quad

    '     Local variables
    Dim As Double z, q, zz
    Dim As quad ac

    z = a - c.hi
    q = a - z
    zz = ((q - c.hi) + (a - (q + z))) - c.lo
    ac.hi = z + zz
    ac.lo = (z - ac.hi) + zz

    Return ac
End Function 'dp_sub_quad



Function negate_quad(Byref a As quad) As quad
    '     Change the sign of a quadruple-precision number.
    '     In many cases, a & b will occupy the same locations.

    Dim As quad b

    b.hi = -a.hi
    b.lo = -a.lo

    Return b
End Function 'negate_quad


'##############################################################################
'##############################################################################

Function nint(Byref x As Double) As Double
    If x<0 Then
        Return Fix(x-0.5)
    Else
        Return Fix(x+0.5)
    End If
End Function

Function pow(Byref x As quad, Byval e As Integer) As quad
    'take x to an integer power
    Dim As quad z, y = x
    Dim As Integer n
    n = Abs(e)
    z=cquad(1.0#, 0.0#)
    While n>0
      While (n Mod 2)=0
        n=n\2
        y=longmul(y,y)
      Wend
      n=n-1
      z=longmul(z,y)
    Wend
    If e<0 Then
        y=cquad(1.0,0.0)
        z=longdiv(y,z)
    End If
    Return z
End Function

'Function qscale(byref x as quad, byval j as integer) as quad
'    dim as quad y
'    Asm
'        mov eax,DWORD PTR [x]
'        lea edx,DWORD PTR [y]
'        fild dword ptr [j]
'        fld qword ptr [eax]   'x
'        fscale
'        fstp qword ptr [edx]  'y
'        fstp st(0) 'clean fpu stack
'        fild dword ptr [j]
'        fld qword ptr [eax+8]   'x
'        fscale
'        fstp qword ptr [edx+8]  'y
'        fstp st(0) 'clean fpu stack
'    End Asm
'    return y
'End Function

Function qscale(Byref a As quad, Byval i As Integer ) As quad
    '     Multiply a by 2^i

    Dim As quad b
    b=cquad(2,0)
    b=pow(b,i)
    b=longmul(a,b)
    'b.hi = a.hi * 2.0^i 'SCALE(a.hi, i)
    'b.lo = a.lo * 2.0^i 'SCALE(a.lo, i)

    Return b
End Function

'FUNCTION qscale(byref a as quad, byval i as integer ) as quad
'    '     Multiply a by 2^i
'
'    dim as quad b
'
'    b.hi = a.hi * 2.0^i 'SCALE(a.hi, i)
'    b.lo = a.lo * 2.0^i 'SCALE(a.lo, i)
'
'    RETURN b
'END FUNCTION

Function longexp(Byref x As quad) As quad
'  Calculate a quadruple-precision exponential
'  Method:
'   x    x.log2(e)    nint[x.log2(e)] + frac[x.log2(e)]
'  e  = 2          = 2
'
'                     iy    fy
'                  = 2   . 2
'  Then
'   fy    y.ln(2)
'  2   = e
'
'  Now y.ln(2) will be less than 0.3466 in absolute value.
'  This is halved and a Pade approximation is used to approximate e^x over
'  the region (-0.1733, +0.1733).   This approximation is then squared.

'  WARNING: No overflow checks'


' Local variables
Dim As quad temp, ysq, sum1, sum2, y
Dim As Integer iy

y = longdiv(x , ln2)
iy = NINT(y.hi)
y = longmul(quad_sub_dp(y , Cdbl(iy)) , ln2)
y = qSCALE(y, -1)

' The Pade series is:
'     p0 + p1.y + p2.y^2 + p3.y^3 + ... + p9.y^9
'     ------------------------------------------
'     p0 - p1.y + p2.y^2 - p3.y^3 + ... - p9.y^9
'
' sum1 is the sum of the odd powers, sum2 is the sum of the even powers

ysq = longmul(y , y)
sum1 = quad_add_dp(ysq , 3960)
sum1 = longmul(sum1,ysq)
sum1 = quad_add_dp(sum1,2162160)
sum1 = longmul(sum1,ysq)
sum1 = quad_add_dp(sum1,302702400)
sum1 = longmul(sum1,ysq)
sum1 = quad_add_dp(sum1,8821612800)
sum1 = longmul(y,sum1)
'sum1 = y * ((((ysq + 3960.)*ysq + 2162160._dp)*ysq + 302702400._dp)*ysq +   _
'               8821612800._dp)
sum2 = mult_int_quad(90,ysq)
sum2 = quad_add_dp(sum2,110880)
sum2 = longmul(sum2,ysq)
sum2 = quad_add_dp(sum2,30270240)
sum2 = longmul(sum2,ysq)
sum2 = quad_add_dp(sum2,2075673600)
sum2 = longmul(sum2,ysq)
sum2 = quad_add_dp(sum2,17643225600)
'sum2 = (((90.*ysq + 110880.)*ysq + 30270240._dp)*ysq + 2075673600._dp)*ysq +  _
'       17643225600._dp

'                     sum2 + sum1         2.sum1
' Now approximation = ----------- = 1 + ----------- = 1 + 2.temp
'                     sum2 - sum1       sum2 - sum1
'
' Then (1 + 2.temp)^2 = 4.temp.(1 + temp) + 1

temp = longsub(sum2,sum1)
temp = longdiv(sum1 ,temp)
y = quad_add_dp(temp,1)
y = longmul(y,temp)
y = qSCALE(y, 2)
y = quad_add_dp(y,1)
y = qSCALE(y, iy)

Return y
End Function 'longexp

Function longlog(Byref x As quad) As quad
    '  Quadruple-precision logarithm to base e
    '  Halley's algorithm using double-precision logarithm as starting value.
    '  Solves:         y
    '          f(y) = e  - x = 0

    'TYPE (quad), INTENT(IN) :: x
    Dim As quad y

    '     Local variables
    Dim As quad expy, f

    y.hi = Log(x.hi)
    y.lo = 0.0
    expy = longEXP(y)
    f = longsub(expy , x)
    f = qSCALE(f, 1)
    expy = longadd(expy , x)
    expy = longdiv(f , expy)
    y = longsub(y , expy)

    Return y
End Function


Function quad_string(Byref value As quad, Byval ier As Integer) As String
    ' Convert a quadruple-precision quantity to a decimal character string.
    ' Error indicator ier = 0 if conversion OK
    '                     = 1 if the length of the string < 36 characters.

    ' Local variables
    Dim As Zstring * 2  sign
    Dim As Zstring * 18 str1, str2
    Dim As String strng
    Dim As quad vl, v
    Dim As Integer dec_expnt, i
    Dim As Double tmp

    'IF (LEN(strng) < 36) THEN
    '  ier = 1
    '  strng = "***"
    '  RETURN strng
    'END IF
    ier = 0

    ' Check if value = zero.
    If (value.hi = zero) Then
      strng = " 0.00"
      Return strng
    End If

    If (value.hi < 0) Then
      sign = "-"
      vl.hi = - value.hi
      vl.lo = - value.lo
    Else
      sign = " "
      vl = value
    End If

    ' Use LOG10 to set the exponent.
    dec_expnt = Int( Log(vl.hi)*0.4342944819032518 )

    ' Get the first 15 decimal digits
    If (dec_expnt <> 14) Then
        v=longLOG(cquad(10.0, zero))
        v=mult_quad_dp(v,Cdbl(14 - dec_expnt))
        v=longexp(v)
'*************
'HACK
v.lo = 0.0
'*************
        vl = longmul(vl , v )
    End If

    str1=format( vl.hi,"################")
    ' Calculate the remainder

    tmp =Val(str1+".0#")
    vl = quad_sub_real(vl , tmp)
    ' If vl is -ve, subtract 1 from the last digit of str1, and add 1 to vl.
    If (vl.hi < -0.5D-16) Then
      tmp = tmp - one
      str1=format( tmp,"################")
      vl = quad_add_dp(vl , one)
    End If

    vl = mult_quad_dp(vl , 1.D15)

    ' write the second 15 digits
    If vl.hi<1 Then
        str2 = "000000000000000"
    Else
        str2=format( vl.hi,"################")
    End If
    strng = sign+Left(str1,1)+"."+Mid(str1,2)+str2
    If dec_expnt>=0 Then
        strng=strng+"e+"+Str(dec_expnt)
    Else
        strng=strng+"e"+Str(dec_expnt)
    End If

'    ' If str2 consists of asterisks, add 1 in the last digit to str1.
'    ' Set str2 to zeroes.
'    IF (str2(2:2) = '*') THEN
'      tmp = tmp + one
'      WRITE(str1, '(f17.0)') tmp
'      IF (str1(1:1) /= ' ') THEN
'        dec_expnt = dec_expnt + 1
'      ELSE
'        str1 = str1(2:17)
'      END IF
'      str2 = '000000000000000.'
'    END IF
'
'    ' Replace leading blanks with zeroes
'    DO i = 1, 15
'      IF (str2(i:i) /= ' ') EXIT
'      str2(i:i) = '0'
'    END DO
'
'    ' Combine str1 & str2, removing decimal points & adding exponent.
'    i = INDEX(str1, '.')
'    str1(i:i) = ' '
'    str2(16:16) = ' '
'    strng = '.' // TRIM(ADJUSTL(str1)) // TRIM(ADJUSTL(str2)) // 'E'
'    WRITE(str1, '(i4.2)') dec_expnt+1
'    strng = TRIM(strng) // ADJUSTL(str1)
'
'    ' Restore the sign.
'    IF (sign = '-') THEN
'      strng = '-' // ADJUSTL(strng)
'    ELSE
'      strng = ADJUSTL(strng)
'    END IF

    Return strng
End Function

'##############################################################################

Function string_quad(Byref value As String) As quad
  Dim As quad qd
  Dim As Integer j, s, d, e, ep, ex, es, i, f, fp, fln
  Dim As String c, f1, f2, f3, vl
  j=1
  s=1
  d=0
  e=0
  ep=0
  ex=0
  es=1
  i=0
  f=0
  fp=0
  f1=""
  f2=""
  f3=""
  vl=Ucase(value) 
  fln=Len(vl) 
  
  While j<=fln 
    c=Mid(vl,j,1) 
    If ep=1 Then
      If c=" " Then
        Goto atof1nxtch 
      Endif 
      If c="-" Then
        es=-es 
        c="" 
      Endif 
      If c="+" Then
        Goto atof1nxtch 
      Endif 
      If (c="0") And (f3="") Then
        Goto atof1nxtch 
      Endif 
      If (c>"/") And (c<":") Then 'c is digit between 0 and 9 
        f3=f3+c 
        ex=10*ex+(Asc(c)-48) 
        Goto atof1nxtch 
      Endif 
    Endif 
    
    If c=" " Then
      Goto atof1nxtch 
    Endif 
    If c="-" Then
      s=-s 
      Goto atof1nxtch 
    Endif 
    If c="+" Then
      Goto atof1nxtch 
    Endif 
    If c="." Then
      If d=1 Then
        Goto atof1nxtch 
      Endif 
      d=1 
    Endif 
    If (c>"/") And (c<":") Then 'c is digit between 0 and 9 
      If ((c="0") And (i=0)) Then
        If d=0 Then
          Goto atof1nxtch 
        Endif 
        If (d=1) And (f=0) Then
          e=e-1 
          Goto atof1nxtch 
        Endif 
      Endif 
      If d=0 Then
        f1=f1+c 
        i=i+1 
      Else 
        If (c>"0") Then
          fp=1 
        Endif 
        f2=f2+c 
        f=f+1 
      Endif 
    Endif 
    If c="E" Then
      ep=1 
    Endif 
    atof1nxtch: 
    j=j+1 
  Wend 
  If fp=0 Then
    f=0 
    f2="" 
  Endif 
  
  ex=es*ex-1+i+e '(es*ex)-(58-i)+e 
  f1=f1+f2
  'f1=Mid(f1,1,1)+"."+Right(f1,Len(f1)-1) 
'  If s<0 then
'    f1="-"+f1
'  Else
'    f1=" "+f1
'  EndIf
  fln=Len(f1) 
  If Len(f1)>30 Then
    f1=Mid(f1,1,30) 
  Endif 
  While Len(f1)<30 
    f1=f1+"0" 
  Wend 
  f2=Str(Abs(ex))
  f2=String(4-Len(f2),"0")+f2
  If ex<0 Then
    f2="E-"+f2 
  Else 
    f2="E+"+f2 
  Endif 
  f2=Left(f1,15)
  f3=Right(f1,15)
  qd.hi=Val(f2)
  qd.lo=0
  qd=mult_quad_dp(qd, 1000000000000000)
  qd=quad_add_dp(qd, Val(f3))
  qd=div_quad_dp(qd,1000000000000000)
  qd=div_quad_dp(qd,1000000000000000)
'  f1=f1+f2
'  fb.fs=f1
'  n.bcd(0) = 0
'  If fb.fc(0)=asc("-") then
'    n.bcd(0)=&h80
'  EndIf
'  If fb.fc(61)=asc("-") then
'    n.bcd(0)=n.bcd(0) or &h40
'  EndIf
'  n.bcd(1)=((fb.fc(62)-asc("0")) Shl 4) or (fb.fc(63)-asc("0"))
'  n.bcd(2)=((fb.fc(64)-asc("0")) Shl 4) or (fb.fc(65)-asc("0"))
'  n.bcd(3)=((fb.fc(1)-asc("0")) Shl 4) or (fb.fc(3)-asc("0"))
'  For i=4 To 31
'    n.bcd(i)=((fb.fc(2*(i-4)+4)-asc("0")) Shl 4) or (fb.fc(2*(i-4)+5)-asc("0"))
'  Next
  Return qd
End Function

'##############################################################################

Dim As quad x, y
Dim As String s
Dim As Integer i, j
Dim As Double t, t1, dy, dx = 3.1415926535
x=cquad(2,0)
'y=pow(x,3)
y=longlog(x)
'y=longexp(ln2)
?quad_string(y, i)
?quad_string(ln2, i)
Sleep
     

in the function exactmul2 I worked around the problem by spliting up a1 = (a - t) + t into a1 = a - t and a1 = a1 + t but am at a loss as to why the first expression does not work and the second one does.
you can see my post at the FreeBasic Forum here http://www.freebasic.net/forum/viewtopic.php?t=16450
Post 22 Sep 2010, 18:38
View user's profile Send private message Reply with quote
jack2



Joined: 06 Jul 2008
Posts: 33
jack2 22 Sep 2010, 19:57
OK, so how do I set the FPU control word to Double Precision?
because if I change
Code:
        mov eax, dword ptr [a] 
        fld qword ptr [eax] 
        fsub qword ptr [t] 
        fadd qword ptr [t] 
        fstp qword ptr [a1] 
    

to
Code:
'##     a1 = (a - t) + t
  mov eax, dword ptr [a]
      fld qword ptr [eax]
 fsub qword ptr [t]
  fstp qword ptr [a1]
 fld qword ptr [a1]
  fadd qword ptr [t]
  fstp qword ptr [a1]
    

it works as expected, this is weird,
Post 22 Sep 2010, 19:57
View user's profile Send private message Reply with quote
baldr



Joined: 19 Mar 2008
Posts: 1651
baldr 22 Sep 2010, 20:13
jack2,

Store/load operations effectively limit precision to 53 bits. fldcw instruction can be used to load x87 Control Word.
Post 22 Sep 2010, 20:13
View user's profile Send private message Reply with quote
jack2



Joined: 06 Jul 2008
Posts: 33
jack2 22 Sep 2010, 20:37
yes, I know that much, but what is the value I need to store in the controlword to force it to double?
Post 22 Sep 2010, 20:37
View user's profile Send private message Reply with quote
baldr



Joined: 19 Mar 2008
Posts: 1651
baldr 22 Sep 2010, 20:47
jack2,

That depends. You may change only precision though: fstcw can be useful.
Post 22 Sep 2010, 20:47
View user's profile Send private message Reply with quote
LocoDelAssembly
Your code has a bug


Joined: 06 May 2005
Posts: 4624
Location: Argentina
LocoDelAssembly 22 Sep 2010, 22:20
0x027F and 0x037F are WinXP and finit defaults respectively.

Here you can see a quick control word reference.

PS: In case isn't obvious, WinXP's default is double precision and finit is extended precision.
Post 22 Sep 2010, 22:20
View user's profile Send private message Reply with quote
jack2



Joined: 06 Jul 2008
Posts: 33
jack2 23 Sep 2010, 00:20
thanks, I examined the control word just before doing the arithmetic and it was 0x37F.
just did a test run, setting the control word to 0x27F gives the expected result. Smile
Post 23 Sep 2010, 00:20
View user's profile Send private message Reply with quote
Display posts from previous:
Post new topic Reply to topic

Jump to:  


< Last Thread | Next Thread >
Forum Rules:
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
You cannot attach files in this forum
You can download files in this forum


Copyright © 1999-2025, Tomasz Grysztar. Also on GitHub, YouTube.

Website powered by rwasa.