Oops, XL did it again

Posted in Computer, Microsoft Excel, multi-precision calculations, Software, Software tools, Spreadsheet, xNumbers on June 22, 2014 by giorgiomcmlx

I really do love software that adheres to the "Law of Least Astonishment (LoLA)" (please see §4.1). Unfortunately, XL definitely does not belong to this group. Quite recently, I came across a violation of the LoLA by XL that is caused by the inconsistent manner of how functions evaluate numerical cells. MS have stated relatively clearly that the binary double-precision floating-point format (DP) is used internally in order to store numerical cell values. On the other hand, in the model-view-controller (MVC) of XL, i.e. a cell, a decimal string comprising at maximum 15 digits is used to display the numerical value of a cell to humans. To my great astonishment, it has turned out that different functions belonging to the same category in the official list of functions of XL2013 may operate on these two very different representations of numerical values. Of course, this can have severe consequences, as a precision of 15 decimal digits roughly corresponds to a precision of only 50 bits, whereas the 53 bits of precision of the IEEE754 DP format roughly correspond to 16 decimal digits.

I suppose virtually nobody regularly working with numbers will expect a cell formula like "= (D5 - FLOOR.PRECISE(D5;1))" to return a negative value in case the value of D5 is positive. In the official description of this function MS states: "Returns a number that is rounded down to the nearest integer or to the nearest multiple of significance. Regardless of the sign of the number, the number is rounded down." Moreover, the function FLOOR.PRECISE is listed in the "math and trigonometry" section of the aforementioned list, so that IMHO people will expect this function to be as accurate as the functions SQR and SIN are, for example. Yet, look what happens if the value of D5 is close to but still smaller than a positive integer. I’ll demonstrate that now in case of the integer 12 (twelve), which in DP big-endian byte order equals 40 28 00 00 00 00 00 00.

The hex representation of a numerical cell value as well as its equivalent exact decimal string were obtained by harnessing the commands "= Dbl2Hex(CellRef)" and "= DCStr(CellRef;100)", resp., which require an appropriate version of the FOSS add-in xNumbers to be installed in XL.

The largest DP value that is smaller than 12, but for which the XL2013x64 functions FLOOR, FLOOR.MATH, FLOOR.PRECISE, INT, ROUNDDOWN, and TRUNC return the mathematically correct result of 11 (eleven), has a hex notation of 40 27 FF FF FF FF FF E3. For all the other 28 DP values from 40 27 FF FF FF FF FF E4 up to 40 27 FF FF FF FF FF FF, the wrong value of 12 (twelve) is returned by all of these functions. Exactly for these 28 values the results of the formula "= (CellRef - FLOOR.PRECISE(CellRef;1))" turn out to be negative: about -4.97E-14 is returned by XL for the input corresponding to 40 27 FF FF FF FF FF E4, and about -1.776E-15 for 40 27 FF FF FF FF FF FF. If the mathematically unnecessary outer "magic parentheses" are omitted in the formula, XL tries to hide the dirty little secrets of its floating-point core from you by artificially zeroing the results for all values from 40 27 FF FF FF FF FF F9 to 40 27 FF FF FF FF FF FF. IMHO, both of these undocumented effects clearly exclude XL from the set of numerical software that might be considered for any serious computation (VBA is different!).

The observed results could be most easily explained by assuming that all these functions do not use the binary DP representation of the cell content in order to calculate their results, but operate on the numerical string comprising 15 decimal digits in the MVC instead. At least we have:
40 27 FF FF FF FF FF E3 hex =
XL: 1.19999999999999E+01
xN: 1.199999999999994848565…E+01
40 27 FF FF FF FF FF E4 hex =
XL: 1.20000000000000E+01
xN: 1.19999999999999502620…E+01
From the given xN values it can be seen that XL2013x64 correctly rounds the binary DP values to strings comprising 15 decimal digits, and that the functions in doubt obviously use this decimal string representation in order to calculate their result. Exactly the same kind of inaccurate results is delivered by the functions that are meant to round up, i.e. CEILING, CEILING.MATH, CEILING.PRECISE (and even ISO.CEILING), and ROUNDUP. In this example, all the 28 results corresponding to input values from 40 28 00 00 00 00 00 01 to 40 28 00 00 00 00 00 1C are flawed. Moreover, the functions ROUND and MROUND are also subject to this kind of flaw if their arguments are close to values at which their result should change.

