## cos in XL14 or how not to do argument reduction

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.

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.