random.choice IndexError due to double-rounding

[stevendaprano]

BPO	24567
Nosy	@tim-one, @rhettinger, @mdickinson, @pitrou, @vstinner, @stevendaprano, @bitdancer, @skrah, @serhiy-storchaka
PRs	bpo-24567: Random subnormal.diff #7954 [3.7] bpo-24567: Random subnormal.diff (GH-7954) #7955 [3.6] bpo-24567: Random subnormal.diff (GH-7954) #7956
Files	choice_timings.txt shuffle_timings.txt handle_double_rounding2.diff handle_double_rounding3.diff sample_timings.txt random_double_round.diff: For Python 3.6 and 3.7

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = 'https://github.com/rhettinger'
closed_at = <Date 2018-06-27.07:07:38.703>
created_at = <Date 2015-07-05.03:04:00.634>
labels = ['3.7', 'type-bug', 'library']
title = 'random.choice IndexError due to double-rounding'
updated_at = <Date 2018-06-27.09:37:20.023>
user = 'https://github.com/stevendaprano'

bugs.python.org fields:

activity = <Date 2018-06-27.09:37:20.023>
actor = 'rhettinger'
assignee = 'rhettinger'
closed = True
closed_date = <Date 2018-06-27.07:07:38.703>
closer = 'rhettinger'
components = ['Library (Lib)']
creation = <Date 2015-07-05.03:04:00.634>
creator = 'steven.daprano'
dependencies = []
files = ['39897', '39898', '39905', '39908', '39937', '47650']
hgrepos = []
issue_num = 24567
keywords = ['patch']
message_count = 61.0
messages = ['246291', '246297', '246303', '246317', '246318', '246336', '246352', '246363', '246366', '246484', '246530', '246531', '246532', '246534', '246536', '246537', '246548', '246550', '246551', '246565', '246566', '246581', '246615', '246616', '246623', '246628', '246630', '246636', '246639', '246641', '246642', '246650', '246656', '246657', '246664', '246668', '246674', '246815', '267785', '267874', '275594', '277199', '277200', '304734', '320392', '320394', '320423', '320434', '320443', '320474', '320479', '320502', '320511', '320516', '320519', '320521', '320526', '320549', '320553', '320555', '320563']
nosy_count = 10.0
nosy_names = ['tim.peters', 'rhettinger', 'mark.dickinson', 'pitrou', 'vstinner', 'steven.daprano', 'r.david.murray', 'skrah', 'serhiy.storchaka', 'Serge Anuchin']
pr_nums = ['7954', '7955', '7956']
priority = 'normal'
resolution = 'wont fix'
stage = 'resolved'
status = 'closed'
superseder = None
type = 'behavior'
url = 'https://bugs.python.org/issue24567'
versions = ['Python 2.7', 'Python 3.6', 'Python 3.7']

[stevendaprano]

While investigating bpo-24546, I discovered that at least some versions of Python on 32-bit Linux have a double-rounding bug in multiplication which may lead to an unexpected IndexError in random.choice.

See http://bugs.python.org/issue24546 for additional discussion, but the short version is that due to double-rounding, if random.random returns (1. - 2.**-53), int(i * random.random()) may return 1 for some i >= 2049, and random.choice will then evaluate seq[len(seq)], raising IndexError.

[skrah]

I've finally gotten hold of a gcc (5.1.0) that produces the unwanted code.

For a *different* test case, even the -O0 and -O2 results differ, since
with -O2 gcc uses mpfr (I think!) for constant folding, which produces
the correct result.

Using -mfpmath=sse seems to fix all test cases that I've tried.

[vstinner]

There is a simple fix. If the result is invalid, drop it and generate a new
number. It should keep the distribution uniform.

Another fix is to only use integers. Why do we use float in choice by the
way?

[mdickinson]

Note that this affects other Random methods, too. There are several places in the source where something of the form int(i * random.random()) is used to select an integer in the half-open range [0, i).

And a nitpick: it's a stretch to describe double-rounding on multiplication (or any other arithmetic operation) as a bug: Python-the-language makes no promises about adhering to IEEE 754 arithmetic rule. (In fact, it makes no promises about using IEEE 745 formats in the first place.)

[mdickinson]

IEEE 745

IEEE 754. Stupid fingers.

[tim-one]

I suppose the simplest "fix" would be to replace relevant instances of

int(random() * N)

with

min(int(random() * N), N-1)

That would be robust against many kinds of arithmetic quirks, and ensure that platforms with and without such quirks would, if forced to the same internal random state, yield the same result when random() happens to return 1. - 2.**-53.

Victor, what - precisely - do you have in mind with "only use integers"? The current method was used because it's simple, efficient, and obvious. Python 3.4.3 appears to (usually) use a different method, though, relying on C code to generate random bytes without going through floats.

[skrah]

Python-the-language makes no promises about adhering to IEEE 754 arithmetic rule.

Still, I think we could switch to -mfpmath=sse, which would at least
guarantee consistent behavior across different gcc versions.

[mdickinson]

Still, I think we could switch to -mfpmath=sse, which would at least
guarantee consistent behavior across different gcc versions.

... which would fix fsum on 32-bit Linux, too. (I thought there was an open issue for this, but now I can't find it.)

[tim-one]

Mark, closest I could find to a substantive SSE-vs-fsum report is here, but it was closed (because the fsum tests were changed to ignore the problem ;-) ):

http://bugs.python.org/issue5593

[skrah]

Qt had a similar initiative regarding -msse2 -mfpmath:

http://lists.qt-project.org/pipermail/development/2013-December/014530.html

They say that also Visual Studio 2012 has switched to sse2 by default.

The only problem are the Linux distributions that are stuck in i386 land
and who would probably disable the flags in their builds.

[pitrou]

I suppose the simplest "fix" would be to replace relevant instances of

int(random() * N)

with

min(int(random() * N), N-1)

That sounds simple and generic. It skews the distribution a tiny little bit, but it doesn't sound significant (perhaps Mark would disagree).

[tim-one]

It skews the distribution a tiny little bit, ...

But it doesn't - that's the point ;-)

If double-rounding doesn't occur at all (which appears to be the case on most platforms), absolutely nothing changes (because min(int(random() * N), N-1) == int(random() * N) on such boxes).

If double-rounding does occur, double-rounding itself may change results "all over the place", and I haven't tried to analyze what effects that has on the distribution. In the comparative handful of cases where int(random() * N) == N on such boxes, clamping that back to N-1 just yields the same result we would have gotten on a box that didn't do double-rounding.

[vstinner]

Again, why not using only integers?

Pseudo-code to only use integers:

def randint(a, b):
    num = b - a
    if not num:
        return a
    nbits = (num + 1).bit_length()
    while True:
        x = random.getrandbits(nbits)
        if x <= num:
            break
    return a + x

(It doesn't handle negative numbers nor empty ranges.)

[tim-one]

Victor, if people want to use getrandbits(), we should backport the Python3 code, not reinvent it from scratch.

Note too Mark's comment: "There are several places in the source where something of the form int(i * random.random()) is used". The min() trick is relatively easy to apply via simple, local editing of such places.

Your "It doesn't handle negative numbers nor empty ranges" should be a hint about what a pain it is to rewrite everything wholesale to use a _fundamentally_ different method.