On the other hand, the functions EVEN, MOD, ODD, and QUOTIENT deliver results that accurately account for the least bit of their arguments. Thus, in order to define an accurate version of the mathematical floor function, you might want to use = QUOTIENT(CellRef;1) + IF(MOD(CellRef;1)=0;0;(SIGN(CellRef)-1)/2), whereas for an accurate version of the mathematical ceiling function, you might want to use = QUOTIENT(CellRef;1) + 1 + IF(MOD(CellRef;1)=0;-1;(SIGN(CellRef)-1)/2). Or, you break free and look out for another numerically much more comprehensive programming language than VBA is, but that can be connected to XL relatively easy. At least, that is what I shall do from now on.

The decimal data sub-type in XL-VBA

Posted in Computer, Microsoft Excel, Programming Languages, Software, Spreadsheet, VBA on June 9, 2014 by giorgiomcmlx

In VBA, the most straightforward way to add some extra precision to calculations with doubles is to try to harness the sub-type decimal of the data type variant. According to the VBA language specification provided by MS quite recently, a variable of sub-type decimal comprises a numerator of 96 bits (or 12 octets) of precision that encode an unsigned integer in binary format, and a denominator comprising 8 bits (or 1 octet) that encode the power of ten by which the numerator has to be divided in order to yield the mathematical value that is represented by the variable. MS further state that there also is a sign bit, and that a decimal is "A rational number represented in a 14 byte data structure". I started to doubt this statement when reading this information on the decimal sub-type (in German). There, a decimal is said to be a data structure of 16 bytes, and moreover, a rather peculiar order of the 12 octets of the numerator is
mentioned. Consequently, I started to check out the details by myself.

For that, in my VBA7.1 I declared a call to the kernel function RtlMoveMemory that can be accessed via the system library kernel32.dll as follows:
Private Declare PtrSafe Sub MoveMem Lib "kernel32" Alias "RtlMoveMemory" (PointerToMemory As Any, PointerToBuffer As Any, ByVal AmountOfOctets As Long)

In order to call this function properly from within a VBA UDF, an array of Bytes of desired length has to be declared. Finally, it will contain an exact copy of the block of computer memory that starts at address PointerToMemory and comprises AmountOfOctets octets. For that, the VBA UDF called HexDump, which is of return type String, is defined like this:
Function HexDump$(ByVal PointerToBufferStart As LongLong, ByVal myOctets As Long)
In HexDump, the amount of octets to be returned is rounded up to the next equal or larger multiple of 16 and
assigned to the variable NumOctets. Then an array of type Byte and of length NumOctets is declared:
ReDim BytArr(1 To NumOctets) As Byte .
When calling MoveMem, the pointer to the start of BytArr in memory is obtained implicitly by simply submitting the reference to its first element:
Call MoveMem(BytArr(1), ByVal PointerToBufferStart, ByVal NumOctets).
The rest of MoveMem is only concerned with creating a proper string from all the elements of BytArr in a form that is easily perceivable by humans.

I then created some decimals according to the following scheme:
Function some_decimals$()

Dim …, d18, d19, d20, d21, d22
d18 = CDec("1234567890123456789012345")
d19 = CDec("123456789012345678901234567")
d20 = CDec("12345678901234567890123456789")
d21 = CDec("1.2345678901234567890123456789")
d22 = CDec("-1.2345678901234567890123456789")
some_decimals = HexDump(VarPtr(d22), 528&)

Quite surprisingly for me at least, the position of the variables in memory is found to be in reversed order of declaration, which means the pointer to the variable declared least has to be submitted to HexDump in order to get back a block of memory that contains all the variables declared. The relevant part of the string returned by HexDump, which encodes the values of all octets in hexadecimal notation, reads as follows:

0E 00 1C 80 32 1B E4 27 15 81 39 6E B1 C9 BE 46
00 00 00 00 00 00 00 00 0E 00 1C 00 32 1B E4 27
15 81 39 6E B1 C9 BE 46 00 00 00 00 00 00 00 00
0E 00 00 00 32 1B E4 27 15 81 39 6E B1 C9 BE 46
00 00 00 00 00 00 00 00 0E 00 00 00 FD 1E 66 00
87 4B 9F 2C A8 F2 58 F1 00 00 00 00 00 00 00 00
0E 00 00 00 6E 05 01 00 79 DF E2 3D 44 A6 36 0F

As I shall show in detail, the first 16 octets indeed encode the value of d22 as a decimal correctly. Yet, the effective amount of memory claimed by each of consecutively declared decimals is neither 14 nor 16 octets, but equals 24 octets, at least on my Win7x64 XL15x64 PC. Consequently, in order to catch the entire block of 22 decimals, an amount of at least 21*24+16 = 520 octets is needed, which is then internally rounded up to 33*16 = 528 octets.

