4

Suppose I have some code such as:

float a, b = ...; // both positive
int s1 = ceil(sqrt(a/b));
int s2 = ceil(sqrt(a/b)) + 0.1;

Is it ever possible that s1 != s2? My concern is when a/b is a perfect square. For example, perhaps a=100.0 and b=4.0, then the output of ceil should be 5.00000 but what if instead it is 4.99999?

Similar question: is there a chance that 100.0/4.0 evaluates to say 5.00001 and then ceil will round it up to 6.00000?

I'd prefer to do this in integer math but the sqrt kinda screws that plan.

EDIT: suggestions on how to better implement this would be appreciated too! The a and b values are integer values, so actual code is more like: ceil(sqrt(float(a)/b))

EDIT: Based on levis501's answer, I think I will do this:

float a, b = ...; // both positive
int s = sqrt(a/b);
while (s*s*b < a) ++s;

Thank you all!

Dan
  • 5,929
  • 6
  • 42
  • 52
  • 2
    What are these calculations for? You can, for example, change `sqrt(x) < y` into `x < y * y`, which keeps your integer operations (as long as `y * y` doesn't get too large). – GManNickG Nov 02 '11 at 17:31
  • 2
    Unfortunately, it's very possible that `s1 != s2`, even if `a/b` is a perfect square. Floating point arithmetic is notoriously inaccurate, so your chances of ever getting `5.00000` (or any number exactly what you want) are rather slim. – ssube Nov 02 '11 at 17:34
  • @peachykeen Even if we are dealing with integer numbers? `5.0` can be perfectly represented as a float, an so can `100.0` and `20.0`. – glglgl Nov 02 '11 at 17:39
  • 3
    @Dan If you have an implementation where `100/4` evaluates to something close to `5`, you definitely have a problem. – glglgl Nov 02 '11 at 17:40
  • @GMan I'll think about that... I have `a` items that I need to distribute into `b` boxes. When `a>b` I want to divide each box into 4, 9, 16, 25, etc. sub-boxes. So the `s` value is the scale (2, 3, 4, 5 respectively), or the number of division of the box along each dimension. For example, 300 units and 20 boxes requires each box be divided into 16 (4x4). – Dan Nov 02 '11 at 17:45
  • If 100.0/4.0 ever evaluates to anything close to 5.0 then you've got bigger problems to worry about... – Paul R Nov 02 '11 at 17:46
  • @glglgl It depends on the integer. 5.0 may be represented, but if those numbers ever change, you might get data-dependent bugs. – ssube Nov 02 '11 at 17:54

4 Answers4

6

I don't think it's possible. Regardless of the value of sqrt(a/b), what it produces is some value N that we use as:

int s1 = ceil(N);
int s2 = ceil(N) + 0.1;

Since ceil always produces an integer value (albeit represented as a double), we will always have some value X, for which the first produces X.0 and the second X.1. Conversion to int will always truncate that .1, so both will result in X.

It might seem like there would be an exception if X was so large that X.1 overflowed the range of double. I don't see where this could be possible though. Except close to 0 (where overflow isn't a concern) the square root of a number will always be smaller than the input number. Therefore, before ceil(N)+0.1 could overflow, the a/b being used as an input in sqrt(a/b) would have to have overflowed already.

