Quad Precision

[srvaldez]

here's the translated Fortran quad precision module, unfortunately I had to use assembly, FB kept giving strange results.
added some operator overloading to make it easier to use.
have split the file so the forum won't choke.
part 1 of 3
"quad.bi"

Code: Select all

'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"

namespace quad_precision

'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
    
    Declare constructor ( Byref rhs As quad )
    Declare constructor ( byref rhs As Integer = 0 )
    Declare constructor ( byref rhs As Long = 0 )
    Declare constructor ( byref rhs As LongInt = 0 )
    Declare constructor ( byref rhs As UInteger = 0 )
    Declare constructor ( byref rhs As ULong = 0 )
    Declare constructor ( byref rhs As ULongInt = 0 )
    Declare constructor ( byref rhs As Single )
    Declare constructor ( byref rhs As Double )
    Declare constructor ( Byref rhs As String ="0" )

    Declare Operator Let( Byref rhs As quad )
    Declare Operator Let( Byref rhs As Integer )
    Declare Operator Let( Byref rhs As Long )
    Declare Operator Let( Byref rhs As Longint )
    Declare Operator Let( Byref rhs As Uinteger )
    Declare Operator Let( Byref rhs As ULong )
    Declare Operator Let( Byref rhs As Ulongint )
    Declare Operator Let( Byref rhs As Single )
    Declare Operator Let( Byref rhs As Double )
    Declare Operator Let( Byref rhs As String )  

    '' implicit step versions
    Declare Operator For ( )
    Declare Operator Step( )
    Declare Operator Next( Byref end_cond As quad ) As Integer
    
    '' explicit step versions
    Declare Operator For ( Byref step_var As quad )
    Declare Operator Step( Byref step_var As quad )
    Declare Operator Next( Byref end_cond As quad, Byref step_var As quad ) As Integer
    Declare Operator += ( Byref rhs As quad )
    Declare Operator += ( Byref rhs As Double )
    Declare Operator += ( Byref rhs As Integer )
    Declare Operator += ( Byref rhs As String )
    Declare Operator -= ( Byref rhs As quad )
    Declare Operator -= ( Byref rhs As Double )
    Declare Operator -= ( Byref rhs As Integer )
    Declare Operator -= ( Byref rhs As String )
    Declare Operator *= ( Byref rhs As quad )
    Declare Operator *= ( Byref rhs As Double )
    Declare Operator *= ( Byref rhs As Integer )
    Declare Operator *= ( Byref rhs As String )
    Declare Operator /= ( Byref rhs As quad )
    Declare Operator /= ( Byref rhs As Double )
    Declare Operator /= ( Byref rhs As Integer )
    Declare Operator /= ( Byref rhs As String )
    Declare Operator ^= ( Byref rhs As quad )
    Declare Operator ^= ( Byref rhs As Double )
    Declare Operator ^= ( Byref rhs As Integer )
    Declare Operator ^= ( Byref rhs As String )
    Declare Operator cast( ) As Integer
    Declare Operator cast( ) As Long
    Declare Operator cast( ) As Longint
    Declare Operator cast( ) As Uinteger
    Declare Operator cast( ) As ULong
    Declare Operator cast( ) As Ulongint
    Declare Operator cast( ) As Single
    Declare Operator cast( ) As Double
    Declare Operator cast( ) As String    
End Type

Type quad_rm
    As Integer n
    As quad rm
end type

Declare Function longadd(Byref a As quad, Byref c As quad) As quad
Declare Function string_quad(Byref value As String) As quad
Declare Function quad_string(Byref value As quad) As String

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

Function exactmul2(Byref a As Double, Byref c As Double) As quad 'Result(ac)
'  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             ' Lahey's optimization removes the brackets
'                                 ' and sets a1 = a which defeats the whole point.
'    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
        Asm
    ' t = constant * a

            movsd     xmm0, qword Ptr [constant]
            mov       eax, dword Ptr [a]                  
            movsd     xmm1, qword Ptr [eax]                 
            mulsd     xmm0, xmm1
            movsd     qword Ptr [t], xmm0

    ' a1 = (a - t) + t ! lahey's optimization removes the brackets
    ' and sets a1 = a which defeats the whole point.

            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             

    ' a2 = a - a1

            mov       eax, dword Ptr [a]                  
            movsd     xmm0, qword Ptr [eax]               
            movsd     xmm1, qword Ptr [a1]                
            subsd     xmm0, xmm1                          
            movsd     qword Ptr [a2], xmm0                

    ' t = constant * c

            movsd     xmm0, qword Ptr [constant]  
            mov       eax, dword Ptr [c]                    
            movsd     xmm1, qword Ptr [eax]                  
            mulsd     xmm0, xmm1                             
            movsd     qword Ptr [t], xmm0                   

    ' c1 = (c - t) + t

            mov       eax, dword Ptr [c]                 
            movsd     xmm0, qword Ptr [eax]              
            movsd     xmm1, qword Ptr [t]                
            subsd     xmm0, xmm1                         
            movsd     xmm1, qword Ptr [t]                
            addsd     xmm0, xmm1                         
            movsd     qword Ptr [c1], xmm0               

    ' c2 = c - c1

            mov       eax, dword Ptr [c]                 
            movsd     xmm0, qword Ptr [eax]              
            movsd     xmm1, qword Ptr [c1]               
            subsd     xmm0, xmm1                         
            movsd     qword Ptr [c2], xmm0               

    ' ac.hi = a * c

            mov       eax, dword Ptr [a]                 
            movsd     xmm0, qword Ptr [eax]              
            mov       eax, dword Ptr [c]                 
            movsd     xmm1, qword Ptr [eax]              
            mulsd     xmm0, xmm1                         
            lea       eax, dword Ptr [ac]                
            movsd     qword Ptr [eax], xmm0                  

    ' ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2

            movsd     xmm0, qword Ptr [a1]               
            movsd     xmm1, qword Ptr [c1]               
            mulsd     xmm0, xmm1                         
            lea       eax, dword Ptr [ac]                
            movsd     xmm1, qword Ptr [eax]              
            subsd     xmm0, xmm1                         
            movsd     xmm1, qword Ptr [a1]               
            movsd     xmm2, qword Ptr [c2]               
            mulsd     xmm1, xmm2                         
            addsd     xmm0, xmm1                         
            movsd     xmm1, qword Ptr [c1]               
            movsd     xmm2, qword Ptr [a2]               
            mulsd     xmm1, xmm2                         
            addsd     xmm0, xmm1                         
            movsd     xmm1, qword Ptr [c2]               
            movsd     xmm2, qword Ptr [a2]               
            mulsd     xmm1, xmm2                         
            addsd     xmm0, xmm1                         
            lea       eax, dword Ptr [ac]                
            movsd     qword Ptr [8+eax], xmm0     'ac.lo      
        End asm
    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.

    Dim As quad ac, z

    '     local variables
    Dim As Double zz, a1, a2, c1, c2, t

'    z = exactmul2(a.hi, c.hi)
    asm
' t = constant * a

        movsd     xmm0, qword Ptr [constant]
        mov       eax, dword Ptr [a]                  
        movsd     xmm1, qword Ptr [eax]                 
        mulsd     xmm0, xmm1
        movsd     qword Ptr [t], xmm0

' a1 = (a - t) + t ! lahey's optimization removes the brackets
' and sets a1 = a which defeats the whole point.

        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             

' a2 = a - a1

        mov       eax, dword Ptr [a]                  
        movsd     xmm0, qword Ptr [eax]               
        movsd     xmm1, qword Ptr [a1]                
        subsd     xmm0, xmm1                          
        movsd     qword Ptr [a2], xmm0                

' t = constant * c

        movsd     xmm0, qword Ptr [constant]  
        mov       eax, dword Ptr [c]                    
        movsd     xmm1, qword Ptr [eax]                  
        mulsd     xmm0, xmm1                             
        movsd     qword Ptr [t], xmm0                   

' c1 = (c - t) + t

        mov       eax, dword Ptr [c]                 
        movsd     xmm0, qword Ptr [eax]              
        movsd     xmm1, qword Ptr [t]                
        subsd     xmm0, xmm1                         
        movsd     xmm1, qword Ptr [t]                
        addsd     xmm0, xmm1                         
        movsd     qword Ptr [c1], xmm0               

' c2 = c - c1

        mov       eax, dword Ptr [c]                 
        movsd     xmm0, qword Ptr [eax]              
        movsd     xmm1, qword Ptr [c1]               
        subsd     xmm0, xmm1                         
        movsd     qword Ptr [c2], xmm0               

' ac.hi = a * c

        mov       eax, dword Ptr [a]                 
        movsd     xmm0, qword Ptr [eax]              
        mov       eax, dword Ptr [c]                 
        movsd     xmm1, qword Ptr [eax]              
        mulsd     xmm0, xmm1                         
        lea       eax, dword Ptr [ac]                
        movsd     qword Ptr [eax], xmm0                  

' ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2

        movsd     xmm0, qword Ptr [a1]               
        movsd     xmm1, qword Ptr [c1]               
        mulsd     xmm0, xmm1                         
        lea       eax, dword Ptr [ac]                
        movsd     xmm1, qword Ptr [eax]              
        subsd     xmm0, xmm1                         
        movsd     xmm1, qword Ptr [a1]               
        movsd     xmm2, qword Ptr [c2]               
        mulsd     xmm1, xmm2                         
        addsd     xmm0, xmm1                         
        movsd     xmm1, qword Ptr [c1]               
        movsd     xmm2, qword Ptr [a2]               
        mulsd     xmm1, xmm2                         
        addsd     xmm0, xmm1                         
        movsd     xmm1, qword Ptr [c2]               
        movsd     xmm2, qword Ptr [a2]               
        mulsd     xmm1, xmm2                         
        addsd     xmm0, xmm1                         
        lea       eax, dword Ptr [ac]                
        movsd     qword Ptr [8+eax], xmm0     'ac.lo      
    End asm
    z=ac
    asm
    'zz = ((a.hi + a.lo) * c.lo + a.lo * c.hi) + z.lo
        mov       eax, dword Ptr [a]
        mov       edx, dword Ptr [a]  'a
        movsd     xmm0, qword Ptr [eax]   'a.hi
        movsd     xmm1, qword Ptr [8+edx] 'a.lo
        addsd     xmm0, xmm1
        mov       eax, dword Ptr [c]
        movsd     xmm1, qword Ptr [8+eax] 'c.lo
        mulsd     xmm0, xmm1
        mov       eax, dword Ptr [a]
        movsd     xmm1, qword Ptr [8+eax] 'a.lo
        mov       eax, dword Ptr [c] 'c
        movsd     xmm2, qword Ptr [eax] 'c.hi
        mulsd     xmm1, xmm2
        addsd     xmm0, xmm1 
        movsd     xmm1, qword Ptr [z+8] 'z.lo
        addsd     xmm0, xmm1
        movsd     qword Ptr [zz], xmm0 
    End asm

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

    Return ac
End Function


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

    Dim As quad c = cquad(zero, zero)
    Dim As Double di = i

    If (i = 1) Then
      c = a
    Elseif (i = -1) Then
      c.hi = -a.hi
      c.lo = -a.lo
    Else
      c = longadd(exactmul2(a.hi, di) , exactmul2(a.lo, di))
    End If

    Return c
End Function


Function mult_int_quad(Byref i As Integer, Byref a As quad) As quad'Result(c)
'     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).

    Dim As quad c, z

    '     local variables
    Dim As Double zz
    z = exactmul2(a.hi, b)

    asm
' zz = a.lo * b + z.lo
        mov       eax, dword Ptr [a]
        movsd     xmm0, qword Ptr [8+eax] 'a.lo
        mov       eax, dword Ptr [b]
        movsd     xmm1, qword Ptr [eax] 'b
        mulsd     xmm0, xmm1 'a.lo * b
        movsd     xmm1, qword Ptr [z+8]'z.lo
        addsd     xmm0, xmm1 'a.lo * b + z.lo
        movsd     qword Ptr [zz], xmm0 '-> zz
' c.hi = z.hi + zz
        movsd     xmm0, qword Ptr [z] 'z.hi
        movsd     xmm1, qword Ptr [zz]
        addsd     xmm0, xmm1  'z.hi + zz
        lea       eax, dword Ptr [c] 
        movsd     qword Ptr [eax], xmm0 '-> c.hi
' c.lo = (z.hi - c.hi) + zz
        lea       eax, dword Ptr [c] 
        movsd     xmm0, qword Ptr [z]
        movsd     xmm1, qword Ptr [eax]
        subsd     xmm0, xmm1 'z.hi - c.hi
        movsd     xmm1, qword Ptr [zz]
        addsd     xmm0, xmm1 '(z.hi - c.hi) + zz
        lea       eax, dword Ptr [c] 
        movsd     qword Ptr [8+eax], xmm0 '-> c.lo
    End asm

    Return c
End Function



Function mult_quad_Real(Byref a As quad, Byref b As Single) As quad 'Result(c)
'  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, Cdbl(b))
    asm
        ' zz = a.lo * b + z.lo

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax] 'a.lo
            mov       eax, dword Ptr [b]
            movss     xmm1, dword Ptr [eax]
            cvtss2sd  xmm1, xmm1
            mulsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z+8] 'z.lo
            addsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0
    ' c.hi = z.hi + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0

    ' c.lo = (z.hi - c.hi) + zz

            lea       eax, dword Ptr [c]
            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 [c]
            movsd     qword Ptr [8+eax], xmm0    
    End asm

Return c
End Function 'mult_quad_Real



Function mult_dp_quad(Byref b As Double, Byref a As quad) As quad 'Result(c)
'  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)
    asm
    ' zz = a.lo * b + z.lo

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax]
            mov       eax, dword Ptr [b]
            movsd     xmm1, qword Ptr [eax]
            mulsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z+8]
            addsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' c.hi = z.hi + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0

    ' c.lo = (z.hi - c.hi) + zz

            lea       eax, dword Ptr [c]
            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 [c]
            movsd     qword Ptr [8+eax], xmm0
    End asm

    Return c
End Function 'mult_dp_quad



Function mult_Real_quad(Byref a As Single, Byref b As quad) As quad 'Result(c)
'  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(b.hi, Cdbl(a))
    asm
    ' zz = b.lo * a + z.lo

            mov       eax, dword Ptr [b]
            movsd     xmm0, qword Ptr [8+eax]
            mov       eax, dword Ptr [a]
            movss     xmm1, dword Ptr [eax]
            cvtss2sd  xmm1, xmm1
            mulsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z+8]
            addsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' c.hi = z.hi + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0

    ' c.lo = (z.hi - c.hi) + zz

            lea       eax, dword Ptr [c]
            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 [c]
            movsd     qword Ptr [8+eax], xmm0
    End asm

    Return c
End Function 'mult_Real_quad



Function longdiv(Byref a As quad, Byref c As quad) As quad 'Result(ac)
'  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)

    asm
    ' zz = ((((a.hi - q.hi) - q.lo) + a.lo) - z*c.lo) / (c.hi + c.lo)

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [eax]
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [a]
            movsd     xmm1, qword Ptr [8+eax]
            addsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z]
            mov       eax, dword Ptr [c]
            movsd     xmm2, qword Ptr [8+eax]
            mulsd     xmm1, xmm2
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [c]
            mov       edx, dword Ptr [c]
            movsd     xmm1, qword Ptr [eax]
            movsd     xmm2, qword Ptr [8+edx]
            addsd     xmm1, xmm2
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' 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 'longdiv



Function div_quad_dp(Byref a As quad, Byref b As Double) As quad 'Result(c)
'     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)
    asm
    ' zz = (((a.hi - q.hi) - q.lo) + a.lo) / b

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [eax]
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [a]
            movsd     xmm1, qword Ptr [8+eax]
            addsd     xmm0, xmm1
            mov       eax, dword Ptr [b]
            movsd     xmm1, qword Ptr [eax]
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' c.hi = z + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0

    ' c.lo = (z - c.hi) + zz

            lea       eax, dword Ptr [c]
            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 [c]
            movsd     qword Ptr [8+eax], xmm0
    End asm

    Return c
End Function 'div_quad_dp


Function div_quad_Real(Byref a As quad, Byref b As Single) As quad'Result(c)
'     Divide a quadruple-precision number (a) by a Real number (b)

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

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

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [eax]
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [a]
            movsd     xmm1, qword Ptr [8+eax]
            addsd     xmm0, xmm1
            mov       eax, dword Ptr [b]
            movss     xmm1, dword Ptr [eax]
            cvtss2sd  xmm1, xmm1
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' c.hi = z + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0
    End asm
    Return c
End Function 'div_quad_Real


Function div_quad_int(Byref a As quad, Byref b As Integer) As quad 'Result(c)
'     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)
    asm

    ' zz = (((a.hi - q.hi) - q.lo) + a.lo) / b

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [eax]
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [a]
            movsd     xmm1, qword Ptr [8+eax]
            addsd     xmm0, xmm1
            mov       eax, dword Ptr [b]
            mov       eax, dword Ptr [eax]
            cvtsi2sd  xmm1, eax
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' c.hi = z + zz

            movsd     xmm0, qword Ptr [z]
            movsd     xmm1, qword Ptr [zz]
            addsd     xmm0, xmm1
            lea       eax, dword Ptr [c]
            movsd     qword Ptr [eax], xmm0

    ' c.lo = (z - c.hi) + zz

            lea       eax, dword Ptr [c]
            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 [c]
            movsd     qword Ptr [8+eax], xmm0
    End asm

    Return c
End Function 'div_quad_int



