Wer ist online?
Buchempfehlung
FreeBASIC-Chat
FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!fb:porticula NoPaste
quad_1x.bi
| Uploader: |  Volta |
| Datum/Zeit: | 30.09.2010 20:38:35 |
'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
'http://www.cmis.csiro.au/Alan_Miller/quad.html (quad.f90)
#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.
Dim As quad ac
Dim As Double a1, a2, c1, c2, t
Dim As UShort ocw, ncw = &h27f
Asm
fstcw Word [ocw]
fldcw Word [ncw]
't = constant * a
mov eax, Dword Ptr [a]
fld qword Ptr [eax]
fmul qword Ptr [constant]
fstp qword Ptr [t]
'a1 = (a - t) + t ' Lahey's optimization removes the brackets
' ' and sets a1 = a which defeats the whole point.
'mov eax, dword ptr [a]
fld qword Ptr [eax]
fsub qword Ptr [t]
fadd qword Ptr [t]
fstp qword Ptr [a1]
'a2 = a - a1
'mov eax, dword ptr [a]
fld qword Ptr [eax]
fsub qword Ptr [a1]
fstp qword Ptr [a2]
't = constant * c
mov eax, Dword Ptr [c]
fld qword Ptr [eax]
fmul qword Ptr [constant]
fstp qword Ptr [t]
'c1 = (c - t) + t
'mov eax, dword ptr [c]
fld qword Ptr [eax]
fsub qword Ptr [t]
fadd qword Ptr [t]
fstp qword Ptr [c1]
'c2 = c - c1
'mov eax, dword ptr [c]
fld qword Ptr [eax]
fsub qword Ptr [c1]
fstp qword Ptr [c2]
'ac.hi = a * c
'mov eax, dword ptr [c]
fld qword Ptr [eax]
mov eax, Dword Ptr [a]
fmul qword Ptr [eax]
fstp qword Ptr [ac]
'ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2
fld qword Ptr [c1]
fmul qword Ptr [a1]
fsub qword Ptr [ac]
fld qword Ptr [c2]
fmul qword Ptr [a1]
fxch st(1)
faddp
fld qword Ptr [a2]
fmul qword Ptr [c1]
fxch st(1)
faddp
fld qword Ptr [a2]
fmul qword Ptr [c2]
fxch st(1)
faddp
fstp qword Ptr [ac+8]
fldcw Word [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
Dim As Double zz
z = exactmul2(a.hi, c.hi)
zz = ((a.hi + a.lo) * c.lo + a.lo * c.hi) + z.lo
ac.hi = z.hi + zz
ac.lo = (z.hi - ac.hi) + zz
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
Dim As Double zz
z = exactmul2(a.hi, b)
zz = a.lo * b + z.lo
c.hi = z.hi + zz
c.lo = (z.hi - c.hi) + zz
Return c
End Function
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).
Dim As quad z, c
Dim As Double zz
z = exactmul2(a.hi, CDbl(b))
zz = a.lo * b + z.lo
c.hi = z.hi + zz
c.lo = (z.hi - c.hi) + zz
Return c
End Function 'mult_quad_Real
Function mult_dp_quad(ByRef b As Double, ByRef a As quad) As quad 'Result(c)
' Multiply a quadruple-precision number (a) by a double-precision number (b).
Dim As quad z, c
Dim As Double zz
z = exactmul2(a.hi, b)
zz = a.lo * b + z.lo
c.hi = z.hi + zz
c.lo = (z.hi - c.hi) + zz
Return c
End Function 'mult_dp_quad
Function mult_Real_quad(ByRef a As Single, ByRef b As quad) As quad 'Result(c)
' Multiply a quadruple-precision number (a) by a double-precision number (b).
Dim As quad z, c
Dim As Double zz
z = exactmul2(b.hi, CDbl(a))
zz = b.lo * a + z.lo
c.hi = z.hi + zz
c.lo = (z.hi - c.hi) + zz
Return c
End Function 'mult_Real_quad
Function longdiv(ByRef a As quad, ByRef c As quad) As quad '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.
Dim As Double z, zz
Dim As quad q, ac
z = a.hi / c.hi
q = exactmul2(c.hi, z)
zz = ((((a.hi - q.hi) - q.lo) + a.lo) - z*c.lo) / (c.hi + c.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'longdiv
Function div_quad_dp(ByRef a As quad, ByRef b As Double) As quad 'Result(c)
' Divide a quadruple-precision number (a) by a double-precision number (b)
Dim As Double z, zz
Dim As quad q, c
z = a.hi / b
q = exactmul2(b, z)
zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
c.hi = z + zz
c.lo = (z - c.hi) + zz
Return c
End Function 'div_quad_dp
Function div_quad_Real(ByRef a As quad, ByRef b As Single) As quad'Result(c)
' Divide a quadruple-precision number (a) by a Real number (b)
Dim As Double z, zz
Dim As quad q, c
z = a.hi / b
q = exactmul2(CDbl(b), z)
zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
c.hi = z + zz
c.lo = (z - c.hi) + zz
Return c
End Function 'div_quad_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)
Dim As Double z, zz
Dim As quad q, c
z = a.hi / b
q = exactmul2(CDbl(b), z)
zz = (((a.hi - q.hi) - q.lo) + a.lo) / b
c.hi = z + zz
c.lo = (z - c.hi) + zz
Return c
End Function 'div_quad_int
Function div_int_quad(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.
Dim As Double z, zz
Dim As quad q, ac
z = CDbl(a) / c.hi
q = exactmul2(c.hi, z)
zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'div_int_quad
Function div_Real_quad(ByRef a As 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.
Dim As Double z, zz
Dim As quad q, ac
z = CDbl(a) / c.hi
q = exactmul2(c.hi, z)
zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'div_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.
Dim As Double z, zz
Dim As quad q, ac
z = a / c.hi
q = exactmul2(c.hi, z)
zz = (((a - q.hi) - q.lo) - z*c.lo) / (c.hi + c.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'div_dp_quad
Function longadd(ByRef a As quad, ByRef c As quad) As quad
' procedure longadd, translated from pascal, from:
' linnainmaa, seppo (1981). software for doubled-precision floating-point
' computations. acm trans. on math. software (toms), 7, 272-283.
Dim As quad ac
Dim As Double z, q, zz
z = a.hi + c.hi
q = a.hi - z
zz = (((q + c.hi) + (a.hi - (q + z))) + a.lo) + c.lo
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function
Function quad_add_int(ByRef a As quad, ByRef c As Integer) As quad 'Result(ac)
Dim As Double z, q, zz
Dim As quad ac
z = a.hi + c
q = a.hi - z
zz = (((q + c) + (a.hi - (q + z))) + a.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'quad_add_int
Function quad_add_Real(ByRef a As quad, ByRef c As Single) As quad 'Result(ac)
Dim As Double z, q, zz
Dim As quad ac
z = a.hi + c
q = a.hi - z
zz = (((q + c) + (a.hi - (q + z))) + a.lo)
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'quad_add_Real
#Include "quad_2.bi"
#Include "quad_3.bi"
#Include "quad_4.bi"