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