Function div_int_quad(Byref a As Integer, Byref c As quad) As quad 'Result(ac)
'  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)
    asm
    ' zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)

            mov       eax, dword Ptr [a]
            mov       eax, dword Ptr [eax]
            cvtsi2sd  xmm0, eax
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z]
            mov       eax, dword Ptr [c]
            movsd     xmm2, qword Ptr [8+eax]
            mulsd     xmm1, xmm2
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [c]
            mov       edx, dword Ptr [c]
            movsd     xmm1, qword Ptr [eax]
            movsd     xmm2, qword Ptr [8+edx]
            addsd     xmm1, xmm2
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' 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 'div_int_quad



Function div_Real_quad(Byref a As Single, Byref c As quad) As quad 'Result(ac)
'  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)
    asm
    ' zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)

            mov       eax, dword Ptr [a]
            movss     xmm0, dword Ptr [eax]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z]
            mov       eax, dword Ptr [c]
            movsd     xmm2, qword Ptr [8+eax]
            mulsd     xmm1, xmm2
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [c]
            mov       edx, dword Ptr [c]
            movsd     xmm1, qword Ptr [eax]
            movsd     xmm2, qword Ptr [8+edx]
            addsd     xmm1, xmm2
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' 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 'div_Real_quad



Function div_dp_quad(Byref a As Double, Byref c As quad) As quad 'Result(ac)
'  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)
    asm
    ' zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)

            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [eax]
            movsd     xmm1, qword Ptr [q]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [q+8]
            subsd     xmm0, xmm1
            movsd     xmm1, qword Ptr [z]
            mov       eax, dword Ptr [c]
            movsd     xmm2, qword Ptr [8+eax]
            mulsd     xmm1, xmm2
            subsd     xmm0, xmm1
            mov       eax, dword Ptr [c]
            mov       edx, dword Ptr [c]
            movsd     xmm1, qword Ptr [eax]
            movsd     xmm2, qword Ptr [8+edx]
            addsd     xmm1, xmm2
            divsd     xmm0, xmm1
            movsd     qword Ptr [zz], xmm0

    ' 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 '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.

    Dim As quad ac

    '     local variables
    Dim As Double z, q, zz
    
    asm
' z = a.hi + c.hi        
        mov       eax, dword Ptr [a]
        mov       edx, dword Ptr [c]
        movsd     xmm0, qword Ptr [eax] 'a
        movsd     xmm1, qword Ptr [edx] 'c
        addsd     xmm0, xmm1            'a.hi + c.hi
        movsd     qword Ptr [z], xmm0   '-> z
' q = a.hi - z
        mov       eax, dword Ptr [a]
        movsd     xmm0, qword Ptr [eax]
        movsd     xmm1, qword Ptr [z]
        subsd     xmm0, xmm1            'a.hi - z
        movsd     qword Ptr [q], xmm0   '-> q
' zz = (((q + c.hi) + (a.hi - (q + z))) + a.lo) + c.lo
        mov       eax, dword Ptr [c]
        movsd     xmm0, qword Ptr [q]
        movsd     xmm1, qword Ptr [eax]
        addsd     xmm0, xmm1            'q + c.hi
        mov       eax, dword Ptr [a]
        movsd     xmm1, qword Ptr [q]
        movsd     xmm2, qword Ptr [z]
        addsd     xmm1, xmm2            'q + z
        movsd     xmm2, qword Ptr [eax] 'a.hi
        subsd     xmm2, xmm1            'a.hi - (q + z)
        addsd     xmm0, xmm2    '(q + c.hi) + (a.hi - (q + z))
        mov       eax, dword Ptr [a]
        movsd     xmm1, qword Ptr [8+eax] 'a.lo
        addsd     xmm0, xmm1 '(((q + c.hi) + (a.hi - (q + z))) + a.lo)
        mov       eax, dword Ptr [c]
        movsd     xmm1, qword Ptr [8+eax] 'c.lo
        addsd     xmm0, xmm1 '(((q + c.hi) + (a.hi - (q + z))) + a.lo) + c.lo
        movsd     qword Ptr [zz], xmm0 '-> zz
' ac.hi = z + zz
        movsd     xmm0, qword Ptr [z]
        movsd     xmm1, qword Ptr [zz]
        addsd     xmm0, xmm1    'z + zz
        lea       eax, dword Ptr [ac]
        movsd     qword Ptr [eax], xmm0 '-> ac.hi
' ac.lo = (z - ac.hi) + zz
        lea       eax, dword Ptr [ac]
        movsd     xmm0, qword Ptr [z]
        movsd     xmm1, qword Ptr [eax] 'ac.hi
        subsd     xmm0, xmm1    'z - ac.hi
        movsd     xmm1, qword Ptr [zz]
        addsd     xmm0, xmm1    '(z - ac.hi) + zz
        lea       eax, dword Ptr [ac]
        movsd     qword Ptr [8+eax], xmm0 '-> ac.lo
    End asm

    Return ac
End Function


Function quad_add_int(Byref a As quad, Byref c As Integer) As quad 'Result(ac)

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

            mov       eax, dword Ptr [c]
            mov       eax, dword Ptr [eax]
            cvtsi2sd  xmm0, eax
            movsd     xmm1, qword Ptr [q]
            addsd     xmm1, xmm0
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [q]
            movsd     xmm2, qword Ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword Ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword Ptr [zz], xmm1

    ' 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 'quad_add_int



Function quad_add_Real(Byref a As quad, Byref c As Single) As quad 'Result(ac)

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

            mov       eax, dword Ptr [c]
            movss     xmm0, dword Ptr [eax]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword Ptr [q]
            addsd     xmm1, xmm0
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [q]
            movsd     xmm2, qword Ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword Ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword Ptr [zz], xmm1

    ' 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 'quad_add_Real



Function quad_add_dp(Byref a As quad, Byref c As Double) As quad 'Result(ac)

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

            mov       eax, dword Ptr [c]
            movsd     xmm0, qword Ptr [q]
            movsd     xmm1, qword Ptr [eax]
            addsd     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
            movsd     qword Ptr [zz], xmm0

    ' 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 'quad_add_dp



Function int_add_quad(Byref c As Integer, Byref a As quad) As quad 'Result(ac)

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

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

            mov       eax, dword Ptr [c]
            mov       eax, dword Ptr [eax]
            cvtsi2sd  xmm0, eax
            movsd     xmm1, qword Ptr [q]
            addsd     xmm1, xmm0
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [q]
            movsd     xmm2, qword Ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword Ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword Ptr [zz], xmm1

    ' 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 'int_add_quad



Function Real_add_quad(Byref c As Single, Byref a As quad) As quad 'Result(ac)

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

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

            mov       eax, dword Ptr [c]
            movss     xmm0, dword Ptr [eax]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword Ptr [q]
            addsd     xmm1, xmm0
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [q]
            movsd     xmm2, qword Ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword Ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword Ptr [a]
            movsd     xmm0, qword Ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword Ptr [zz], xmm1

    ' 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 'Real_add_quad

[srvaldez]

part 2 of 3 "quad.bi"

Code: Select all

Function dp_add_quad(byref c as double, byref a as quad) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

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

            mov       eax, dword ptr [c]
            movsd     xmm0, qword ptr [q]
            movsd     xmm1, qword ptr [eax]
            addsd     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
            movsd     qword ptr [zz], xmm0

    ' 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 'dp_add_quad



Function longsub(byref a as quad, byref c as quad) as quad 'Result(ac)
'  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

    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

            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

    ' 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 'longsub



Function quad_sub_int(byref a as quad, byref c as integer) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = a.hi - c

            mov       eax, dword ptr [a]
            mov       edx, dword ptr [c]
            mov       edx, dword ptr [edx]
            cvtsi2sd  xmm0, edx
            movsd     xmm1, qword ptr [eax]
            subsd     xmm1, xmm0
            movsd     qword ptr [z], xmm1

    ' 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) + (a.hi - (q + z))) + a.lo)

            mov       eax, dword ptr [c]
            mov       eax, dword ptr [eax]
            cvtsi2sd  xmm0, eax
            movsd     xmm1, qword ptr [q]
            subsd     xmm1, xmm0
            mov       eax, dword ptr [a]
            movsd     xmm0, qword ptr [q]
            movsd     xmm2, qword ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword ptr [a]
            movsd     xmm0, qword ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword ptr [zz], xmm1

    ' 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 'quad_sub_int



Function quad_sub_Real(byref a as quad, byref c as single) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = a.hi - c

            mov       eax, dword ptr [a]
            mov       edx, dword ptr [c]
            movss     xmm0, dword ptr [edx]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword ptr [eax]
            subsd     xmm1, xmm0
            movsd     qword ptr [z], xmm1

    ' 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) + (a.hi - (q + z))) + a.lo)

            mov       eax, dword ptr [c]
            movss     xmm0, dword ptr [eax]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword ptr [q]
            subsd     xmm1, xmm0
            mov       eax, dword ptr [a]
            movsd     xmm0, qword ptr [q]
            movsd     xmm2, qword ptr [z]
            addsd     xmm0, xmm2
            movsd     xmm2, qword ptr [eax]
            subsd     xmm2, xmm0
            addsd     xmm1, xmm2
            mov       eax, dword ptr [a]
            movsd     xmm0, qword ptr [8+eax]
            addsd     xmm1, xmm0
            movsd     qword ptr [zz], xmm1

    ' 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 'quad_sub_Real



