Implement floordiv reduction in parfors.

[sklam]

assume
x_init // x0 // x1 // x2 // ... // xn
can be computed as
floor((max_i64 // x0 // x1 // x2 // ... // xn) * (x_init / max_i64))

However, there may be numerical problems:

from numba import njit, prange
import numpy as np


@njit(parallel=True)
def foo(x, ys):
    for i in prange(ys.size):
        x //= ys[i]
    return x


ys = np.array([5, 3, 2, 1, 5, 1, 8, 1, 1, 5, 4, 1, 8, 3, 1, 1, 2, 2, 2, 7, 6, 7, 1, 7, 4, 7, 1])
init = 9999 * np.prod(ys)
res = foo(init, ys)
print("res = ", res)
exp = foo.py_func(init, ys)
print("exp = ", exp)
assert res == exp

The above will fail with:

res =  9998
exp =  9999
Traceback (most recent call last):
  File "/Users/siu/dev/numba/chk.py", line 36, in <module>
    assert res == exp
AssertionError

[sklam]

The error has to do with the precision loss in the float64 division during the finalizing reduction for all the threads intermediate. I think it is more stable if the reduction is done in descending values, for example:

import numpy as np
from math import floor

max_i64 = np.iinfo(np.int64).max

def simulated_reduce_finalization(arr, init):
    c = 1.0

    for x in reversed(sorted(list(arr))):
        c *= x  / max_i64


    print("before floor", init * c)
    print("     floored", floor(init * c))
    return floor(init * c)

# The below values are obtained from debug print after the gufunc
intermediates = [
307445734561825860,
1844674407370955161,
1152921504606846975,
461168601842738790,
1152921504606846975,
3074457345618258602,
4611686018427387903,
2305843009213693951,
219604096115589900,
1317624576693539401,
329406144173384850,
1317624576693539401,

]

data = np.array(intermediates, dtype=np.int64)
init = 2655048388608000
res = simulated_reduce_finalization(data, init)
print("res", res)

will print:

before floor 9999.0
     floored 9999
res 9999

[sklam]

Saving a test script here:

from numba import njit, prange
import numpy as np

stats = dict(ok=0, bad=0)


@njit(parallel=True)
def foo(x, ys):
    for i in prange(ys.size):
        x //= ys[i]
    return x


def test_one(val, a):
    expected = foo.py_func(val, a)
    got = foo(val, a)
    if expected != got:
        stats['bad'] += 1

    else:
        stats['ok'] += 1


def test_divs(val):
    for i in range(4,25):
        a = np.random.randint(1, 10, size=i)
        test_one(val, a)


def test():
    np.random.seed(0)

    for i in range(10000):
        for j in range(1):
            test_divs(i+1)

    print(stats)

test()

[sklam]

The numeric issues with this approach will make it impractical