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?