Function quad_sub_dp(byref a as quad, byref c as double) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = a.hi - c

            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) + (a.hi - (q + z))) + a.lo)

            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
            movsd     qword ptr [zz], xmm0

    ' 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 'quad_sub_dp



Function int_sub_quad(byref a as integer, byref c as quad) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = dble(a) - c.hi

            mov       eax, dword ptr [a]
            mov       eax, dword ptr [eax]
            cvtsi2sd  xmm0, eax
            mov       eax, dword ptr [c]
            movsd     xmm1, qword ptr [eax]
            subsd     xmm0, xmm1
            movsd     qword ptr [z], xmm0

    ' q = dble(a) - z

            mov       eax, dword ptr [a]
            mov       eax, dword ptr [eax]
            cvtsi2sd  xmm0, eax
            movsd     xmm1, qword ptr [z]
            subsd     xmm0, xmm1
            movsd     qword ptr [q], xmm0

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

            mov       eax, dword ptr [c]
            movsd     xmm0, qword ptr [q]
            movsd     xmm1, qword ptr [eax]
            subsd     xmm0, xmm1
            mov       eax, dword ptr [a]
            mov       eax, dword ptr [eax]
            cvtsi2sd  xmm1, eax
            movsd     xmm2, qword ptr [q]
            movsd     xmm3, qword ptr [z]
            addsd     xmm2, xmm3
            subsd     xmm1, xmm2
            addsd     xmm0, xmm1
            mov       eax, dword ptr [c]
            movsd     xmm1, qword ptr [8+eax]
            subsd     xmm0, xmm1
            movsd     qword ptr [zz], xmm0

    ' 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 'int_sub_quad



Function Real_sub_quad(byref a as single, byref c as quad) as quad'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = a - c.hi

            mov       eax, dword ptr [a]
            movss     xmm0, dword ptr [eax]
            cvtss2sd  xmm0, xmm0
            mov       eax, dword ptr [c]
            movsd     xmm1, qword ptr [eax]
            subsd     xmm0, xmm1
            movsd     qword ptr [z], xmm0

    ' q = a - z

            mov       eax, dword ptr [a]
            movss     xmm0, dword ptr [eax]
            cvtss2sd  xmm0, xmm0
            movsd     xmm1, qword ptr [z]
            subsd     xmm0, xmm1
            movsd     qword ptr [q], xmm0

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

            mov       eax, dword ptr [c]
            movsd     xmm0, qword ptr [q]
            movsd     xmm1, qword ptr [eax]
            subsd     xmm0, xmm1
            mov       eax, dword ptr [a]
            movss     xmm1, dword ptr [eax]
            cvtss2sd  xmm1, xmm1
            movsd     xmm2, qword ptr [q]
            movsd     xmm3, qword ptr [z]
            addsd     xmm2, xmm3
            subsd     xmm1, xmm2
            addsd     xmm0, xmm1
            mov       eax, dword ptr [c]
            movsd     xmm1, qword ptr [8+eax]
            subsd     xmm0, xmm1
            movsd     qword ptr [zz], xmm0

    ' 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 'Real_sub_quad




Function dp_sub_quad(byref a as double, byref c as quad) as quad 'Result(ac)

    '     Local variables
    dim as double z, q, zz
    dim as quad ac

    asm
    ' z = a - 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 - 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 - (q + z))) - c.lo

            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 [c]
            movsd     xmm1, qword ptr [8+eax]
            subsd     xmm0, xmm1
            movsd     qword ptr [zz], xmm0

    ' 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 'dp_sub_quad



Function negate_quad(byref a as quad) as quad 'Result(b)
'     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



Sub quad_eq_int(byref a as quad, byref i as integer)
'     Assignment

    a.hi = i
    a.lo = 0

    Return
End Sub 'quad_eq_int



Sub quad_eq_Real(byref a as quad, byref r as single)
    '     Assignment

    a.hi = r
    a.lo = zero

    Return
End Sub 'quad_eq_Real



Sub quad_eq_dp(byref a as quad, byref d as double)
    '     Assignment

    a.hi = d
    a.lo = zero

    Return
End Sub 'quad_eq_dp



Sub int_eq_quad(byref i as integer, byref a as quad)
    '     Assignment

    i = a.hi

    Return
End Sub 'int_eq_quad



Sub Real_eq_quad(byref r as single, byref a as quad)
    '     Assignment

    r = a.hi

    Return
End Sub 'Real_eq_quad



Sub dp_eq_quad(byref d as double, byref a as quad)
    '     Assignment

    d = a.hi

    Return
End Sub 'dp_eq_quad



Function quad_lt(byref x as quad, byref y as quad) as integer 'Result(is_it)
'     Comparison of 2 logical numbers

    dim as integer is_it

    ' Local variable
    dim as quad dIff

    dIff = longsub(x , y)
    is_it = (dIff.hi < zero)

    Return is_it
End Function 'quad_lt



Function quad_le(byref x as quad, byref y as quad) as integer 'Result(is_it)
'     Comparison of 2 logical numbers

    dim as integer is_it

    ' Local variable
    dim as quad dIff

    dIff = longsub(x , y)
    is_it = (not (dIff.hi > zero))

    Return is_it
End Function 'quad_le



Function quad_eq(byref x as quad, byref y as quad) as integer'Result(is_it)
    '     Comparison of 2 logical numbers

    dim as integer is_it

    ' Local variable
    dim as quad dIff

    dIff = longsub(x , y)
    is_it = (dIff.hi = zero)

    Return is_it
End Function 'quad_eq



Function quad_ge(byref x as quad, byref y as quad) as integer 'Result(is_it)
'     Comparison of 2 logical numbers

    dim as integer is_it
    ' Local variable
    dim as quad dIff

    dIff = longsub(x , y)
    is_it = (not (dIff.hi < zero))

    Return is_it
End Function 'quad_ge



Function quad_gt(byref x as quad, byref y as quad) as integer 'Result(is_it)
'     Comparison of 2 logical numbers

    dim as integer is_it

    ' Local variable
    dim as quad dIff

    dIff = longsub(x , y)
    is_it = (dIff.hi > zero)

    Return is_it
End Function 'quad_gt



