Wer ist online?
Buchempfehlung
FreeBASIC-Chat
FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!fb:porticula NoPaste
quad_2.bi
| Uploader: |  Volta |
| Datum/Zeit: | 30.09.2010 20:37:40 |
Function quad_add_dp(Byref a As quad, Byref c As Double) 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_dp
Function int_add_quad(Byref c As Integer, Byref a As quad) 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 'int_add_quad
Function Real_add_quad(Byref c As Single, Byref a As quad) 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 'Real_add_quad
Function dp_add_quad(Byref c As Double, Byref a As quad) 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 '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.
Dim As Double z, q, zz
Dim As quad ac
z = a.hi - c.hi
q = a.hi - z
zz = (((q - c.hi) + (a.hi - (q + z))) + a.lo) - c.lo
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'longsub
Function quad_sub_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_sub_int
Function quad_sub_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_sub_Real
Function quad_sub_dp(Byref a As quad, Byref c As Double) 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_sub_dp
Function int_sub_quad(Byref a As Integer, Byref c As quad) As quad 'Result(ac)
Dim As Double z, q, zz
Dim As quad ac
z = a - c.hi
q = a - z
zz = ((q - c.hi) + (a - (q + z))) - c.lo
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'int_sub_quad
Function Real_sub_quad(Byref a As Single, Byref c As quad) As quad'Result(ac)
Dim As Double z, q, zz
Dim As quad ac
z = a - c.hi
q = a - z
zz = ((q - c.hi) + (a - (q + z))) - c.lo
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'Real_sub_quad
Function dp_sub_quad(Byref a As Double, Byref c As quad) As quad 'Result(ac)
Dim As Double z, q, zz
Dim As quad ac
z = a - c.hi
q = a - z
zz = ((q - c.hi) + (a - (q + z))) - c.lo
ac.hi = z + zz
ac.lo = (z - ac.hi) + zz
Return ac
End Function 'dp_sub_quad
Function negate_quad(Byref a As quad) As quad '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
End Sub 'quad_eq_int
Sub quad_eq_Real(Byref a As quad, Byref r As Single)
' Assignment
a.hi = r
a.lo = zero
End Sub 'quad_eq_Real
Sub quad_eq_dp(Byref a As quad, Byref d As Double)
' Assignment
a.hi = d
a.lo = zero
End Sub 'quad_eq_dp
Sub int_eq_quad(Byref i As Integer, Byref a As quad)
' Assignment
i = a.hi
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
End Sub 'Real_eq_quad
Sub dp_eq_quad(Byref d As Double, Byref a As quad)
' Assignment
d = a.hi
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
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
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
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
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
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(one, zero)
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(one,zero)
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,zero)
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, one)
y = longmul(temp , y)
y = qscale(y, 2)
y = quad_add_dp(y , one)
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 = zero
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.
Dim As Double t, res
Dim As quad b, tt
' Check that ahi >= 0.
If (a.hi < zero) Then
Print " *** Negative argument for longsqrt ***"
Return b
Elseif (a.hi = zero) Then
b.hi = zero
b.lo = zero
Return b
End If
' First approximation is t = Sqr(a).
t = Sqr(a.hi)
tt = exactmul2(t, t)
res = t + (((a.hi - tt.hi) - tt.lo) + a.lo) * 0.5 / t
b.lo = (t - res) + (((a.hi - tt.hi) - tt.lo) + a.lo) * 0.5 / t
b.hi = res
Return b
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