flat assembler
Message board for the users of flat assembler.
![]() |
Author |
|
jack2 22 Sep 2010, 18:08
Ok just for fun I was translating Alan Miller's Quad Module ( http://jblevins.org/mirror/amiller/ ) to FreeBasic and came across a problem,
in the function exactmul2 there is this line of code Code: a1 = (a - t) + t which FreeBasic trnslates to Code: mov eax, dword ptr [a] fld qword ptr [eax] fsub qword ptr [t] fadd qword ptr [t] fstp qword ptr [a1] which look OK to me, however it gives the wrong result, the following code works Code: a1 = a - t a1 = a1 + t which is translated into Code: ;## a1 = a - t mov eax, dword ptr [a] fld qword ptr [eax] fsub qword ptr [t] fstp qword ptr [a1] ;## a1 = a1 + t fld qword ptr [t] fadd qword ptr [a1] fstp qword ptr [a1] my queston is: why do they give different results? |
|||
![]() |
|
jack2 22 Sep 2010, 18:38
LocoDelAssembly I can give you an example but it's in FreeBasic, nonetheless here's the code
Code: '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" '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 End Type Declare Function longadd(Byref a As quad, Byref c As quad) As quad Function cquad(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 = ( 0.3141592653589793D+01, 0.1224646799147353D-15 ), _ piby2 = ( 0.1570796326794897D+01, -.3828568698926950D-15 ), _ piby3 = ( 0.1047197551196598D+01, -.3292527815701405D-15 ), _ piby4 = ( 0.7853981633974484D+00, -.8040613248383182D-16 ), _ piby6 = ( 0.5235987755982990D+00, -.1646263907850702D-15 ), _ twopi = ( 0.6283185307179586D+01, 0.2449293598294707D-15 ), _ ln_pi = ( 0.1144729885849400D+01, 0.2323105560877391D-15 ), _ sqrtpi = ( 0.1772453850905516D+01, -.7666586499825800D-16 ), _ fact_pt5 = ( 0.8862269254527582D+00, -.1493552349616447D-15 ), _ sqrt2pi = ( 0.2506628274631001D+01, -.6273750096546544D-15 ), _ lnsqrt2pi = ( 0.9189385332046728D+00, -.3878294158067242D-16 ), _ one_on2pi = ( 0.1591549430918953D+00, 0.4567181289366658D-16 ), _ two_on_rtpi = ( 0.1128379167095513D+01, -.4287537502368968D-15 ), _ deg2rad = ( 0.1745329251994330D-01, -.3174581724866598D-17 ), _ rad2deg = ( 0.5729577951308232D+02, -.1987849567057628D-14 ), _ ln2 = ( 0.6931471805599454D+00, -.8783183432405266D-16 ), _ ln10 = ( 0.2302585092994046D+01, -.2170756223382249D-15 ), _ log2e = ( 0.1442695040888964D+01, -.6457785410341630D-15 ), _ log10e = ( 0.4342944819032519D+00, -.5552037773430574D-16 ), _ log2_10 = ( 0.3321928094887362D+01, 0.1661617516973592D-15 ), _ log10_2 = ( 0.3010299956639812D+00, -.2803728127785171D-17 ), _ euler = ( 0.5772156649015330D+00, -.1159652176149463D-15 ), _ e = ( 0.2718281828459045D+01, 0.1445646891729250D-15 ), _ sqrt2 = ( 0.1414213562373095D+01, 0.1253716717905022D-15 ), _ sqrt3 = ( 0.1732050807568877D+01, 0.3223954471431004D-15 ), _ sqrt10 = ( 0.3162277660168380D+01, -.6348773795572286D-15 ) Function exactmul2(Byref a As Double, Byref c As Double) As quad ' 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 t = constant * a ' a1 = (a - t) + t ' FreeBasic's optimization removes the brackets ' and sets a1 = a which defeats the whole point. ' splitting the expression makes it work a1 = a - t a1 = a1 + t 'asm ''## a1 = a - t ' mov eax, dword ptr [a] ' fld qword ptr [eax] ' fsub qword ptr [t] ' fstp qword ptr [a1] ' ''## a1 = a1 + t ' fld qword ptr [t] ' fadd qword ptr [a1] ' fstp qword ptr [a1] 'end asm 'asm ''## a1 = (a - t) + t ' mov eax, dword ptr [a] ' fld qword ptr [eax] ' fsub qword ptr [t] ' fadd qword ptr [t] ' fstp qword ptr [a1] 'end asm ' asm ' mov eax, DWORD PTR [a] ' movsd xmm0, QWORD PTR [eax] ' movsd xmm1, QWORD PTR [t] ' subsd xmm0, xmm1 ' movsd xmm1, QWORD PTR [t] ' addsd xmm0, xmm1 ' movsd QWORD PTR [a1], xmm0 ' end asm 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 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. ' Local variables Dim As Double zz Dim As quad z, ac 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 Return ac End Function 'longmul Function mult_quad_int(Byref a As quad, Byval i As Integer) As quad ' 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_quad_int Function mult_int_quad(Byval i As Integer, Byref a As quad) As quad ' 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). ' Local variables 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_quad_dp Function mult_quad_real(Byref a As quad, Byref b As Double) As quad ' 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, 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 ' 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) 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 Double, Byref b As quad) As quad ' Multiply a quadruple-precision number (a) by a single-precision number (b). ' Local variables Dim As quad z, c Dim As Double zz z = exactmul2(b.hi, 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 ' 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) 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 ' 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) 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 Double) As quad 'Divide a quadruple-precision number (a) by a single-precision number (b) ' Local variables 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_real Function div_quad_int(Byref a As quad, Byval b As Integer) As quad '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) 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(Byval a As Integer, Byref c As quad) As quad ' 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) 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 Double, Byref c As quad) As quad ' 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) 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 ' 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) 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. ' Local variables 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 'longadd Function quad_add_int(Byref a As quad, Byval c As Integer) As quad ' Local variables 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 Double) As quad ' Local variables 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 Function quad_add_dp(Byref a As quad, Byref c As Double) As quad ' Local variables 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(Byval c As Integer, Byref a As quad) As quad ' Local variables 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 Double, Byref a As quad) As quad ' Local variables 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 ' Local variables 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 ' ' 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 ' ' 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 longsub(Byref a As quad, Byref c As quad) As quad ' 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 quad ac ' Local variables Dim As Double z, q, zz, temp ' asm z = a.hi - c.hi ' mov eax, DWORD PTR [a] ' mov edx, DWORD PTR [c] ' movsd xmm0, QWORD PTR [eax] ' movsd xmm1, QWORD PTR [edx] ' subsd xmm0, xmm1 ' movsd QWORD PTR [z], xmm0 q = a.hi - z ' mov eax, DWORD PTR [a] ' movsd xmm0, QWORD PTR [eax] ' movsd xmm1, QWORD PTR [z] ' subsd xmm0, xmm1 ' movsd QWORD PTR [q], xmm0 ' zz = (((q - c.hi) + (a.hi - (q + z))) + a.lo) - c.lo 'line above rewritten in simple sub-expressions to make it work in FB temp = q + z temp = a.hi - temp zz = q - c.hi zz = zz + temp zz = zz + a.lo zz = zz - c.lo 'asm ' mov eax, DWORD Ptr [c] ' movsd xmm0, QWORD Ptr [q] ' movsd xmm1, QWORD Ptr [eax] ' subsd xmm0, xmm1 ' mov eax, DWORD Ptr [a] ' movsd xmm1, QWORD Ptr [q] ' movsd xmm2, QWORD Ptr [z] ' addsd xmm1, xmm2 ' movsd xmm2, QWORD Ptr [eax] ' subsd xmm2, xmm1 ' addsd xmm0, xmm2 ' mov eax, DWORD Ptr [a] ' movsd xmm1, QWORD Ptr [8+eax] ' addsd xmm0, xmm1 ' mov eax, DWORD Ptr [c] ' movsd xmm1, QWORD Ptr [8+eax] ' subsd xmm0, xmm1 ' movsd QWORD Ptr [zz], xmm0 'End asm ac.hi = z + zz ' movsd xmm0, QWORD PTR [z] ' movsd xmm1, QWORD PTR [zz] ' addsd xmm0, xmm1 ' lea eax, DWORD PTR [ac] ' movsd QWORD PTR [eax], xmm0 ac.lo = (z - ac.hi) + zz ' lea eax, DWORD PTR [ac] ' movsd xmm0, QWORD PTR [z] ' movsd xmm1, QWORD PTR [eax] ' subsd xmm0, xmm1 ' movsd xmm1, QWORD PTR [zz] ' addsd xmm0, xmm1 ' lea eax, DWORD PTR [ac] ' movsd QWORD PTR [8+eax], xmm0 ' ' end asm Return ac End Function Function quad_sub_int(Byref a As quad, Byval c As Integer) As quad ' Local variables 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, Byval c As Double) As quad ' Local variables 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 ' Local variables 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(Byval a As Integer, Byref c As quad) As quad ' Local variables Dim As Double z, q, zz Dim As quad ac z = Cdbl(a) - c.hi q = Cdbl(a) - z zz = ((q - c.hi) + (Cdbl(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 Double, Byref c As quad) As quad ' Local variables 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 ' Local variables 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 ' 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 '############################################################################## '############################################################################## 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 pow(Byref x As quad, Byval 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 (n Mod 2)=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 x as quad, byval j as integer) as quad ' dim as quad y ' Asm ' mov eax,DWORD PTR [x] ' lea edx,DWORD PTR [y] ' fild dword ptr [j] ' fld qword ptr [eax] 'x ' fscale ' fstp qword ptr [edx] 'y ' fstp st(0) 'clean fpu stack ' fild dword ptr [j] ' fld qword ptr [eax+8] 'x ' fscale ' fstp qword ptr [edx+8] 'y ' fstp st(0) 'clean fpu stack ' End Asm ' return y 'End Function Function qscale(Byref a As quad, Byval i As Integer ) As quad ' Multiply a by 2^i Dim As quad b b=cquad(2,0) b=pow(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 qscale(byref a as quad, byval i as integer ) as quad ' ' Multiply a by 2^i ' ' dim as quad 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 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' ' Local variables Dim As quad temp, ysq, sum1, sum2, y Dim As Integer iy y = longdiv(x , ln2) iy = NINT(y.hi) y = longmul(quad_sub_dp(y , Cdbl(iy)) , 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) sum1 = longmul(sum1,ysq) sum1 = quad_add_dp(sum1,2162160) sum1 = longmul(sum1,ysq) sum1 = quad_add_dp(sum1,302702400) sum1 = longmul(sum1,ysq) sum1 = quad_add_dp(sum1,8821612800) sum1 = longmul(y,sum1) 'sum1 = y * ((((ysq + 3960.)*ysq + 2162160._dp)*ysq + 302702400._dp)*ysq + _ ' 8821612800._dp) sum2 = mult_int_quad(90,ysq) sum2 = quad_add_dp(sum2,110880) sum2 = longmul(sum2,ysq) sum2 = quad_add_dp(sum2,30270240) sum2 = longmul(sum2,ysq) sum2 = quad_add_dp(sum2,2075673600) sum2 = longmul(sum2,ysq) sum2 = quad_add_dp(sum2,17643225600) 'sum2 = (((90.*ysq + 110880.)*ysq + 30270240._dp)*ysq + 2075673600._dp)*ysq + _ ' 17643225600._dp ' 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) y = longmul(y,temp) y = qSCALE(y, 2) y = quad_add_dp(y,1) y = qSCALE(y, iy) Return y End Function 'longexp 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_string(Byref value As quad, Byval ier As Integer) 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 Dim As Double tmp 'IF (LEN(strng) < 36) THEN ' ier = 1 ' strng = "***" ' RETURN strng 'END IF 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 = Int( 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,Cdbl(14 - dec_expnt)) v=longexp(v) '************* 'HACK v.lo = 0.0 '************* vl = longmul(vl , v ) End If str1=format( vl.hi,"################") ' Calculate the remainder tmp =Val(str1+".0#") vl = quad_sub_real(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 If vl.hi<1 Then str2 = "000000000000000" Else str2=format( vl.hi,"################") 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 ' ' If str2 consists of asterisks, add 1 in the last digit to str1. ' ' Set str2 to zeroes. ' IF (str2(2:2) = '*') THEN ' tmp = tmp + one ' WRITE(str1, '(f17.0)') tmp ' IF (str1(1:1) /= ' ') THEN ' dec_expnt = dec_expnt + 1 ' ELSE ' str1 = str1(2:17) ' END IF ' str2 = '000000000000000.' ' 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 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 '(es*ex)-(58-i)+e f1=f1+f2 'f1=Mid(f1,1,1)+"."+Right(f1,Len(f1)-1) ' If s<0 then ' f1="-"+f1 ' Else ' f1=" "+f1 ' EndIf 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)) qd=div_quad_dp(qd,1000000000000000) qd=div_quad_dp(qd,1000000000000000) ' f1=f1+f2 ' fb.fs=f1 ' n.bcd(0) = 0 ' If fb.fc(0)=asc("-") then ' n.bcd(0)=&h80 ' EndIf ' If fb.fc(61)=asc("-") then ' n.bcd(0)=n.bcd(0) or &h40 ' EndIf ' n.bcd(1)=((fb.fc(62)-asc("0")) Shl 4) or (fb.fc(63)-asc("0")) ' n.bcd(2)=((fb.fc(64)-asc("0")) Shl 4) or (fb.fc(65)-asc("0")) ' n.bcd(3)=((fb.fc(1)-asc("0")) Shl 4) or (fb.fc(3)-asc("0")) ' For i=4 To 31 ' n.bcd(i)=((fb.fc(2*(i-4)+4)-asc("0")) Shl 4) or (fb.fc(2*(i-4)+5)-asc("0")) ' Next Return qd End Function '############################################################################## Dim As quad x, y Dim As String s Dim As Integer i, j Dim As Double t, t1, dy, dx = 3.1415926535 x=cquad(2,0) 'y=pow(x,3) y=longlog(x) 'y=longexp(ln2) ?quad_string(y, i) ?quad_string(ln2, i) Sleep in the function exactmul2 I worked around the problem by spliting up a1 = (a - t) + t into a1 = a - t and a1 = a1 + t but am at a loss as to why the first expression does not work and the second one does. you can see my post at the FreeBasic Forum here http://www.freebasic.net/forum/viewtopic.php?t=16450 |
|||
![]() |
|
jack2 22 Sep 2010, 19:57
OK, so how do I set the FPU control word to Double Precision?
because if I change Code: mov eax, dword ptr [a] fld qword ptr [eax] fsub qword ptr [t] fadd qword ptr [t] fstp qword ptr [a1] to Code: '## a1 = (a - t) + t mov eax, dword ptr [a] fld qword ptr [eax] fsub qword ptr [t] fstp qword ptr [a1] fld qword ptr [a1] fadd qword ptr [t] fstp qword ptr [a1] it works as expected, this is weird, |
|||
![]() |
|
baldr 22 Sep 2010, 20:13
jack2,
Store/load operations effectively limit precision to 53 bits. fldcw instruction can be used to load x87 Control Word. |
|||
![]() |
|
jack2 22 Sep 2010, 20:37
yes, I know that much, but what is the value I need to store in the controlword to force it to double?
|
|||
![]() |
|
baldr 22 Sep 2010, 20:47
jack2,
That depends. You may change only precision though: fstcw can be useful. |
|||
![]() |
|
LocoDelAssembly 22 Sep 2010, 22:20
0x027F and 0x037F are WinXP and finit defaults respectively.
Here you can see a quick control word reference. PS: In case isn't obvious, WinXP's default is double precision and finit is extended precision. |
|||
![]() |
|
jack2 23 Sep 2010, 00:20
thanks, I examined the control word just before doing the arithmetic and it was 0x37F.
just did a test run, setting the control word to 0x27F gives the expected result. ![]() |
|||
![]() |
|
< Last Thread | Next Thread > |
Forum Rules:
|
Copyright © 1999-2025, Tomasz Grysztar. Also on GitHub, YouTube.
Website powered by rwasa.