Associativity of floating-point multiplication in a special case

Question

It's well known that IEEE floating-point multiplication is not associative. However, consider the special case where a, b, and c are 32-bit signed integers (in C):

double da = (double) a;
double db = (double) b;
double dc = (double) c;

Now, does da*(db*dc) == (da*db)*dc always return 1? Here, doubles are 64-bit double-precision IEEE floating point numbers. Round-to-even is being used.

I tried several examples by hand, such as a = pow(2, 30) + 1, b = pow(2, 30)+1, c = pow(2, 30) + pow(2, 6), numbers large enough to ensure there would be some rounding involved, because the exact mathematical answer has no representation.

Unfortunately, I couldn't find a counterexample. This expression looks like it may return 0. Could it possibly always return 1? Why?

Answer 1

I managed to find a counterexample. Consider a = pow(2,30) + 1, b = pow(2,30) + 1, c = pow(2,30) + pow(2, 6) + 1.

Then if one computes this by hand using round-to-even, we get:

(da*db)*dc == pow(2, 90) + pow(2, 66) + pow(2, 61) + pow(2, 60) + pow(2, 38)

On the other hand:

da*(db*dc) == pow(2, 90) + pow(2, 66) + pow(2, 61) + pow(2, 60)

Notice that the binary representations of these two results differ only in the very last bit. Indeed, 90-38=52 , the least significant fractional bit in the binary representation of a 64-bit double precision floating point.

Answer 2

For questions like this, my first impulse is to try a brute-force search with random numbers. In this case, generating random 32-bit integers with uniform distribution demonstrates that da*(db*dc) != (da*db)*dc for about one third of the cases.

In comments, Pascal Cuoq made a convincing argument that the largest operand must be > 2²⁶ in magnitude for inequality to occur. This means we can further restrict a search to operands < 2²⁷ in magnitude:

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static uint32_t z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

int main (void)
{
    long long count = 0LL, fail = 0LL;
    double minmax = 1.0 / 0.0;
    do {
        int32_t a = KISS & (((uint32_t)1 << 27) - 1);
        int32_t b = KISS & (((uint32_t)1 << 27) - 1);
        int32_t c = KISS & (((uint32_t)1 << 27) - 1);
        volatile double da = (double) a;
        volatile double db = (double) b;
        volatile double dc = (double) c;
        volatile double dab = da * db;
        volatile double dbc = db * dc;
        volatile double p1 = dab * dc;
        volatile double p2 = da * dbc;

        count++;
        if (p1 != p2) {
            double max = fmax (fabs(c), fmax (fabs(da), fabs(db)));
            fail++;
            if (max < minmax) {
                minmax = max;
                printf ("a=%08x  b=%08x  c=%08x  p1=% 22.13a  p2=% 22.13a  minmax=%08x  count=%13lld  fail=%13lld\n",
                        a, b, c, p1, p2, (int32_t)minmax, count, fail);
            }
        }
    } while (1);
    return EXIT_SUCCESS;
}

The output of the above code strongly suggests that inequality can occur if the magnitude of either the product da*db or the product db*dc exceeds 2⁵³, and that the largest factor therefore must exceed √2⁵³ in magnitude, that is, it must be greater than 0x05A82799.

