Feedback-ID: ic6394509:Fastmail
User-Agent: Cyrus-JMAP/3.11.0-alpha0-515-g87b2bad5a-fm-20240604.001-g87b2bad5
MIME-Version: 1.0
Message-Id: <d87545fa-5b98-4474-a1ab-5af1ed78fb4d@app.fastmail.com>
In-Reply-To: <f2f123a6-4b63-4932-9ccb-1474625e6285@posteo.net>
References: <7f95163f-2019-4416-a042-6e2141619e5d@app.fastmail.com>
 <f2f123a6-4b63-4932-9ccb-1474625e6285@posteo.net>
Date: Thu, 13 Jun 2024 11:08:25 +0200
From: "Joel Jacobson" <joel@compiler.org>
To: "Aaron Altman" <aaronaltman@posteo.net>,
 pgsql-hackers <pgsql-hackers@postgresql.org>
Subject: Re: Optimize numeric.c mul_var() using the Karatsuba algorithm
Content-Type: text/plain;charset=utf-8
Content-Transfer-Encoding: quoted-printable
Archived-At: <https://www.postgresql.org/message-id/d87545fa-5b98-4474-a1ab-5af1ed78fb4d%40app.fastmail.com>
Precedence: bulk

On Tue, Jun 11, 2024, at 19:16, Aaron Altman wrote:
> Hi Joel, thanks for posting this.=C2=A0 Although I have only a cursory=20
> familiarity with fast multiplication algorithms, I'd like to try and=20
> give it a review.=C2=A0 To start with, can you help me understand the =
choice=20
> of this algorithm versus others like Toom?=C2=A0 If this looks correct=
 on a=20
> closer view I'll propose it for inclusion. Along the way though I'd li=
ke=20
> to have it explicitly called out whether this is superior in general t=
o=20
> other choices, better for more realistic use cases, simpler, clearer t=
o=20
> license or something similar.=C2=A0 It would be nice for future dicuss=
ions to=20
> have some context around whether it would make sense to have condition=
s=20
> to choose other algorithms as well, or if this one is generally the be=
st=20
> for what Postgres users are usually doing.
>
> Continuing with code review in any case.=C2=A0 Interested to hear more.

Hi Aaron, thanks for looking at this!

The choice of best algorithm depends on the factor sizes.

The larger factor sizes, the more complicated algorithm can be "afforded=
".

List of fast multiplication algorithms,
ordered by factor sizes they are suitable for:

- Long multiplication, aka "schoolbook" multiplication.
- Karatsuba
- Toom-3
- Sch=C3=B6nhage=E2=80=93Strassen algorithm (Fast Fourier Transform)

The Toom-3 algorithm can be modified to split the smaller and larger fac=
tors
into different number of parts. The notation used at Wikipedia is e.g. T=
oom-2.5
which I think means splitting the larger into three parts, and the small=
er into
two parts, while GMP uses Toom32 to mean the same thing.
Personally, I think GMPs notation is easier to understand as the number =
of parts
can be directly derived from the name.

I experimented with implementing Toom-3 as well, but there was only a ma=
rginal
win, at very large factor sizes, and since numeric's max ndigits
(number of base-digits) is capped at 32768, I didn't think it was worth =
it,
since it adds quite a lot of complexity.

The Karatsuba algorithm is the next step in the hierarchy of fast multip=
lication
algorithms, and all other bigint libs I've looked at implement Karatsuba,
even if they also implement Toom-3, since Karatsuba is faster than Toom-=
3 for
sufficiently small factors, but that are at the same time sufficiently l=
arge for
Karatsuba to be faster than schoolbook.

I was initially surprised by the quite large threshold, where Karatsuba =
started
to be a win over schoolbook.

I think the explanation why mul_var() stays fast up to quite remarkably =
high
factor sizes, could be a combination of several things, such as:

- mul_var() is already heavily optimized, with clever tricks,
  such as deferred carry propagation.

- numeric uses NBASE=3D10000, while other bigint libs usually use a powe=
r of two.

In the Karatsuba implementation, I tried to keep the KARATSUBA_CONDITION=
()
quite simple, but it's way more complex than what most bigint libs use,
that usually just check if the smaller factor is smaller than some thres=
hold,
and if so, use schoolbook. For instance, this is what Rust's num-bigint =
does:

    if x.len() <=3D 32 {
        // Long multiplication
    } else if x.len() * 2 <=3D y.len() {
        // Half-Karatsuba, for factors with significant length disparity
    } else if x.len() <=3D 256 {
        // Karatsuba multiplication
    } else {
        // Toom-3 multiplication
    }

Source: https://github.com/rust-num/num-bigint/blob/master/src/biguint/m=
ultiplication.rs#L101

Side note: When working on Karatsuba in mul_var(), I looked at some othe=
r bigint
implementations, to try to understand their threshold functions.
I noticed that Rust's num-bigint didn't optimise for factors with signif=
icant
length disparity, so I contributed a patch based on my "Half-Karatsuba" =
idea,
that I got when working with mul_var(), which has now been merged:
https://github.com/rust-num/num-bigint/commit/06b61c8138ad8a9959ac54d977=
3d0a9ebe25b346

In mul_var(), if we don't like the complexity of KARATSUBA_CONDITION(),
we could go for a more traditional threshold approach, i.e. just checking
the smaller factor size. However, I believe that would be at the expense
of missing out of some performance gains.

I've tried quite hard to find the best KARATSUBA_CONDITION(), but I foun=
d it to
be a really hard problem, the differences between different CPU architec=
tures,
in combination with wanting a simple expression, means there is no obvio=
us
perfect threshold function, there will always be a trade-off.

I eventually stopped trying to improve it, and just settled on the versi=
on in
the patch, and thought that I'll leave it up to the community to give fe=
edback
on what complexity for the threshold function is motivated. If we absolu=
tely
just want to check the smallest factor size, like Rust, then it's super =
simple,
then the threshold can easily be found just by testing different values.
It's when both factor sizes are input to the threshold function that mak=
es it
complicated.

/Joel