Feat/mul div

[aritkulova]

• This PR suggests a bug fix and I've added the necessary tests.
• This PR introduces a new feature and I've discussed the update in an Issue or with the team.
• This PR is just a minor change like a typo fix.

---

[Hrom131]

Shouldn't we panic if we divide by 0?

[aritkulova]

I agree that panicking on division by zero would be less error-prone, but the jet's division functions don't panic on division by zero, so here it feels more like the default behavior

[schoen]

We had a discussion on a call with @roconnor-blockstream about this behavior, which other people found unusual or surprising. I think it has to do with some functional programming conventions related to wanting functions to be total when possible. I don't know that he necessarily felt that all library functions written using the jet needed to replicate the behavior, but I don't remember his exact position about this.

I think @apoelstra is also familiar with this discussion and might be able to weigh in.

[LesterEvSe]

Also let's add more tests to cover new u512 division functionality.

1. Force entry into algorithm_d_512_256:

a = generate_u256(2^128, U256::MAX)
b = generate_u256(2^128, U256::MAX)

result_high = high 256 bits of a.full_mul(b)
c = generate_u256(max(2^128, result_high + 1), U256::MAX)

Guarantees result_high != 0 and denom_high != 0.

2. Add two normalize_to_threshold_512_127 tests.
   2.1. norm == 1. Forces the divisor's top bit to already be set:

a = generate_u256(2^128, U256::MAX)
b = generate_u256(2^128, U256::MAX)

result_high = high 256 bits of a.full_mul(b)
c = generate_u256(max(2^255, result_high + 1), U256::MAX)

2.2. norm > 1. Forces the divisor's top bit to be clear:

a = generate_u256(2^128, U256::MAX)
b = generate_u256(2^128, 2^192)

result_high = high 256 bits of a.full_mul(b)   // < 2^192, well under the 2^255 ceiling
c = generate_u256(max(2^128, result_high + 1), 2^255 - 1)

3. Same a, b as (1), but c = result_high + 1.

4. Exact division with remainder == 0. Set a = c for some c in range [2^128, U256::MAX], pick b as a*b >= 2^256.

5. Minimal denom_high:

a = generate_u256(2^128, 2^129)
b = generate_u256(2^128, 2^129)

result_high = high 256 bits of a.full_mul(b)   // in {1, 2, 3}
c = generate_u256(max(2^128, result_high + 1), 2^129 - 1)