In the German resource mentioned above, it is claimed that the 16 octets of a decimal start with an identifier octet having a value of 14 dec = 0E hex, and that the following octet remains unused. The 3rd octet is said to encode the power of 10 of the denominator (here: 28 dec = 1C hex or 0 dec = 0 hex), and the 4th octet is claimed to encode a positive sign by a value of 0, and a negative one by 128 dec = 80 hex. As can easily be seen from the list given above, all of these statements hold true. Therefore, the trailing 12 octets are needed to encode the 96 bits of the numerator. In order to understand the values returned by HexDump, I harnessed xNumbers for converting “12345678901234567890123456789” into a binary string. This string is found to comprise 94 bits, so that it is presented here on three lines each of which gives 4 octets, starting with the most significant one. We have:

OOIO-OIII IIIO-OIOO OOOI-IOII OOII-OOIO
2-7 E-4 1-B 3-2
OIOO-OIIO IOII-IIIO IIOO-IOOI IOII-OOOI
4-6 B-E C-9 B-1
OIIO-IIIO OOII-IOOI IOOO-OOOI OOOI-OIOI
6-E 3-9 8-1 1-5

so that the 94 bits in big endian byte order and in hexadecimal notation read:
27 E4 1B 32 46 BE C9 B1 6E 39 81 15
which in little endian byte order reads:
15 81 39 6E B1 C9 BE 46 32 1B E4 27 .
When this sequence is compared to that returned by HexDump, the peculiar byte order by which variable values of sub-type decimal store the numerator in computer memory becomes obvious:
32 1B E4 27 15 81 39 6E B1 C9 BE 46
The first block in memory comprises the 4 most significant consecutive octets (green), and the second one the 8 least significant (blue) of them, which is just the opposite of what is said in the mentioned German resource.

When comparing the 16 net octets of the values assigned to d20, d21, and d22, we have:

d22 = CDec("-1.2345678901234567890123456789")
d22 = 0E 00 1C 80 32 1B E4 27 15 81 39 6E B1 C9 BE 46
d21 = CDec("+1.2345678901234567890123456789")
d21 = 0E 00 1C 00 32 1B E4 27 15 81 39 6E B1 C9 BE 46
d20 = CDec( "+12345678901234567890123456789")
d20 = 0E 00 00 00 32 1B E4 27 15 81 39 6E B1 C9 BE 46

This means the power of 10 of the denominator, which equals 28 dec = 1C hex for both d22 and d21, and 0 dec/hex for d20, and that is encoded in the 3rd octet, can also be interpreted as the position of where to insert the decimal separator when counting positions of the unsigned numerator from right to left. Finally, there is a sign bit, but effectively, a decimal uses a sign octet.

Summarising, variables of sub-type decimal encode their values very in-efficiently. On the other hand, they allow for exact calculations with decimal fractions, which makes them the optimal data type in the area of finance. Due to the fact that only the basic algebraic operators support decimals, any function that is to be evaluated needs to be implemented from scratch using only the four basic operators +,-,·, and /, or calling other functions supporting decimals that have already been implemented. Not every double can be converted directly into a decimal due to either the limited range of the latter of at maximum ±(2^96-1) = ±7.9228162514264337593543950335E28 or their minimum resolution of 1E-28 or a combination of both effects. Moreover, only doubles having no more than 29 significant decimal figures can be exactly casted into decimals. When scaling doubles in order to make them fit the limits of decimals, only factors equalling integral powers of 2 should be used in order to leave the floating-point mantissa unchanged. This will be demonstrated in a forthcoming post.

XL2013x64: Extremely slow VBA in copied cells

Posted in Computer, Microsoft Excel, Software, Spreadsheet, Uncategorized, VBA on May 25, 2014 by giorgiomcmlx

As mentioned before, the execution speed of VBA UDFs in copied cells is painfully slow on my XL2013x64 PC. It takes a second or so to calculate the UDFs in a copied block of only about 10 cells. Yet, after being copied, the re-calculation of exactly these UDFs happens instantly; for example, if an input value is changed.

Switching XL into manual re-calculation mode (thanks to John Beyers for this hint) completely removed this slow-down, but I do not like to develop sheets this way. Then, I happened to come across this post by Charles Williams on a VBE window refresh bug in XL. Though it had obviously been there, I had not noticed the flickering of the borderline around the block of selected cells during UDF calculations until I finished reading this post. So I started to try to copy a small block of cells to a location that made the selected block of cells completely disappear from the view pane. And voila, VBA executes at full speed in automatic calculation mode.

Thus, MS seem to have built a refresh bug similar to that in the VBE window now into the selection indication process of XL2013x64, too.

VBA speed in XL2013x64 on Win7x64

Posted in Computer, Microsoft Excel, Software, Software tools, VBA, xNumbers on April 13, 2014 by giorgiomcmlx

