In my edit of 1/28/09 I added the function Str96A() to do much the same as Str$(), added more options to String96A(), and added Tweak96A() to make small changes to the significand. Also one-step routines Subtract96A() and Divide96A() were added to simplify those calculations. My edit of 1/29/09 corrected a fatal error in the Type96 structure and simplified member names. The look forward reference capability of PBCC5 does not apply to UDT structures so Type/End Type statements must still appear before their reference. The equates and format description are also moved here to a more logical place. I also added a sub CheckOut96A() for internal use to check the output array dimensions, made a minor revision to Case 3 in the String96A(), and added missing array checks to the start of Tweak96A() and Val96A() functions. Sign96A() had a couple typos left over from when it initially only checked for equality. Some rearranging of statements was done to put the checks first. My edit of 2/5/09 corrected an error in Val96A() when evaluating denormals. Also the entire text was reread to make numerous revisions. My 2/9/09 edit modified Val96A() to allow optional use of PowerOfTen96A() funtion (included in my post of the same date). See the statements immediately after the IncludePowerSign label. Also added are PowerOfTwo96T() and PowerOfTwo96A() functions that return an integer power of two.
This suite of functions has passed all my tests. But this must be considered as a work in progress and is subject to revision. The various blocks of code should be combined and put in an include file; I have separated them here so they can be more easily viewed. If required, e.g. in PBCC4, either add Delare statements or rearrange.
I had been posting some notes I have made on various PowerBASIC functions on the PB Console Compiler forum more or less as an academic exercise. To have baseline data with which to work, I needed arrays of a data that represented the most accurate extended-precision representations of numbers, rather like an enhanced Val() function. To do this initially I wrote a routine that doubled or halved a string of decimal digits to convert them to binary and then put that into the extended-precision format. Without any restriction on string length, this made the process very slow. The original strings would come from web, the Microsoft calculator included in Windows, or from an
expression such as "1E"+Format$(n). While preparing a post on the Exp() function, I recognized that a routine to multiply floating point numbers would be needed, and that is how these functions were conceived.
One of the posts mentioned above provides information on the 80-bit extended precision format. These functions use an enhanced version of that format called by Intel double extended-precision. This enhanced format adds (on the low-memory end) 32 bits to the 64-bit significand and adds (on the high-memory end) 16 bits to the 16-bit sign/exponent. The 96-bit significand (hence "96" in the naming) is equivalent to
almost 29 significant decimal digits. The 32 bits add ample extra precision to avoid accumulation of rounding errors and the other 16 bits contain a few flags to simplify the detection of zero, infinity and denormals. Denormals are those numbers very near zero that have reduced precision and use, depending on how one interprets it, either
an exponent bias of &H3FFE instead of the standard &H3FFF or else an implied binary point one bit further into the significand. Denormals have very limited support in PowerBASIC (they are, for example, printed as zero) but can be created and manipulated in the Floating Point Unit (FPU) that does extended-precision arithmetic.
This gives a total size of 128 bits (16 bytes) to this enhanced format with a 96-bit significand. All the original development of the functions was done dividing these 16 bytes into arrays of four double words dimensioned 0 to 3. The functions that use these arrays have "96A" as a suffix to their names. One can equally well use a user-defined type (UDT) and a Type structure and a set of functions with suffix "96T" are given below. Also provided are the flag locations in equate statements. There flag locations in bits 112 to 127 can be changed to others at will; those below 112 are tied to the
extended-precision format. The Type structure and equate definitions must precede their actual usage. The code was written for PowerBASIC Console Compiler (PBCC) Version 5.01. The functions should work with earlier PBCC versions with Declare statements or rearranging. The functions should also work with PBWin and, with a few modifications, with PBDOS. The range of available options can be seen by scanning the "96T" functions below.
Arguments may be repeated. I.e.:
Add96A(a(), a(), a()) doubles and Multiply96A(a(), a(), a()) squares
so reset output arrays before reuse. Output arrays have their bounds checked by CheckOut(). Input arrays pass through Create96A() and can have the forms:
(a(3),a(2),a(1),a(0)) in 96-bit significand format
a(3)=a(2)=a(1)=0 with a(0) in long integer format
Hi(a(3))=0 with (Lo(a(3)),a(2),a(1)) in extended format (a(0) is ignored)
Most functions return the rounded extended-precision value. Functions using double word arrays have suffix "96A"; those using Type96 UDT's have suffix "96T".
For the rest of the description, the use of the double word array format will be assumed and, other than the equates, none of the above code is needed. There are four ways to input data into a form that can be used, and the function Convert96A() is a general purpose routine to do that with numbers. Convert96A() has the additional capability of converting an optional string using the Val96A() function. Like almost all functions Convert96A() also returns the rounded best extended-precision value corresponding to the array. These routines normally follow PowerBASIC usage of bankers rounding and may decrement the 96-bit significand insignificantly to simplify rounding. Elements are called 0 to 3. The four ways are:
1. Use the optional string to input a positive or negative decimal value as one would with Val() except that spaces and commas are removed automatically. Like Val(), the characters E, e, D, d introduce a power of ten but will degrade accuracy but will be faster than a longer string.
2. Load the array with a valid 96-bit significand value with correct flags. Incorrect flag assignment will result in a warning, and, if the array is located in protected memory, a General Protection Fault will result.
3. Reset the array and load element 0 with a long integer value. The following will work correctly even if MyLong& is negative: MyArray???(0)=MyLong&
4. Zero element 3 (clearing flag bits) and then load the array beginning with element 1 with an extended-precision value. Element 0 will be ignored and zeroed automatically. This can be done by direct addressing and a statement pair:
Dim MyExt(0) As Ext At VarPtr(MyArray(1))
MyArray(3)=0: MyExt(0)=ExtVariable##
Obviously, if ways 2, 3 or 4 are used, the string input must be omitted or null. A value of +zero would be loaded identically by either 3 or 4 with a reset array. For output, like almost all the functions, Convert96A() also returns a rounded extended-precision value. Other functions doing I/O are String96A() and Val96A(). There is also the function String96A() that can be easily expanded to put the output into different string formats. The function Val96A() does the string conversion. PowerOfTwo96A() returns an integer power of two. Prior to PBCC5 with its ability to resolve forward references, these, like almost all these functions, need Declare statements.
The next functions deal with addition and multiplication on which many of the other functions rely. They are appropriately named Add96A() and Multiply96A(). Each has two input arrays as the last two parameters and an output array, any or all of which can be the same array. As multiplying two 96-bit significands requires up to 192 bits for the
product, an intermediate array of 24 bytes is used for each function. The addition function aligns the terms by shifting bits and adding (with the appropriate sign) the smaller to the larger using a quad array to hold the sum. The multiplication methodology is standard elementary school long multiplication of double words. It should be noted that as there is no quad-word datatype in PowerBASIC, multiplying two double words is best done in assembly using the opcode for unsigned multiply, the EDA:EAX register pair, and then adding the high and low double words to the appropriate quad array element.
The next set of functions do some simple chores. Sign96A() returns a long integer equivalent to the sign of the difference (or, if just one array passed, the sign of that value). Copy96A(), Negative96A() and Absolute96A() respectively copies, negates and takes the absolute value of an 96-bit significand format array. Negative96A() can be used before Add96A() for the purpose of subtraction. The remaining three are a bit more complicated but show how these functions fit together. They are Power96A(), Reciprocal96A(), and SquareRoot96A(). The Power96A() function raises a 96-significand value to an integer power. The Reciprocal96A() function can be used before Multiply96A() to permit division. The functions Subtract96A() and Divide96A() were added to simplify those operations. Reciprocal96A() and SquareRoot96A() make use of an approximate value calculated in the normal PB fashion to create extended-precision approximations to the desired result and an additive correction term to refine that value. This is possible because the expression used to invert the expression ((1/x)*x=1 and Sqr(x)*Sqr(x)=x) can be computed to near 96-bit accuracy. One important note, a modification is made to SquareRoot96() so it actually calculates the equivalent to the PB expression:
Sgn(x)*Sqr(Abs(x)).
In this last batch of code are three subroutines that are used internally. The sub BankersRounding96A() tweaks element 0 to simplify bankers' rounding and addresses a peculiarity in rounding denormals. The sub CheckOut96A() checks the output array dimensions. The sub Finish96A() tidies up the longer 6-element arrays used in addition and multiplication. The sub RoundUp96A() rounds up a 96-bit significand and, by setting element 0 of a dummy copy to &HFFFFFFFF, can be used to round up a 64-bit significand.
There is quite a bit of redundancy in these functions. This is intentional as they are still being developed. It also helps to track coding errors. Also since this suite is intended to be supplemental and used in R&D work, it should not form part of any final program. One use for these functions is when it is desirable to create a table of mathematical constants embedded in the code as !dw or !dd data. The coding would be something like:
If one wants a set of data for testing, I would suggest checking out the Factorial96A() function that appears in a later post or there is also a program that will give the best extended-precision values for all integer powers of ten on my post in the PBCC forum. Note however the correction to the first line of data later in that thread. To find it do a search for the function's name: BestPowerOfTen.
________________
Mark Longley-Cook
Code:
' As used here, my "96-bit significand" format is as follows: ' Bit(s) Meaning ' 127 1 (Signature bit indicating valid format) ' 126-112 Flag bits (see below) and reserved _______ ' 111 1 if negative (or NaN); otherwise 0 | Equal to ' 110-96 Biased two's exponent (all 1's if Inf/NaN)| extended- ' 95 Always one (MSB of significand)* | precision ' -- Implicit binary point | (truncated) ' 94-32 Next significant 63 bits of significand __| 80-bit format ' 31-0 Least significant 32 bits of significand ' * Bit 95 is zero if denormal or zero as in extended-precision. ' Exponent bias is &H3FFF for normals (&H3FFE for denormals). ' In addition to signature, reserved, and flag bits (112-127): ' Zeros always have bits 0-110 reset. ' Denormals always have bits 95-110 reset. ' Infinities always have bits 0-94 reset and bits 95-110 set. ' NaN's always have bits 0-93 reset and bits 94-111 set. ' Flag bit positions between 112 and 127 are arbitrary. ' They are in separate nibbles for readability. ' If changed, modify in-code constants, e.g., Infinity in Factorial96A() %FlagSig96=127 ' Set if 96-bit significand format %FlagNaN96=124 ' Set if NaN FFFF C0000000 00000000 %FlagInf96=120 ' Set if Infinity {7FFF|FFFF} 80000000 00000000 %FlagDen96=116 ' Set if Denormal {0000|8000} binary 0xxxxxx... %FlagZer96=112 ' Set if Zero {0000|8000} 00000000 00000000 ' Following flag bits determined by extended-precision format: %FlagNeg96=111 ' =79+32 Set if Negative (or NaN) %FlagNor96= 95 ' =63+32 Set if finite normal (or infinity or NaN) %FlagInc96= 31 ' =-1+32 Set if round up for 64-bit significand ' UDT below required if using the 96T style. ' Recommended for arrays of 96-bit significand values. Type Type96 Dword L As Long ' Can use to load a long integer value X As Ext ' Can use to load an extended-precision value F As Word ' Flags and reserved bits End Type
I had been posting some notes I have made on various PowerBASIC functions on the PB Console Compiler forum more or less as an academic exercise. To have baseline data with which to work, I needed arrays of a data that represented the most accurate extended-precision representations of numbers, rather like an enhanced Val() function. To do this initially I wrote a routine that doubled or halved a string of decimal digits to convert them to binary and then put that into the extended-precision format. Without any restriction on string length, this made the process very slow. The original strings would come from web, the Microsoft calculator included in Windows, or from an
expression such as "1E"+Format$(n). While preparing a post on the Exp() function, I recognized that a routine to multiply floating point numbers would be needed, and that is how these functions were conceived.
One of the posts mentioned above provides information on the 80-bit extended precision format. These functions use an enhanced version of that format called by Intel double extended-precision. This enhanced format adds (on the low-memory end) 32 bits to the 64-bit significand and adds (on the high-memory end) 16 bits to the 16-bit sign/exponent. The 96-bit significand (hence "96" in the naming) is equivalent to
almost 29 significant decimal digits. The 32 bits add ample extra precision to avoid accumulation of rounding errors and the other 16 bits contain a few flags to simplify the detection of zero, infinity and denormals. Denormals are those numbers very near zero that have reduced precision and use, depending on how one interprets it, either
an exponent bias of &H3FFE instead of the standard &H3FFF or else an implied binary point one bit further into the significand. Denormals have very limited support in PowerBASIC (they are, for example, printed as zero) but can be created and manipulated in the Floating Point Unit (FPU) that does extended-precision arithmetic.
This gives a total size of 128 bits (16 bytes) to this enhanced format with a 96-bit significand. All the original development of the functions was done dividing these 16 bytes into arrays of four double words dimensioned 0 to 3. The functions that use these arrays have "96A" as a suffix to their names. One can equally well use a user-defined type (UDT) and a Type structure and a set of functions with suffix "96T" are given below. Also provided are the flag locations in equate statements. There flag locations in bits 112 to 127 can be changed to others at will; those below 112 are tied to the
extended-precision format. The Type structure and equate definitions must precede their actual usage. The code was written for PowerBASIC Console Compiler (PBCC) Version 5.01. The functions should work with earlier PBCC versions with Declare statements or rearranging. The functions should also work with PBWin and, with a few modifications, with PBDOS. The range of available options can be seen by scanning the "96T" functions below.
Arguments may be repeated. I.e.:
Add96A(a(), a(), a()) doubles and Multiply96A(a(), a(), a()) squares
so reset output arrays before reuse. Output arrays have their bounds checked by CheckOut(). Input arrays pass through Create96A() and can have the forms:
(a(3),a(2),a(1),a(0)) in 96-bit significand format
a(3)=a(2)=a(1)=0 with a(0) in long integer format
Hi(a(3))=0 with (Lo(a(3)),a(2),a(1)) in extended format (a(0) is ignored)
Most functions return the rounded extended-precision value. Functions using double word arrays have suffix "96A"; those using Type96 UDT's have suffix "96T".
Code:
' General conversion and tweaking significand: Function Convert96T(a As Type96, Opt s As String) As Ext ' Convert in-situ extended-precision or long-integer or string. ' This does not have to be called as is done internally. Dim a96(3) As Dword At VarPtr(a) Function=Convert96A(a96(), s) End Function Function PowerOfTwo96T(y As Type96, n As Long) As Ext ' y=2^n Integer power of two Dim y96(3) As Dword At VarPtr(y) Function=PowerOfTwo96A(y96(), n) End Function Function Str96T(a As Type96, n As Long) As String ' Output=Str$(a,n) (n<=29). Dim a96(3) As Dword At VarPtr(a) Function=Str96A(a96(), n) End Function Function String96T(a As Type96, Opt o As Long) As String ' Convert a to string depending on option code. Dim a96(3) As Dword At VarPtr(a) Function=String96A(a96(), o) End Function Function Tweak96T(a As Type96, n As Quad) As Ext ' Significand(a)=Significand(a)+n in-situ Dim a96(3) As Dword At VarPtr(a) Function=Tweak96A(a96(), n) End Function Function Val96T(y As Type96, s As String) As Ext ' y=Val(s) Convert from string Dim y96(3) As Dword At VarPtr(y) Function=Val96A(y96(), s) End Function ' Arithmetic functions: ' Absolute, Add, Copy, Multiply, Negative, Power, Reciprocal, Sign, and SquareRoot. ' Results appear as first argument and extended-precision value as function return. ' The exception is Sign with a long-integer result as the function return value. Function Absolute96T(y As Type96, a As Type96) As Ext ' y=Abs(a) Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=Absolute96A(y96(), a96()) End Function Function Add96T(y As Type96, a As Type96, b As Type96) As Ext ' y=a+b Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Dim b96(3) As Dword At VarPtr(b) Function=Add96A(y96(), a96(), b96()) End Function Function Copy96T(y As Type96, a As Type96) As Ext ' y=a Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=Copy96A(y96(), a96()) End Function Function Divide96T(y As Type96, a As Type96, b As Type96) As Ext ' y=a/b Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Dim b96(3) As Dword At VarPtr(b) Function=Divide96A(y96(), a96(), b96()) End Function Function Multiply96T(y As Type96, a As Type96, b As Type96) As Ext ' y=a*b Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Dim b96(3) As Dword At VarPtr(b) Function=Multiply96A(y96(), a96(), b96()) End Function Function Negative96T(y As Type96, a As Type96) As Ext ' y=-a Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=Negative96A(y96(), a96()) End Function Function Power96T(y As Type96, a As Type96, n As Long) As Ext ' y=a^n (n must be an integer) Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=Power96A(y96(), a96(), n) End Function Function Reciprocal96T(y As Type96, a As Type96) As Ext ' y=1/a Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=Reciprocal96A(y96(), a96()) End Function Function Sign96T(a As Type96, Opt b As Type96) As Long ' y=Sgn(a-b) (y=Sgn(a) if b is not passed) Dim a96(3) As Dword At VarPtr(a) If VarPtr(b) Then ' b was passed y=Sgn(a-b) Dim b96(3) As Dword At VarPtr(b) Function=sign96A(a96(), b96()) Else ' b was not passed y=Sgn(a) Function=Sign96A(a96()) End If End Function Function SquareRoot96T(y As Type96, a As Type96) As Ext ' y=Sqr(Abs(a))*Sgn(a) (Note treatment for x<0) Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Function=SquareRoot96A(y96(), a96()) End Function Function Subtract96T(y As Type96, a As Type96, b As Type96) As Ext ' y=a-b Dim y96(3) As Dword At VarPtr(y) Dim a96(3) As Dword At VarPtr(a) Dim b96(3) As Dword At VarPtr(b) Function=Subtract96A(y96(), a96(), b96()) End Function
1. Use the optional string to input a positive or negative decimal value as one would with Val() except that spaces and commas are removed automatically. Like Val(), the characters E, e, D, d introduce a power of ten but will degrade accuracy but will be faster than a longer string.
2. Load the array with a valid 96-bit significand value with correct flags. Incorrect flag assignment will result in a warning, and, if the array is located in protected memory, a General Protection Fault will result.
3. Reset the array and load element 0 with a long integer value. The following will work correctly even if MyLong& is negative: MyArray???(0)=MyLong&
4. Zero element 3 (clearing flag bits) and then load the array beginning with element 1 with an extended-precision value. Element 0 will be ignored and zeroed automatically. This can be done by direct addressing and a statement pair:
Dim MyExt(0) As Ext At VarPtr(MyArray(1))
MyArray(3)=0: MyExt(0)=ExtVariable##
Obviously, if ways 2, 3 or 4 are used, the string input must be omitted or null. A value of +zero would be loaded identically by either 3 or 4 with a reset array. For output, like almost all the functions, Convert96A() also returns a rounded extended-precision value. Other functions doing I/O are String96A() and Val96A(). There is also the function String96A() that can be easily expanded to put the output into different string formats. The function Val96A() does the string conversion. PowerOfTwo96A() returns an integer power of two. Prior to PBCC5 with its ability to resolve forward references, these, like almost all these functions, need Declare statements.
Code:
Function Convert96A(a() As Dword, _ ' Input array Opt s As String) As Ext ' Input string does Val96A() ' Called with all input arrays at start of functions. ' Called with output array at end of functions. ' Convert in-situ extended-precision in (a(3),a(2),a(1)) or ' long integer in a(0) to 96-bit significand format. ' Also checks bounds of input arrays. ' Second parameter (IncrSig) intended for internal use. ' ' Methodology: ' If s is provided, immediately calls Val96A() to convert it. ' Otherwise, a() is checked for following possibilities in this order: ' 1. If already has 96-bit significand signature bit set: ' a. Set flag bits as required ' b. If %FlagInc96 is set round up copy of a() ' c. Return rounded extended-precision value and exit ' 2. If a(3)=a(2)=a(1)=0 (includes extended of value +zero): ' a. Treat a(0) as long integer value ' b. Set signature bit ' c. Set %FlagZer96 if zero ' d. Zero a(0) (no rounding required or possible) ' e. Return extended-precision value and exit ' 3. Otherwise ' a. Treat (Lo(a(3)),a(2),a(1)) as extended-precision value ' b. Set signature bit ' c. Set any flag bits required ' d. Zero a(0) (no rounding required or possible) ' e. Return extended-precision value and exit ' Call CheckOut96A(a()) Local i As Long Local sg As Long ' Sign Dim wa(7) As Word At VarPtr(a(0)) ' w(7)=flags; w(6)=sign & exponent Dim la(0) As Long At VarPtr(a(0)) Dim aa(0) As Ext At VarPtr(a(1)) Dim b(3) As Dword ' Working copy used if signature flag set Dim wb(7) As Word At VarPtr(b(0)) ' wb(7)=flags; wb(6)=sign/exp. Dim bb(0) As Ext At VarPtr(b(1)) Dim q64(0) As Quad At VarPtr(b(1)) ' Extended significand Dim w64(0) As Word At VarPtr(b(3)) ' Extended sign/exponent ' ' Check if string input and use Val96A() If s<>"" Then Call Val96A(a(), s) If Bit(a(0),%FlagSig96)=0 Then ?"Error Val return": WaitKey$ End If ' ' 1. Check if already in 96-bit significand format. ' Since a() may be in protected memory, work on a copy. If Bit(a(0),%FlagSig96) Then ' General protection error will occur if ' a() in protected momory and flags set incorrectly. ' Make a working copy in unprotected memory. For i=0 To 3: b(i)=a(i): Next i ' Copy96A() calls Convert96A() If wb(6)=&HFFFF?? And wb(5)=>&HC000?? Then ' NaN Reset b() Bit Set b(0),%FlagSig96 Bit Set b(0),%FlagNaN96 For i=%FlagNor96-1 To %FlagNeg96: Bit Set b(0),i: Next i Function=(0##/0##) ' All NaN's have sign bit set GoSub CheckFlags Exit Function End If If Bit(b(0),%FlagNeg96) Then sg=-1 Else sg=+1 b(3)=b(3) And &H00007FFF??? ' Clear flags and sign If b(2)=0 And b(1)=0 And b(0)=0 Then If b(3)<>0 Then ?"Error zero significand":WaitKey$ Reset b() Bit Set b(0),%FlagSig96 Bit Set b(0),%FlagZer96 If sg=-1 Then Bit Set b(0),%FlagNeg96 Function=sg*0## GoSub CheckFlags Exit Function End If If b(3)=&H00007FFF??? Then Reset b() Bit Set b(0),%FlagSig96 Bit Set b(0),%FlagInf96 If sg=-1 Then Bit Set b(0),%FlagNeg96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set b(0),i: Next i Function=sg/0## GoSub CheckFlags Exit Function End If Bit Set b(0),%FlagSig96 If sg=-1 Then Bit Set b(0),%FlagNeg96 If Bit(b(0),%FlagNor96)=0 Then Bit Set b(0),%FlagDen96 Call BankersRounding96A(b()) GoSub CheckFlags ' Use the working copy as scratch for rounding return value. If Bit(b(0),%FlagInc96) Then b(0)=&HFFFFFFFF??? ' Forces rounding of 64-bit significand Call RoundUp96A(b()) End If Function=bb(0) ' Altered b() is not a part of a() Exit Function CheckFlags: ' Check if flags correct: i=0 If a(0)<>b(0) Or a(1)<>b(1) Or a(2)<>b(2) Then Incr i ? "Warning significands don't match":WaitKey$ End If If Lo(Word,a(3))<>Lo(Word,a(3)) Then Incr i ? "Warning sign/exponents don't match":WaitKey$ End If If Hi(Word,a(3))<>Hi(Word,a(3)) Then Incr i ? "Warning flags don't match":WaitKey$ End If If i Then For i=0 To 3: a(i)=b(i): Next i End If Return End If ' ' 2. Check if long integer in a(0) (or +zero extended-precision). If a(3)=0 And a(2)=0 And a(1)=0 Then Bit Set a(0),%FlagSig96 Function=la(0) Select Case a(0) Case &H00000000& ' +Zero Bit Set a(0),%FlagZer96 ' Zero Exit Function Case &H80000000& ' -2^31 (cannot simply negate) wa(6)=&HC01E??: Swap a(2),a(0) ' So significand OK Function=-2^31 Exit Function End Select wa(6)=&H401E?? ' &H3FFF+31=&H401E If Bit(a(0),%FlagInc96) Then ' Bit 31 set if negative long Bit Set a(0),%FlagNeg96: la(0)=-la(0) End If Do ' As Bit(a(0),%FlagInc96)=0, do loop done at least once. Shift Left a(0),1 Decr wa(6) Loop Until Bit(a(0),%FlagInc96) ' %FlagInc96-31 Swap a(2),a(0) ' Significand has MSB set and a(0)=0 Exit Function End If ' ' 3. Must be extended-precision (other than +zero). wa(7)=0: a(0)=0 ' Clear flags and a(0) Bit Set a(0),%FlagSig96 If wa(6)=&H8000?? And a(2)=0 And a(1)=0 Then Bit Set a(0),%FlagZer96 ' -Zero ElseIf a(2)=&H800000000??? And a(1)=0 Then Select Case wa(6) Case &H7FFF??,&HFFFF??: Bit Set a(0),%FlagInf96 ' Infinity End Select ElseIf Bit(a(0),%FlagNor96)=0 Then Bit Set a(0),%FlagDen96 ' Denormal ElseIf wa(6)=&HFFFF?? And wa(5)=>&HC000?? Then wa(5)=&HC000??: Bit Set a(0),%FlagNaN96 ' NaN For i=0 To 4: wa(i)=0: Next i End If Function=aa(0) Exit Function End Function Function PowerOfTwo96A(y() As Dword, _ ' Output array n As Long) As Ext ' Input integer power of two ' y <-- 2^n Call CheckOut96A(y()) Local i As Long Reset y() Bit Set y(0),%FlagSig96 Select Case n Case <-16477 ' Zero Bit Set y(0),%FlagZer96 Case -16477 To -16383 ' Denormal ' -16477 To -16446 are zero in extended-precision Bit Set y(0),%FlagDen96 Bit Set y(0),n+16477 ' 16477=16382+95 Case -16382 To 16383 ' Normal y(3)=y(3)+&H00003FFF???+n Bit Set y(0),%FlagNor96 Case >16383 ' Infinity Bit Set y(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i End Select Function=Convert96A(y()) End Function Function Str96A(a() As Dword, _ ' Input array ByVal n As Long) As String ' Input length (29 max) ' Output <-- Str$(a,n) Convert to n-digit string (n<=29). ' The 29th digit is subject to error as 2^95=0.396E29<1E29. ' If n<=0 then defaults to 29. ' Methodology: Use Str$() and Val96A() to get correction term. Call Convert96A(a()) Local a0,i,j,k As Long: a0=Asc("0") Dim d(3) As Dword ' |a()|, possibly biased by 1E64 Dim t(3) As Dword ' Temporary array Dim s(2) As String ' String components of result Dim e(2) As Long ' ...E... exponent If Bit(a(0),%FlagNaN96) Then Function="Not-a-Number":Exit Function If Bit(a(0),%FlagInf96) Then i=Bit(a(0),%FlagNeg96) If i Then Function="-Infinity" Else Function="+Infinity" Exit Function End If If Bit(a(0),%FlagZer96) Then If Bit(a(0),%FlagNeg96) Then Function="-0" Else Function=" 0" Exit Function End If If n<=0 Or n>29 Then n=29 If Bit(a(0),%FlagDen96) Then t(0)=10 Call Power96A(t(), t(),64) ' 10^64 Call Multiply96A(t(), a(),t()) ' Bias by 10^64 s(0)=Trim$(Str$(Absolute96A(d(), t()),18)) Else s(0)=Trim$(Str$(Absolute96A(d(), a()),18)) End If Call Val96A(t(), s(0)) Call Subtract96A(d(), d(),t()) ' d() = correction to |a()| s(1)=Trim$(Str$(Absolute96A(t(), d()),18)) ' d() retains sign of correction term For i=0 To 1 If InStr(s(i),".") Then If InStr(s(i),"E")=0 Then s(i)=s(i)+"E0" ElseIf InStr(s(i),"E") Then Replace "E" With ".E" In s(i) Else s(i)=s(i)+".E0" End If j=InStr(s(i),"E") e(i)=Val(Mid$(s(i),j+1)) If Bit(a(0),%FlagDen96) Then e(i)=e(i)-64 s(i)=Left$(s(i),j-1) e(i)=e(i)+InStr(s(i),".")-2 ' Convert to x.xxxxx... Replace "." With "" In s(i) Next i If e(1)=>e(0) Then ? "Error exponent":WaitKey$ s(1)=String$(e(0)-e(1),"0")+s(1) e(1)=e(0): e(2)=e(0) k=Max&(Len(s(0)),Len(s(1))) s(0)=s(0)+String$(k-Len(s(0)),"0") s(1)=s(1)+String$(k-Len(s(1)),"0") If Bit(d(0),%FlagNeg96) Then j=0 For i=k To 1 Step -1 j=Asc(s(0),i)-Asc(s(1),i)-j If j<0 Then s(2)=Chr$(j+10+a0)+s(2):j=1 Else s(2)=Chr$(j+a0)+s(2):j=0 End If Next i If Trim$(s(2),"0")="" Then Function=" 0": Exit Function Do While Asc(s(2))=a0 s(2)=Mid$(s(2),2) Decr e(2) Loop Else j=0 For i=k To 1 Step -1 j=Asc(s(0),i)-a0+Asc(s(1),i)-a0+j s(2)=Chr$((j Mod 10)+a0)+s(2) j=j\10 Next i If j Then s(2)="1"+s(2): Incr e(2) End If s(2)=RTrim$(s(2),"0") Select Case Len(s(2)) Case <=n: s(2)=s(2)+String$(n-Len(s(2)),a0):j=0 Case Else Select Case Asc(s(2),n+1)-a0 Case 0 To 4: j=0 Case 5 Select Case Asc(s(2),n)-a0 Case 0,2,4,6,8: j=Min&(1,Len(s(2))-(n+1)) Case 1,3,5,7,9: j=1 End Select Case 6 To 9: j=1 End Select s(2)=Left$(s(2),n) End Select Do While j For i=n To 1 Step -1 j=Asc(s(2),i)-a0+j Asc(s(2),i)=(j Mod 10)+a0 j=j\10 Next i If j Then s(2)="1"+s(2): Incr e(2) Select Case Asc(s(2),-1)-a0 Case 0 To 5:j=0 Case 6 To 9:j=1 End Select s(2)=Left$(s(2),n) Else j=0 End If Loop If Len(s(2))>n+1 Then j=Asc(s(2),n+1)=>a0+5:s(2)=Left$(s(2),n) If Bit(a(0),%FlagNeg96) Then s(2)="-"+s(2) Else s(2)=" "+s(2) Function=Left$(s(2),2)+"."+Mid$(s(2),3)+"E"+Format$(e(2)) End Function Function String96A(a() As Dword, _ ' Input array Opt ByVal o As Long) _ ' Input optional action to be taken As String ' String output depending on option o ' o<0 adds a $CrLf prefix for Clipboard statements ' This function is subject to expansion and modification so output: ' Generally used with Print for testing results, and ' For inbedded constants (e.g. see the factorial post). Local x As Ext x=Convert86A(a()) Local i As Long Local d As Dword Local s As String Local t As String Dim b(3) As Dword Dim c(3) As Dword Dim ic(0) As Long At VarPtr(c(0)) Dim wx(5) As Word At VarPtr(x) If Bit(a(0),%FlagNaN96) Then ? "Warning NaN":WaitKey$ ' Clipboard needs $CrLf between lines. If o<0 Then o=Abs(o):s=$CrLf Select Case o Case 0 ' Each DWord in hex format separated by space ' Full 29-digit decimal value appended For i=3 To 0 Step -1: s=s+Hex$(a(i),8)+$Spc: Next i s=s+"="+$Spc+Str96A(a(),29) Case 1 ' Suitable string to put 80 bits on Clipboard (PBCC5+) s=s+" !dw " ' Clipboard needs $CrLf between lines. For i=0 To 4 s=s+"&H"+Hex$(wx(i),4)+Choose$(i+1,",",",",",",", ","") Next i Case 2 ' Suitable string to put 128-bits on Clipboard (PBCC5+) s=s+"!dd " For i=0 To 3 s=s+"&H"+Hex$(a(i),8)+Choose$(i+1,", ",",",", ","") Next i Case 3 ' Str96A(...,29) with space separators t=Str96A(a(),29) ' String to avoid conflict with o=-3 i=InStr(t,".") s=s+Left$(t,i+5)+ _ $Spc+Mid$(t,i+06,5)+$Spc+Mid$(t,i+11,5)+ _ $Spc+Mid$(t,i+16,5)+$Spc+Mid$(t,i+21,5)+ _ $Spc+Mid$(t,i+26,3)+$Spc+Mid$(t,i+29) Case 8 ' Octal output for comparing with Donald E. Knuth's ' "The Art of Computer Programming" (any vol.) appendix ' This was written for values near 1.0 -- modify if testing others. If Bit(a(0),%FlagNeg96) Then s="-" Else s="" Call Absolute96A(b(),a()) x=Abs(x) If x=>Exp2(31) Then ? "Error Too large" ic(0)=Fix(x) s=s+Oct$(ic(0))+"." Decr ic(0) Call Subtract96A(b(), b(),c()) ' 1<b<2 If b(3)<>&H80003FFF??? Then ? "Error in octal":WaitKey$ For i=92 To 2 Step -3 s=s+Oct$(4*Bit(b(0),i+2)+2*Bit(b(0),i+1)+Bit(b(0),i)) If (i Mod 15)=5 Then s=s+$Spc Next i s=s+Oct$(4*Bit(b(0),1)+2*Bit(b(0),0)) Case Else: ?"Error Unsupported option ";o:WaitKey$ End Select Function=s End Function Function Tweak96A(a() As Dword, _ ' Input array ByVal n As Quad) As Ext _ ' Input tweak ' Significand(a) <-- Significand(a)+n in-situ ' Will not tweak NaN, Infinity, or Zero. Call Convert96A(a()) Local i As Long Dim q(0) As Quad At VarPtr(a(0)) If n=0 Then Function=Convert96A(a()): Exit Function If Bit(a(0),%FlagNaN96) Then ? "Warning Tweak NaN";:n=0 If Bit(a(0),%FlagInf96) Then ? "Warning Tweak Infinity";:n=0 If Bit(a(0),%FlagZer96) Then ? "Warning Tweak Zero";:n=0 If Bit(a(0),%FlagDen96) And a(2)=0 And q(0)+n=0 Then ? "Warning Tweak denormal to zero";:n=0 End If If q(0)<0 Xor q(0)+n<0 Then ? "Warning Tweak negative quad";:n=0 If n=0 Then ? " cannot be done":WaitKey$ q(0)=q(0)+n Function=Convert96A(a()) End Function Function Val96A(y() As Dword, _ ' Output array ByVal s As String) As Ext ' Input string ' Converts decimal string to 96-bit significand array ' Change numeric string to a() ' Allowed characters: ' Commas and spaces (removed at start) ' Initial "-" (also after E) ' One decimal point and one E (or e,D,d) ' Digits (0123456789) ' Most accurate if E avoided, but if s is very long will be slow. ' ' Methodology for decimal significand: ' 1. If 1.0 <= |value| < 2.0 goto 4 ' 2. Halve or double value changing exponent ' 3. Goto 1 ' 4. If |value| => 1.0 then ' 4a. value <-- value - 1 ' 4b. Set significand bit ' 5. Double |value| and move bit pointer towards LSB ' 6. While bit pointer within significand goto 4 ' 7. Round final value ' Call CheckOut96A(y()) Local i,j As Long Local e As Long ' Decimal exponent Local n As Long ' -1 if negative, else 0 Local dp As Long ' Decimal point is just before this character Local bn As Long ' Bit number to set Local t As String ' Temporary string Dim b(0) As Integer At VarPtr(y(3)) ' Binary exponent Dim p(3) As Dword ' Power-of-ten factor Reset y() Bit Set y(0),%FlagSig96 s=Remove$(s, Any " ,+") i=InStr(s, Any "EeDd") If i Then e=Val(Mid$(s,i+1)) ' If i=0, e=0 s=Left$(s,i-1) End If n=(Asc(s)=Asc("-")): If n Then s=Mid$(s,2) If InStr(s,".")=0 Then s=s+"." If Tally(s,".")<>1 Then ? "Error multiple decimals":WaitKey$ s=Trim$(s,"0") If Remove$(s, Any "01234456789.")<>"" Then ? "Error character(s) ";Remove$(s,Any"01234456789."):WaitKey$ End If If s="." Then ' Zero Bit Set y(0),%FlagZer96 If n Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) Exit Function ' Power makes no difference End If ' Initial values: b(0)=&H3FFF%: bn=3*32-1: dp=InStr(s,".") s=Remove$(s,".") If dp<>1 Then ' Halve until <2.0 Do Until dp=2 And Asc(s)=Asc("1") j=0 For i=1 To Len(s) j=Asc(s,i)-Asc("0")+10*j Mid$(s,i,1)=Chr$(j\2+Asc("0")) j=Bit(j,0) Next i Incr b(0) If b(0)=&H7FFF% Then ? "Warning: Decimal significand too long":WaitKey$ Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1:Bit Set y(0),i:Next i GoTo IncludePowerSign End If If j Then s=s+"5" j=Len(s): s=LTrim$(s,"0"): dp=dp-(j-Len(s)) s=RTrim$(s,"0") If Len(s)<dp Then s=s+String$(dp-Len(s),"0") Loop Else ' Double until =>1.0 Do Until dp=2 And Asc(s)=Asc("1") Decr b(0) If b(0)=&H0000% Then Bit Set y(0),%FlagDen96 Exit Do End If GoSub DoubleString Loop End If ' Create significand Do If dp=2 And Asc(s)=Asc("1") Then If bn=>0 Then Bit Set y(0),bn Mid$(s,1,1)="0" Else Call RoundUp96A(y()) Exit Do End If End If If bn<0 Then Exit Do If Trim$(s,"0")="" Then Call BankersRounding96A(y()) Exit Do End If Decr bn GoSub DoubleString Loop IncludePowerSign: #If -1 ' Use True (-1) if PowerOfTen96A() exists, else use 0. ' More accurate, but PowerOfTen96A() is almost 600KB. Call PowerOfTen96A(p(),e) Call Multiply96A(y(), y(),p()) #Else ' Less accurate but shorter. ' Using 5 instead of 10 avoids denormals in Power96A(). p(0)=5: Call Power96A(p(), p(),e) Call Multiply96A(y(), y(),p()) b(0)=b(0)+e #EndIf Select Case b(0) Case =>&H7FFF ' Infinity Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1:Bit Set y(0),i:Next i Case <=0 ' Denormal Do While b(0) Shift Right y(0),1 If Bit(y(1),0) Then Bit Set y(0),31 Shift Right y(1),1 If Bit(y(2),0) Then Bit Set y(1),31 Shift Right y(2),1 Incr b(0) Loop If y(0)=0 And y(1)=0 And y(2)=0 Then Bit Set y(0),%FlagZer96 Else Bit Set y(0),%FlagDen96 End If End Select If n Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) Exit Function DoubleString: ' GoSub routine called when decrementing b(0) or bn j=0 For i=Len(s) To 1 Step -1 j=2*(Asc(s,i)-Asc("0"))+j Mid$(s,i,1)=Chr$((j Mod 10)+Asc("0")) j=j\10 Next i If j Then s="1"+s Incr dp End If s=RTrim$(s,"0") If Len(s)<dp Then s=s+String$(dp-Len(s),"0") Return End Function
product, an intermediate array of 24 bytes is used for each function. The addition function aligns the terms by shifting bits and adding (with the appropriate sign) the smaller to the larger using a quad array to hold the sum. The multiplication methodology is standard elementary school long multiplication of double words. It should be noted that as there is no quad-word datatype in PowerBASIC, multiplying two double words is best done in assembly using the opcode for unsigned multiply, the EDA:EAX register pair, and then adding the high and low double words to the appropriate quad array element.
Code:
Function Add96A(y() As Dword, _ ' Output array a() As Dword, b() As Dword) As Ext ' Input arrays ' y() <-- a() + b() ' Precede with Negative96A() to subtract. ' a(), b(), and y() are DWord arrays dimensioned (0:3). ' a() and/or b() can be in extended-precision or long format as, e.g.: ' Dim aa(0) As Ext At VarPtr(a(1)):aa(0)=xxx## or a(0)=xxx& ' For compatability with Finish96A() (called by Multiply96A()), ' a working DWord array dimensioned (0:5) is used. ' Methodology: shift to align and add. ' Call CheckOut96A(y()) Call Convert96A(a()): Call Convert96A(b()) Local i,j As Long Local xbd As Long ' True (-1) if excess bits discarded Local x As Ext Dim aa(0) As Ext At VarPtr(a(1)) ' Extended approximation Dim bb(0) As Ext At VarPtr(b(1)) ' Extended approximation Dim mn(3) As Dword: Dim mx(3) As Dword ' Min and max abs values Dim eab(1) As Long ' Exponents of a() and b() Dim sab(1) As Long ' Signs of a() and b() Dim ab(5,1) As Dword ' Matrix of a() & b() extended significands Dim ab3(1) As Dword ' a(3) and b(3) in array format Dim q(5) As Quad ' q(5) = temporary significand array (quad) ' Check if NaN If Bit(a(0),%FlagNaN96) Or _ Bit(b(0),%FlagNaN96) Or _ (Bit(a(0),%FlagInf96) And _ Bit(b(0),%FlagInf96) And _ (Bit(a(0),%FlagNeg96) Xor _ Bit(b(0),%FlagNeg96) ) ) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagNaN96 For i=%FlagNor96-1 To %FlagNeg96: Bit Set y(0),i: Next i ' NaN Function=Convert96A(y()) Exit Function End If ' Check if |a()|=|b()| If IsFalse((a(3) Xor b(3)) And &H00007FFF???) And _ a(2)=b(2) And a(1)=b(1) And a(0)=b(0) Then If Bit(a(0),%FlagZer96) And Bit(b(0),%FlagZer96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 If Bit(a(0),%FlagNeg96) And Bit(b(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 End If ElseIf Bit(a(0),%FlagNeg96) Xor Bit(b(0),%FlagNeg96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 Else ' Must be double Call Copy96A(y(), a()) Select Case y(3) And &H00007FFF??? Case &H00007FFF??? Case &H00007FFE??? Reset y() y(3)=a(3)+1 Bit Set y(0),%FlagInf96 Bit Set y(0),%FlagNor96 Case &H00000000??? If Bit(y(0),%FlagNor96) Then ? "Error normal":WaitKey$ Shift Left y(2),1 If Bit(y(1),31) Then Bit Set y(2),0 Shift Left y(1),1 If Bit(y(0),31) Then Bit Set y(1),0 Shift Left y(0),1 If Bit(y(0),%FlagNor96) Then Bit Reset y(0),%FlagDen96 Incr y(3) End If Case Else Incr y(3) End Select End If Function=Convert96A(y()) Exit Function End If ' Check if y()=a() If Bit(a(0),%FlagInf96) Or Bit(b(0),%FlagZer96) Then Call Copy96A(y(), a()) Function=Convert96A(y()) Exit Function End If ' Check if y()=b() If Bit(a(0),%FlagZer96) Or Bit(b(0),%FlagInf96) Then Call Copy96A(y(), b()) Function=Convert96A(y()) Exit Function End If ' Done with special cases Call Absolute96A(mn(),a()): Call Absolute96A(mx(),b()) Select Case Sign96A(mn(),mx()) Case -1: j=0 ' |a{}|<|b{}| Case +1: j=1 ' |a{}|>|b{}| Case Else: ? "Error Unexpected equal arrays":WaitKey$ End Select For i=3 To 5: ab(i,j)=a(i-3): ab(i,1-j)=b(i-3): Next i ab3(j)=a(3): ab3(1-j)=b(3) ' eab(j)=exponent treating (ab(5,j),ab(4,j),ab(3,j)} as integer For j=0 To 1 eab(j)=(ab3(j) And &H00007FFF???)-&H00003FFF???-95 If eab(j)=-16478 Then Incr eab(j) ' Denormal (&H3FFF+95=16478) eab(j)=eab(j)+32 If (ab3(j) And &H00008000???) Then sab(j)=-1 Else sab(j)=+1 Next j ' ab3() no longer needed. ' Normalize ab(,j) For j=0 To 1 Do Until Bit(ab(5,j),31) For i=5 To 3 Step -1 Shift Left ab(i,j),1 If Bit(ab(i-1,j),31) Then Bit Set ab(i,j),0 Next i Decr eab(j) Loop Next j ' Make eab(0)=>eab(1) If eab(0)<eab(1) Then For i=3 To 5: Swap ab(i,0),ab(i,1): Next i Swap eab(0),eab(1) Swap sab(0),sab(1) End If If eab(0)-eab(1)<6*32 Then Do Until eab(0)=eab(1) xbd=xbd Or IsTrue(Bit(ab(0,1),0)) For i=0 To 5 Shift Right ab(i,1),1 If i<5 Then If Bit(ab(i+1,1),0) Then Bit Set ab(i,1),31 Next i Incr eab(1) Loop Else ' Smaller term is too small to matter other than in rounding xbd=-1 For i=0 To 5 ab(i,1)=&H00000000??? Next i eab(1)=eab(0) sab(1)=sab(0) End If ' Add For i=0 To 5 q(i)=ab(i,0)*sab(0)+ab(i,1)*sab(1) Next i ' Make temporary significand (q()) positive x=0 For i=0 To 5 x=x*2^-32+q(i) Next i If Sgn(x)=-1 Then For i=0 To 5: q(i)=-q(i): Next i Call Finish96A(y(), q(),eab(0),Sgn(x),xbd) Function=Convert96A(y()) Exit Function End Function Function Multiply96A(y() As Dword, _ ' Output array a() As Dword, b() As Dword) As Ext ' Input arrays ' y() <-- a() * b() ' Precede with Reciprocal96A() to divide. ' a(), b(), and y() are DWord arrays dimensioned (0:3) ' a() and/or b() can be in extended-precision or long format as, e.g., ' Dim aa(0) As Ext At VarPtr(a(1)):aa(0)=xxx## or a(0)=xxx& ' ' Methodology: ' As extended: <-Unrounded extended-precision-> Maxima: ' DWord: [Lo(Word,a(3))] a(2) a(1) a(0) 2^32-1 ' DWord: [Lo(Word,b(3))] b(2) b(1) b(0) 2^32-1 ' -------------------- ' Multiply DWords: ' QWord*: p(2,0) p(1,0) p(0,0) 2^64-2^33+1 ' QWord*: p(2,1) p(1,1) p(0,1) 2^64-2^33+1 ' QWord*: p(2,2) p(1,2) p(0,2) 2^64-2^33+1 ' ======================================= ' * As QWord datatype doesn't exist in PB, will use Quad and assembly. ' Split into DWords: ' DWord: ph(2,0) ph(1,0) ph(0,0) 2^32-2 ' DWord: pl(2,0) pl(1,0) pl(0,0) 2^32-1 ' DWord: ph(2,1) ph(1,1) ph(0,1) 2^32-2 ' DWord: pl(2,1) pl(1,1) pl(0,1) 2^32-1 ' DWord: ph(2,2) ph(1,2) ph(0,2) 2^32-2 ' DWord: pl(2,2) pl(1,2) pl(0,2) 2^32-1 ' ----------------------------------------------- ' Sum DWords: ' Quad: s(5) s(4) s(3) s(2) s(1) s(0) ~5*2^32 ' Max. s(): n-2 3*n-5 5*n-8 5*n-6 3*n-4 n-1 (n=2^32) ' ============================================= ' Split into DWords (sh(5)=sh(0)=0: ' DWord: sh(4) sh(3) sh(2) sh(1) ' Max. sh():2 4 4 2 ' DWord: sl(5) sl(4) sl(3) sl(2) sl(1) sl(0) ' Max. sl():n-2 n-1 n-1 n-1 n-1 n-1 (n=2^32) ' --------------------------------------------- ' Sum DWords: ' Quad: s(5) s(4) s(3) s(2) s(1) s(0) ' DWord: t(5) t(4) t(3) t(2) t(1) t(0) 2^32-1 ' ============================================== ' Normalize: ' DWord: t(5) t(4) t(3) ' Round: <-Unrounded extended-precision-> ' DWord: [Lo(Word,y(3))] y(2) y(1) y(0) 2^32-1 ' ======================================== ' Call CheckOut96A(y()) Call Convert96A(a()): Call Convert96A(b()) Local i,j,sg As Long Local ea,eb As Long ' Exponents unbiased Local d,d0,d1 As Dword Local g As Ext ' Significand terms only (no exponent nor sign): Dim wa(4) As Word At VarPtr(a(3)) Dim p(2,2) As Local Quad Dim ph(2,2) As Local Dword Dim pl(2,2) As Local Dword Dim s(5) As Local Quad Dim sh(5) As Local Dword Dim sl(5) As Local Dword d=Bit(a(0),%FlagNeg96) Xor Bit(b(0),%FlagNeg96) If IsTrue(d) Then sg=-1 Else sg=+1 ' Check if NaN If Bit(a(0),%FlagNaN96) Or _ Bit(b(0),%FlagNaN96) Or _ (Bit(a(0),%FlagInf96) And _ Bit(b(0),%FlagZer96) ) Or _ (Bit(b(0),%FlagInf96) And _ Bit(a(0),%FlagZer96) ) Then Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagNaN96 For i=%FlagNor96-1 To %FlagNeg96: Bit Set y(0),i: Next i ' NaN Function=Convert96A(y()) Exit Function End If ' Check if Zero If Bit(a(0),%FlagZer96) Or Bit(b(0),%FlagZer96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 If d Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) Exit Function End If ' Check if Infinity If Bit(a(0),%FlagInf96) Or Bit(b(0),%FlagInf96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 If d Then Bit Set y(0),%FlagNeg96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i ' Inf Function=Convert96A(y()) Exit Function End If ' Done with special cases ' Neither factor equals zero so not all t() will be zero. For i=0 To 2:For j=0 To 2 d0=a(i):d1=b(j) ' Do this in assembler since no QWord datatype in PB. !mov eax,d0 !mov edx,d1 !mul edx ' Unsigned multiply edx:eax <-- eax * edx !mov d0,edx !mov d1,eax ph(i,j)=d0 pl(i,j)=d1 s(i+j+1)=s(i+j+1)+ph(i,j) s(i+j+0)=s(i+j+0)+pl(i,j) Next j:Next i ' ea=exponent treating (a(2),a(1),a(0)} as integer ea=(a(3) And &H00007FFF???)-&H00003FFF???-95 If ea=-16478 Then Incr ea ' Denormal (&H3FFF+95=16478) ea=ea+32 ' eb=exponent treating (b(2),b(1),b(0)} as integer eb=(b(3) And &H00007FFF???)-&H00003FFF???-95 If eb=-16478 Then Incr eb ' Denormal (&H3FFF+95=16478) eb=eb+32 ' ea+eb=exponent treating (t(5),t(4),t(3),t(2),t(1),t(0)} as integer Call Finish96A(y(), s(),ea+eb+64,sg,0) Function=Convert96A(y()) Exit Function End Function
Sgn(x)*Sqr(Abs(x)).
Code:
Function Sign96A(a() As Dword, Opt b() As Dword) _ ' Input arrays As Long ' Output long ' Output <-- Sgn(a()-b()) (or Sgn(a()) if b() not passed) Call Convert96A(a()) Local i,j As Long If Bit(a(0),%FlagNaN96) Then ? "Error NaN in Sign":WaitKey$ If ArrayAttr(b(),0) Then ' b() was passed. Return Sgn(a()-b()) Call Convert96A(b()) If Bit(b(0),%FlagNaN96) Then ? "Error NaN in Sign":WaitKey$ If Bit(a(0),%FlagZer96) And Bit(b(0),%FlagZer96) Then Function=0: Exit Function End If j=IsTrue Bit(a(0),%FlagNeg96) i=j Xor IsTrue(Bit(b(0),%FlagNeg96)) If i Then i=1+2*j Else j=1+2*j If (a(3)And&H00007FFF???)>(b(3)And&H00007FFF???) Then i=+j If (a(3)And&H00007FFF???)<(b(3)And&H00007FFF???) Then i=-j If i Then Exit If If a(2)>b(2) Then i=+j If a(2)<b(2) Then i=-j If i Then Exit If If a(1)>b(1) Then i=+j If a(1)<b(1) Then i=-j If i Then Exit If If a(0)>b(0) Then i=+j If a(0)<b(0) Then i=-j End If Else ' b() was not passed. Return Sgn(a()) If Bit(a(0),%FlagZer96) Then i=0 Else i=1-2*Bit(a(0),%FlagNeg96) End If End If Function=i End Function Function Copy96A(y() As Dword, _ ' Output array a() As Dword) As Ext ' Input array ' y() <-- a() Call CheckOut96A(y()) Call Convert96A(a()) Local i As Long For i=0 To 3: y(i)=a(i): Next i Function=Convert96A(y()) End Function Function Negative96A(y() As Dword, _ ' Output array a() As Dword) As Ext ' Input array ' y() <-- -a() Call CheckOut96A(y()) Call Convert96A(a()) Call Copy96A(y(), a()) If Bit(y(0),%FlagNaN96)=0 Then Bit Toggle y(0),%FlagNeg96 Function=Convert96A(y()) End Function Function Absolute96A(y() As Dword, _ ' Output array a() As Dword) As Ext ' Input array ' y() <-- |a()| Call CheckOut96A(y()) Call Convert96A(a()) Call Copy96A(y(), a()) If Bit(y(0),%FlagNaN96)=0 Then Bit Reset y(0),%FlagNeg96 Function=Convert96A(y()) End Function Function Power96A(y() As Dword, _ ' Output array a() As Dword, _ ' Input array ByVal n As Long) As Ext ' Integer power ' y() <-- a()^n ' n is restricted to long integer ' Although 0^0=1 in PB, here 0^0=NaN Call CheckOut96A(y()) Call Convert96A(a()) Local i As Long Local sg As Long ' Sign of n Dim b(3) As Dword ' b()=|a{}| Dim t(3) As Dword ' Tempoary in case y=a If Bit(a(0),%FlagNaN96) Then Call Copy96A(y(), a()) Function=Convert96A(y()) Exit Function End If If Bit(a(0),%FlagInf96) Then Reset y() Bit Set y(0),%FlagSig96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 Select Case Sgn(n) Case -1 Bit Set y(0),%FlagZer96 Case 0 Bit Set y(0),%FlagNaN96 For i=%FlagNor96-1 To %FlagNeg96: Bit Set y(0),i: Next i Case +1 Bit Set y(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i End Select Function=Convert96A(y()) Exit Function End If If Bit(a(0),%FlagZer96) Then Reset y() Bit Set y(0),%FlagSig96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 Select Case Sgn(n) Case -1 Bit Set y(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i Case 0 Bit Set y(0),%FlagNaN96 For i=%FlagNor96-1 To %FlagNeg96: Bit Set y(0),i: Next i Case +1 Bit Set y(0),%FlagZer96 End Select Function=Convert96A(y()) Exit Function End If ' From here use t() in case y() overlays a() Bit Set t(0),%FlagSig96 For i=%FlagNor96 To %FlagNeg96-2: Bit Set t(0),i: Next i ' One i=Abs(n) Call Absolute96A(b(), a()) ' Working array (>0) If n<0 Then Call Reciprocal96A(b(), b()) Do Until i=0 Or _ IsTrue(Bit(t(0),%FlagInf96)) Or _ IsTrue(Bit(t(0),%FlagZer96)) If Bit(i,0) Then Call Multiply96A(t(), b(),t()) End If Shift Right i,1 If i>0 Then Call Multiply96A(b(), b(),b()) If Bit(b(0),%FlagInf96) Then Reset t() Bit Set t(0),%FlagSig96 Bit Set t(0),%FlagInf96 For i=%FlagNor96 To %FlagNeg96-1:Bit Set t(0),i:Next i Exit Do End If If Bit(b(0),%FlagZer96) Then Reset t() Bit Set t(0),%FlagSig96 Bit Set t(0),%FlagZer96 Exit Do End If End If Loop Copy96A(y(),t()) ' Only way answer is negative is if a() is negative and n is odd. If Bit(a(0),%FlagNeg96)+Bit(n,0)=2 Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) End Function Function Reciprocal96A(y() As Dword, _ ' Output array a() As Dword) As Ext _ ' Input array ' y() <-- 1/a() ' Methodology: If a*y=1, y=r+x, where r~1/a then ' a*r+a*x=1, x=(1-a*r)/a, y=r+(1-a*r)*r ' ' Examples of significands for values near 2^n: ' Significand: &B1000...0001 &B1111...1111 &B1000...0000 ' a: 1 + 2^-63 1 - 2^-64 1 ' y: 1 - 2^-63 1 + 2^-64 1 ' r: &B1111...1110 &B1000...0000 &B1000...0000 ' r: 1 - 2^-63 1 1 ' a*r: 1 - 2^126 = 1 1 - 2^-64 1 ' 1-a*r: 0 2^64 0 ' (1-a*r)*r=x: 0 2^64 0 ' r+x=y: 1 - 2^-63 1 + 2^-64 1 (all OK) ' ' 3.6451995318824746025284059336194e-4951 = smallest denormal ' 8.4052578577802337661123444847897e-4933 = 1/largest normal ' 3.3621031431120935058981578641335e-4932 = largest denormal ' 3.3621031431120935062626778173218e-4932 = smallest normal ' 2.9743287383930794127143983165700e+4931 = 1/smallest normal ' 1.1897314953572317650212638530310e+4932 = largest normal ' Call CheckOut96A(y()) Call Convert96A(a()) Local i As Long Dim aa(0) As Ext At VarPtr(a(1)) Dim r(3) As Dword ' Approximate reciprical Dim rr(0) As Ext At VarPtr(r(1)) Dim o(3) As Dword ' One Dim t(3) As Dword ' Temporary scratch Dim x(3) As Dword ' Correction If Bit(a(0),%FlagNaN96) Then Call Copy96A(y(), a()) Function=Convert96A(y()) Exit Function End If If Bit(a(0),%FlagInf96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) Exit Function End If rr(0)=1/aa(0) ' r= approximately 1/a ' rr(0) can be denormal (OK) or infinite (not OK) Call Convert96A(r()) If Bit(a(0),%FlagZer96) Or Bit(r(0),%FlagInf96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i Function=Convert96A(y()) Exit Function End If Call Multiply96A(t(),a(),r()) ' t=a*r Call Negative96A(t(),t()) ' t=-a*r o(0)=1: Call Convert96A(o()) ' o= exactly 1 Call Add96A(t(),o(),t()) ' t=1-a*r Call Multiply96A(x(),r(),t()) ' x=r*(1-a*r) Call Add96A(y(),r(),x()) ' y=r+r*(1-a*r) Function=Convert96A(y()) End Function Function SquareRoot96A(y() As Dword, _ ' Output array a() As Dword) As Ext ' Input array ' y() <-- Sgn(a())*Sqr(Abs(a())) ' This is differant than PB for negative arguments. ' ' Methodology: ' If y*y=+a, y=r+x, where r*r~a then ' r*r+2*r*x~a, x=(a-r*r)/(2*r), y=r+(a/r-r)*(1/2) ' Call CheckOut96A(y()) Call Convert96A(a()) Local i As Long Dim b(3) As Dword Dim bb(0) As Ext At VarPtr(b(1)) Dim r(3) As Dword ' Approximate root Dim rr(0) As Ext At VarPtr(r(1)) Dim t(3) As Dword ' Temporary scratch Dim tt(0) As Ext At VarPtr(t(1)) Dim x(3) As Dword ' Correction If Bit(a(0),%FlagNaN96) Then Call Copy96A(y(), a()) Function=Convert96A(y()) Exit Function End If If Bit(a(0),%FlagZer96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 Function=Convert96A(y()) Exit Function End If If Bit(a(0),%FlagInf96) Then Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 If Bit(a(0),%FlagNeg96) Then Bit Set y(0),%FlagNeg96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i Function=Convert96A(y()) Exit Function End If ' If 0<x<1, x<Sqr*(x)<1; if 1<x, 1<Sqr(x)<x so function is always ' nearer to Sgn(x) than x unless x=(-1,0,+1). There is no ' migration towards zero nor infinity. Call Absolute96A(b(),a()) ' b=|a| rr(0)=Sqr(Abs(bb(0))) ' r= approximately Sqr(|a|) Call Multiply96A(t(),r(),r()) ' t=r*r Call Negative96A(t(),t()) ' t=-r*r Call Add96A(x(),b(),t()) ' x=|a|-r*r Call Add96A(t(),r(),r()) ' r=2*r Call Reciprocal96A(t(),t()) ' t=1/(2*r) Call Multiply96A(x(),x(),t()) ' x=(|a|-r*r)/(2*r) Call Add96A(y(),r(),x()) ' y=r+(|a|-r*r)/(2*r) If Bit(a(0),%FlagNeg96) Then Negative96A(y(),y()) ' Xfer sign Function=Convert96A(y()) ' y=Sqr(Abs(a))*Sgn(a) End Function Function Subtract96A(y() As Dword, _ ' Output array a() As Dword, b() As Dword) As Ext ' Input arrays ' y <-- a-b Call CheckOut96A(y()) Call Convert96A(a()): Call Convert96A(b()) Dim t(3) As Dword ' Leave it to other functions to catch NaN's etc. Call Negative96A(t(), b()) Call Add96A(y(), a(),t()) Function=Convert96A(y()) End Function Function Divide96A(y() As Dword, _ ' Output array a() As Dword, b() As Dword) As Ext ' Input arrays ' y <-- a/b Call CheckOut96A(y()) Call Convert96A(a()): Call Convert96A(b()) Dim t(3) As Dword ' Leave it to other functions to catch NaN's etc. Call Reciprocal96A(t(), b()) Call Multiply96A(y(), a(),t()) Function=Convert96A(y()) End Function
Code:
Sub BankersRounding96A(y() As Dword) ' For internal use only ' Handles resetting %FlagInc96 when y(0)=&H80000000??? ' Tweaking y(0) to &H7FFFFFFF??? conforms to bankers' rounding ' without introducing any significant error in computation. ' Calling Convert96a() might be unwise since this is internal. ' so just check bounds and %FlagSig96. Call CheckOut96A(y()) If Bit(y(0),%FlagSig96)=0 Then ? "Error Not 96-bit sig.":WaitKey$ If y(0)=&H80000000??? Then ' For denormals, PB never rounds up exact halves ' but here just check if y(1) is odd and if so ' make y(0)=&H7FFFFFFF??? resetting &FlagInc96. If Bit(y(1),0) Then Decr y(0) ' Next If branch reverses the above Decr if denormal. ' PB (or likely the Pentium 4 FPU) never rounds up ' a denormal if the excess is exactly one half. ' Since this is inconsistent with PB documentation ' it is commented for now. '----------------------------------------' ' If Bit(y(0),%FlagDen96) Then Incr y(0) ' '----------------------------------------' End If End If End Sub Sub CheckOut96A(y() As DWord) ' Check dimensions of output array If ArrayAttr(y(),3)<>1 Then ? "Error Output dimensions":WaitKey$ If LBound(y)<>0 Or UBound(y)<>3 Then ? "Error Bounds 0:3":WaitKey$ End Sub Sub Finish96A( _ ' For internal use only y() As Dword, _ ' y(3) = result in 96-bit significand format q() As Quad, _ ' q(5) = temporary significand array (quad) ec As Long, _ ' ec = exponent of y(), the result sg As Long, _ ' sg = sign of result = -1/+1 xbd As Long) ' xbd = -1 if excess set bits were discarded ' Terminal sub for Multiply96A() and Add96A() If ArrayAttr(q(),3)<>1 Then ? "Error Quad dimensions":WaitKey$ If LBound(q)<>0 Or UBound(q)<>5 Then ? "Error Bounds 0:5":WaitKey$ Call CheckOut96A(y()) Local i As Long Local d As Dword Local f As Ext ' Result of Round_192_to_96 sub Dim t(5) As Dword ' q(5) = temporary significand array (dword) ' Variables and arrays only used after y(3) established ' in last section rounding (0:5) array to (0:3) array: Local IncrSig As Long ' True (-1) if rounding incr. significand Local Two32 As Quad: Two32=&H0000000100000000&& ' constant 2^32 Dim yy(0) As Ext At VarPtr(y(1)) ' Approx. result Local yy As Ext ' Refined yy(0) for returned extended result Dim s64(0) As Quad At VarPtr(yy) ' 64-bit significand of f Dim w(4) As Word At VarPtr(yy) ' w(4)=sign & exponent of f If sg=0 Then ? "Unexpected zero":WaitKey$ ' NB: In what follows, Mak(Quad,0,1)=&H100000000&& If q(5)<0 Then ?"Error q(5) sign":WaitKey$ For i=0 To 4 Do While q(i)<0 ' Negative only only if called from Add Decr q(i+1) q(i)=q(i)+Two32 Loop Do While q(i)=>Two32 Incr q(i+1) q(i)=q(i)-Two32 Loop Next i Do While q(5)=>Two32 For i=5 To 0 Step -1 If i Then If Bit(q(i),0) Then q(i-1)=q(i-1)+Two32 Shift Right q(i),1 Next i Incr ec Loop For i=0 To 5 If q(i)<0 Or q(i)=>Two32 Then ?"Error range":WaitKey$ Next i ' As 3*32=96 and 5*32=160, the maximum product has ' q(5)=(2^96-1)*(2^96-1)\2^160=(2^192-2^97+1)\2^160=2^32-2<2^32 For i=1 To 5 If Hi(Dword,q(i)) Then ? "Overflow":WaitKey$ t(i)=q(i) Next i ' Normalize t() ' Maximum t(5)=t(4)=&HFFFFFFFF???; t(3)=&HFFFFFFFE???; ' t(2)=t(1)=&H00000000???; t(0)=&H00000001??? ' Normalize t() Do Until t(5) For i=5 To 1 Step -1: t(i)=t(i-1): Next i Reset t(0) ec=ec-32 Loop Do Until Bit(t(5),31) For i=5 To 0 Step -1 Shift Left t(i),1 If i Then If Bit(t(i-1),31) Then Bit Set t(i),0 Next i Decr ec Loop ec=ec+63+&H3FFF?? ' Apply biases ' Now ec=>-16445 then normal result (=> 2^-16382=2^-&H3FFE) ' and ec<=-16446 then denormal result (< 2^-16382) Select Case ec Case <=-96 ' Zero ec=-96 Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagZer96 If sg=-1 Then Bit Set y(0),%FlagNeg96 Exit Sub Case -95 To 0 ' Denormal ' (&H0001 8000 ...)/2=&H0000 4000 ... (not &H0000 8000 ...) ' so do one extra shift right. Decr ec Do For i=0 To 5 Shift Right t(i),1 If i<5 Then If Bit(t(i+1),0) Then Bit Set t(i),31 Next i Incr ec Loop Until ec=0 y(3)=0 Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagDen96 If sg=-1 Then Bit Set y(0),%FlagNeg96 Case =>&H00007FFF& ' Infinity ec=&H00007FFF& Reset y() Bit Set y(0),%FlagSig96 Bit Set y(0),%FlagInf96 If sg=-1 Then Bit Set y(0),%FlagNeg96 For i=%FlagNor96 To %FlagNeg96-1: Bit Set y(0),i: Next i Exit Sub Case Else y(3)=ec ' Exponent Bit Set y(0),%FlagSig96 If sg=-1 Then Bit Set y(0),%FlagNeg96 End Select ' Rounds (0:5) array to (0:3) array For i=0 To 2: y(i)=t(i+3): Next i yy=yy(0) IncrSig=t(2)>&H80000000??? If t(2)=&H80000000??? Then If t(1) Then IncrSig=-1 Else If t(0) Then IncrSig=-1 Else IncrSig=xbd Or IsTrue(Bit(y(0),0)) End If End If End If If IncrSig Then Call RoundUp96A(y()) Call BankersRounding96A(y()) End Sub Sub RoundUp96A( _ ' For internal use only y() As Dword) ' Input array ' Rounds up 96-bit significand Call CheckOut96A(y()) Local d As Dword Incr y(0) If y(0)=0 Then ' If y(0) initially &HFFFFFFFF???, rounding starts here Incr y(1) If y(1)=0 Then Incr y(2) If y(2)=0 Then ' Pass 2^n boundary Incr y(3) Bit Set y(0),%FlagNor96 If (y(3) And &H00007FFF???)=&H00007FFF??? Then Bit Set y(0),%FlagInf96 End If ElseIf y(2)=&H80000000??? Then ' Denormal becomes normal Bit Reset y(0),%FlagDen96 Incr y(3) End If End If End If End Sub
Code:
%MaxN=... Function MyFunction(N As DWord) As Ext Dim MyArray(%MaxN) As Ext At CodePtr(MyLabel) If N<=%MaxN Then Function=MyArray(N) Exit Function #Align 2 ' Optional PBCC version 5+ MyLabel: !dw &H...,&H...,&H...,&H..., &H... ' MyArray(0) | Output96A() !dw &H...,&H...,&H...,&H..., &H... ' MyArray(1) | can create !dw &H...,&H...,&H...,&H..., &H... ' MyArray(2) | these with !dw &H...,&H...,&H...,&H..., &H... ' MyArray(3) | the use of ... | Clipboard !dw &H...,&H...,&H...,&H..., &H... ' MyArray(%MaxN) | statements. End Function
________________
Mark Longley-Cook
Comment