ni...@lysator.liu.se (Niels Möller) writes:
> Below is one variant, that seems to pass tests. I haven't done any
> benchmarks yet.
Just gave it a try with tune/speed, together with all the latest div1
variants. This one on my laptop, and not looking that great:
$ ./tune/speed -c -s1 -p100000 mpn_hgcd2_1 mpn_hgcd2_2 mpn_hgcd2_3 mpn_hgcd2_4
mpn_hgcd2_5 mpn_hgcd2_binary
overhead 6.02 cycles, precision 100000 units of 8.33e-10 secs, CPU freq 1200.00
MHz
mpn_hgcd2_1 mpn_hgcd2_2 mpn_hgcd2_3 mpn_hgcd2_4 mpn_hgcd2_5
mpn_hgcd2_binary
1 #1668.90 1863.72 1670.73 1757.54 1738.50
2044.25
Had a look at the disassembly for the binary algorithm. The
double-precision loop needs, 20 instructions for just the conditional
swap logic, 23 for the clz + shift + subtract, 8 for the shift+add
updates of the u matrix.
By doing another conditional subtraction in each iteration, I get down
to 1865 cycles. See below version. In the case acnt == bcnt, the
condition for this is always false. But when acnt < bcnt, we have
2^{k-1} <= a < 2^k
2^{k-2} <= (2^s b) < 2^{k-1}
so floor (a / (b 2^s) can be 1, 2, or 3. So one could consider yet
another conditional subtraction (but then q = 3 2^s is no longer just a
single shift and add).
If I also try the SIMD trick, I get down to 1790 cycles, but I'd need a
special middle iteration to make it correct. As long as both a,b >= 2^96
(more than one and a half limb), matrix coefficients fit in a half limb.
But we don't switch to single precision until a,b < 2^96. So we'll have
one or more iterations with a >= 2^96 > b, and there u01 and u11 may
exceeed half a limb. May still be worth it to add that middle loop, for
5% speedup.
Seems hard to beat method 1 on this hardware, but binary may well be a
winner on machines with slow multiply, slow division, but fast
count_leading_zeros.
Regards,
/Niels
#define SWAP_MASK(mask, x, y) do { \
mp_limb_t swap_xor = ((x) ^ (y)) & mask; \
(x) ^= swap_xor; \
(y) ^= swap_xor; \
} while (0);
int
mpn_hgcd2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
struct hgcd_matrix1 *M)
{
mp_limb_t u00, u01, u10, u11, swapped;
if (ah < 2 || bh < 2)
return 0;
if (ah > bh || (ah == bh && al > bl))
{
sub_ddmmss (ah, al, ah, al, bh, bl);
if (ah < 2)
return 0;
u00 = u01 = u11 = 1;
u10 = 0;
}
else
{
sub_ddmmss (bh, bl, bh, bl, ah, al);
if (bh < 2)
return 0;
u00 = u10 = u11 = 1;
u01 = 0;
}
swapped = 0;
for (;;)
{
mp_limb_t mask, dh, dl;
int acnt, bcnt, shift;
if (ah == bh)
goto done;
mask = -(mp_limb_t) (ah < bh);
swapped ^= mask;
SWAP_MASK (mask, ah, bh);
SWAP_MASK (mask, al, bl);
SWAP_MASK (mask, u00, u01);
SWAP_MASK (mask, u10, u11);
ASSERT (ah > bh);
if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
{
ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
break;
}
count_leading_zeros (acnt, ah);
count_leading_zeros (bcnt, bh);
shift = bcnt - acnt;
ASSERT (shift >= 0);
/* Avoid subtraction underflow */
shift -= (shift > 0);
dh = (bh << shift) | ((bl >> 1) >> (GMP_LIMB_BITS - 1 - shift));
dl = bl << shift;
sub_ddmmss (ah, al, ah, al, dh, dl);
mask = -(mp_limb_t) (ah > dh);
shift -= mask;
sub_ddmmss (ah, al, ah, al, mask & dh, mask & dl);
if (ah < 2)
{
/* A is too small, so decrement q. */
mp_limb_t q = ((mp_limb_t) 1 << shift) - 1;
u01 += q * u00;
u11 += q * u10;
goto done;
}
u01 += (u00 << shift);
u11 += (u10 << shift);
}
/* Single-precision loop */
for (;;)
{
mp_limb_t mask, dh;
int acnt, bcnt, shift;
if (ah == bh)
break;
mask = -(mp_limb_t) (ah < bh);
swapped ^= mask;
SWAP_MASK (mask, ah, bh);
SWAP_MASK (mask, u00, u01);
SWAP_MASK (mask, u10, u11);
ASSERT (ah > bh);
count_leading_zeros (acnt, ah);
count_leading_zeros (bcnt, bh);
shift = bcnt - acnt;
ASSERT (shift >= 0);
/* Avoid subtraction underflow */
shift -= (shift > 0);
dh = bh << shift;
ah -= dh;
mask = - (mp_limb_t) (ah > dh);
shift -= mask;
ah -= mask & dh;
if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
{
/* A is too small, so decrement q. */
mp_limb_t q = ((mp_limb_t) 1 << shift) - 1;
u01 += q * u00;
u11 += q * u10;
break;
}
u01 += (u00 << shift);
u11 += (u10 << shift);
}
done:
SWAP_MASK (swapped, u00, u01);
SWAP_MASK (swapped, u10, u11);
ASSERT (u00 * u11 - u01 * u10 == 1);
M->u[0][0] = u00; M->u[0][1] = u01;
M->u[1][0] = u10; M->u[1][1] = u11;
return 1;
}
--
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