When playing with my newly acquired XL2013x64 and xNumbers on a i7-16GB-RAM PC running Win7x64prof , I experienced a dramatic loss of speed as compared to that I used to observe with XL2010x64 on my old i5-4GB-RAM PC running an identical OS. Some search on the Web proved that this phenomenon is wide-spread, see e.g. here.

For example, when copying the formula =xSqr(xAdd(xMult("2";B4;200);"";200);200) (in cell c4; b4:b105 contain 1,2,3,…,101) down a hundered times over existing data in the range c5:c105, it takes an unbelievable amount of 7 seconds before the results appear. Moreover, this time span is nearly independent from the multi-threading settings. On the contrary, erasing all the data in cells c5:c105 manually before copying down the formula makes the results appear instantly. This phenomenon is reproducible, and it does not depend sensitively on the kind of a formula, at least on my PCs.

So quite obviously, writing to a non-emtpy cell is the primary cause for the observed speed penalty of VBA in XL2013. At the moment, I guess XL2013 saves any content it finds in a cell that is to be overriden due to some built-in web collaboration features. I’ll have to figure out how to switch them off in case my guess is correct.

cos in XL14 or how not to do argument reduction

Posted in Computer, Microsoft Excel, multi-precision calculations, Programming Languages, Software, Software tools, Spreadsheet, Uncategorized, VBA, xNumbers on March 12, 2014 by giorgiomcmlx

While playing around with formulae in the area of computational geometry, I’ve noticed a severe flaw in the way XL14 implements the cosine of numbers close to π/2. Let XL14 and VBA7 denote results obtained from calculation in XL14 and VBA7, resp.. All final results are double-precision floating-point numbers that are uniquely characterised by their hex notation, which can be obtained by the "dbl2hex" function of the XL-add-in xNumbers. The "dcstr" function of this add-in allows for obtaining the complete sequence of decimal figures corresponding to any binary floating-point number. XN60 refers to multi-precision results returned by the XL-add-in xNumbers version 6.0.5. The letter c indicates cosine values. For example, I found:

x_XL14_wks = 1.57079632679461E+00

x_XL14_dec = 1.57079632679460701183415949344635009765625
x_XL14_hex = 3F F9 21 FB 54 44 28 00
c_XL14_wks = 2.89607422244986E-13

c_XL14_dec = 2.89607422244986256743004560121335089206695556640625E-13
c_XL14_hex = 3D 54 61 1A 80 00 00 00
c_XL14_err = 496815802 ULPs
c_VBA7_wks = 2.89607397162198E-13

c_VBA7_dec = 2.89607397162198205938334649267199577700980517303008809903985820710659027099609375E-13
c_VBA7_hex = 3D 54 61 1A 62 63 31 46
c_VBA7_err = 0 ULPs
c_XN60_dec = 2.89607397162198193401344438286352963771616699886035390348378210024734360417199…E-13

Applying the "x2dbl3" function of xNumbers, which implements the to-nearest-ties-to-even rounding scheme, to c_XN60_dec yields a result identical to that obtained in VBA7 directly. So I had to find a reasonable explanation for the really huge error of the XL14 cosine function around π/2.

A first clue came from the many zeroes at the end of the hex representation of c_XL14, which usually indicates heavy cancellation. Yet, there was no explicit subtraction involved, and so a hidden subtraction had to be responsible. The most likely reason for such a hidden subtraction is a procedure called argument reduction. Nowadays, most complicated functions are calculated by an algebraically simple approximation that is valid over only a small range of all possible input parameters, so that parameter values not covered by the range of such an approximation have to be mapped using some kind of relation. In case of the cosine function, there is the mathematical identity cos(x) = sin(π/2-x), which may be used for π/4 ≤ x ≤ π/2 in order to map some range of cosine input parameters to the input range 0 ≤ x ≤ π/4 of the sine function.

A first try did not explain the result:

pio2hi_XL14_wks = 1.57079632679490E+00
pio2hi_XL14_dec = 1.5707963267948965579989817342720925807952880859375
pio2hi_XN60_dec = 1.57079632679489661923132169163975144209858469968755291…
pio2hi_XL14_hex = 3F F9 21 FB 54 44 2D 18
pio2lo_XL14_wks = 6.12323399573677E-17
pio2lo_XL14_dec = 6.12323399573676603586882014729198302312846062338790031…E-17
pio2lo_XN60_dec = 6.12323399573676588613032966137500529104874722961539082…E-17
pio2lo_XL14_hex = 3C 91 A6 26 33 14 5C 07

