# Computing a conservatively rounded square of a double on broken compilers

I am trying to compute $x^2$ where $x$ is a double precision floating point number. I need conservative rounding, which means I need $x \cdot x$ rounded up and $x \cdot x$ rounded down. If the rounding mode is set to FE_UPWARDS, the code is

double lo = -(x * -x);
double hi = x * x;


However, this does not work for me on clang 3.4: the compiler notices the "common subexpression" and does only one multiply. Does anyone know how to convince clang to not do this, ideally without harming other optimizations? In particular, I would like to avoid noninline functions.

• "I need conservative rounding, which means I need $x \cdot x$ rounded up and $x \cdot x$ rounded down." I'm a little confused; what's the difference between those two expressions? – Geoff Oxberry Mar 26 '14 at 19:59
• If the floating point number $x$ is considered as a real number, $x \cdot x$ is a real number that usually cannot be represented as a float. I need float point numbers on either side. – Geoffrey Irving Mar 26 '14 at 20:04
• I believe GCC has a flag that forces the compiler to do IEEE conforming floating point arithmetic. This would presumably eliminate the "optimization" of the double negation. – Wolfgang Bangerth Mar 29 '14 at 1:14
• Yes, gcc falls into the class of "nonbroken compilers" for purposes of this question. Clang lacks -frounding-math. – Geoffrey Irving Mar 29 '14 at 1:20
• llvm (and thus presumably clang) has an enable-unsafe-fp-math option - what does clang -mno-enable-unsafe-fp-path do? – mabraham Apr 7 '14 at 22:36

I believe this can be fixed by packing the interval bounds into SSE registers and performing all the interval operations with SSE as described in

http://hal.inria.fr/docs/00/28/84/56/PDF/intervals-sse2-long-paper.pdf

This should be faster than my current code anyways, and clang shouldn't apply SSE within the same instruction.

If your compiler has nextafter, then you can do something like this, which works regardless of which rounding mode is set:

double x2 = x*x;
double lo = nextafter(x2,-1);
double hi = nextafter(x2,x2+1);


but I don't know whether nextafter is inlined (probably not).

• Is there a reason for the $x2+1$ argument to nextafter, instead of just $+1$? – Geoffrey Irving Apr 15 '14 at 1:25
• @GeoffreyIrving, nextafter(x,y) returns the next representable floating-point value after x in the direction of y. – lhf Apr 15 '14 at 10:53
• If you want to go a specific direction, your best bet is to just pass + or - infinity to nextafter. If x2 is so big that x2+1 == x2, you'll just get x2 back from nextafter. – Peter Cordes Feb 22 '16 at 22:56

This is a sad answer, so I am not going to accept it in hopes that someone else has a better solution:

double hi = x * x;
double lo = (2*epsilon-1)*hi;


This is safe because if we're rounding up, \begin{aligned} -\left((2\epsilon-1) \odot (x \odot x)\right) &\le (1 + \epsilon)^2(1-2\epsilon) x^2 \\ &\le (1+2\epsilon+\epsilon^2)(1-2\epsilon)x^2 \\ &\le (1-\epsilon^2-\epsilon^3)x^2 \\ &\le x^2 \end{aligned} so our code correctly computes the lower bound.