Jerry Coffin
  • 476,176
  • 80
  • 629
  • 1,111
  • I agree that we should expect to trust `ceil`, but I can't find where (if anywhere) in the C standard it actually talks about the accuracy of the functions in `math.h`. Assuming 0.5ulp accuracy, `ceil` *has* to return an integer, and it would be pretty unreasonable for an implementation to fail to do that. But we know that, say `cos` isn't guaranteed to be 0.5ulp accurate (or `sqrt`, which I think which answers the second question), so why should `ceil` be? Other than because the obvious implementation is to mask off some bits of mantissa (adding 1 if those bits weren't already 0), I mean. – Steve Jessop Nov 02 '11 at 17:51
  • I guess it depends whether we read "compute the smallest integer value not less than x" to mean, "it is guaranteed that the value returned is an integer, equal to ..." or to mean, "the exact mathematical result is ..., but the value returned is subject to usual inaccuracy". – Steve Jessop Nov 02 '11 at 17:55
  • @SteveJessop: Realistically, about the only time you stand a chance of running into a problem is with a 64-bit int. With the typical situation of a 32-bit int and a double with a 53-bit mantissa, you're doing to overflow the int long before there's any difficulty converting the double to a precise integer value. – Jerry Coffin Nov 02 '11 at 18:19
  • As for why we'd expect a difference, I think it's pretty simple. Using your example (`sqrt`) all non-integer results are irrational, so they can't *possibly* be represented accurately. The only place there's room for any question at all about it producing a truly accurate answer is when we present it with a perfect square. – Jerry Coffin Nov 02 '11 at 18:51
  • Well, I expect `ceil` to be accurate simply because I can imagine a simple implementation which is. All that I'm agonising over is whether it's *guaranteed* to be accurate, and if so where (C99 and IEEE 754 being candidates). Simultaneously, the reason that I don't expect `sqrt` to necessarily be accurate within 0.5ulp is not (as is the case for you) simply because it can't be perfectly accurate. It's because I happen to know that it's difficult to implement it to 0.5ulp accuracy without prohibitive cost. Also nothing to do with the standard, which is why I'm not entirely happy there either... – Steve Jessop Nov 03 '11 at 00:28
  • `sqrt` is guaranteed to be correctly rounded on any platform that adheres to IEEE-754. As for `ceil`: "The ceil functions return x, expressed as a floating-point number." (§7.12.9.1). x is always exactly representable, and if it were not returned that would constitute a serious bug in a platform's math library. – Stephen Canon Nov 03 '11 at 10:07
  • (Note also that the cost of implementing `sqrt` with correct rounding (actually a stronger statement than "0.5 ulp") is not at all prohibitive. Most modern systems have a correctly-rounded square root operation in hardware. Even on those that don't, even a reference integer implementation generally only requires a couple hundred cycles at worst (often less). Note that (in particular) correctly-rounded square root is *easier* than correctly-rounded divide. – Stephen Canon Nov 03 '11 at 10:13
  • @Stephen: I know the standard *says* it returns the ceiling, but then again `cos` *says* that it returns the cosine. In the latter case, we understand this to mean, "at any rate some nearby value, not necessarily the nearest". In the former case, we understand it to mean, "exactly". All I ask is whether the C99 standard specifies that there's a difference, or if it's left to common-sense QOI and/or IEEE 754. And OK, `sqrt` may have been a bad example of a function permitted by C99 to be less accurate than 0.5ulp, it's easier than most. – Steve Jessop Nov 04 '11 at 12:54
  • QOI and IEEE-754 (clause 5.3.1 in the 2008 standard, clause 5.5 in the 1985 standard). – Stephen Canon Nov 04 '11 at 13:05
3

You may want to write an explicit function for your case. e.g.:

/* return the smallest positive integer whose square is at least x */
int isqrt(double x) {
  int y1 = ceil(sqrt(x));
  int y2 = y1 - 1;
  if ((y2 * y2) >= x) return y2;
  return y1;
}

This will handle the odd case where the square root of your ratio a/b is within the precision of double.

levis501
  • 4,117
  • 23
  • 25
  • I like this approach. It's provably correct, without having to know about or rely on floating point implementation details. – Dan Nov 02 '11 at 17:58
1

Equality of floating point numbers is indeed an issue, but IMHO not if we deal with integer numbers.

If you have the case of 100.0/4.0, it should perfectly evaluate to 25.0, as 25.0 is exactly representable as a float, as opposite to e.g. 25.1.

glglgl
  • 89,107
  • 13
  • 149
  • 217
-2

s1 will always equal s2.

The C and C++ standards do not say much about the accuracy of math routines. Taken literally, it is impossible for the standard to be implemented, since the C standard says sqrt(x) returns the square root of x, but the square root of two cannot be exactly represented in floating point.

Implementing routines with good performance that always return a correctly rounded result (in round-to-nearest mode, this means the result is the representable floating-point number that is nearest to the exact result, with ties resolved in favor of a low zero bit) is a difficult research problem. Good math libraries target accuracy less than 1 ULP (so one of the two nearest representable numbers is returned), perhaps something slightly more than .5 ULP. (An ULP is the Unit of Least Precision, the value of the low bit given a particular value in the exponent field.) Some math libraries may be significantly worse than this. You would have to ask your vendor or check the documentation for more information.

So sqrt may be slightly off. If the exact square root is an integer (within the range in which integers are exactly representable in floating-point) and the library guarantees errors are less than 1 ULP, then the result of sqrt must be exactly correct, because any result other than the exact result is at least 1 ULP away.

Similarly, if the library guarantees errors are less than 1 ULP, then ceil must return the exact result, again because the exact result is representable and any other result would be at least 1 ULP away. Additionally, the nature of ceil is such that I would expect any reasonable math library to always return an integer, even if the rest of the library were not high quality.

As for overflow cases, if ceil(x) were beyond the range where all integers are exactly representable, then ceil(x)+.1 is closer to ceil(x) than it is to any other representable number, so the rounded result of adding .1 to ceil(x) should be ceil(x) in any system implementing the floating-point standard (IEEE 754). That is provided you are in the default rounding mode, which is round-to-nearest. It is possible to change the rounding mode to something like round-toward-infinity, which could cause ceil(x)+.1 to be an integer higher than ceil(x).

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312