The difference between x_XL14_wks and pio2hi_XL14 is only about 3E-13, so that the sine of the mapped argument is exactly equal to the mapped argument. The most simple reduction formula (pio2hi – x_XL14) yields a result (hex: 3D 54 60 00 00 00 00 00) that is inaccurate and does not agree withXL14 as well. The more complicated reduction formula ((pio2hi – x_XL14) + pio2hi) yields the correct result (hex: 3D 54 61 1A 62 63 31 46), but of course does not agree with XL14, too. So I used my recent finding that XL14 harnesses the extDP format when dealing with two numbers (unfortunately, its formula parser does not support inlining, as I will explain in a forthcoming post) for a last try. The binary values of both π/2 and x_XL14 are given by (V: border between half-bytes, X: end of DP mantissa, Y: end of extDP mantissa):
first 9 half-bytes of π/2 :

+IV.IOOIVOOIOVOOOIVIIIIVIOIIVOIOIVOIOOVOIOOVOIOOV
+IV.IOOIVOOIOVOOOIVIIIIVIOIIVOIOIVOIOOVOIOOVOIOOV

first 9 half-bytes of x_XL14
The first 9 half-bytes of both numbers cancel out completely.
trailing 4 half-bytes of π/2 plus extDP section:

OOIOVIIOIVOOOIVIOOOXOIOOOIIOIOOYIIOO…

which will be rounded up to extDP

+OOIOVIIOIVOOOIVIOOOXOIOOOIIOIOIYO

from which the trailing 4 half-bytes of x_XL14 will be subtracted

-OOIOVIOOOVOOOOVOOOOXOOOOOOOOOOOYO
=OOOOVOIOIVOOOIVIOOOXOIOOOIIOIOIYO

Due to cancellation of many leading bits, the border indicators are misplaced now, but the bit sequence exactly matches that of the result obtained in XL14:

OOOOOIOIOOOIIOOOOIOOOIIOIOIO
IOIOOOIIOOOOIOOOIIOIOIOOOOOOO…
I.VOIOOVOIIOVOOOIVOOOIVIOIOVIOOOVOOOO…
    4    6    1    1    A    8    0…
c_XL14_hex = 3D 54 61 1A 80 00 00 00

This works for all examples I’ve checked. Summarising, the cosine function of XL14 is severely flawed around π/2 (at least). "Forensic doublelogy" indicates that this flaw is caused by an argument reduction process performed using the 64 mantissa bits of the extDP floating-point format, whereas the obviously exploited mathematical identity required a full double-double approach. Fortunately, no flaws around π/2 were found in the cosine function of VBA7. The following graph illustrates the situation.

graph of errors of cosine in XL14 around pi half

The value of pio2_extDP is slightly larger than the mathematical value of π/2 (contrary, the DP value pio2hi is too small, as can be seen from the fact that pio2lo is positive!!). Thus, in combination with the structure of the exploited mathematical identity, the errors of the cosine in XL14 around π/2 are not random, but show a strong systematic component: all the cosine values for x just below π/2 are too large, and all the cosine values for x just above π/2 are too small.

To SSE, or not to SSE, that is the question

Posted in binary, decimal, gnumeric, Maxima, Microsoft Excel, multi-precision calculations, Number processing, Number systems, Software tools, Spreadsheet, VBA, xNumbers on February 25, 2014 by giorgiomcmlx

Some time ago, I discovered that in a thirty year old floating-point (FP) library the value of the mathematical expression 1 – x2 was determined by calculating 2·(1-x)-(1-x)2 for 0.5 ≤ x ≤ 1. At a first glance, this method appeared rather peculiar to me, but then I realised that the Lemma of Sterbenz (see, e.g. page 45 of this slide collection on FP calculations) guarantees that 1-x is calculated in FP without any round-off error in the specified range. The multiplication by 2 does not introduce errors as well. Yet, there are still two sources of errors: firstly, in case (1-x) has more than 26 significant bits, (1-x)2 is not exactly representable in FP, and secondly, finally taking the difference might cause a round-off error as well.

In my tests though, this old formula outperformed all other versions on average, and so I shared it with developer friends. Surprisingly, some of them responded with a statement like "The results produced by this formula do not differ from those obtained using (1-x)·(1+x)." In the end, the reason for this paradox could be identified: all my tests were executed on SSE registers, calculations on which strictly adhere to the double-precision floating-point format, whereas the calculations not showing any improvement had harnessed the intermediate extended double-precision format of the x87 stack. Modern x87-compatible FPUs provide both of these two very different kinds of intermediate FP number storage.

