Quad Precision
you can download the new binutils from http://sourceforge.net/downloads/mingw/ ... .tar.lzma/
change to your FreeBasic bin\win32 folder and rename your current version of as.exe to something else then copy the new version from the download.
change to your FreeBasic bin\win32 folder and rename your current version of as.exe to something else then copy the new version from the download.
dodicat wrote: Unfortunately I get a crash at the first for y loop.
Using a later version of as will correct the assembler errors, but when the code is run on a processor without SSE2 support it will trigger an illegal-instruction exception and crash. SSE2 support was introduced with the Intel P4 and the AMD Opteron/Athlon 64. The critical difference here between the original SSE and SSE2 is support for 64-bit floating-point values.I use pentium 3's, but I'll try it on a xeon machine
@dodicat
here a version that should work for you (I hope)
part 1 of 4 "quad.bi"
here a version that should work for you (I hope)
part 1 of 4 "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
Dim as Ushort ocw, ncw = &h27f
' 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
fstcw word ptr [ocw]
fldcw word ptr [ncw]
fld qword ptr [constant]
mov eax, dword ptr [a]
fld qword ptr [eax]
fmulp st(1), st
fstp qword ptr [t]
' a1 = (a - t) + t
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [t]
fsubp st(1), st
fld qword ptr [t]
faddp st(1), st
fstp qword ptr [a1]
' a2 = a - a1
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [a1]
fsubp st(1), st
fstp qword ptr [a2]
' t = constant * c
fld qword ptr [constant]
mov eax, dword ptr [c]
fld qword ptr [eax]
fmulp st(1), st
fstp qword ptr [t]
' c1 = (c - t) + t
mov eax, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [t]
fsubp st(1), st
fld qword ptr [t]
faddp st(1), st
fstp qword ptr [c1]
' c2 = c - c1
mov eax, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [c1]
fsubp st(1), st
fstp qword ptr [c2]
' ac.hi = a * c
mov eax, dword ptr [a]
fld qword ptr [eax]
mov eax, dword ptr [c]
fld qword ptr [eax]
fmulp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2
fld qword ptr [a1]
fld qword ptr [c1]
fmulp st(1), st
lea eax, dword ptr [ac]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [a1]
fld qword ptr [c2]
fmulp st(1), st
faddp st(1), st
fld qword ptr [c1]
fld qword ptr [a2]
fmulp st(1), st
faddp st(1), st
fld qword ptr [c2]
fld qword ptr [a2]
fmulp st(1), st
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
fldcw word ptr [ocw]
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
' zz = ((a.hi + a.lo) * c.lo + a.lo * c.hi) + z.lo
mov eax, dword ptr [a]
mov edx, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [8+edx]
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fmulp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
mov eax, dword ptr [c]
fld qword ptr [eax]
fmulp st(1), st
faddp st(1), st
fld qword ptr [z+8]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z.hi + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z.hi - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [8+eax]
mov eax, dword ptr [b]
fld qword ptr [eax]
fmulp st(1), st
fld qword ptr [z+8]
faddp st(1), st
fstp qword ptr [zz]
' c.hi = z.hi + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z.hi - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [8+eax]
mov eax, dword ptr [b]
fld dword ptr [eax]
fmulp st(1), st
fld qword ptr [z+8]
faddp st(1), st
fstp qword ptr [zz]
' c.hi = z.hi + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z.hi - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [8+eax]
mov eax, dword ptr [b]
fld qword ptr [eax]
fmulp st(1), st
fld qword ptr [z+8]
faddp st(1), st
fstp qword ptr [zz]
' c.hi = z.hi + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z.hi - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [8+eax]
mov eax, dword ptr [a]
fld dword ptr [eax]
fmulp st(1), st
fld qword ptr [z+8]
faddp st(1), st
fstp qword ptr [zz]
' c.hi = z.hi + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z.hi - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fld qword ptr [z]
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fmulp st(1), st
fsubp st(1), st
mov eax, dword ptr [c]
mov edx, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [8+edx]
faddp st(1), st
fdivp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
mov eax, dword ptr [b]
fld qword ptr [eax]
fdivp st(1), st
fstp qword ptr [zz]
' c.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
mov eax, dword ptr [b]
fld dword ptr [eax]
fdivp st(1), st
fstp qword ptr [zz]
' c.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
mov eax, dword ptr [b]
mov eax, dword ptr [eax]
push eax
fild dword ptr [esp]
pop eax
fdivp st(1), st
fstp qword ptr [zz]
' c.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [eax]
' c.lo = (z - c.hi) + zz
lea eax, dword ptr [c]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [c]
fstp qword ptr [8+eax]
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]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
fld qword ptr [z]
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fmulp st(1), st
fsubp st(1), st
mov eax, dword ptr [c]
mov edx, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [8+edx]
faddp st(1), st
fdivp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld dword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
fld qword ptr [z]
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fmulp st(1), st
fsubp st(1), st
mov eax, dword ptr [c]
mov edx, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [8+edx]
faddp st(1), st
fdivp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [q]
fsubp st(1), st
fld qword ptr [q+8]
fsubp st(1), st
fld qword ptr [z]
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fmulp st(1), st
fsubp st(1), st
mov eax, dword ptr [c]
mov edx, dword ptr [c]
fld qword ptr [eax]
fld qword ptr [8+edx]
faddp st(1), st
fdivp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [edx]
faddp st(1), st
fstp qword ptr [z]
' q = a.hi - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = (((q + c.hi) + (a.hi - (q + z))) + a.lo) + c.lo
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [q]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld dword ptr [eax]
fld qword ptr [q]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
End asm
Return ac
End Function 'quad_add_Real
Last edited by srvaldez on Oct 02, 2010 3:21, edited 4 times in total.
part 2 of 4 "quad.bi"
Code: Select all
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]
fld qword ptr [q]
fld qword ptr [eax]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [q]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld dword ptr [eax]
fld qword ptr [q]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
End asm
Return ac
End Function 'Real_add_quad
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]
fld qword ptr [q]
fld qword ptr [eax]
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [edx]
fsubp st(1), st
fstp qword ptr [z]
' q = a.hi - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = (((q - c.hi) + (a.hi - (q + z))) + a.lo) - c.lo
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fsubp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
push edx
fild dword ptr [esp]
pop edx
fld qword ptr [eax]
fsubrp st(1), st
fstp qword ptr [z]
' q = a.hi - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = (((q - c) + (a.hi - (q + z))) + a.lo)
mov eax, dword ptr [c]
mov eax, dword ptr [eax]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [q]
fsubrp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld dword ptr [edx]
fld qword ptr [eax]
fsubrp st(1), st
fstp qword ptr [z]
' q = a.hi - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = (((q - c) + (a.hi - (q + z))) + a.lo)
mov eax, dword ptr [c]
fld dword ptr [eax]
fld qword ptr [q]
fsubrp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [edx]
fsubp st(1), st
fstp qword ptr [z]
' q = a.hi - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = (((q - c) + (a.hi - (q + z))) + a.lo)
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
mov dword ptr [-12+ebp], eax
fild dword ptr [-12+ebp]
mov eax, dword ptr [c]
fld qword ptr [eax]
fsubp st(1), st
fstp qword ptr [z]
' q = dble(a) - z
mov eax, dword ptr [a]
mov eax, dword ptr [eax]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = ((q - c.hi) + (dble(a) - (q + z))) - c.lo
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
fsubp st(1), st
mov eax, dword ptr [a]
mov eax, dword ptr [eax]
push eax
fild dword ptr [esp]
pop eax
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fsubp st(1), st
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fsubp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld dword ptr [eax]
mov eax, dword ptr [c]
fld qword ptr [eax]
fsubp st(1), st
fstp qword ptr [z]
' q = a - z
mov eax, dword ptr [a]
fld dword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = ((q - c.hi) + (a - (q + z))) - c.lo
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
fsubp st(1), st
mov eax, dword ptr [a]
fld dword ptr [eax]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fsubp st(1), st
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fsubp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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]
fld qword ptr [eax]
fld qword ptr [edx]
fsubp st(1), st
fstp qword ptr [z]
' q = a - z
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [z]
fsubp st(1), st
fstp qword ptr [q]
' zz = ((q - c.hi) + (a - (q + z))) - c.lo
mov eax, dword ptr [c]
fld qword ptr [q]
fld qword ptr [eax]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [q]
fld qword ptr [z]
faddp st(1), st
fld qword ptr [eax]
fsubrp st(1), st
faddp st(1), st
mov eax, dword ptr [c]
fld qword ptr [8+eax]
fsubp st(1), st
fstp qword ptr [zz]
' ac.hi = z + zz
fld qword ptr [z]
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [eax]
' ac.lo = (z - ac.hi) + zz
lea eax, dword ptr [ac]
fld qword ptr [z]
fld qword ptr [eax]
fsubp st(1), st
fld qword ptr [zz]
faddp st(1), st
lea eax, dword ptr [ac]
fstp qword ptr [8+eax]
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
function qsgn(byref a as quad) as integer
return sgn(a.hi)
end function
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]
fld qword ptr [eax]
fld qword ptr [tt]
fsubp st(1), st
fld qword ptr [tt+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fld qword ptr [_longSqr_half_dp]
fmulp st(1), st
fld qword ptr [t]
fdivp st(1), st
fld qword ptr [t]
faddp st(1), st
fstp qword ptr [res]
' b.lo = (t - res) + (((a.hi - tt.hi) - tt.lo) + a.lo) * 0.5_dp / t
fld qword ptr [t]
fld qword ptr [res]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [eax]
fld qword ptr [tt]
fsubp st(1), st
fld qword ptr [tt+8]
fsubp st(1), st
mov eax, dword ptr [a]
fld qword ptr [8+eax]
faddp st(1), st
fld qword ptr [_longSqr_half_dp]
fmulp st(1), st
fld qword ptr [t]
fdivp st(1), st
faddp st(1), st
lea eax, dword ptr [b]
fstp qword ptr [8+eax]
' b.hi = res
lea eax, dword ptr [b]
fld qword ptr [res]
fstp qword ptr [eax]
End asm
Return b
asm
_longSqr_half_dp: .int &h000000000,&h03fe00000
End asm
End Function 'longsqrt
Function FBfloor(Byref x as double) as integer
'Fast Rounding of Floating Point Numbers
'in C/C++ on Wintel Platform
'Laurent de Soras
'Updated on 2008.03.08
'web: http://ldesoras.free.fr
dim as integer N
dim as single round_towards_m_i = -0.5
dim as longint tmp
Asm
mov eax, dword ptr [x]
fld qword ptr [eax]
fadd st(0),st(0)
fadd dword ptr [round_towards_m_i]
fistp qword ptr [tmp]
mov eax, dword ptr [tmp]
sar dword ptr [tmp+4],1
rcr eax,1
mov dword ptr [N], eax
End Asm
return N
end function
Last edited by srvaldez on Oct 02, 2010 3:23, edited 3 times in total.
part 3 of 4 "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
Last edited by srvaldez on Oct 02, 2010 3:27, edited 4 times in total.
part 4 of 4 "quad.bi"
Code: Select all
'************************************************************
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
Last edited by srvaldez on Oct 02, 2010 3:29, edited 3 times in total.
And for me, as the most recent system I have here is also a P3. Thank you for going to the trouble.srvaldez wrote:@dodicat
here a version that should work for you (I hope)
It still will not assemble with as 2.15.94 (although looking at the assembly code I cannot see any problem):
Code: Select all
test.asm: Assembler messages:
test.asm:64: Error: suffix or operands invalid for `fstcw'
test.asm:65: Error: suffix or operands invalid for `fldcw'
test.asm:127: Error: suffix or operands invalid for `fldcw'
But with as 2.19.1 (from FB version 0.21.0 07-16-2010) and 2.20.51 (from here) it assembles OK and appears to run correctly. I added another simple test to the code:
Code: Select all
print
print "pi from atn function"
y=1.0
print longatn(y)*4
print " 3.14159265358979323846264338327950..."
Code: Select all
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.73205080756887729352744634151e+0
2.00000000000000000000000000000e+0
2.23606797749978969640917366873e+0
log function
0.00
6.93147180559945309417232121458e-1
1.09861228866810969139524523692e+0
1.38629436111989061883446424292e+0
1.60943791243410037460075933323e+0
pi from atn function
3.14159265358979323846264338328e+0
3.14159265358979323846264338327950...
Both versions of quad.bi have a problem working with 10^n when n is even:
Results:
1.000000000000000000000000000000e+1 1.000000000000000000000000000000e+0
9.99999999999999999999999999999e+2 3.16227766016837933199889354443e+1
1.000000000000000000000000000000e+3 1.000000000000000000000000000000e+1
1.000000000000000000000000000000e+4 3.16227766016837933199889354443e+2
1.000000000000000000000000000000e+5 9.99999999999999999999999999999e+2
Code: Select all
Print cquad(Cast(Double,100)),longSqr(cquad(Cast(Double,100)))
Print cquad(Cast(Double,1000)),longSqr(cquad(Cast(Double,1000)))
Print cquad(Cast(Double,10000)),longSqr(cquad(Cast(Double,10000)))
Print cquad(Cast(Double,100000)),longSqr(cquad(Cast(Double,100000)))
Print cquad(Cast(Double,1000000)),longSqr(cquad(Cast(Double,1000000)))
1.000000000000000000000000000000e+1 1.000000000000000000000000000000e+0
9.99999999999999999999999999999e+2 3.16227766016837933199889354443e+1
1.000000000000000000000000000000e+3 1.000000000000000000000000000000e+1
1.000000000000000000000000000000e+4 3.16227766016837933199889354443e+2
1.000000000000000000000000000000e+5 9.99999999999999999999999999999e+2
Thanks srvaldez, for your trouble.srvaldez wrote:the problem was in the function quad_string , due to the fact that the FB function Int does not return the floor of the value as is stated in the manual.
I updated the code to use the crt function floor and it seems to be working OK, the posts are updated.
It's working fine now on my p3.
I get identical results to Michaelw.
Job well done.
I'll test it out with my matrices, inverses and determinants now.
I am surprised because I am convinced that the FreeBasic function "INT(number)" returns the floor of number.srvaldez wrote:the problem was in the function quad_string , due to the fact that the FB function Int does not return the floor of the value as is stated in the manual.
I updated the code to use the crt function floor and it seems to be working OK, the posts are updated.
Remark : "-INT(-number)" returns the ceil of number.