[Libre-soc-isa] [Bug 817] Big Integer Math (sv.adde, sv.subfe, sv.madded, 128 by 64-bit -> 64-bit div/rem, maybe more...)

[bugzilla-daemon at libre-soc.org]

https://bugs.libre-soc.org/show_bug.cgi?id=817

--- Comment #32 from Jacob Lifshay <programmerjake at gmail.com> ---
(In reply to Luke Kenneth Casson Leighton from comment #29)
>  128         // Compute estimate qhat of q[j] from top 2 digits
>  129         uint64_t dig2 = ((uint64_t)un[j + n] << 32) | un[j + n - 1];
>  130         qhat = dig2 / vn[n - 1];
>  131         rhat = dig2 % vn[n - 1];
> 
> ehm.... ehmehmehm... qhat and rhat are uint64_t (and un/vn 32 bit*)
> therefore when vn/un
> are uint64_t[] qhat and rhat must be uint128_t which is NOT happening,
> we are NOT going to return QTY 4 64 bit registers from divmod2du.
> 
> thus even on scalar usage if un[j+n] > vn[n-1] an overflow on the 128/64
> divide will occur.
> 
> any clues?

yes...knuth's original algorithm is written in mix assembly where the operation
used is a 2-word by 1-word div/rem where both quotient and remainder are 1
word. the only case where that can overflow is where un[j + n] == vn[n - 1] in
which case it sets qhat to the max value that fits in a single word (in the
64-bit case, 2**64-1).

the initial normalization step creates the invariant that un[j + n] <= vn[n -
1] and that invariant is maintained by the loop.

the c code just has qhat and rhat be 2-word values to prevent UB when the
overflow occurs, and also for convenience to avoid needing to constantly cast
them to 2-word values. they can totally be 1-word values at the cost of putting
the initial qhat computation inside of an if that checks for overflow.

original algorithm step:
D3. [Calculate q̂.] If uⱼ = v₁, set q̂ ← b - 1; otherwise set q̂ ← ⌊(uⱼb +
uⱼ₊₁)/v₁⌋. Now test if v₂q̂ > (uⱼb + uⱼ₊₁ - q̂v₁)b + uⱼ₊₂; if so, decrease q̂ by 1
and repeat this test. (The latter test determines at high speed most of the
cases in which the trial value q̂ is one too large, and it eliminates all cases
where q̂ is two too large; see exercises 19, 20, 21.)

in ascii:
D3. [Calculate qhat.] If u[j] = v[1], set qhat <- b - 1; otherwise set qhat <-
floor((u[j]*b + u[j+1])/v[1]). Now test if v[2]*qhat > (u[j]*b + u[j+1] -
qhat*v[1])*b + u[j+2]; if so, decrease qhat by 1 and repeat this test. (The
latter test determines at high speed most of the cases in which the trial value
qhat is one too large, and it eliminates all cases where qhat is two too large;
see exercises 19, 20, 21.)

-- 
You are receiving this mail because:
You are on the CC list for the bug.