The SSE have been invented in favour of fast parallel data processing, whereas the x87 double-extended stack provides both 11 additional mantissa bits and 4 additional exponent bits to intermediate calculations in order to try to maintain an overall accuracy of 53 bits and to reduce the probability of intermediate under-/over-flows as well. The decision whether a piece of software uses the x87 or the SSE instruction sets has to be made up at compile time, because the generated byte code is different. Contrary, both the rounding-mode and the exception-handling properties may be influenced at runtime by changing the so called floating-point control word (FPUCW). I could not find an article on Wikipedia, so please see e.g. this collection of or this introduction to the FPUCW world. I also have come across these imho nice summaries that are written from C-programmers points of view: The pitfalls of verifying floating-point computations as well as SciLab is not naive. Unfortunately, values of the FPUCW having an identical meaning differ across platforms. For example, a value of the MS-FPUCW of hex 8001F corresponds to a GNU-GCC-FPUCW value of hex 37F …  Fortunately, the ability to change the rounding mode during runtime enables programmers to harness interval arithmetic in order to compute reliable FP error estimates of their code for both the x87 and SSE instructions sets. The rounding mode can be chosen from the set UP, DOWN, CHOP, and TO NEAREST – TIES TO EVEN, the latter usually being the default setting. Only if the x87 instruction has been selected at compile time, the precision of the calculations may be changed during runtime as well.

The 11 additional significand bits of the x87 extDP format correspond to roughly 3 decimal figures, and thus are a valuable tool in order to minimise the effects of cancellation in intermediate calculations. The 4 additional exponent bits allow for an intermediate product of up to 16 huge (tiny) doubles without causing an over-(under-)flow. Yet, harnessing the extDP format bears one dirty little secret: when the result of an intermediate calculation is transferred from the x87 stack firstly into an extDP x87 register and finally into the DP main memory, the intermediate value is subject to a double-rounding procedure. In combination with the standard rounding mode "TO NEAREST – TIES TO EVEN" this can lead to a slightly wrong result in two situations: an intermediate calculation rounds up/down to an extDP value representing an exact DP tie that is again rounded up/down to a DP value. The following tables give an example:

decimal number    DP hex    comment
+18446744073709551616.0 43F0000000000000 exact DP Num1
+6143.0 40B7FF0000000000 exact DP Num2
+18446744073709557759.0 Num1 + Num2,
not representable in DP
+18446744073709555712.0 43F0000000000001 DP_lo: nearest DP
+18446744073709559808.0 43F0000000000002 DP_hi: next larger DP
-2047.0 error of DP_lo
+2049.0 error of DP_hi

When adding Num1 to Num2 on the x87 extDP stack, the hex representation of the final DP sum in memory will be 43F0000000000002 instead of the closest DP 43F0000000000001. This can be understood by analysing the binary representation of the numbers involved. In order to clearly show the rounding effects, some separators are used: a "V" indicates the border between any two DP half-bytes, each of which corresponds to a single character in the hex representation, an "X" shows where the DP bits end and the extDP bits begin, and finally a "Y" marks the end of the extDP bit sequence.

Num1
+IOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO.O
Num2
+IOIIIIIIIIIII.O
Num1 + Num2
+IOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOIOIIIIIIIIIII.O
same with separators
+IVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOIXOIIIIIIIIIIYI.O
1st step: round (up) to extDP
+IVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOIXIOOOOOOOOOOYO.O
2nd step: round DP tie (up) to even
+IVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOOOVOOIOXOOOOOOOOOOOYO.O

In this example, the addition creates a bit sequence, the mantissa of which exceeds the 64 bits of the extDP format. Thus, in a first step, the result is rounded to 64 bits. These extDP bits happen to represent a so-called tie (a number exactly half-way between the two adjacent DP numbers) of the DP format. In this example, the tie-breaking rule TO EVEN finally causes another rounding up procedure, so that the value stored in memory does not represent the DP number closest to the mathematical result of the summation, but the next larger one. Yet, it is still one of the two DP numbers bracketing the true value.

Only two different bit sequences among the 2^13 possible ones lead to this kind of error:
leadingbitsIXOIIIIIIIIIIYItrailingbits ++
leadingbitsOXIOOOOOOOOOOYOtrailingbits —

so that the odds ratio of this error should be expected to be close to 1/2^12 = 1/4096 ≈ 2.44E-4 . Of course, if the bit sequence of an intermediate result does not exceed the 64 bits of the extDP format, there will be only a single rounding step yielding a correct result. When analysing the results of 5.75E5 renormalisation steps of double-double numbers in VBA7, the sum of which is determined to exceed the 64 bit extDP limit, and that were introduced here at the end of this post, I found 1480 errors and thus an odds ratio of about 2.57E-4. Due to the relatively large sample size, I guess the deviation from the theoretical value is caused by non-randomness of the bit sequences involved.

