[perl #37003] [PATCH] Change GCD and LCM opcode to use Binary GCD algorithm

[via RT]

# New Ticket Created by  Imran Ghory 
# Please include the string:  [perl #37003]
# in the subject line of all future correspondence about this issue. 
# <URL: https://rt.perl.org/rt3/Ticket/Display.html?id=37003 >


Presently GCD/LCM use the Euclidlean algorithm to calculate the GCD,
the attached path replaces it with the Binary GCD algorithm. The
primary advantage of the Binary GCD algorithm is that it doesn't use
division like the Euclidlean algorithm (it uses bitwise shifts
instead) and so is more efficient. On a PC it is roughly 25% faster,
on embedded platforms which don't have processor based division
support it can be around 80% faster.

The following code can be used to test both timings and correctness by
running against pre-patch and post-patch code and comparing the
results:

.sub _main
     .local int i,j
     i = 0
     j = 0
loop:
loop2:             
     gcd $I3, i, j
     lcm $I4, i, j
     print i
     print ","
     print j
     print "->("
     print $I3
     print ","
     print $I4
     print ")\n"
     inc j
     if j <= 1000 goto loop2
     j = 0
     inc i 
     if i <= 1000 goto loop
     end
.end

--- clean/parrot/ops/math.ops	2005-08-25 00:30:16.000000000 +0100
+++ parrot/ops/math.ops	2005-08-25 00:37:13.786736920 +0100
@@ -974,20 +974,30 @@
 =cut
 
 inline op gcd(out INT, in INT, in INT) :advanced_math {
-
-  UINTVAL q = 0;
-  UINTVAL c = 0;
-  INTVAL temp2 = $2 < 0 ? -$2 : $2;
-  INTVAL temp3 = $3 < 0 ? -$3 : $3;
-  while (temp3 != 0) {
-    q = (UINTVAL)floor( (FLOATVAL)temp2/temp3 );
-    c = temp2 - temp3*q;
-    temp2 = temp3;
-    temp3 = c;
-  }
-  $1 = temp2;
-  goto NEXT();
+ INTVAL p = 0;
+ INTVAL a = $2 < 0 ? -$2 : $2;
+ INTVAL b = $3 < 0 ? -$3 : $3;
+ 
+ if(a==0) { $1=b; goto NEXT(); }
+ if(b==0) { $1=a; goto NEXT(); }
+ 
+ while( !((a | b) & 1) ) {
+   a>>=1;
+   b>>=1;
+   p++;
+ }
+
+ while(a>0) {
+   if      (!(a & 1)) a>>=1;
+   else if (!(b & 1)) b>>=1;
+   else if (a<b)      b = (b-a)>>1;
+   else               a = (a-b)>>1;
+ }
+ 
+ $1 = b<<p;
+ goto NEXT();
 }
+ 
 
 ########################################
 
@@ -998,26 +1008,31 @@
 =cut
 
 inline op lcm(out INT, in INT, in INT) :advanced_math {
-  UINTVAL q = 0;
-  UINTVAL c = 0;
-  INTVAL temp2 = $2;
-  INTVAL temp3 = $3;
-  INTVAL saved_var2 = temp2;
-  INTVAL saved_var3 = temp3;
-  while (temp3 != 0) {
-    q = (UINTVAL)floor( (FLOATVAL)temp2/temp3 );
-    c = temp2 - temp3*q;
-    temp2 = temp3;
-    temp3 = c;
-  }
-  if (saved_var2 == 0) {
-    $1 = 0;
-  }
-  else {
-    saved_var2 = saved_var2 / temp2;
-    $1 = saved_var2*saved_var3;
-  }
-  goto NEXT();
+ INTVAL gcd = 0;
+ INTVAL p = 0;
+ INTVAL a = $2 < 0 ? -$2 : $2;
+ INTVAL b = $3 < 0 ? -$3 : $3;
+ INTVAL saved_var1 = a, saved_var2 = b; 
+
+ if(a==0 || b==0) { $1=0; goto NEXT(); }
+ 
+ while( !((a | b) & 1) ) {
+   a>>=1;
+   b>>=1;
+   p++;
+ }
+
+ while(a>0) {
+   if      (!(a & 1)) a>>=1;
+   else if (!(b & 1)) b>>=1;
+   else if (a<b)      b = (b-a)>>1;
+   else               a = (a-b)>>1;
+ }
+ 
+ gcd = b<<p;
+ saved_var1 /= gcd;
+ $1 = saved_var1*saved_var2;
+ goto NEXT();
 }
 
 ########################################