# Multi-limb inverse for mpn_divrem_1

Torbjorn Granlund tege at swox.com
Mon Sep 15 01:41:07 CEST 2003

```We are planning to radically speed mpn_divrem_1, mpn_mod_1, and
friends by allowing for a multi-limb divisor reciprocal.

The problem with division by a single-limb divisor is that the
limb multiply instructions lie on the recurrency path.  Pipelined
multiply units are largely wasted.

By using a n-limb inverse, we can develop n quotient limbs per
iteration, at the cost of O(n^2) multiplies, but without
significantly increased recurrency depth between iterations.  It
will be practical to use n of up to about 4, bigger n result in
too many limb multiplies.

Here is a proof-of-concept implementation for n = 2.
(Karl Hasselström and I wrote this.)

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"

/* # of limbs in inverse	# of full limb mults	# of limb half mults
1                             1+1                     1+2
2                             4+2                     7+3
3                             8+3                    14+5
4                            13+4                    23+7
*/

mp_limb_t
xmpn_divrem_1 (mp_ptr qp, mp_ptr up, mp_size_t un, mp_limb_t vl)
{
mp_limb_t vip;
mp_limb_t tp;
mp_limb_t rp;
mp_limb_t pp;
int i;
mp_limb_t r;

#if GMP_NAIL_BITS != 0
#error not yet for anils
#endif

/* Compute 2-limb inverse.  */
if ((vl << 1) == 0)
{
vip = ~ (mp_limb_t) 0;
vip = ~ (mp_limb_t) 0;
}
else
{
/* For now use mpn_divrem_1.  */
mp_limb_t pow2;
pow2 = 0;  pow2 = 0;  pow2 = -vl;
mpn_divrem_1 (tp, 0, pow2, 3, vl);
vip = tp;
vip = tp;
}

qp += un - 1;
up += un - 1;

if (up >= vl)
{
qp = 1;
rp = up - vl;
}
else
{
qp = 0;
rp = up;
}
qp -= 2;
up -= 2;

for (i = un - 3; i >= 0; i -= 2)
{
rp = rp;
rp = up;
mpn_mul_basecase (tp, rp, 2, vip, 2);
add_ssaaaa (qp, qp, tp, tp, rp, rp);

mpn_mul_1 (pp, qp, 2, vl);
rp = rp;
rp = up;
sub_ddmmss (rp, rp, rp, rp, pp, pp);

if (rp != 0)
{
mpn_incr_u (qp, 1);
mpn_decr_u (rp, vl);
}
if (rp != 0)
{
mpn_incr_u (qp, 1);
mpn_decr_u (rp, vl);
}
if (rp >= vl)
{
mpn_incr_u (qp, 1);
mpn_decr_u (rp, vl);
}
ASSERT (rp == 0);
qp -= 2;
up -= 2;
}
if ((un & 1) == 0)
{
udiv_qrnnd_preinv (qp, r, rp, up, vl, vip);
}
else
{
r = rp;
}
return r;
}
```
{"ft_lang_label":"__label__en","ft_lang_prob":0.5937903,"math_prob":0.99361354,"size":2206,"snap":"2022-27-2022-33","text_gpt3_token_len":924,"char_repetition_ratio":0.123069935,"word_repetition_ratio":0.0691358,"special_character_ratio":0.44424298,"punctuation_ratio":0.20194174,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9816965,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-08-10T21:28:01Z\",\"WARC-Record-ID\":\"<urn:uuid:ff61d3d7-7b4b-4c1e-bf07-d490a79a3ea4>\",\"Content-Length\":\"4962\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:87cbb16e-432f-4960-9030-3fa77c0dd8e8>\",\"WARC-Concurrent-To\":\"<urn:uuid:f638e682-6504-47da-bb2f-f02cd3ef9f12>\",\"WARC-IP-Address\":\"130.242.124.102\",\"WARC-Target-URI\":\"https://gmplib.org/list-archives/gmp-devel/2003-September/000220.html\",\"WARC-Payload-Digest\":\"sha1:N5QDYZXATXTRI5ZFLS6Y6X7COLLXBVJY\",\"WARC-Block-Digest\":\"sha1:JCU5SFCL5NZSFBHYOW7PSNYE4Z5TRJCR\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-33/CC-MAIN-2022-33_segments_1659882571210.98_warc_CC-MAIN-20220810191850-20220810221850-00130.warc.gz\"}"}