Double-rounding is unavoidable when using the extDP registers of the x87 FPU, but it only slightly increases the representation error of about 0.5 ULPs by at maximum 0.05&percent;. Moreover, this only happens with an odds ratio of about 1/4096. Thus, this disadvantage of rarely occurring false double rounding is usually more than counterbalanced by the gain of 11 bits of intermediate precision in every calculation. However, there is a special kind of multi-precision software that cannot be run on the x87 stack due to the double rounding, i.e. all the double-double as well as quad-double software that use chained DP variables in order to store multi-precision values (see, e.g. this implementation). The double rounding on the x87 stack can lead to violations of the basic assumption of this kind of software, i.e., that the leading DP value always represents the DP value closest to the multi-precision value stored in a chain of doubles. Therefore, arbitrary precision tools being determined to use the x87 stack, for example Xnumbers, need to harness completely different algorithms for internal number storage. The numerical trick using a double-double representation of π/2 introduced in the aforementioned post works reliably, because the firstly calculated intermediate difference does not involve more than 64 significant bits.

The FP instruction set actually used by third-party numerical software is often unknown and needs to be determined by using test cases. In order to avoid problems with the digits chopping feature of most spreadsheet software, examples need to be created that use exact DP values corresponding to less than 14 significant decimal digits. The previously given example can be reformulated as sum of +1.8446744E+19 ≡ 43EFFFFFFDDAD230 and +73709557759 ≡ 4231296E97FF0000. In order to check the DP result for correctness without any debugging capabilities like those provided by Xnumbers for MS XL, e.g., both DP numbers that were summed up before now need to be successively subtracted from the sum, the bigger one at first, and the smaller one afterwards. It is important to assure that all intermediate differences are rounded to DP, which can most easily be achieved by assigning them to different DP variables, or different cells in a spreadsheet. Contrary to the first addition, the two final subtractions do not introduce any FP error, because they are subject to the lemma of Sterbenz. A final result of -2047 corresponds to calculations that have not been executed on the x87, whereas a value of +2049 indicates the opposite. Here are my PY3- and CAS-Maxima scripts:

Python 3     CAS-Maxima
num1 = +1.8446744E+19 +1.8446744E+19;
num2 = +73709557759.0 +73709557759.0;
sum12 = (num1+num2) %o1+%o2;
sum12sub1 = (sum12-num1) %-%o1;
sum12sub1sub2=(sum12sub1-num2) %-%o2;
print(sum12sub1sub2)

The following table summarises the results I found in some spreadsheet and scripting software running on MS WIN7x64prof. On Linux, all SW listed in the table that is available on that OS as well returns -2047.

spreadsheet SW   result        scripting SW   result
GNUmeric 1.12.9 +2049 PY2.7.3×64 -2047
GScalc 11.4.2 -2047 PY3.2.3×64 -2047
KSspread 9.1.0 +2049 R2.15.2×64 -2047
LOcalc 4.1.4.2 -2047 R3.0.0×32 +2049
LT SYM 3.0.1 -2047 R3.0.0×64 -2047
MS XL14/VBA7 +2049 FreeMat 4.2 -2047
SM PM 2012 -2047 Maxima x32 +2049

As you can see from the table, the situation is completely messed up. It seems to me as if on the Win32 platform the x87 with its extDP format was used by default, whereas on the Win64 platform the SSE part of the FPU was involved. Of course, any software can override the default setting, and MS XL64 as well as VBA7 indeed do so. Generally speaking, specifying the FPUCW as well as the kind of FPU during compilation is a very good idea in case of numerical libraries that were developed and optimised for only one of the two scenarios. In case these two setting are omitted and the unknown defaults are used at runtime, small deviations are to be expected only in case the software was designed not to harness the extDP precision of the x87 stack. Otherwise, SSE lacking the additional 11 bits of precision will very likely cause things like polynomial approximations to special functions to fail at runtime. In particular, source code of libraries originally designed for the x87 must not be re-compiled for SSE without closely inspecting the results. Libraries that were built before the advent of the x87 can be a very valuable source of inspiration for changes in the code that are likely to become necessary.

How to call a MSVCpp2010x64-DLL from XL14x64

Posted in C++, Compiler, Computer, Microsoft Excel, Programming Languages, Software, Software tools, Spreadsheet, VBA on December 20, 2012 by giorgiomcmlx

Belonging to the group of people that sociologists call "early adopters" usually is a burden with respect to newly invented technologies, and unfortunately, I cannot avoid this situation all of the time… Fortunately, nowadays it is possible to search for suitable keywords on the Internet, and quite often, this strategy is much more effective in order to solve a problem than reading through manuals. Rarely, the result of such a search will be a thick searchable manual that covers all aspects one is interested in. A recent example of this situation occurred when I wanted to call functions in a self-compiled C++ DLL from VBA7, i.e., the macro language than comes with the 64bit version of Microsoft Office 2010.

