Anton Ertl <anton@mips.complang.tuwien.ac.at> schrieb:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>Since then every compiler where I have checked the asm/machine code
>>output will use effectively the same approach, sometimes simplifying it
>>slightly for smaller divisors.
>
> For u/7 on AMD64 (latency measured on Skylake) there are at least the
> following sequences:
>
> gcc Robison/fish ertl
> latency=8c latency=6.25c latency=6c
> %Cl=0x4924924924924925
> %C=0x2492492492492493 %C=0x9249249249249249 %Ch=0x2492492492492492
> movabs $C,%rdx movabs $C,%rax movabs $Cl,%rax
> mov %rdi,%rax mov %rax,%rcx mul %rdi
> mul %rdx mul %rdi mov %rdx,%rcx
> sub %rdx,%rdi add %rcx,%rax movabs $Ch,%rax
> shr %rdi adc $0x0,%rdx mul %rdi
> lea (%rdx,%rdi),%rax shr $0x2,%rdx add %rcx,%rax
> shr $0x2,%rax mov %rdx,%rax adc $0x0,%rdx
> mov %rdx,%rax

Makes me wonder... what do you think of the 128-bit division and
remainder by some constants with gcc11, like

unsigned int divrem_10 (__uint128_t x, __uint128_t *div)
{ *div = x / 10;
return x % 10;
}