a=0658c99a  b=066c9dd5  c=0679e799  p1= 0x1.080f1700c8935p+80  p2= 0x1.080f1700c8934p+80  minmax=0679e799  count=           48  fail=            1
a=05b06599  b=062ad373  c=0303ee05  p1= 0x1.a72fc3446c46ep+78  p2= 0x1.a72fc3446c46dp+78  minmax=062ad373  count=          248  fail=            6
a=05b3b91a  b=05c6607f  c=061761fd  p1= 0x1.9129296aa189ap+79  p2= 0x1.9129296aa189bp+79  minmax=061761fd  count=          552  fail=           14
a=05ffbe73  b=055f71db  c=0413771e  p1= 0x1.06c12351727fbp+79  p2= 0x1.06c12351727fcp+79  minmax=05ffbe73  count=         1029  fail=           33
a=05ae5163  b=05e25901  c=017bd0b4  p1= 0x1.8cc25e308318ap+77  p2= 0x1.8cc25e3083189p+77  minmax=05e25901  count=         1160  fail=           41
a=0349eff2  b=05a40dd7  c=05b2b3c7  p1= 0x1.a6d583d48b807p+78  p2= 0x1.a6d583d48b806p+78  minmax=05b2b3c7  count=         3187  fail=          129
a=05ae99c5  b=05a4e5a1  c=03d88fe6  p1= 0x1.ed5c24a47c416p+78  p2= 0x1.ed5c24a47c417p+78  minmax=05ae99c5  count=        43479  fail=         1722
a=05626e68  b=05ad8275  c=05acdbe3  p1= 0x1.5b017380396c4p+79  p2= 0x1.5b017380396c5p+79  minmax=05ad8275  count=        49985  fail=         1952
a=01f7f82a  b=05a9f16b  c=05ad0265  p1= 0x1.fa48865afc8abp+77  p2= 0x1.fa48865afc8acp+77  minmax=05ad0265  count=       362414  fail=        14190
a=016047d9  b=05a5d671  c=05ab6f59  p1= 0x1.608381a720becp+77  p2= 0x1.608381a720bebp+77  minmax=05ab6f59  count=      1701498  fail=        66487
a=059ecf08  b=05a83be9  c=05aa9cf9  p1= 0x1.685523897b927p+79  p2= 0x1.685523897b926p+79  minmax=05aa9cf9  count=      3412539  fail=       132973
a=05a9203d  b=05a8b04b  c=050c6609  p1= 0x1.436f803e2c24bp+79  p2= 0x1.436f803e2c24ap+79  minmax=05a9203d  count=      4074531  fail=       158854
a=00e22912  b=05a8f039  c=05a846ad  p1= 0x1.c49a83e32465fp+76  p2= 0x1.c49a83e32465ep+76  minmax=05a8f039  count=     17021780  fail=       665646
a=03bc399a  b=05a8bdf9  c=05a7ca6d  p1= 0x1.de2fac272a096p+78  p2= 0x1.de2fac272a095p+78  minmax=05a8bdf9  count=     17716380  fail=       692903
a=054d26c1  b=05a84b09  c=05a8a8c7  p1= 0x1.53704223556b3p+79  p2= 0x1.53704223556b4p+79  minmax=05a8a8c7  count=     45074477  fail=      1762937
a=03afcaa4  b=05a88a8d  c=05a7f119  p1= 0x1.d7f3ccaacca2dp+78  p2= 0x1.d7f3ccaacca2cp+78  minmax=05a88a8d  count=    118913730  fail=      4649640
a=024efebb  b=05a7e34d  c=05a881b5  p1= 0x1.2783cfcc9de75p+78  p2= 0x1.2783cfcc9de74p+78  minmax=05a881b5  count=    254043151  fail=      9930225
a=05a7d8e7  b=05a87f8b  c=050f538e  p1= 0x1.43d6f3ad2f956p+79  p2= 0x1.43d6f3ad2f957p+79  minmax=05a87f8b  count=    401539580  fail=     15695528
a=007790b6  b=05a859c9  c=05a852bb  p1= 0x1.de61a9e2da63ap+75  p2= 0x1.de61a9e2da63bp+75  minmax=05a859c9  count=    455973582  fail=     17824719
a=01505b17  b=05a8270d  c=05a830d1  p1= 0x1.505d1a5677dfbp+77  p2= 0x1.505d1a5677dfap+77  minmax=05a830d1  count=   1109123125  fail=     43353113
a=0161a988  b=05a82bed  c=05a82917  p1= 0x1.61aaf393e420bp+77  p2= 0x1.61aaf393e420cp+77  minmax=05a82bed  count=   6209738039  fail=    242736587
a=05a82a7d  b=05a8267b  c=01c962b9  p1= 0x1.c96347ff18b85p+77  p2= 0x1.c96347ff18b84p+77  minmax=05a82a7d  count=  43616243136  fail=   1704950604
a=05a8293b  b=05a828a3  c=01a44ac0  p1= 0x1.a44b86270a7d5p+77  p2= 0x1.a44b86270a7d6p+77  minmax=05a8293b  count=  51943039028  fail=   2030420245
a=05a82769  b=05a8292d  c=048dbc75  p1= 0x1.236f64a4577fap+79  p2= 0x1.236f64a4577fbp+79  minmax=05a8292d  count= 476560049469  fail=  18628158958
a=05a82871  b=05a828fd  c=01c044b1  p1= 0x1.c04561ac59666p+77  p2= 0x1.c04561ac59667p+77  minmax=05a828fd  count= 846496019034  fail=  33088575026
a=003f4e0f  b=05a828d7  c=05a828f7  p1= 0x1.fa71612c258eap+74  p2= 0x1.fa71612c258e9p+74  minmax=05a828f7  count= 864908706983  fail=  33808303937
a=043c7359  b=05a827db  c=05a828d9  p1= 0x1.0f1d1e47c37e2p+79  p2= 0x1.0f1d1e47c37e3p+79  minmax=05a828d9  count=1147388053277  fail=  44850164685
a=05a82895  b=05a82887  c=02e5fc6a  p1= 0x1.72feb235ca1e7p+78  p2= 0x1.72feb235ca1e6p+78  minmax=05a82895  count=1383932059721  fail=  54096517230
a=005b76e8  b=05a8284f  c=05a82793  p1= 0x1.6ddbcc2614301p+75  p2= 0x1.6ddbcc2614300p+75  minmax=05a8284f  count=2212385094246  fail=  86479620757
a=034f873b  b=05a827f1  c=05a82755  p1= 0x1.a7c3a2fcb346bp+78  p2= 0x1.a7c3a2fcb346ap+78  minmax=05a827f1  count=2360223870382  fail=  92258517174
a=02b2ee37  b=05a827e7  c=05a82789  p1= 0x1.597729fe1b0f2p+78  p2= 0x1.597729fe1b0f3p+78  minmax=05a827e7  count=4588948516457  fail= 179376767823
a=0023898a  b=05a827db  c=05a82779  p1= 0x1.1c4c566dfffcbp+74  p2= 0x1.1c4c566dfffccp+74  minmax=05a827db  count=6159639135691  fail= 240773437792
a=05a827a1  b=05a827bd  c=0066ef2d  p1= 0x1.9bbcc027b0c20p+75  p2= 0x1.9bbcc027b0c21p+75  minmax=05a827bd  count=7528083873615  fail= 294264463804

