CODEWITHSUNDEEP
algorithmmathfloating-point
relatively_random
relatively_random

Reputation: 5126

Calculating a floating-point value to multiply an integer with to produce largest number less than 1

For a given 32-bit integer y, pick a double-precision floating point number x so that x*(double)y is as large as possible, but strictly less than 1.

Is there an algorithm that can calculate this?

I'm mostly just curious, I don't have a real problem which can't be solved in another way.

Update:

Question is a bit imprecise so here's an update:

Thanks to @DavidEisenstat and @EricPostpischil for noticing these.

Upvotes: 2

Views: 251

Answers (2)

Rory Daulton
Rory Daulton

Reputation: 22544

Here is an algorithm that solves your problem but depends on some bit-twiddling.

First calculate and store

x_attempt = 1.0 / (double) y

This value will be very close to your desired value of x but it may be a bit off due to the approximations in floating-point math.

Now check the value of x_attempt * (double) y. If that value is not strictly less than 1.0, "nudge" the value of x_attempt down to the next smaller value that can be stored in a double variable. Check x_attempt * (double) y again and keep nudging down until the value is strictly less than 1.0. This uses a loop but it will execute very few times, assuming that floating-point arithmetic on your platform is any good at all.

If that value is strictly less than 1.0, "nudge" the value of x_attempt up until that product is 1.0 or more. Then set x to the previous value of x_attempt.

I implemented my own "nudge" up and down routines in Borland Delphi back when that was my preferred language. Ask if you need help in coding such routines for your language and environment.

I was motivated by your question to code the "nudge" up and down routines in my current preferred language, Python 3.7. Here they are, without my unit testing routines.

import math

MIN_DOUBLE = 2 ** -1074  # 4.940656458412465442e-324
MIN_DOUBLE_DENORMAL = MIN_DOUBLE
MIN_DOUBLE_NORMAL = 2 ** -1022  # 2.225073858507201383e-308
MAX_DOUBLE = 1.7976931348623157082e308  # just under 2 ** 1024
EPSILON_DOUBLE_HALF = 2 ** -53  # changes most significand values
EPSILON_DOUBLE_FOURTH = 2 ** -54  # changes significand of -0.5


def nudged_up(x: float) -> float:
    """Return the next larger float value.

    NOTES:  1.  This routine expects that `float` values are implemented
                as  IEEE-754 binary64 and includes denormal values. No
                check is done on these assumptions.
            2.  This raises an OverflowError for MAX_DOUBLE. All other
                float values do not raise an error.
            3.  For +INF, -INF, or NAN, the same value is returned.
    """
    if x == 0.0:
        return MIN_DOUBLE_DENORMAL
    significand, exponent = math.frexp(x)
    if exponent < -1021:  # denormal
        return x + MIN_DOUBLE_DENORMAL
    if significand == -0.5:  # nudging will change order of magnitude
        diff = EPSILON_DOUBLE_FOURTH
    else:  # usual case
        diff = EPSILON_DOUBLE_HALF
    return math.ldexp(significand + diff, exponent)


def nudged_down(x: float) -> float:
    """Return the next smaller float value.

    NOTES:  1.  This routine expects that `float` values are implemented
                as  IEEE-754 binary64 and includes denormal values. No
                check is done on these assumptions.
            2.  This raises an OverflowError for -MAX_DOUBLE. All other
                float values do not raise an error.
            3.  For +INF, -INF, or NAN, the same value is returned.
    """
    return -nudged_up(-x)

And here is Python 3.7 code that answers your problem. Note that an error is raised if the input parameter is zero--you may want to change that.

def lower_reciprocal(y: int) -> float:
    """Given a (32-bit) integer y, return a float x for which
    x * float(y) is as large as possible but strictly less than 1.

    NOTES:  1.  If y is zero, a ZeroDivisionError exception is raised.
    """
    if y < 0:
        return -lower_reciprocal(-y)
    float_y = float(y)
    x = 1.0 / float_y
    while x * float_y < 1.0:
        x = nudged_up(x)
    while x * float_y >= 1.0:
        x = nudged_down(x)
    return x

Upvotes: 4

alias
alias

Reputation: 30428

This is probably not what you had in mind, but a great alternative for answering such questions is to use off-the-shelf SMT solvers, which can solve constraint satisfaction problems over many domains; including IEEE floating point. See http://smtlib.cs.uiowa.edu/ for details.

A good (and free!) SMT solver is Microsoft's Z3: https://github.com/Z3Prover/z3

SMT-solvers are scripted using the so called SMTLib language, which is more machine oriented. But they also provide programmatic APIs for many languages (C, C++, Java, Python) and there are many high-level language bindings that make them easy to use, including those from Scala, O'Caml, Go, Haskell, to name a few.

For instance, here's how one can code your query in Haskell's SBV binding to Z3. If you don't read Haskell, do not worry; the idea is to illustrate how to code such problems quickly and at a very high level:

import Data.SBV

q :: SInt32 -> IO Double
q y = do LexicographicResult m <- optimize Lexicographic $ do
                                      x <- sDouble "x"
                                      let res = (sFromIntegral (y::SInt32) * x)
                                      constrain $ res .< 1
                                      maximize "mx" res
         case getModelValue "x" m of
           Just x  -> return x
           Nothing -> error "No such value exists!" -- shouldn't happen!

With his piece of code, we can use the Haskell interpreter ghci (https://www.haskell.org/ghc/) to query various values:

*Main> q 12
8.333333333333331e-2
*Main> q 821
1.2180267965895247e-3

These answers are quick to compute, and they can be further used and combined with other constraints to find such values of interest easily.

If you're interested also check out Python bindings for z3, which are probably the easiest to get started with, though not as type-safe as other bindings, as you can imagine. See here: https://ericpony.github.io/z3py-tutorial/guide-examples.htm

Upvotes: 1

Related Questions