The first thing I realised was that it is impossible to call a 32bit DLL from a 64bit application directly (and vice versa). Then it turned out that the so-called express edition of Visual Studio 2010 does not contain 64bit versions of the compilers. Therefore, I had to follow the advice provided by MS and install the WindowsSDK version 7.1 in order to enable building x64 DLLs. Unfortunately, I got this finally working only after having previously de-installed all of the VS2010express stuff present on my PC, because otherwise, the installation of the WinSDK7.1 failed at an early stage (may be due to some updates…). After having re-installed all the desired languages in VS2010Express and afterwards the WinSDK_7p1 as well, I was happy to find VCpp2010Exp working as expected.

VCpp2010Exp provides a project template meant for creating DLLs: select file → new → project → Win32 → Win32-project, enter a name, and click ok. On the following welcome screen, do not click the finish but the continue button. A properties window will then be displayed in which you select the DLL radio button, and check the box labelled export symbols as well. Clicking the finish button will make VCpp2010Exp create a DLL project template that contains two files: the source code in *.cpp, and the header information in *.h . You’ll need to edit both of them later. Now the target platform has to be switched from Win32 to x64 in two steps. Firstly, select project → properties, change the content of the configuration button to all configurations, then select configuration properties → general in the list displayed on the left, change the entry for platform toolset in the general section of the list displayed on the right from v100 (or v90 or …) to Windows7.1SDK, and click the apply button. Secondly, after VCpp2010Exp has finished applying this change, click the configuration manger button and set the x64 platform (it does not matter whether the debug or release configuration is active) by clicking on platform, selecting new, changing the entry of new platform to x64, clicking the ok button, and finally clicking the close button of the configuration manager.

I decided to start with a very simple x64 DLL that only provides two different functions, i.e., a function of type integer (in VBA7, this corresponds to a function of type long), and another function of type void (corresponding to a sub in VBA7). Both functions calculate the product of two numbers provided as input variables; the integer function returns the result as its own value, whereas the void function returns it as value of a third variable, which consequently must be transferred by reference (in the declaration, an ampersand has to be noted directly in front of the variable’s name in Cpp, and the ByRef keyword in VBA7).

The C++ file named DLL_Test.cpp thus reads:

#include "stdafx.h"
#include "DLL_Test.h"

DLL_TEST_API int DLL_Mult_FNC(int a, int b)
{
return a*b;
}

DLL_TEST_API void DLL_Mult_SUB(int a, int b, int &c)
{
c = a*b;
return;
}

The meaning of DLL_TEST_API is defined in the header file DLL_Test.h together with some other very important settings:

#ifdef DLL_TEST_EXPORTS
#define DLL_TEST_API __declspec(dllexport)
#else
#define DLL_TEST_API __declspec(dllimport)
#endif

extern "C"
{
DLL_TEST_API int DLL_Mult_FNC(int a, int b);

DLL_TEST_API void DLL_Mult_SUB(int a, int b, int &c);
}

The block of macro-defining statements at the beginning has been added by VCpp2010Exp automatically, and saves a bit of typing if many functions are to be exported. The magic extern "C" block around the function declarations has been added manually after finding out that otherwise the names of the exported functions are so badly mangled that they cannot be called from VBA7. Imho, the tool Dependency Walker is quite valuable in such a situation. The hint that name mangling could be the reason for my DLL functions not yielding any result on the first try was mentioned in a post on Stack overflow, in which user arsane posted a link to a very comprehensive survey of calling conventions.

In VBA7, the functions need to be declared in a standard module; the name of a function in VBA7, which directly follows the function or sub keywords, has to be different from the name exported by the (C++) DLL, which has to be given directly after the Alias keyword.

Option Explicit

Declare PtrSafe Sub SUB_AmultB Lib "insert-correct-path-to-dll-here\DLL_Test.dll" Alias "DLL_Mult_SUB" (ByVal ia&, ByVal ib&, ByRef ic&)
Declare PtrSafe Function FNC_AmultB& Lib "insert-correct-path-to-dll-here\DLL_Test.dll" Alias "DLL_Mult_FNC" (ByVal ia&, ByVal ib&)

Function WS_DLL_Mult_SUB&(ByVal ia&, ByVal ib&)
Dim ic&
Call SUB_AmultB(ia, ib, ic)
WS_DLL_Mult_SUB = ic
End Function

Function WS_DLL_Mult_FNC&(ByVal ia&, ByVal ib&)
WS_DLL_Mult_FNC = FNC_AmultB(ia, ib)
End Function

Do not get confused by the completely different meanings of ampersands in C++ and VBA7, resp.: in C++, a leading ampersand denotes a call by reference, whereas in VBA7, a trailing ampersand defines the variable to be of type long.

In case you find this recipe is not working for you, please let me know…