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);

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?