Answer 3

does da*(db*dc) == (da*db)*dc always return 1?

Inspired by user394438 who used pairs of powers-of-2 for a counterexample, I tried the below to find a small triple a,b,c counterexample. This meets the requirements of Pascal Cuoq's comment: at least 1 of a,b,c needs to exceed 2²⁶. I am confident that a tighter condition is needed: at least 1 of a,b,c needs to exceed 2²⁷ to exercise the typical 53 binary digit double. In this case 0x8000001 = 2²⁷ + 1.

3, 0x4000001, 0x8000001

A key coding issue is to make certain intermediate products do not use long double. C allows (a*b)*c to perform as long double even with double a,b,c. Check FLT_EVAL_METHOD or use volatile. Else with long double in play with its typical 64+ bit precision, finding a counterexample may be moot.

void test(volatile double a, volatile double b, volatile double c) {
  volatile double bc = b * c;
  volatile double ab = a * b;
  volatile double p1 = ab * c;
  volatile double p2 = a * bc;

  // da*(db*dc) == (da*db)*dc ?
  if (p1 != p2) {
    static double sum_small = 1.0 / 0.0;
    double sum = a + b + c;
    if (sum < sum_small) {
      sum_small = sum;
      printf("%a\n%a\n", p1, p2);
      printf("%x %x %x\n", (unsigned) a, (unsigned) b, (unsigned) c);
      fflush(stdout);
    }
  }
}

double r() {
  int p2a = rand()%31;
  int p2b = rand()%31;
  return pow(2, p2a) + pow(2, p2b);
}

int main() {
  test(pow(2, 30) + 1, pow(2, 30) + 1, pow(2, 30) + pow(2, 6) + 1);
  for (unsigned long long i = 0; i < 1000ULL * 1000 * 1000 * 1000; i++) {
    test(r(), r(), r());
    //test(rand(), rand(), rand());
  }
  return 0;
}