? (see https://godbolt.org/z/czqWEn8jo ).

Terje Mathisen <terje.mathisen@tmsw.no> writes:
>I suspect that for divisors with less than max significant bits, your
>version will also generate a sufficiently accurate factional part which,
>after multiplication, can generate the remainder from the division.

Interesting idea, but I did not investigate it; the classical approach
of multiplying the quotient by the divisor and subtracting the result
from the dividend looked good enough (and it's not clear if it's
slower).

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

Thomas Koenig wrote:
> Anton Ertl <anton@mips.complang.tuwien.ac.at> schrieb:
>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>> Since then every compiler where I have checked the asm/machine code
>>> output will use effectively the same approach, sometimes simplifying it
>>> slightly for smaller divisors.
>>
>> For u/7 on AMD64 (latency measured on Skylake) there are at least the
>> following sequences:
>>
>> gcc Robison/fish ertl
>> latency=8c latency=6.25c latency=6c
>> %Cl=0x4924924924924925
>> %C=0x2492492492492493 %C=0x9249249249249249 %Ch=0x2492492492492492
>> movabs $C,%rdx movabs $C,%rax movabs $Cl,%rax
>> mov %rdi,%rax mov %rax,%rcx mul %rdi
>> mul %rdx mul %rdi mov %rdx,%rcx
>> sub %rdx,%rdi add %rcx,%rax movabs $Ch,%rax
>> shr %rdi adc $0x0,%rdx mul %rdi
>> lea (%rdx,%rdi),%rax shr $0x2,%rdx add %rcx,%rax
>> shr $0x2,%rax mov %rdx,%rax adc $0x0,%rdx
>> mov %rdx,%rax
>
> Makes me wonder... what do you think of the 128-bit division and
> remainder by some constants with gcc11, like
>
> unsigned int divrem_10 (__uint128_t x, __uint128_t *div)
> {
> *div = x / 10;
> return x % 10;
> }
>
> ? (see https://godbolt.org/z/czqWEn8jo ).
>

It looks like 25-30 cycles of latency on a 3+ wide core?

It isn't horrible, but it also isn't the best we can do:

What I think is that this particular code will pretty much always occur
in the context of transforming multiple digits to decimal, in which case
I am sure I could outperform a loop of this code by a big factor using
my divide & conquer approach:

An uint128 can require up to 39 decimal digits, so we start by splitting
the input into 13-digit chunks with a constant div-mod by 10^13, then we
can extract all 39 digits in parallel:

Start by scaling each chunk by 2^60/1e12, leaving the top digit in the
top 4 bits of a 64-bit register, then repeatedly mask away those 4 bits,
multiply the remainder by 5 (a single-cycle LEA) and move the mask down
one bit position, so effectively a mul-by-10.

Since all these operations are fast and unrolled by the three parallel
chunks, we'll get close to a digit/cycle, possibly staying under 100
cycles for the full conversion?

An easier approach will do the same using divrem_1e13() then convert
each chunk using standard 64-bit operations:

void u128_to_ascii(char *buf, int buflen, __uint128_t n)
{ if (buflen < 40) return null;
__uint128_t t;
uint64_t c0, c1, c2;

t = n / (__uint128_t) 10000000000000LL;
c0 = (uint64_t) n % 10000000000000LL;

c2 = (uint64_t) t / (__uint128_t) 10000000000000LL;
c0 = (uint64_t) t % 10000000000000LL;

sprintf(buf,"%013d%013d%013d", c2, c1, c0);

return buf;
}

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Anton Ertl wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>> I suspect that for divisors with less than max significant bits, your
>> version will also generate a sufficiently accurate factional part which,
>> after multiplication, can generate the remainder from the division.
>
> Interesting idea, but I did not investigate it; the classical approach
> of multiplying the quotient by the divisor and subtracting the result
> from the dividend looked good enough (and it's not clear if it's
> slower).

It is probably just marginally slower, by needing a couple more
housekeeping instructions, but nothing significant.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Terje Mathisen <terje.mathisen@tmsw.no> schrieb:
> Thomas Koenig wrote:
>> Anton Ertl <anton@mips.complang.tuwien.ac.at> schrieb:
>>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>>> Since then every compiler where I have checked the asm/machine code
>>>> output will use effectively the same approach, sometimes simplifying it
>>>> slightly for smaller divisors.
>>>
>>> For u/7 on AMD64 (latency measured on Skylake) there are at least the
>>> following sequences:
>>>
>>> gcc Robison/fish ertl
>>> latency=8c latency=6.25c latency=6c
>>> %Cl=0x4924924924924925
>>> %C=0x2492492492492493 %C=0x9249249249249249 %Ch=0x2492492492492492
>>> movabs $C,%rdx movabs $C,%rax movabs $Cl,%rax
>>> mov %rdi,%rax mov %rax,%rcx mul %rdi
>>> mul %rdx mul %rdi mov %rdx,%rcx
>>> sub %rdx,%rdi add %rcx,%rax movabs $Ch,%rax
>>> shr %rdi adc $0x0,%rdx mul %rdi
>>> lea (%rdx,%rdi),%rax shr $0x2,%rdx add %rcx,%rax
>>> shr $0x2,%rax mov %rdx,%rax adc $0x0,%rdx
>>> mov %rdx,%rax
>>
>> Makes me wonder... what do you think of the 128-bit division and
>> remainder by some constants with gcc11, like
>>
>> unsigned int divrem_10 (__uint128_t x, __uint128_t *div)
>> {
>> *div = x / 10;
>> return x % 10;
>> }
>>
>> ? (see https://godbolt.org/z/czqWEn8jo ).
>>
>
> It looks like 25-30 cycles of latency on a 3+ wide core?
>
> It isn't horrible, but it also isn't the best we can do:
>
> What I think is that this particular code will pretty much always occur
> in the context of transforming multiple digits to decimal, in which case
> I am sure I could outperform a loop of this code by a big factor using
> my divide & conquer approach:

Actually, truth be told, this was partially inspired by my search
to find non-trivial numbers whose sum of digits is the same for
every prime number up to 17 (you may remember the discussion about
this some time ago).

I've found 70911040973874056146188543 and 77332999599545910254098143
with an algorithm that used the approach above, which gained me
an OEIS entry (even if it is one of the "less interesting" ones).

> An uint128 can require up to 39 decimal digits, so we start by splitting
> the input into 13-digit chunks with a constant div-mod by 10^13, then we
> can extract all 39 digits in parallel:

Hm, yes, that approach could have sped up my calculations somewhat :-)

> Start by scaling each chunk by 2^60/1e12, leaving the top digit in the
> top 4 bits of a 64-bit register, then repeatedly mask away those 4 bits,
> multiply the remainder by 5 (a single-cycle LEA) and move the mask down
> one bit position, so effectively a mul-by-10.
>
> Since all these operations are fast and unrolled by the three parallel
> chunks, we'll get close to a digit/cycle, possibly staying under 100
> cycles for the full conversion?

I am tempted to put this approach libgfortran, but unfortunately
the general overhad of Fortran I/O is very high.

> An easier approach will do the same using divrem_1e13() then convert
> each chunk using standard 64-bit operations:
>
> void u128_to_ascii(char *buf, int buflen, __uint128_t n)
> {
> if (buflen < 40) return null;
> __uint128_t t;
> uint64_t c0, c1, c2;
>
> t = n / (__uint128_t) 10000000000000LL;
> c0 = (uint64_t) n % 10000000000000LL;
>
> c2 = (uint64_t) t / (__uint128_t) 10000000000000LL;
> c0 = (uint64_t) t % 10000000000000LL;
>
> sprintf(buf,"%013d%013d%013d", c2, c1, c0);
>
> return buf;
> }

That, as well.