src/integer/gcd/cl_low_gcd.cc

   1 // gcd().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/integer.h"
   8
   9
  10 // Implementation.
  11
  12 namespace cln {
  13
  14 // Liefert den ggT zweier Integers.
  15 // gcd(a,b)
  16 // > a,b: zwei Integers
  17 // < ergebnis: (gcd a b), ein Integer >=0
  18   uintV gcd (uintV a, uintV b)
  19 // binäre Methode:
  20 // (gcd a b) :==
  21 //   (prog ((j 0))
  22 //     1 {a,b >0}
  23 //       (when (oddp a) (if (oddp b) (go 2) (go 4)))
  24 //       (when (oddp b) (go 3))
  25 //       (incf j) (setq a (/ a 2)) (setq b (/ b 2))
  26 //       (go 1)
  27 //     2 {a,b >0, beide ungerade}
  28 //       (cond ((> a b) (setq a (- a b)) (go 3))
  29 //             ((= a b) (go 5))
  30 //             ((< a b) (setq b (- b a)) (go 4))
  31 //       )
  32 //     3 {a,b >0, a gerade, b ungerade}
  33 //       (repeat (setq a (/ a 2)) (until (oddp a)))
  34 //       (go 2)
  35 //     4 {a,b >0, a ungerade, b gerade}
  36 //       (repeat (setq b (/ b 2)) (until (oddp b)))
  37 //       (go 2)
  38 //     5 {a=b>0}
  39 //       (return (ash a j))
  40 //   )
  41 // Statt j zu erhöhen und immer Bit 0 von a und b abfragen,
  42 // fragen wir stattdessen immer Bit j von a und b ab; Bits j-1..0 sind =0.
  43 {
  44       #ifdef DUMMER_GGT // so macht's ein Mathematiker:
  45       if (a==0) { return b; }
  46       if (b==0) { return a; }
  47       var uintV bit_j = bit(0);
  48       loop
  49         { // a,b >0
  50           if (!((a & bit_j) ==0))
  51             { if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
  52           if (!((b & bit_j) ==0)) goto even_odd;
  53           // a,b >0 gerade
  54           bit_j = bit_j<<1;
  55         }
  56       #else // Trick von B. Degel:
  57       var uintV bit_j = (a | b); // endet mit einer 1 und j Nullen
  58       bit_j = bit_j ^ (bit_j - 1); // Maske = bit(j) | bit(j-1) | ... | bit(0)
  59       if (!((a & bit_j) ==0))
  60         { if (!((b & bit_j) ==0))
  61             goto odd_odd;
  62           else
  63             if (b==0)
  64               return a;
  65             else
  66               goto odd_even;
  67         }
  68       if (!((b & bit_j) ==0))
  69         { if (a==0)
  70             return b;
  71           else
  72             goto even_odd;
  73         }
  74       return 0; // a=b=0 -> Ergebnis 0
  75       #endif
  76       loop
  77         { odd_odd: // a,b >0, beide ungerade
  78           // Vergleiche a und b:
  79           if (a == b) break; // a=b>0 -> fertig
  80           if (a > b) // a>b ?
  81             { a = a-b;
  82               even_odd: // a,b >0, a gerade, b ungerade
  83               do { a = a>>1; } while ((a & bit_j) ==0);
  84             }
  85             else // a<b
  86             { b = b-a;
  87               odd_even: // a,b >0, a ungerade, b gerade
  88               do { b = b>>1; } while ((b & bit_j) ==0);
  89             }
  90         }
  91       // a=b>0
  92       return a;
  93 }
  94
  95 }  // namespace cln