Computing desk
< April 20	<< Mar \| April \| May >>	April 22 >

Welcome to the Wikipedia Computing Reference Desk Archives
The page you are currently viewing is a transcluded archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

April 21

Dividing floating point numbers down to zero?

I understand that there are many different floating point formats and hardware. Assuming that one were to select some random real number and then repeated divide that number in half. Is there some guarantee that, on any given platform, the number would eventually resolve to zero? It seems to work fine on IEEE 754 types, I'm just not so sure about the others... Earl of Arundel (talk) 17:08, 21 April 2021 (UTC)[reply]

For fixed format floating point numbers, I would expect that once the minimum representable (closest to zero) value is reached, the next division would yield an arithmetic underflow and a zero result. However, some arbitrary-precision arithmetic implementations may not have such limitations and would allow you to continue to halve the value indefinitely, limited only by the memory available to represent the mantissa and exponent of each value. I don't know if such implementations actually exist, and haven't researched the available implementations. I did try one online arbitrary precision arithmetic page, [1], and was able to calculate 0.5^(10^18) = 6e-301029995663981196 - well beyond what you could achieve with fixed format floating point. At some point after that, this tool yields "Complete loss of accurate digits". Other implementations will likely differ. -- Tom N talk/contrib 02:37, 22 April 2021 (UTC)[reply]

Another possibility would be a fixed format, but with a division algorithm such that when the smallest representable non-zero number is divided by 2, it results in the same number. Then repeated halving would never reach zero. CodeTalker (talk) 04:00, 22 April 2021 (UTC)[reply]

So the short answer is: no, there is no such guarantee. There seem to be several possibilities. One is that indeed zero is always soon reached (possibly with a raised condition like underflow, but this causing an exception would in practice be unacceptably obnoxious). A second is that you can keep halving until the available space for the exponent is exhausted, which can take ages – the value 0.5^(10^18) mentioned above suggests that 60 bits were available for representing the exponent, not counting a sign bit, while the error message as well as the printed value show that mantissa space was sacrificed for the exponent – a more precise value being 6.113094441455699e-301029995663981196. Finally, one might reach a non-zero fixed point. The iterated transformation x ↦ (x+1)/(x+2), starting from x = 1, results mathematically in a strictly descending sequence, but in floating-point arithmetic you'll probably reach a fixed point in less than 20 iterations. --Lambiam 10:04, 22 April 2021 (UTC)[reply]

Ok so subdividing floating point numbers in such a way may not be the best approach after all. I appreciate the input everyone. Thanks! Earl of Arundel (talk) 18:27, 22 April 2021 (UTC)[reply]

So it sounds like you're trying to accomplish something. Just tell us what it is and what you want advice on. See XY problem. --47.155.96.47 (talk) 18:45, 22 April 2021 (UTC)[reply]

It sounds to me like he wants to learn about a certain property of different floating-point representations, and asked exactly what he meant to ask. --184.147.181.129 (talk) 20:50, 22 April 2021 (UTC)[reply]

However, the formulation that some choice "may not be the best approach after all" suggests an unstated objective with regard to which different approaches can be ranked, going beyond properties shared by all representations. Adriaan van Wijngaarden has attempted to axiomatize these shared properties (in Section 2 of "Numerical analysis as an independent science", BIT 6 (1966) 66–81). I don't know if IEEE floating-point arithmetic actually satisfies these axioms. --Lambiam 13:07, 23 April 2021 (UTC)[reply]

I'm looking for some way to generate random floats "naturally" (ignoring implementation specifics of various floating point formats, that is). My first thought was to repeatedly "partition" the number until some "exit condition" was met. In this case, it was to divide the subrange in half until it went to zero (I've since switched to checking for an unchanged value). One glaring problem with that approach however is that it is quite inefficient. 1000+ iterations for IEEE 754 types and something like 25,000+ iterations for a "long double" C type. So I don't know. It would be nice to have a generic method for this, but at THAT cost it is a bit hard to justify. Earl of Arundel (talk) 14:13, 23 April 2021 (UTC)[reply]

Without knowing which OS you run it's hard to advise. random() and rand() should be present on all systems and generate pseudo-random numbers. If you generate the seed from a good source they are adequate for most purposes. If you need truly random numbers, and are on a Linux system, then look at man 4 random and user either /dev/random or /dev/urandom as appropriate.

I was actually referring to the problem in more generic terms. Which is to say computing a floating point number within some range (0 to 1 for example) using any given random bitstream, be it from a PRNG or some stochastic source (such as /dev/urandom). Bits are extracted from one end and then "somehow" applied in such a way as to compute a real number. So ideally I would like to avoid ALL platform specific assumptions if at all possible. Maybe there is some way to "naturally" detect the number of bits in the mantissa of a given type? Earl of Arundel (talk) 16:05, 23 April 2021 (UTC)[reply]

The following algorithm, expressed in pseudocode, should give the number of mantissa bits.

Set x = 1.0, n = 0.

While 0.5 * (x + 0.5) < x, set x = 0.5 * (x + 0.5), n = n + 1.

If x > 0.5, set n = n + 1.

Assume you have a sequence b₀, ..., b_n−1 of random bits. To convert this sequence to a number r in the range 0.0 <= r < 1.0, execute the following:

Set r = 0.0.

For i from 0 to n − 1, set r = 0.5 * (r + b_i).

--Lambiam 10:15, 24 April 2021 (UTC)[reply]

Absolutely brilliant, Lambiam! How does your algorithm actually work though, if I may ask? Earl of Arundel (talk) 13:35, 24 April 2021 (UTC)[reply]

For the first part, it is best understood by looking at the binary representation of the course of values of x in successive transformations. Assuming a mantissa of 8 bits, we get:

x n

1.0000000 0

0.11000000 1

0.10100000 2

0.10010000 3

0.10001000 4

0.10000100 5

0.10000010 6

0.10000001 7

0.10000000 8

In the initial rounds, the invariant holds that x = (1+2⁻ⁿ)/2, which implies that 0.5 < x ≤ 1. Mathematically, this should keep holding indefinitely, but because of the limited precision of floating point arithmetic there cannot be an indefinitely descending sequence; it breaks down when the 2⁻⁽ⁿ⁺¹⁾ bit either gets shifted out – which happens in IEEE floating point – or gets stuck in the least significant position. The part generating a random value should also be obvious if you look at the successive binary representations of r: the bits in the mantissa of the values are those of the random bitstream in reverse order. --Lambiam 23:47, 24 April 2021 (UTC)[reply]

I see. Well it's a very elegant solution. Much simpler than anything I could come up with, that's for sure! Thanks again...Earl of Arundel (talk) 01:08, 25 April 2021 (UTC)[reply]

If you're programming in C, you can do some interesting low-level things with floating-point numbers, in a machine-independent way, using the library functions ldexp, frexp, and modf. —Steve Summit (talk) 16:13, 28 April 2021 (UTC)[reply]