# Modular multiplication without overflow

**URL:** <https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421>\
**Category:** Performance\
**Created:** [November 18, 2022, 1:47am UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421 "2022-11-18T01:47:32Z")\
**Posts on this page:** 13\
**Page:** 2

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2022, 6:54pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/21 "2022-11-21T18:54:23Z")

</div>

I was going to say that perhaps I’m being dense but doesn’t this definition work well for all bit sizes less that 128:

```julia
mulmod(a::Integer, b::Integer, c::Integer) =
    mod(widemul(a, b), c) % typeof(c)

```

The tricky case is how to do this for 128-bit arguments where `c` actually is too big to fit in a 64-bit integer.

Also note: regrettably we followed the C definition of `%` and it actually does `rem` whereas what everyone wants here is `mod`.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [November 21, 2022, 7:29pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/22 "2022-11-21T19:29:43Z")

</div>

I agree that we should use native arithmetic (with widening) whenever possible. Although I’m not sure anyone is disputing that.

At the widest native types (eg `Int128`) I imagine it’s still worthwhile to _try_ to complete the operation in native arithmetic. A run-time overflow check can slow-path to the algorithm discussed here when required, but the significant speed difference makes the “try native first” approach enticing. I don’t imagine that very large arguments are common in most domains.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2022, 7:38pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/23 "2022-11-21T19:38:38Z")

</div>

Very closely related: factorials and binomial coefficients! Both of those are useful to have modulo some other number, because both overflow very quickly.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2022, 7:58pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/24 "2022-11-21T19:58:38Z")

</div>

Modular variants of those would certainly be useful to have! Not sure if they belong in base or a package though. Maybe `ModularMethods.jl` or something like that?

---

<div class="post-metadata">

**Author:** ![cjdoris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cjdoris/32/213133_2.png) [@cjdoris](https://discourse.julialang.org/u/cjdoris)\
**Post date:** [November 21, 2022, 9:06pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/25 "2022-11-21T21:06:07Z")

</div>

The OP was specifically about UInt128.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2022, 9:13pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/26 "2022-11-21T21:13:05Z")

</div>

Ah yeah, missed that!

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2022, 9:52pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/27 "2022-11-21T21:52:34Z")

</div>

I do think there should be a package specifically for modular arithmetic, but I think these two are basic enough they could go in base alongside the current modular power.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2022, 9:58pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/28 "2022-11-21T21:58:13Z")

</div>

You say regrettably–are you thinking of changing it in Julia 2.0? I think it’d be nice for a lot of people coming from Python, at the very least.

---

<div class="post-metadata">

