7

Is there a value of type double (IEEE 64-bit float / binary64), K, such that K * K == 3.0? (The irrational number is of course "square root of 3")

I tried:

static constexpr double Sqrt3 = 1.732050807568877293527446341505872366942805253810380628055806;
static_assert(Sqrt3 * Sqrt3 == 3.0);

but the static assert fails.

(I'm guessing neither the next higher nor next lower floating-point representable number square to 3.0 after rounding? Or is the parser of the floating point literal being stupid? Or is it doable in IEEE standard but fast math optimizations are messing it up?)

I think the digits are right:

$ python

>>> N = 1732050807568877293527446341505872366942805253810380628055806
>>> N * N
2999999999999999999999999999999999999999999999999999999999996\
607078976886330406910974461358291614910225958586655450309636

Update

I've discovered that:

static_assert(Sqrt3 * Sqrt3 < 3.0); // pass
static_assert(Sqrt3 * Sqrt3 > 2.999999999999999); // pass
static_assert(Sqrt3 * Sqrt3 > 2.9999999999999999); // fail

So the literal must produce the next lower value.

I guess I need to check the next higher value. Could bit-dump the representation maybe and then increment the last bit of the mantissa.

Update 2

For posterity: I wound up going with this for the Sqrt3 constant and the test:

static constexpr double Sqrt3 = 1.7320508075688772;
static_assert(0x1.BB67AE8584CAAP+0 == 1.7320508075688772);
static_assert(Sqrt3 * Sqrt3 == 2.9999999999999996);
Andrew Tomazos
  • 66,139
  • 40
  • 186
  • 319
  • 3
    The compiler will not use as many digits as you have given, it will round to a valid double number. So the check you did in Python doesn't tell you anything. – Mark Ransom Jun 06 '22 at 00:25
  • @MarkRansom: Yeah I know that, I'm just checking/showing that the literal is entered correctly. – Andrew Tomazos Jun 06 '22 at 00:26
  • 1
    To give a mathematical answer: Every floating-point number is scaled by a power of two, and the mantissa is in binary, hence even something like sqrt(3) must be expressed as a fraction where the denominator is a power of two. But the square root of a non-square number can never be represented as a fraction, thus what you get via sqrt(3) is not "the" square root of 3, but the closest approximation to that value (some might say it corresponds to a *range* of values near that approximation). A similar example is sin(pi), which is *not* 0, as neither pi nor sin(x) are accurate enough. – IS4 Jun 06 '22 at 00:27
  • @IS4: Yes, so it's the first case? ie *neither the next higher nor next lower floating-point representable number square to exactly 3.0* ? – Andrew Tomazos Jun 06 '22 at 00:28
  • @IS4: I don't think that answers the question, though. It's true that `sqrt(3)` won't be exactly √3, but that doesn't tell us whether `sqrt(3) * sqrt(3)` will be exactly 3 (since the multiplication will *also* involve rounding). – ruakh Jun 06 '22 at 00:29
  • @IS4. You're missing a step. When you square, you lose the lowest bit, so it's still possible to get 3 from your argument. You have to show that the error flows into the lowest bit after squaring – Mad Physicist Jun 06 '22 at 00:31
  • 1
    Indeed. As real numbers, it cannot represent anything that would produce 3 when squared since it must be rational, but as @ruakh points out (and what I was going to add), rounding *could* occur during the multiplication, and so the question is why it doesn't happen here. That's why this remains a comment and not an answer. – IS4 Jun 06 '22 at 00:32
  • @IS4: To be clear, when I say square to exactly 3.0, I mean that after rounding the result of the multiplication into a floating point 64-bit double, the real number that corresponds to that double is 3 (which is exactly representable in a double). – Andrew Tomazos Jun 06 '22 at 00:35
  • If I was going to guess, the reason it is not rounded is that *both* operands to multiplication have some error, which also gets amplified in the multiplication. The error in the result is no longer smaller than what could be rounded off. In contrast, something like 0.3333... * 3 should be equal to 1 as the error is only in the first operand. – IS4 Jun 06 '22 at 00:38
  • Does this answer your question? [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) – duffymo Jun 06 '22 at 00:43
  • @duffymo: Thanks but no, I know how floating point math works. – Andrew Tomazos Jun 06 '22 at 00:47
  • Square root of 3 is an irrational number. "feasible"? No. Your code tells you so. You seem to suggest that the IEEE-754 standard is either wrong or implemented incorrectly in Python. – duffymo Jun 06 '22 at 01:19
  • I know what you're asking. – duffymo Jun 06 '22 at 01:32
  • 1
    FWIW, the smallest non-square numbers for which `sqrt(n) * sqrt(n) == n` as an IEEE-754 `double` under the default rounding mode seem to be 11, 14, and 17. – Steve Summit Jun 06 '22 at 15:47
  • Perhaps worth noting that assuming IEEE 754 binary floating-point, round-ties-to-even, and a correctly rounded `sqrt` function, it's sufficient to test whether `sqrt(n) * sqrt(n) == n`: if the equality fails, then no other `x` can possibly square to `n`. This follows from the fact that (under the same conditions) `sqrt(x*x) == x` for any non-huge non-tiny positive float `x`. – Mark Dickinson Jun 06 '22 at 15:50

6 Answers6

11

The answer is no; there is no such K.

The closest binary64 value to the actual square root of 3 is equal to 7800463371553962 × 2-52. Its square is:

60847228810955004221158677897444 × 2-104

This value is not exactly representable. It falls between (3 - 2-51) and 3, which are respectively equal to

60847228810955002264642499117056 × 2-104

and

60847228810955011271841753858048 × 2-104

As you can see, K * K is much closer to 3 - 2-51 than it is to 3. So IEEE 754 requires the result of the operation K * K to yield 3 - 2-51, not 3. (The compiler might convert K to an extended-precision format for the calculation, but the result will still be 3 - 2-51 after conversion back to binary64.)

Furthermore, if we go to the next representable value after K in the binary64 format, we will find that its square is closest to 3 + 2-51, which is the next representable value after 3.

This result should not be too surprising; in general, incrementing a number by 1 ulp will increment its square by roughly 2 ulps, so you have about a 50% chance, given some value x, that there is a K with the same precision as x such that K * K == x.

Brian Bi
  • 111,498
  • 10
  • 176
  • 312
  • Very interesting, so there was a 50% chance assuming it was a uniform distribution. I find myself wondering if about half of irrational roots do hit their square exactly, but I guess once they get longer the square will take up more space in the mantissa, leading to even less chance you'll hit the interval. – Andrew Tomazos Jun 06 '22 at 01:02
  • 2
    Another way to see that there is a 50% chance that `sqrt(x)*sqrt(x)==x` (for positive normal numbers) is that `sqrt` maps the interval `[1,4]` to the interval `[1, 2]` which only has half the numbers, so squaring the result must miss the original `x` in half the cases. – chtz Jun 06 '22 at 12:00
8

The C standard does not dictate the default rounding mode. While it is typically round-to-nearest, ties-to-even, it could be round-upward, and some implementations support changing the mode. In such case, squaring 1.732050807568877193176604123436845839023590087890625 while rounding upward produces exactly 3.

#include <fenv.h>
#include <math.h>
#include <stdio.h>

#pragma STDC FENV_ACCESS ON


int main(void)
{
    volatile double x = 1.732050807568877193176604123436845839023590087890625;
    fesetround(FE_UPWARD);
    printf("%.99g\n", x*x);  // Prints “3”.
}

x is declared volatile to prevent the compiler from computing x*x at compile-time with a different rounding mode. Some compilers do not support #pragma STDC FENV_ACCESS but may support fesetround once the #pragma line is removed.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • You wouldn't want to change the rounding for your whole program just to change one inconsequential result. – Mark Ransom Jun 06 '22 at 01:16
  • 1
    Very cool. If we refer to Brian Bis answer then `FE_UPWARD` rounds upward to the further away `60847228810955011271841753858048 × 2-104` (3) rather than rounding down to the nearer `60847228810955002264642499117056 × 2^-104` (3-2^51). (But of course Mark is right, we wouldn't change global rounding mode as a practical solution - I think you're just pointing this out for fun) – Andrew Tomazos Jun 06 '22 at 01:16
  • 1
    For completeness: There are *two* solutions, as I noted in comment on the question: `K=0x1.bb67ae8584caap+0`, multiply rounded toward pos. inf, and `K=0x1.bb67ae8584cabp+0`, multiply rounded towards neg. inf, or towards zero. – njuffa Jun 06 '22 at 01:25
  • @njuffa: Yes, there are the "next lower" and "next higher" values I was refering to in the question (in hexadecimal floating-point literal format no less). Thanks! – Andrew Tomazos Jun 06 '22 at 01:28
  • So it is possible - nice. – chux - Reinstate Monica Jun 06 '22 at 03:44
3

Testing with Python is valid I think, since both use the IEEE-754 representation for doubles along with the rules for operations on same.

The closest possible double to the square root of 3 is slightly low.

>>> Sqrt3 = 3**0.5
>>> Sqrt3*Sqrt3
2.9999999999999996

The next available value is too high.

>>> import numpy as np
>>> Sqrt3p = np.nextafter(Sqrt3,999)
>>> Sqrt3p*Sqrt3p
3.0000000000000004

If you could split the difference, you'd have it.

>>> Sqrt3*Sqrt3p
3.0
Mark Ransom
  • 299,747
  • 42
  • 398
  • 622
  • `numpy.nextafter` for the win! Thanks! https://numpy.org/doc/stable/reference/generated/numpy.nextafter.html – Andrew Tomazos Jun 06 '22 at 00:48
  • @AndrewTomazos I thought they put `nextafter` into Python itself somewhere, but I couldn't find it in a hurry. – Mark Ransom Jun 06 '22 at 00:49
  • Might be interesting to see how its implemented and what underlying C functions it uses. There is a C header that has IEEE float manipulation routines / intrinsics. Maybe it uses them. – Andrew Tomazos Jun 06 '22 at 00:51
  • 1
    @MarkRansom [`math.nextafter`](https://docs.python.org/3/library/math.html#math.nextafter) (In the relatively recent Python 3.9) – Artyer Jun 06 '22 at 01:24
  • @Artyer thanks. That explains why I couldn't find it, I'm stuck on Python 3.8. – Mark Ransom Jun 06 '22 at 01:26
  • "since both use the IEEE-754 representation for doubles " --> "both" seems to imply C specifies IEEE-754 format and functional details. C does not specify that, even though OP did tag question as such. – chux - Reinstate Monica Jun 06 '22 at 14:13
  • @chux-ReinstateMonica of course you're right, neither C or C++ specify floating point details to the best of my knowledge. It was only in the question itself where IEEE was specified; my wording was sloppy. – Mark Ransom Jun 06 '22 at 14:39
1

In the Ruby language, the Float class uses "the native architecture's double-precision floating point representation" and it has methods named prev_float and next_float that let you iterate through different possible floats using the smallest possible steps. Using this, I was able to do a simple test and see that there is no double (at least on x86_64 Linux) that meets your criterion. The Ruby interpreter is written in C, so I think my results should be applicable to the C double type.

Here is the Ruby code:

x = Math.sqrt(3)
4.times { x = x.prev_float }
9.times do
  puts "%.20f squared is %.20f" % [x, x * x]
  puts "Success!" if x * x == 3
  x = x.next_float
end

And the output:

1.73205080756887630500 squared is 2.99999999999999644729
1.73205080756887652704 squared is 2.99999999999999733546
1.73205080756887674909 squared is 2.99999999999999822364
1.73205080756887697113 squared is 2.99999999999999866773
1.73205080756887719318 squared is 2.99999999999999955591
1.73205080756887741522 squared is 3.00000000000000044409
1.73205080756887763727 squared is 3.00000000000000133227
1.73205080756887785931 squared is 3.00000000000000177636
1.73205080756887808136 squared is 3.00000000000000266454
David Grayson
  • 84,103
  • 24
  • 152
  • 189
0

Is there a value of type double, K, such that K * K == 3.0?

Yes.


K = sqrt(n); and K * K == n may be true, even when √n is irrational.

Note that K, the result of sqrt(n), as a double, is a rational number.


  • Various rounding modes: @Eric

  • K * K rounds to n

Example: Roots n: 11, 14 and 17 when squared are n.

for (int i = 10; i < 20; i++) {
  double x = sqrt(i);
  double y = x * x;
  printf("%2d %.25g\n", i, y);
}

10 10.00000000000000177635684
11 11
12 11.99999999999999822364316
13 12.99999999999999822364316
14 14
15 15.00000000000000177635684
16 16
17 17
18 17.99999999999999644728632
19 19.00000000000000355271368
  • Different precision

Rather than 53 bits with common double, say the FP math was done with 24. Roots n: 3, 5 and 10 when squared are n.

for (int i = 2; i < 11; i++) {
  float x = sqrtf(i);
  printf("%2d  %.25g\n", i, x*x);
}

 2  1.99999988079071044921875
 3  3
 4  4
 5  5
 6  6.000000476837158203125
 7  6.999999523162841796875
 8  7.999999523162841796875
 9  9
10  10

or say the FP math was done with 64 bits. Roots n: 5, 6 and 10 when squared are n.

for (int i = 2; i < 11; i++) {
  long double x = sqrtl(i);
  printf("%2d  %.35Lg\n", i, x*x);
}

 2  1.9999999999999999998915797827514496
 3  3.0000000000000000002168404344971009
 4  4
 5  5
 6  6
 7  6.9999999999999999995663191310057982
 8  7.9999999999999999995663191310057982
 9  9
10  10

With various precisions, (note C does not specify a fixed precision), K * K == 3.0 is possible.

  • FLT_EVAL_METHOD == 2

When FLT_EVAL_METHOD == 2, intermediate calculations may be done at higher precession, thus affecting the product of k*k.

(Have yet to come up with a good simple example.)

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
-1

sqrt(3) is irrational, which means that there is no rational number k such that k*k == 3. A double can only represent rational numbers; therefore, there is no double k such that k*k == 3.

If you can accept a number that is close to satisfying k*k == 3, then you can use std::numeric_limits (in <type_traits>, if memory serves) to see if you’re within some minimal interval around 3. It may look like:

assert( abs(k*k - 3.) <= abs(k*k + 3.) * std::numeric_limits<double>::epsilon * X);

Epsilon is the smallest difference from one that double can represent. We scale it by the sum of the two values to compare in order to bring its magnitude in line with the numbers we’re checking. X is a scaling factor that lets you adjust the precision you accept.

If this is a theoretical question: no. If it’s a practical question: yes, up some level of precision.

Seymore Glass
  • 53
  • 1
  • 4
  • √3 is irrational. `sqrt(3)` returns a `double` and all finite `double` are rational. – chux - Reinstate Monica Jun 06 '22 at 03:46
  • 1
    I think the OP understands that √3 is irrational and cannot be represented exactly using any finite representation. But that wasn't the question. Suppose that D(x) is an operator that means "The closest approximation you can get to x, using type `double`." So we can't represent √3 , but we can represent D(√3 ). And we may not be able to represent D(√3 ) × D(√3 ), either. But we can represent D(D(√3 ) × D(√3 )). So the question is, can we make D(D(√3 ) × D(√3 )) == 3? – Steve Summit Jun 06 '22 at 15:40