src/integer/gcd/cl_I_gcd_aux.cc

   1 // partial_gcd().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_I.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_integer.h"
  13 #include "cl_D.h"
  14
  15 void partial_gcd (uintD z1, uintD z2, partial_gcd_result* erg)
  16   { var uintD x1 = 1;
  17     var uintD y1 = 0;
  18     var uintD x2 = 0;
  19     var uintD y2 = 1;
  20     for (;;)
  21       {
  22         // Hier ist z1-y1>=z2+y2.
  23         // Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
  24         { var uintD zaehler = z1-y1;
  25           var uintD nenner = z2+y2; // z2+y2 <= z1-y1 < beta !
  26           if (floor(zaehler,8) >= nenner) // zaehler >= 8*nenner ?
  27             // ja -> Dividieren lohnt sich wohl
  28             { var uintD q = floorD(zaehler,nenner);
  29               x1 += muluD_unchecked(q,x2); // x1 := x1+q*x2
  30               y1 += muluD_unchecked(q,y2); // y1 := y1+q*y2
  31               z1 -= muluD_unchecked(q,z2); // z1 := z1-q*z2
  32             }
  33             else
  34             // nein -> ein paarmal subtrahieren ist wohl schneller
  35             do { x1 += x2; y1 += y2; z1 -= z2; } // (x1,y1,z1) := (x1+x2,y1+y2,z1-z2)
  36                while (z1-y1 >= nenner);
  37         }
  38         if (z2-x2 <= z1+x1-1) break;
  39         // Hier ist z2-x2>=z1+x1.
  40         // Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
  41         { var uintD zaehler = z2-x2;
  42           var uintD nenner = z1+x1; // z1+x1 <= z2-x2 < beta !
  43           if (floor(zaehler,8) >= nenner) // zaehler >= 8*nenner ?
  44             // ja -> Dividieren lohnt sich wohl
  45             { var uintD q = floorD(zaehler,nenner);
  46               x2 += muluD_unchecked(q,x1); // x2 := x2+q*x1
  47               y2 += muluD_unchecked(q,y1); // y2 := y2+q*y1
  48               z2 -= muluD_unchecked(q,z1); // z2 := z2-q*z1
  49             }
  50             else
  51             // nein -> ein paarmal subtrahieren ist wohl schneller
  52             do { x2 += x1; y2 += y1; z2 -= z1; } // (x2,y2,z2) := (x2+x1,y2+y1,z2-z1)
  53                while (z2-x2 >= nenner);
  54         }
  55         if (z1-y1 <= z2+y2-1) break;
  56       }
  57     // Keine Subtraktion mehr möglich.
  58     erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; // Ergebnis
  59   }