Buchempfehlung
Windows-Programmierung. Das Entwicklerhandbuch zur WIN32-API
Windows-Programmierung. Das Entwicklerhandbuch zur WIN32-API
"Der" Petzold, das über 1000 Seiten starke Standardwerk zum Win32-API - besonders nützlich u. a. bei der GUI-Programmierung in FreeBASIC! [Mehr Infos...]
FreeBASIC-Chat
Es sind Benutzer im FreeBASIC-Chat online.
(Stand:  )
FreeBASIC bei Twitter
Twitter FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!

fb:porticula NoPaste

Info
Info / Hilfe
Liste
Übersicht / Liste
Neu
Datei hochladen
Suche
Quellcode suchen
Download
Dateidownload

quad_2.bi

Uploader:RedakteurVolta
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