function quad_pow_int(byref x as quad, byref 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 (bit(n,0)=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 a as quad, byref i as integer ) as quad
    '     Multiply a by 2^i

    dim as quad b
    b=cquad(2,0)
    b=quad_pow_int(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 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 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'

    dim as quad y

    ' local variables
    dim as quad temp, ysq, sum1, sum2
    dim as integer iy


    y = longdiv(x , ln2)

    iy = nint(csng(y.hi))

    y = quad_sub_dp(y, cdbl(iy))
    y = longmul(y, 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.0)
    sum1=longmul(sum1,ysq)
    sum1=quad_add_dp(sum1 , 2162160.0)
    sum1=longmul(sum1,ysq)
    sum1=quad_add_dp(sum1 , 302702400.0)
    sum1=longmul(sum1,ysq)
    sum1=quad_add_dp(sum1 , 8821612800.0)
    sum1=longmul(y,sum1)

    sum2=mult_quad_dp(ysq, 90.0)
    sum2=quad_add_dp(sum2 , 110880.0)
    sum2=longmul(sum2,ysq)
    sum2=quad_add_dp(sum2 , 30270240.0)
    sum2=longmul(sum2,ysq)
    sum2=quad_add_dp(sum2 , 2075673600.0)
    sum2=longmul(sum2,ysq)
    sum2=quad_add_dp(sum2 , 17643225600.0)

    '                     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.0)

    y = longmul(temp , y)
    y = qscale(y, 2)
    y = quad_add_dp(y , 1.0)

    y = qscale(y, iy)

    return y
end function

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_pow_Real(byref a as quad, byref r as single) as quad 'Result(b)
'     Raise a quadruple=precision number (a) to a power.

    dim as quad b

    If (a.hi < zero) Then
      print _
      " *** Error: attempt to raise negative quad. prec. number to a Real power ***"
      b = cquad(zero, zero)
    ElseIf (a.hi > zero) Then
      b = longExp( mult_quad_Real(longLog(a), r) )
    Else
      b = cquad(zero, zero)
    End If

    Return b
End Function 'quad_pow_Real



Function quad_pow_dp(byref a as quad, byref d as double) as quad 'Result(b)
'     Raise a quadruple=precision number (a) to a power.

    dim as quad b

    If (a.hi < zero) Then
      print _
      " *** Error: attempt to raise negative quad. prec. number to a Real power ***"
      b = cquad(zero, zero)
    ElseIf (a.hi > zero) Then
      b = longExp( mult_quad_dp(longLog(a), d) )
    Else
      b = cquad(zero, zero)
    End If

    Return b
End Function 'quad_pow_dp



Function quad_pow_quad(byref a as quad, byref q as quad) as quad 'Result(b)
'     Raise a quadruple=precision number (a) to a power.

    dim as quad b

    If (a.hi < zero) Then
      print _
      " *** Error: attempt to raise negative quad. prec. number to a Real power ***"
      b = cquad(zero, zero)
    ElseIf (a.hi > zero) Then
      b = longExp( longmul(longLog(a), q) )
    Else
      b = cquad(zero, zero)
    End If

    Return b
End Function 'quad_pow_quad



Function qAbs(byref a as quad) as quad'Result(b)
'     Absolute value of a quadruple-precision number

    dim as quad b

    If (a.hi < zero) Then
      b.hi = - a.hi
      b.lo = - a.lo
    Else
      b.hi = a.hi
      b.lo = a.lo
    End If

    Return b
End Function 'qabs



Function longSqr(byref a as quad) as quad 'Result(b)
' This is modIfied from procedure sqrt2 of:
'    Dekker, T.J. (1971). 'A floating-point technique for extending the
'    available precision', Numer. Math., 18, 224-242.

'     Local variables
    dim as double t, res
    dim as quad b, tt

    ' Check that ahi >= 0.

    If (a.hi < 0.0) Then
      print " *** Negative argument for longsqrt ***"
      Return b
    ElseIf (a.hi = 0.0) Then
      b.hi = 0.0
      b.lo = 0.0
      Return b
    End If

    ' First approximation is  t = Sqr(a).

    t = Sqr(a.hi)
    tt = exactmul2(t, t)
    asm
    ' res = t + (((a.hi - tt.hi) - tt.lo) + a.lo) * 0.5_dp / t

            mov       eax, dword ptr [a]
            movsd     xmm0, qword ptr [eax]
            movsd     xmm1, qword ptr [tt]
            subsd     xmm0, xmm1
            movsd     xmm1, qword ptr [tt+8]
            subsd     xmm0, xmm1
            mov       eax, dword ptr [a]
            movsd     xmm1, qword ptr [8+eax]
            addsd     xmm0, xmm1
            movsd     xmm1, qword ptr [_longSqr_half_dp]
            mulsd     xmm0, xmm1
            movsd     xmm1, qword ptr [t]
            divsd     xmm0, xmm1
            movsd     xmm1, qword ptr [t]
            addsd     xmm1, xmm0
            movsd     qword ptr [res], xmm1

    ' b.lo = (t - res) + (((a.hi - tt.hi) - tt.lo) + a.lo) * 0.5_dp / t

            movsd     xmm0, qword ptr [t]
            movsd     xmm1, qword ptr [res]
            subsd     xmm0, xmm1
            mov       eax, dword ptr [a]
            movsd     xmm1, qword ptr [eax]
            movsd     xmm2, qword ptr [tt]
            subsd     xmm1, xmm2
            movsd     xmm2, qword ptr [tt+8]
            subsd     xmm1, xmm2
            mov       eax, dword ptr [a]
            movsd     xmm2, qword ptr [8+eax]
            addsd     xmm1, xmm2
            movsd     xmm2, qword ptr [_longSqr_half_dp]
            mulsd     xmm1, xmm2
            movsd     xmm2, qword ptr [t]
            divsd     xmm1, xmm2
            addsd     xmm0, xmm1
            lea       eax, dword ptr [b]
            movsd     qword ptr [8+eax], xmm0

    ' b.hi = res

            lea       eax, dword ptr [b]
            movsd     xmm0, qword ptr [res]
            movsd     qword ptr [eax], xmm0
    end asm
    Return b
    asm
        _longSqr_half_dp:	.int	&h000000000,&h03fe00000
    end asm
End Function 'longsqrt

[srvaldez]

part 3 of 3 "quad.bi"

Code: Select all

Function quad_string(Byref value As quad) 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, ier
    Dim As Double tmp

    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 = FBfloor( 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,(14 - dec_expnt))
        v=longexp(v)
        vl = longmul(vl , v )
    End If

    str1=format( vl.hi,"################")
    ' calculate the remainder
    tmp =Val(str1+".0#")
    vl = quad_sub_dp(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

    str2=format( vl.hi,"################")
    '    end if
    If Len(str2)<15 Then
        str2=String(15-Len(str2),"0")+str2
    End If



    ' if str2 consists of asterisks, add 1 in the last digit to str1.
    ' set str2 to zeroes.
    If (Len(str2)>15) Then
        tmp = tmp + one
        str1=format( tmp,"#################")
    If (Left(str1,1) <> " ") Then
        'dec_expnt = dec_expnt + 1
    Else
        str1 = Mid(str1,2)
    End If
        str2 = "000000000000000"
    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

    ' 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, pt = cquad(10.0, 0.0)
  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
  f1=f1+f2
  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))
  pt = quad_pow_int(pt, 29-ex)
  qd=longdiv(qd,pt)

  Return qd
End Function

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


Sub longmodr(Byref a As quad, Byref b As quad, Byref n As Integer, Byref rm As quad)

' Extended arithmetic calculation of the 'rounded' modulus:
'  a = n.b + rm
' where all quantities are in quadruple-precision, except the Integer
' number of multiples, n.   The absolute value of the remainder (rm)
' is not greater than b/2.
' The result is exact.   remainder may occupy the same location as either input.

' Programmer: Alan Miller

' Latest revision - 11 September 1986
' Fortran version - 4 December 1996

' Local variables

    Dim As quad temp

    ' Check that b.hi .ne. 0

    If (b.hi = 0.0) Then
      Print " *** Error in longmodr - 3rd argument zero ***"
      Return
    End If

    ' Calculate n.

    temp = longdiv(a , b)
    n = nint(Csng(temp.hi))

    ' Calculate remainder preserving full accuracy.

    temp = exactmul2(Cdbl(n), b.hi)
    rm.hi = a.hi
    rm.lo = zero
    temp = longsub(rm , temp)
    rm.hi = a.lo
    rm.lo = zero
    temp = longadd(rm , temp)
    rm = exactmul2(Cdbl(n), b.lo)
    rm = longsub(temp , rm)

    Return
End Sub 'longmodr



Sub longcst(Byref a As quad, Byref b As quad, Byref sine As Integer,_
            Byref cosine As Integer, Byref tangent As Integer)

    ' Local variables

    Dim As Integer pos1
    Dim As quad d, term, temp, angle, sum1, sum2, sin1
    Dim As Integer npi, ipt, i
    Dim As Double tol15 = 1.E-15, tol30 = 1.E-30

    ' Sin(i.pi/40), i = 0(1)20
    Static As quad table(0 To 20)
        table( 0) = cquad(0.0000000000000000E+00,  0.0000000000000000E+00)
        table( 1) = cquad(0.7845909572784494E-01,  0.1464397249532491E-17)
        table( 2) = cquad(0.1564344650402309E+00,  -.2770509565052586E-16)
        table( 3) = cquad(0.2334453638559054E+00,  0.2058612230858154E-16)
        table( 4) = cquad(0.3090169943749475E+00,  -.8267172724967036E-16)
        table( 5) = cquad(0.3826834323650898E+00,  -.1005077269646159E-16)
        table( 6) = cquad(0.4539904997395468E+00,  -.1292033036231312E-16)
        table( 7) = cquad(0.5224985647159488E+00,  0.6606794454708078E-16)
        table( 8) = cquad(0.5877852522924732E+00,  -.1189570533007057E-15)
        table( 9) = cquad(0.6494480483301838E+00,  -.1134903961116171E-15)
        table(10) = cquad(0.7071067811865476E+00,  -.4833646656726458E-16)
        table(11) = cquad(0.7604059656000310E+00,  -.1036987135483477E-15)
        table(12) = cquad(0.8090169943749476E+00,  -.1381828784809282E-15)
        table(13) = cquad(0.8526401643540922E+00,  0.4331886637554353E-16)
        table(14) = cquad(0.8910065241883680E+00,  -.1474714419679880E-15)
        table(15) = cquad(0.9238795325112868E+00,  -.9337725537817898E-16)
        table(16) = cquad(0.9510565162951536E+00,  -.7008780156242836E-16)
        table(17) = cquad(0.9723699203976766E+00,  0.4478912629332321E-16)
        table(18) = cquad(0.9876883405951378E+00,  -.4416018005989794E-16)
        table(19) = cquad(0.9969173337331280E+00,  0.1235153006196267E-16)
        table(20) = cquad(0.1000000000000000E+01,  0.0000000000000000E+00)



    ' Reduce angle to range (-pi/2, +pi/2) by subtracting an Integer multiple of pi.

    longmodr(a, pi, npi, angle)

    ' Find nearest multiple of pi/40 to angle.

    longmodr(angle, piby40, ipt, d)

    ' Sum 1 = 1 - d**2/2' + d**4/4' - d**6/6' + ...
    ' Sum 2 = d - d**3/3' + d**5/5' - d**7/7' + ...

    sum1.hi = zero
    sum1.lo = zero
    sum2.hi = zero
    sum2.lo = zero
    pos1 = 0
    term = d
    i = 2
    L20: If (Abs(term.hi) > tol15) Then
      term = longmul(term , d)                                ' Use quad. precision
      If (i = 2 Or i = 4 Or i = 8) Then
        term.hi = term.hi / i
        term.lo = term.lo / i
      Else
        term = div_quad_int(term , i)
      End If
      If (pos1) Then
        sum1 = longadd(sum1 , term)
      Else
        sum1 = longsub(sum1 , term)
      End If
    Else
      term.hi = term.hi * d.hi / Cdbl(i)             ' Double prec. adequate
      If (pos1) Then
        sum1.lo = sum1.lo + term.hi
      Else
        sum1.lo = sum1.lo - term.hi
      End If
    End If

    ' Repeat for sum2

    i = i + 1
    If (Abs(term.hi) > tol15) Then
      term = longmul(term, div_quad_int( d, i))                      ' Use quad. precision
      If (pos1) Then
        sum2 = longadd(sum2 , term)
      Else
        sum2 = longsub(sum2 , term)
      End If
    Else
      term.hi = term.hi * d.hi / Cdbl(i)             ' Double prec. adequate
      If (pos1) Then
        sum2.lo = sum2.lo + term.hi
      Else
        sum2.lo = sum2.lo - term.hi
      End If
    End If

    i = i + 1
    pos1 = Not pos1
    If (Abs(term.hi) > tol30) Then Goto L20

    sum1 = quad_add_dp(sum1 , 1.0)                              ' Now add the 1st terms
    sum2 = longadd(sum2 , d)                                  ' for max. accuracy

    ' Construct sine, cosine or tangent.
    ' Sine first.    Sin(angle + d) = Sin(angle).Cos(d) + Cos(angle).Sin(d)

    If (sine Or tangent) Then
      If (ipt >= 0) Then
        temp = table(ipt)
      Else
        temp = negate_quad(table( -ipt))
      End If
      b = longmul(sum1 , temp)
      If (ipt >= 0) Then
        temp = table( 20-ipt)
      Else
        temp = table( 20+ipt)
      End If
      b = longadd(b , longmul(sum2 , temp))
      If (npi <> 2*(npi\2)) Then
        b = negate_quad(b)
      End If
      If (tangent) Then
        sin1 = b
      End If
    End If

    ' Cosine or tangent.

    If (cosine Or tangent) Then
      If (ipt >= 0) Then
        temp = table( ipt)
      Else
        temp = negate_quad(table( -ipt))
      End If
      b = longmul(sum2 , temp)
      If (ipt >= 0) Then
        temp = table( 20-ipt)
      Else
        temp = table( 20+ipt)
      End If
      b = longsub(longmul(sum1 , temp) , b)
      If (npi <> 2*(npi\2)) Then
        b = negate_quad(b)
      End If
    End If

    ' Tangent.

    If (tangent) Then

    ' Check that bhi .ne. 0

      If (b.hi = 0.d0) Then
        Print " *** Infinite tangent - routine longcst ***"
        b.hi = 1.0D+308
        b.lo = 0.d0
        Return
      End If
      b = longdiv(sin1 , b)
    End If

    Return
End Sub 'longcst



' Extended accuracy arithmetic sine, cosine & tangent (about 31 decimals).
' Calculates  b = sin, cos or tan (a), where all quantities are in
' quadruple-precision, using table look-up and a Taylor series expansion.
' The result may occupy the same locations as the input value.
' Much of the code is common to all three Functions, and this is in a
' Sub longcst.


Function longSin(Byref a As quad) As quad 'Result(b)

    Dim As quad b

    ' Local variables

    Dim As Integer sine, cosine, tangent

    ' Set logical variables for sine Function.

    sine = -1
    cosine = 0
    tangent = 0
    longcst(a, b, sine, cosine, tangent)

    Return b
End Function 'longsin



Function longCos(Byref a As quad) As quad 'Result(b)

    Dim As quad b

    ' Local variables

    Dim As Integer sine, cosine, tangent

    ' Set logical variables for sine Function.

    sine = 0
    cosine = -1
    tangent = 0
    longcst(a, b, sine, cosine, tangent)

    Return b
End Function 'longcos



Function longTan(Byref a As quad) As quad 'Result(b)

    Dim As quad b

    ' Local variables

    Dim As Integer sine, cosine, tangent

    ' Set logical variables for sine Function.

    sine = 0
    cosine = 0
    tangent = -1
    longcst(a, b, sine, cosine, tangent)

    Return b
End Function 'longtan



Function longAsin(Byref a As quad) As quad 'Result(b)

' Quadratic-precision arc sine (about 31 decimals).
' One Newton-Raphson iteration to solve:  f(b) = Sin(b) - a = 0,
' except when a close to -1 or +1.
' The result (b) may occupy the same location as the input values (a).
' Use ACOS when |a| is close to 1.

    ' Local variables
    Dim As quad y, b, c

    ' Check that -1 <= a.hi <= +1.

    If (a.hi < -1.0 Or a.hi > 1.0) Then
      Print " *** Argument outside range for longasin ***"
      Return b
    End If

    If (Abs(a.hi) < 0.866) Then
      ' First approximation is  y = Asin(a).
      ' Quadruple-precision result is  y - [Sin(y) - a]/Cos(y).

      y.hi = Asin(a.hi)
      y.lo = zero
      'b = y + (a - Sin(y)) / Cos(y.hi)
      b = longadd(y,div_quad_dp(longsub(a,longsin(y)),Cos(y.hi)))
    Else
      ' Calculate Acos(c) where c = Sqr(1 - a^2)
      c = longSqr(dp_sub_quad(one , longmul(a,a)))
      y.hi = Acos(c.hi)
      y.lo = zero
      'b = y + (Cos(y) - c) / Sin(y.hi)
      b = longadd(y,div_quad_dp(longsub(longcos(y),c),Sin(y.hi)))
      If (a.hi < zero) Then b = negate_quad(b)
    End If


    Return b
End Function 'longasin



Function longAcos(Byref a As quad) As quad 'Result(b)

' Quadratic-precision arc cosine (about 31 decimals).
' Newton-Raphson iteration to solve: f(b) = Cos(b) - a = 0.
' The result (b) may occupy the same location as the input values (a).
' When |a| is near 1, use formula from p.175 of
' `Software Manual for the Elementary Functions' by W.J. Cody, Jr. &
' W. Waite, Prentice-Hall, 1980.

    ' Local variables
    Dim As quad y, b, c

    ' Check that -1 <= a.hi <= +1.

    If (a.hi < -1.0 Or a.hi > 1.0) Then
      Print "*** Argument outside range for longacos ***"
      Return b
    End If

    If (Abs(a.hi) < 0.866) Then
      ' First approximation is  y = Acos(a).
      ' Quadruple-precision result is  y + [Cos(y) - a]/Sin(y).

      y.hi = Acos(a.hi)
      y.lo = 0.0
      'b = y + (Cos(y) - a) / Sin(y.hi)
      b = longadd(y, div_quad_dp(longsub(longcos(y),a),Sin(y.hi)))
    Else
      ' Calculate 2.Asin(c) where c = Sqr([1 - |a|]/2)
      'c = Sqr((one - Abs(a))/2)
      c = longsqr(div_quad_int(dp_sub_quad(one, qabs(a)),2))
      y.hi = Asin(c.hi)
      y.lo = zero
      'b = (y - (Sin(y) - c) / Cos(y.hi))*2
      b = mult_quad_int(longsub(y,div_quad_dp(longsub(longsin(y),c),Cos(y.hi))),2)
      If (a.hi < zero) Then b = longsub(pi , b)
    End If

    Return b
End Function 'longacos



Function longAtn(Byref a As quad) As quad 'Result(b)

' Quadratic-precision arc tangent (about 31 decimals).
' Newton-Raphson iteration to solve: f(b) = Tan(b) - a = 0.
' The result (b) may occupy the same location as the input values (a).

    ' Local variables
    Dim As quad b, y
    Dim As Double t

    ' First approximation is  y = Atn(a).
    ' Quadruple-precision result is  y - [Tan(y) - a] * Cos(y)**2.

    y.hi = Atn(a.hi)
    y.lo = 0.0
    'b = y - (Tan(y) - a) * (Cos(y.hi))**2
    t = Cos(y.hi)
    t = t*t
    b = longsub(y,mult_quad_dp(longsub(longtan(y),a),t))
    Return b
End Function 'longatan



Function qAtan2(Byref y As quad, Byref x As quad) As quad 'Result(b)

' Quadratic-precision arc tangent (about 31 decimals).
' As for arc tangent (y/x) except that the result is in the range
'       -pi < ATAN2 <= pi.
' The signs of x and y determine the quadrant.

    ' Local variables
    Dim As quad b, z
    Dim As Double t

    ' First approximation is  z = Atan2(y, x).
    ' Quadruple-precision result is  z - [Tan(z) - (y/x)] * Cos(z)**2.

    z.hi = Atan2(y.hi, x.hi)
    z.lo = 0.0
    If (x.hi = zero) Then
        b = z
    Else
        t = Cos(z.hi)
        t = t*t
        'b = z - (Tan(z) - y/x) * (Cos(z.hi))**2
        b = longsub(z, mult_quad_dp(longsub(longtan(z),longdiv(y, x)),t))
    End If

    Return b
End Function 'qatan2



Function quad_sum(a() As quad) As quad 'Result(s)

    ' Quadruple-precision SUM

    ' Local variables

    Dim As Integer i
    Dim As quad s
    
    s = cquad(zero, zero)

    For i=Lbound(a) To Ubound(a)
        s = longadd(s , a(i))
    Next

    Return s
End Function 'quad_sum



Function quad_dot_product(a() As quad, b() As quad) As quad 'Result(ab)

' Quadruple-precision DOT_PRODUCT

    ' Local variables

    Dim As Integer i, n
    Dim As quad ab

    ab = cquad(zero, zero)
    n = Ubound(a)
    If (n <> Ubound(b)) Or (Lbound(a)<>lbound(b)) Then
      Print "** Error invoking DOT_PRODUCT - dIfferent argument sizes **"
      Print " Size of 1st argument = "; n,   _
                                   "   Size of 2nd argument = "; Ubound(b)
      Return ab
    End If

    For i = Lbound(a) To n
      ab = longadd(ab, longmul( a(i),b(i)))
    Next

    Return ab
End Function 'quad_dot_product


Function quad_int(Byref a As quad) As Integer
    Dim as integer i

    i=int(a.hi)

    Return i
End Function


'Declare Function cquad Overload ( Byref lhs As quad ) As quad
Declare Function cquad ( Byref lhs As Integer ) As quad
Declare Function cquad ( Byref lhs As Long ) As quad
Declare Function cquad ( Byref lhs As Longint ) As quad
Declare Function cquad ( Byref lhs As Uinteger ) As quad
Declare Function cquad ( Byref lhs As ULong ) As quad
Declare Function cquad ( Byref lhs As Ulongint ) As quad
Declare Function cquad ( Byref lhs As Single )  As quad
Declare Function cquad ( Byref lhs As Double )  As quad
Declare Function cquad ( Byref lhs As String )  As quad

Function cquad ( Byref lhs As quad ) As quad
    Return lhs
End Function

Function cquad ( Byref lhs As Integer ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As Long ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As Longint ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As Uinteger ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As ULong ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0 
    Return retval
End Function

Function cquad ( Byref lhs As Ulongint ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As Single ) As quad
    Dim As quad retval
    retval.hi = Cdbl(lhs)
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As Double ) As quad
    Dim As quad retval
    retval.hi = lhs
    retval.lo = 0.0
    Return retval
End Function

Function cquad ( Byref lhs As String ) As quad
    Dim As quad retval
    retval = string_quad ( lhs )
    Return retval
End Function

'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Operator + ( Byref lhs As quad, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = longadd ( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Integer ) As quad
    Dim As quad retval
    retval = quad_add_int(  lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As Integer, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_add_quad ( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Long ) As quad
    Dim As quad retval
    retval = quad_add_int( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As Long, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_add_quad(lhs , rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Longint ) As quad
    Dim As quad retval
    retval = quad_add_dp(lhs, Cdbl(rhs) )
    Return retval
End Operator

Operator + ( Byref lhs As Longint, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = dp_add_quad(Cdbl(lhs), rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Uinteger ) As quad
    Dim As quad retval
    retval = quad_add_int( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As Uinteger, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_add_quad( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As ULong ) As quad
    Dim As quad retval
    retval = quad_add_int( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As ULong, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_add_quad( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Single ) As quad
    Dim As quad retval
    retval = quad_add_Real(lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As Single, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = Real_add_quad( lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As quad, Byref rhs As Double ) As quad
    Dim As quad retval
    retval = quad_add_dp(lhs, rhs )
    Return retval
End Operator

Operator + ( Byref lhs As Double, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = dp_add_quad( lhs, rhs )
    Return retval
End Operator

Operator quad.+= ( Byref rhs As quad )
    Dim As quad retval
    this = longadd(this, rhs )
End Operator

Operator quad.+= ( Byref rhs As Double )
    Dim As quad retval
    this = quad_add_dp(this, rhs )
End Operator

Operator quad.+= ( Byref rhs As Integer )
    Dim As quad retval
    this = quad_add_int(this, rhs )
End Operator

Operator quad.+= ( Byref rhs As String )
    Dim As quad retval
    retval = string_quad( rhs )
    this = longadd(this, retval )
End Operator
'-------------------------------------------------------------

Operator - ( Byref lhs As quad, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = longsub(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad, Byref rhs As Integer ) As quad
    Dim As quad retval
    retval = quad_sub_int(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As Integer, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_sub_quad(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad, Byref rhs As Long ) As quad
    Dim As quad retval
    retval = quad_sub_int(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As Long, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = int_sub_quad(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad, Byref rhs As Longint ) As quad
    Dim As quad retval
    retval = quad_sub_dp(lhs, Cdbl(rhs) )
    Return retval
End Operator

Operator - ( Byref lhs As Longint, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = dp_sub_quad(Cdbl(lhs), rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad, Byref rhs As Single ) As quad
    Dim As quad retval
    retval = quad_sub_Real(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As Single, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = Real_sub_quad(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad, Byref rhs As Double ) As quad
    Dim As quad retval
    retval = quad_sub_dp(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As Double, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = dp_sub_quad(lhs, rhs )
    Return retval
End Operator

Operator - ( Byref lhs As quad ) As quad
    Dim As quad retval
    retval = negate_quad(lhs )
    Return retval
End Operator

Operator quad.-= ( Byref rhs As quad )
    Dim As quad retval
    this = longsub(this, rhs )
End Operator

Operator quad.-= ( Byref rhs As Double )
    Dim As quad retval
    this = quad_sub_dp(this, rhs )
End Operator

Operator quad.-= ( Byref rhs As Integer )
    Dim As quad retval
    this = quad_sub_int(this, rhs )
End Operator

Operator quad.-= ( Byref rhs As String )
    Dim As quad retval
    retval = string_quad( rhs )
    this = longsub(this, retval )
End Operator

'************************************************************

Operator * ( Byref lhs As quad, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = longmul(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As quad, Byref rhs As Integer ) As quad
    Dim As quad retval
    retval = mult_quad_int(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As Integer, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = mult_int_quad(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As quad, Byref rhs As Long ) As quad
    Dim As quad retval
    retval = mult_quad_int(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As Long, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = mult_int_quad(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As quad, Byref rhs As Longint ) As quad
    Dim As quad retval
    retval = mult_quad_dp(lhs, Cdbl(rhs) )
    Return retval
End Operator

Operator * ( Byref lhs As Longint, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = mult_dp_quad(Cdbl(lhs), rhs )
    Return retval
End Operator

Operator * ( Byref lhs As quad, Byref rhs As Single ) As quad
    Dim As quad retval
    retval = mult_quad_Real(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As Single, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = mult_Real_quad(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As quad, Byref rhs As Double ) As quad
    Dim As quad retval
    retval = mult_quad_dp(lhs, rhs )
    Return retval
End Operator

Operator * ( Byref lhs As Double, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = mult_dp_quad(lhs, rhs )
    Return retval
End Operator

Operator quad.*= ( Byref rhs As quad )
    Dim As quad retval
    this = longmul(this, rhs )
End Operator

Operator quad.*= ( Byref rhs As Double )
    Dim As quad retval
    this = mult_quad_dp(this, rhs )
End Operator

Operator quad.*= ( Byref rhs As Integer )
    Dim As quad retval
    this = mult_quad_int(this, rhs )
End Operator

Operator quad.*= ( Byref rhs As String )
    Dim As quad retval
    retval = string_quad( rhs )
    this = longmul(this, retval )
End Operator

'//////////////////////////////////////////////////////////////////

Operator / ( Byref lhs As quad, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = longdiv(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As quad, Byref rhs As Integer ) As quad
    Dim As quad retval
    retval = div_quad_int(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As Integer, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = div_int_quad(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As quad, Byref rhs As Long ) As quad
    Dim As quad retval
    retval = div_quad_int(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As Long, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = div_int_quad(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As quad, Byref rhs As Longint ) As quad
    Dim As quad retval
    retval = div_quad_dp(lhs, Cdbl(rhs) )
    Return retval
End Operator

Operator / ( Byref lhs As Longint, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = div_dp_quad(Cdbl(lhs), rhs )
    Return retval
End Operator

Operator / ( Byref lhs As quad, Byref rhs As Single ) As quad
    Dim As quad retval
    retval = div_quad_Real(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As Single, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = div_Real_quad(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As quad, Byref rhs As Double ) As quad
    Dim As quad retval
    retval = div_quad_dp(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As Double, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = div_dp_quad(lhs, rhs )
    Return retval
End Operator

Operator quad./= ( Byref rhs As quad )
    Dim As quad retval
    this = longdiv(this, rhs )
End Operator

Operator quad./= ( Byref rhs As Double )
    Dim As quad retval
    this = div_quad_dp(this, rhs )
End Operator

Operator quad./= ( Byref rhs As Integer )
    Dim As quad retval
    this = div_quad_int(this, rhs )
End Operator

Operator quad./= ( Byref rhs As String )
    Dim As quad retval
    retval = string_quad( rhs )
    this = longdiv(this, retval )
End Operator

'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Operator ^ ( Byref lhs As quad, Byref rhs As quad ) As quad
    Dim As quad retval
    retval = quad_pow_quad(lhs, rhs )
    Return retval
End Operator

Operator ^ ( Byref lhs As quad, Byref rhs As Integer ) As quad
    Dim As quad retval
    retval = quad_pow_int(lhs, rhs )
    Return retval
End Operator

Operator ^ ( Byref lhs As quad, Byref rhs As Longint ) As quad
    Dim As quad retval
    retval = quad_pow_dp(lhs, Cdbl(rhs) )
    Return retval
End Operator

Operator ^ ( Byref lhs As quad, Byref rhs As Double ) As quad
    Dim As quad retval
    retval = quad_pow_dp(lhs, rhs )
    Return retval
End Operator

Operator ^ ( Byref lhs As Quad, Byref rhs As Single ) As quad
    Dim As quad retval
    retval = quad_pow_Real(lhs, rhs )
    Return retval
End Operator

Operator mod ( Byref lhs As quad, Byref rhs As quad ) As quad_rm
    Dim As quad_rm retval
    Dim As quad rm
    Dim as integer n
    longmodr ( lhs, rhs, n, rm )
    retval.n = n
    retval.rm = rm
    Return retval
End Operator

Operator quad.^= ( Byref rhs As quad )
    Dim As quad retval
    this = quad_pow_quad(this, rhs )
End Operator

Operator quad.^= ( Byref rhs As Double )
    Dim As quad retval
    this = quad_pow_dp(this, Cdbl(rhs) )
End Operator

Operator quad.^= ( Byref rhs As Integer )
    Dim As quad retval
    this = quad_pow_int(this, rhs )
End Operator

Operator quad.^= ( Byref rhs As String )
    Dim As quad retval
    retval = string_quad( rhs )
    this = quad_pow_quad(this, retval )
End Operator

'<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
Operator <> ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = Not quad_eq( lhs, rhs )
    Return relop
End Operator

Operator = ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = quad_eq( lhs, rhs )
    Return relop
End Operator

Operator > ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = quad_gt( lhs, rhs )
    Return relop
End Operator

Operator >= ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = quad_ge( lhs, rhs )
    Return relop
End Operator

Operator <= ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = quad_le( lhs, rhs )
    Return relop
End Operator

Operator < ( Byref lhs As quad, Byref rhs As quad ) As Integer
    Dim As Integer relop
    relop = quad_lt( lhs, rhs )
    Return relop
End Operator


Operator quad.cast ( ) As Integer
    Operator = cint(this.hi)
End Operator

Operator quad.cast ( ) As Long
    Operator = clng(this.hi)
End Operator

Operator quad.cast ( ) As Longint
    Operator = clngint(this.hi)
End Operator

Operator quad.cast ( ) As Uinteger
    Operator = cuint(this.hi)
End Operator

Operator quad.cast ( ) As ULong
    Operator = culng(this.hi)
End Operator

Operator quad.cast ( ) As Ulongint
    Operator = culngint(this.hi)
End Operator

Operator quad.cast ( ) As Single
    Operator = csng(this.hi)
End Operator

Operator quad.cast ( ) As Double
    Operator = this.hi
End Operator

Operator quad.cast ( ) As String
    Operator = quad_string( this )
End Operator

constructor quad ( Byref rhs As quad )
	this = rhs
end constructor

constructor quad ( Byref rhs As Integer )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As Long )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As LongInt )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As UInteger )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As ULong )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As ULongInt )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As Single )
    this.hi = cdbl(rhs )
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As Double )
    this.hi = rhs
    this.lo = 0.0
end constructor

constructor quad ( Byref rhs As String )
    this = string_quad ( rhs )
end constructor

Operator quad.let ( Byref rhs As quad )
    this.hi = rhs.hi
    this.lo = rhs.lo
End Operator

Operator quad.let ( Byref rhs As Integer )
    this.hi = Cdbl(rhs)
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Long )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Longint )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Uinteger )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As ULong )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Ulongint )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Single )
    this.hi = Cdbl(rhs )
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As Double )
    this.hi = rhs
    this.lo = 0.0
End Operator

Operator quad.let ( Byref rhs As String )
    this = string_quad ( rhs )
End Operator

'=========================================================
'' For next for quadended type
''
'' implicit step versions
'' 
'' In this example, we interpret implicit step
'' to mean 1
Operator quad.for( )
End Operator
 
Operator quad.step( )
        this += 1 'this = this+1 '
End Operator 
 
Operator quad.next( Byref end_cond As quad ) As Integer
        Return this <= end_cond
End Operator
 
 
'' explicit step versions
'' 
Operator quad.for( Byref step_var As quad )
End Operator
 
Operator quad.step( Byref step_var As quad )
        this += step_var 'this = this + step_var '    
End Operator 
 
Operator quad.next( Byref end_cond As quad, Byref step_var As quad ) As Integer
        If step_var < 0 Then
                Return this >= end_cond
        Else
                Return this <= end_cond
        End If
End Operator


End Namespace

[srvaldez]

here's small test

Code: Select all

#include "quad.bi"
dim as quad y
print "sin function"
for y=1 to 5
    print longsin(1/y)
next
print
print "square root function"
for y=1 to 5
    print longsqr(y)
next
print
print "log function"
for y=1 to 5
    print longlog(y)
next

print
print "pi from atn function"
y=1.0
print longatn(y)*4
print " 3.14159265358979323846264338327950..."
 
sleep

[dodicat]

here's small test

Hi srvaldez
A power of work gone into this I see.
Unfortunately I get a crash at the first for y loop.
Win XP pro. sp3
fb 21
I've downloaded twice and tried, same each time

[srvaldez]

sorry to hear it crashing on your PC, I will stitch the files that I posted and see if something got messed up.

[srvaldez]

I seem to remember when dkl released version 21 that he had included new binary tools, but apparently some people had problems with it so he uploaded FB again but with old tools, if I remember right, but since I had not experienced any problems with the new tools I kept them, maybe that's why it works my machine without problems.
by the way, here's what the test looks like.

sin function
8.41470984807896506652502321631e-1
4.79425538604203000273287935216e-1
3.27194696796152244173344085268e-1
2.47403959254522929596848704849e-1
1.98669330795061215459412627118e-1

square root function
1.00000000000000000000000000000e+0
1.41421356237309504880168872421e+0
1.73205080756887729352744634150e+0
2.00000000000000000000000000000e+0
2.23606797749978969640917366873e+0

log function
0.00
6.93147180559945309417232121458e-1
1.09861228866810969139524523692e+0
1.38629436111989061883446424292e+0
1.60943791243410037460075933323e+0

[dodicat]

Hi srvaldez
Pity it won't work here, I would like to have tried it on my matrix inverses etc. to compare rounding errors for large matrices.
Never mind.

[anonymous1337]

This is great. Good code like this shouldn't have to suffer the problems of some other code or tools.

[jdebord]

When I try to compile under Windows (with FB 0.21.1) I obtain a series of assembler errors ot the type:

Code: Select all

test.asm:4132: Error: suffix or operands invalid for `movsd'

I did not test with Linux yet.

[rvdm]

I've got de same problem as jdebord with fb 0.21.1.

[MichaelW]

When I try to compile under Windows (with FB 0.21.1) I obtain a series of assembler errors of the type…

The problem is apparently the 2.15.94 version of as (from December 01, 2005). The 2.19.1 version included with FBC 0.21.0 works, as does the 2.20 version I tested.

[jdebord]

It works in Linux with FB 0.21.1 and AS 2.20.1

Congratulations for the translation job. It's really impressive!

[srvaldez]

thanks jdebord :)
@MichaelW
your post is right on but a bit hard to understand, as = gnu assembler, so for people wanting to try out this module they need version 2.19.1 or greater of the gnu assembler.

[dodicat]

thanks jdebord :)
@MichaelW
your post it right on but a bit hard to understand, as = gnu assembler, so for people wanting to try out this module they need version 2.19.1 or greater of the gnu assembler.

Hi srvaldez
I tried your quad on a different machine
Win 2000
FreeBASIC Compiler - Version 0.21.0 (07-16-2010) for win32 (target:win32)

Still no luck.

a mass of assembler errors.
I use pentium 3's, but I'll try it on a xeon machine, on XP and Linux, but I have to rig it up, and it's heavy.