I think this should have an attribution, something like

Based on the VCL library, which is (c) Copyright 2012-2020 Agner Fog and licensed under
the Apache License version 2.0.

@rgommers is that correct?

I think this should have an attribution, something like

Based on the VCL library, which is (c) Copyright 2012-2020 Agner Fog and licensed under
the Apache License version 2.0.

@rgommers is that correct?

No that sounds wrong. Reminder, we don't do Apache2, see for example: #13447 (comment).

We could change that at some point, but it's a significant change and needs a strong motivation and decision on the mailing list.

VCL is a little hard to find, so here is a link: https://github.com/vectorclass/version2. If this PR took code from there, that seems like a blocker for accepting it.

Now read the PR description, the whole PR is based on that, and it seems valuable. It may be worth considering, especially if there are no good alternatives. We'd be deciding to give up on GPLv2 compatibility.

@rgommers, we can ask for re-licensing, since we're not taking the same exact code plus the original code itself is based on T. Granlund and P. L. Montgomery work.

Most of the code was rewritten, while some of the code is very similar to Apache 2.0 based VSL, which is not compatible with current 3-clause BSD License(Apache 2.0 is more restrictive), so we should either rewrite the similar code or ask Agner Fog modify the License.

@ganesh-k13, I added intrinsics for the 64-bit division so we can drop libdiv.
@rgommers, @mattip, I changed the term "based" to "adopted" since we're not using the same code.
I don't think we need to leave any license notice because the original work belong to "T. Granlund and P. L. Montgomery." in the first place.

@mattip,

Coverage is missing for some of the code paths. Is there a way to increase the testing?

I added a missing case when the divisor equals 2, but still not enough if the CI worker doesn't support avx512.
anyway, the current worker seems to support it and codecov/patch — 97.08% of diff hit (target 82.50%).

It seems there is not a test for divisor == 0? Or maybe that case is not needed, so maybe add a comment where coverage is complaining that the code path should never be reached

@mattip, calling npyv_divisor_{dtype}(0) from python or c level will raise a hardware arithmetic exception. I added LCOV_EXCL_LINE instead to suppress the warning.

Thanks @seiko2plus

Hey @seiko2plus , I was trying to add dispatch for signed and in the case of signed 64 bits, the "round down" is not happening. Here is the code used to test it: diff.

s8, s16, and s32 seem to work fine.
Simple reproducer

Just confirming from bt:

(gdb) n
simd_divide_by_scalar_contig_s64 (args=0x7fffffffd4b0, len=5) at numpy/core/src/umath/loops_arithmetic.dispatch.c.src:40
40          for (; len >= vstep; len -= vstep, src += vstep, dst += vstep) {
(gdb) n
41              npyv_@sfx@ a = npyv_load_@sfx@(src);
(gdb) 
42              npyv_@sfx@ c = npyv_divc_@sfx@(a, divisor);
(gdb) 
43              npyv_store_@sfx@(dst, c);
(gdb) p a
$1 = {-9223372036854775808, -4611686018427387904, 0, 4611686018427387903}
(gdb) p scalar
$7 = -9223372036854775808
(gdb) p c
$8 = {1, 0, 0, 0}

The last 0 should be -1 as it is 4611686018427387903//-9223372036854775808. This is true for all dividends > 0 with divisors < 0.

@ganesh-k13, this patch follows the C compilers behavior for the signed division which mean:
1- Quotient rounded towards zero(truncate) rather than negative infinity(floor)
2- Wrap around overflow rather than returning zero and raising arithmetic exception.
This case isn't standard but most compilers follow this behavior at least when optimization is enabled

Here an example in Python for implementing floor division on the top truncate division:

"""
the method that been used in this implementation is based on T. Granlund and P. L. Montgomery
−dsign = SRL(d, N − 1);
−nsign = (n < −dsign );
/* SLT, signed */
−qsign = EOR(−nsign , −dsign );
q = TRUNC((n − (−dsign ) + (−nsign ))/d) − (−qsign );
"""
from numpy.core._simd import baseline as v

for d in [-100, -1, -50, 100, 1, 100]:
    nsign_d = v.setall_s16(int(d < 0))
    divisor = v.divisor_s16(d)
    noverflow = v.cvt_b16_s16(v.setall_s16(-1))
    for n in [-0x8000, -255, 255, -100, 1, 100]:
        a         = v.load_s16(range(n, n + v.nlanes_s16))
        nsign_a   = v.cvt_s16_b16(v.cmplt_s16(a, nsign_d))
        nsign_a   = v.and_s16(nsign_a, v.setall_s16(1))
        diff_sign = v.sub_s16(nsign_a, nsign_d)
        to_ninf   = v.xor_s16(nsign_a, nsign_d)
        trunc     = v.divc_s16(v.add_s16(a, diff_sign), divisor)
        if d == -1:
            greater_min = v.cmpgt_s16(a, v.setall_s16(-0x8000))
            noverflow = v.and_b16(noverflow, greater_min)
            floor = v.ifsub_s16(greater_min, trunc, to_ninf, v.zero_s16())
        else:
            floor = v.sub_s16(trunc, to_ninf)
        assert floor == [x//d for x in a]

    if v.tobits_b16(noverflow) != (1 << v.nlanes_s16)-1:
        0//0

Thanks a lot for the info @seiko2plus. As python prefers round to -inf, shall I replace/add the npyv_divc_s* with round to -inf as shown in paper? Will be happy to raise a PR with round to -inf. [EDIT] Or do we build on top of trunc as shown above?

[EDIT 2] My main question is, do we edit npyv_divc_s* or add code in the dispatch to modify the divisor.

@ganesh-k13, just handle it directly within the dispatch-able source file, note we will still need to use truncate division for timedelta so we can completely erase libdivide.