| Dirk Laurie on Fri, 24 Aug 2012 13:31:35 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| How to do fast integer division |
I have an application that does lots of integer divisions with
large operands but small quotients. As a model for that,
I handcoded the extended Euclidean algorithm as follows:
{mbezout(a,b)= my(X); X=[1,0,a;0,1,b]; if(a<b, i=2;j=1, j=2;i=1);
while(X[j,3]>0, X[i,] -= truncate(X[i,3]/X[j,3])*X[j,];
i=j; j=3-i);
X[i,]
}
This code is pretty standard, so I was astounded to see how
slow it is. Timing on 10000-digit numbers:
? b1=bezout(a1,a2);
time = 16 ms.
? b3=mybezout(a1,a2);
time = 38,622 ms.
That's a factor of over 2000 slower, not to be explained
by overhead.
Suspecting that "truncate(X[i,3]/X[j,3])" hides a multitude
of sins (starting with reduction to lowest terms), I changed
it to "truncate(X[i,3]*1./X[j,3])" and the time dropped down
to 200ms. There is a tiny chance that this is not exact, so
I did "divrem(X[i,3],X[j,3])[1]" instead and got 221ms.
Is that the closest one can get to a fast integer quotient
in Pari/GP?