**Author:** ![Lilith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lilith/32/27492_2.png) [@Lilith](https://discourse.julialang.org/u/Lilith)\
**Post date:** [November 22, 2022, 1:35am UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/29 "2022-11-22T01:35:24Z")

</div>

> [@mikmoore](#):
>
> you’d only need to resort to the complicated (and presumably much slower) version when all the above contingencies fail.

> [@StefanKarpinski](#):
>
> doesn’t this definition work well for all bit sizes less that 128:

> [@mikmoore](#):
>
> I agree that we should use native arithmetic (with widening) whenever possible. Although I’m not sure anyone is disputing that.

Surprisingly, the more complex method I defined [above](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/14) is actually 4x faster on UInt64s than promotion to UInt128 via the widemul method in microbenchmarks. Corroborating this, replacing the widemul approach in Base’s `powmod(::Int, ::Int, ::Int)` method with my ugly version also results in a significant speedup in that slightly-less-micro benchmark.

This is because while 128-bit integers are native in Julia they are not supported natively on my hardware. YMMV.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 22, 2022, 2:38am UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/30 "2022-11-22T02:38:03Z")

</div>

It would certainly be on a list of changes to consider.

---

<div class="post-metadata">

**Author:** ![cjdoris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cjdoris/32/213133_2.png) [@cjdoris](https://discourse.julialang.org/u/cjdoris)\
**Post date:** [November 22, 2022, 11:38am UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/31 "2022-11-22T11:38:55Z")

</div>

Here’s a version for `UInt128` which runs in 100ns on my computer (vs 122ns for lilith’s version and 350ns for the wiki version). Looking back, I think it is roughly the same algorithm as lilith’s.

Profiling reveals that almost all the time is spent evenly in the 8 division operations (`mod`, `fldmod` (counts for 2) and `sub_mod`). I’m not sure if you can do it in fewer divisions.

> **Code**
>
> ```julia
> using BitIntegers: UInt256
> using BenchmarkTools: @benchmark
> 
> function mul_full(x::UInt128, y::UInt128)
> z = UInt256(x) * UInt256(y)
> return (mod(z, UInt128), mod(z >> 128, UInt128))
> end
> 
> function sub_mod(x::UInt128, y::UInt128, m::UInt128)
> if x > y
> return mod(x - y, m)
> else
> return m - mod1(y - x, m)
> end
> end
> 
> function mul_mod(x::UInt128, y::UInt128, m::UInt128)
> # fast branch: no overflow
> if iszero(x >> 64) && iszero(y >> 64)
> return mod(x * y, m)
> end
> # fast branch: no overflow after reducing mod m
> if iszero(m >> 64)
> return mod(mod(x, m) * mod(y, m), m)
> end
> # m in base 2^64
> # m = m0 + m1 * u
> # where u = 2^64
> m0 = mod(mod(m, UInt64), UInt128)
> m1 = m >> 64
> # z=x*y in base 2^128
> # z == z0 + u^2 * z1
> z0, z1 = mul_full(x, y)
> # require z1 < m
> # so that a1 <= 2^64
> z1 = mod(z1, m)
> # reduce u * z1 by dividing by m1
> # and use the fact that u * m1 == -m0 (mod m)
> # so z == z0 + u * z2 (mod m)
> a1, a0 = fldmod(mod(z1, m), m1)
> z2 = sub_mod(a0 << 64, a1 * m0, m)
> # reduce u * z2 in the same way
> # so z == z0 - z3 (mod m)
> b1, b0 = fldmod(z2, m1)
> z3 = sub_mod(b1 * m0, b0 << 64, m)
> # final sub_mod gets the answer
> return sub_mod(z0, z3, m)
> end
> 
> function mul_mod_big(x::UInt128, y::UInt128, m::UInt128)
> return convert(UInt128, mod(big(x) * big(y), big(m)))
> end
> 
> function mul_mod_bad(x::UInt128, y::UInt128, m::UInt128)
> return mod(x * y, m)
> end
> 
> function mul_mod_256(x::UInt128, y::UInt128, m::UInt128)
> return convert(UInt128, mod(UInt256(x) * UInt256(y), UInt256(m)))
> end
> 
> function check(x, y, m)
> z0 = mul_mod(x, y, m)
> z1 = mul_mod_big(x, y, m)
> @assert z0 == z1
> end
> 
> function check_rand(n=1_000_000)
> for _ in 1:n
> check(rand(UInt128), rand(UInt128), rand(UInt128))
> end
> end
> 
> function benchmark(f=mul_mod)
> @benchmark $f(x, y, m) setup=(x=rand(UInt128); y=rand(UInt128); m=rand(UInt128))
> end
> 
> ```

---

<div class="post-metadata">

**Author:** ![cjdoris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cjdoris/32/213133_2.png) [@cjdoris](https://discourse.julialang.org/u/cjdoris)\
**Post date:** [November 22, 2022, 4:10pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/32 "2022-11-22T16:10:25Z")

</div>

Here’s a 81ns version (and it’s short!).

This version still divides out 64 bits at a time, but directly accumulates the answer in a `UInt256`.

The while loops are O(1) in the average case, but I’m not sure what the worst case is. Maybe they can be improved (I didn’t think long enough, but they do need to be \>1 sometimes).

> **Code**
>
> ```julia
> using BitIntegers: UInt256
> 
> function mul_mod(x::UInt128, y::UInt128, m::UInt128)
> # M is a multiple of m such that the top bit is set
> M = m << leading_zeros(m)
> @assert leading_zeros(M) == 0
> # express M = (m1 << 64) + m0
> m0 = (M % UInt64) % UInt128
> m1 = M >> 64
> @assert M == (m1 << 64) + m0
> # take the product
> z0 = UInt256(x) * UInt256(y)
> # eliminate the top 64 bits of z0
> q1 = fld(UInt128(z0 >> 128), m1)
> z1 = z0 - (UInt256(M) * UInt256(q1) << 64)
> # correct for small overflow
> while z1 > z0
> z1 += (UInt256(M) << 64)
> end
> # eliminate the next 64 bits of z0
> q2 = fld(UInt128(z1 >> 64), m1)
> z2 = z1 - (UInt256(M) * UInt256(q2))
> # correct for small overflow
> while z2 > z1
> z2 += (UInt256(M) << 0)
> end
> # final answer
> return UInt128(z2) % m
> end
> 
> ```

---

<div class="post-metadata">

**Author:** ![Lilith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lilith/32/27492_2.png) [@Lilith](https://discourse.julialang.org/u/Lilith)\
**Post date:** [November 23, 2022, 12:08pm UTC](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421/33 "2022-11-23T12:08:16Z")

</div>

Nice, thank you!

> [@cjdoris](#):
>
> I’m not sure if you can do it in fewer divisions.

It is possible to do this without any divisions beyond the final % operation. `q` should be calculated using floating point multiplication, not integer division. For UInt128s, it should take 3 rounds consisting of Float64 multiplication of `z` by `inv(m1)` to compute `q`, rounding `q` to `Int64`, and multiplying of `q` by `M` to compute what should be subtracted from `z`. Each round should be capable of shaving off at least 50 bits each without any corrective while loops.

[Previous page](https://discourse.julialang.org/t/modular-multiplication-without-overflow/90421.md?page=1)