Archive for the Uncategorized Category

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.

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.

Mind the gap when using IEEE754 DP numbers

Posted in binary, Computer, Mathematics, Microsoft Excel, Number processing, Number systems, Programming Languages, Software, Spreadsheet, Uncategorized, VBA on September 11, 2011 by giorgiomcmlx

The gap between two adjacent IEEE754 DP numbers is more or less proportional to the magnitude of these numbers. The overwhelming majority of DP numbers do not exactly equal a power of 2. In this case, the gaps between such a DP number and its direct neighbours on both the smaller as well as the bigger sides have a magnitude of about 1.11E-16 times the smallest power of 2 that is larger than the number. For example, consider the DP number 1.5 (hex 3F F8 00 00 00 00 00 00). The smallest power of 2 which is bigger than 1.5 is 2^1=2. Thus, the biggest DP number being smaller than 1.5 is about 1.5 – 2⋅1.11E-16 = 1,499999999999999778 (exactly hex 3F F7 FF FF FF FF FF FF), and the smallest DP number being bigger than 1.5 is about 1.5 + 2⋅1.11E-16 = 1,500000000000000222 (exactly hex 3F F8 00 00 00 00 00 01).

At all exact powers of 2 the gap between adjacent DP numbers doubles in size. This means the question “does x equal 1.0” will yield the answer “true” for any real x in the range 1.0-5.55E-17 <= x < 1.0+1.11E-16, if rounding to the nearest is used in order to convert an arbitrary real number x into an IEEE754 DP number. Please note that in this special case, the interval limits are located asymmetrically around 1.0! That happens for all exact powers of 2. The size of the gap between two adjacent DP numbers in the interval from 1 to 2, i.e. 2.22E-16, is called machine epsilon.

The following figures show the discretisation effect occurring when using DP numbers. In particular, the region around 1 is explored using cell formulae in the spreadsheet application Microsoft®™ Excel®™ 2002 as well as using the VBA 6®™ built-in programming language. The FOSS add-in xNumbers is used to show how the results look like if the calculations are carried out with considerably higher numerical accuracy, i.e., 30 instead of the 17 decimal places in case of the IEEE754 DP numbers. The vertical scaling had to be shifted downwards by 1 in the figures, because XL cannot deal with axis limits of 1+/-2E-15, for example. The discretisation effect occurs in the first step of all shown calculations, i.e., when the relatively small quantity x is added to 1, and is thus not influenced by the final shift process.

machine epsilon in MS Excel and VBA

There is a remarkable difference between the results obtained using XL directly on the one hand and those obtained by invoking the programming language VBA 6 from XL on the other hand. MS seems to have implemented some “politically correct mathematics” into XL here in order to maximise the interval for which the answer to “does x equal 1” yields true… The results of the XL cell formulae are clearly not IEEE754 compliant.

details of machine epsilon in MS Excel and VBA

DP numbers are most dense around 0. Neglecting the special concept of “denormalised numbers”, the absolute value of the smallest numbers different from 0 is 2^(-1022), which is about 2.225E-308. The absolute value of the biggest DP number is 2^1023+2^1022+…+2^972+2^971 (hex
7F EF FF FF FF FF FF FF), which nearly equals 1.7976931348623157E+308. The adjacent DP number on its lower side, i.e. 1.7976931348623155E+308, is smaller by an amount of 2^971, which is about 2.00E+292. In the biggest dyad of the DP number system, the gap of 2.00E+292 between two adjacent DP numbers is huge in absolute terms, but still small in relative terms: 2.00E+292/1,79769E+308 = 1.11E-16, just as it was to be expected.

I have “translated” a routine from FORTRAN into a VBA array function that gives you all relevant information about the numerical characteristics of your computer. You may download the VBA6 code and use it if you follow the GPL v3 (included).

The size of the DP gap is usually small even in absolute terms for all numbers humans deal with in their every day or business life. However, it should be kept in mind that the DP gap around a number is roughly proportional to its magnitude. For example, in the range between 2^52 and 2^53, which equals about 4.5036E+15 to 9.0072E+15, only the integers (gap is 2^0=1) are exactly representable by DP numbers. In the next bigger dyad from 9.0072E+15 to 1.80144E+16, only the even integers (gap is 2^1=2) are exactly representable by DP numbers, and so on. Of course, this is only a consequence of limiting the mantissa of IEEE DP numbers 53 bits or roughly 17 significant decimal digits.

Even in science and engineering, numbers are very rarely known with an accuracy of more than 10 significant decimal digits. In fact, even most scientific measurements typically yield an accuracy of less than 5 significant decimal digits. This is why usually virtually nobody is bothered by neither the discreteness of the set of DP numbers in the continuum of the real numbers nor their limited accuracy of about 17 decimal digits. Consequently, people are quite surprised if even very simple calculations fail due to an effect called catastrophic cancellation. This will be the subject of a follow-up post.

What is this???

Posted in Uncategorized on August 6, 2011 by giorgiomcmlx

I hope “this” will develop into a lively virtual conversation about some aspects of how to apply the laws of mathematics and physics to solve problems using computers. In order to address many people, I will try to use open source software in the examples wherever it is possible. I’m interested in the history of science and the ideas of people having pushed it as well.

Because most of my examples will be inspired by geodesy, geophysics or astronomy, resp., I decided to summarise these subjects by the name “Carolo” in the title of this blog. This is meant as a tribute to the great former scientist Johann Carl Friedrich Gauß. In fact, “Carolo” appears as Latin version of his second German Christian name on the title page of his “Disquisitiones arithmeticae”.

Many Europeans and probably wine connoisseurs all around the world will know that “Barolo” is a famous red wine produced in Italy. According to the region of Italy where most of my practical examples will originate from, I should have chosen “Montepulciano” instead, but that does not rhyme on “Carolo”.