src/integer/gcd/cl_I_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 #include "cl_I.h"
  13 #include "cl_DS.h"
  14 #include "cl_D.h"
  15 #include "cl_xmacros.h"
  16
  17 namespace cln {
  18
  19 #define GCD_ALGO 3  // 1: binär, 2: Schulmethode, 3: Lehmer
  20
  21
  22 #if (GCD_ALGO == 1)
  23
  24 // binäre Methode:
  25 // (gcd a b) :==
  26 // b=0 --> (abs a)
  27 // a=0 --> (abs b)
  28 // sonst:
  29 //   (abs a) und (abs b) in zwei Buffer packen, als Unsigned Digit Sequences.
  30 //   [Schreibe oBdA wieder a,b]
  31 //   (prog ((j 0))
  32 //     1 {a,b >0}
  33 //       (if (evenp a)
  34 //         (if (evenp b)
  35 //           (progn (incf j) (setq a (/ a 2)) (setq b (/ b 2)) (go 1))
  36 //           (go 4)
  37 //         )
  38 //         (while (evenp b) (setq b (/ b 2)))
  39 //       )
  40 //     2 {a,b >0, beide ungerade}
  41 //       (cond ((> a b))
  42 //             ((= a b) (go 5))
  43 //             ((< a b) (rotatef a b))
  44 //       )
  45 //     3 {a,b >0, beide ungerade, a>b}
  46 //       (setq a (- a b))
  47 //     4 {a,b >0, a gerade, b ungerade}
  48 //       (repeat (setq a (/ a 2)) (until (oddp a)))
  49 //       (go 2)
  50 //     5 {a=b>0}
  51 //       (return (ash a j))
  52 //   )
  53 // weil es oft auftritt (insbesondere bei GCD's mehrerer Zahlen):
  54 // a=1 oder b=1 --> 1
  55
  56   const cl_I gcd (const cl_I& a, const cl_I& b)
  57     { if (eq(a,1)) { return 1; } // a=1 -> 1
  58       if (eq(b,1)) { return 1; } // b=1 -> 1
  59       if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
  60       if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
  61       CL_ALLOCA_STACK;
  62       var uintD* a_MSDptr;
  63       var uintC a_len;
  64       var uintD* a_LSDptr;
  65       var uintD* b_MSDptr;
  66       var uintC b_len;
  67       var uintD* b_LSDptr;
  68       // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
  69       #define I_abs_to_NUDS(x)  \
  70         I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
  71         if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren:                        */\
  72           { neg_loop_lsp(x##_LSDptr,x##_len); }                                              \
  73         if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen)       */\
  74           { msshrink(x##_MSDptr); x##_len--; }
  75       I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
  76       I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
  77       // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
  78       // beides NUDS, und a_len>0, b_len>0.
  79       // Macro: Halbiere x.
  80       #define halb(x)  \
  81         { shift1right_loop_msp(x##_MSDptr,x##_len,0); /* um 1 Bit rechts schieben */             \
  82           if (mspref(x##_MSDptr,0) == 0) { msshrink(x##_MSDptr); x##_len--; } /* normalisieren */\
  83         }
  84       // Macro: Ob x gerade ist.
  85       #define evenp(x)  \
  86         ((lspref(x##_LSDptr,0) & bit(0)) ==0)
  87       { var uintL j = 0;
  88         label_1: // a,b >0
  89           if (evenp(a))
  90             { if (evenp(b))
  91                 { j++; halb(a); halb(b); goto label_1; }
  92                 else
  93                 goto label_4;
  94             }
  95           while (evenp(b)) { halb(b); }
  96         label_2: // a,b >0, beide ungerade
  97           // Vergleiche a und b:
  98           if (a_len > b_len) goto label_3; // a>b ?
  99           if (a_len == b_len)
 100             { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
 101               if (vergleich > 0) goto label_3; // a>b ?
 102               if (vergleich == 0) goto label_5; // a=b ?
 103             }
 104           // a<b -> a,b vertauschen:
 105           swap(uintD*, a_MSDptr,b_MSDptr);
 106           swap(uintC, a_len,b_len);
 107           swap(uintD*, a_LSDptr,b_LSDptr);
 108         label_3: // a,b >0, beide ungerade, a>b
 109           // subtrahiere a := a - b
 110           if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
 111             { dec_loop_lsp(a_LSDptr lspop b_len,a_len-b_len); }
 112           while (mspref(a_MSDptr,0) == 0) { msshrink(a_MSDptr); a_len--; } // normalisieren
 113         label_4: // a,b >0, a gerade, b ungerade
 114           do { halb(a); } while (evenp(a));
 115           goto label_2;
 116         label_5: // a=b>0
 117           // a zu einer NDS machen:
 118           return ash(NUDS_to_I(a_MSDptr,a_len),j); // ggT der ungeraden Anteile als Integer, mal 2^j
 119       }
 120       #undef evenp
 121       #undef halb
 122       #undef I_abs_to_NUDS
 123     }
 124
 125 #endif /* GCD_ALGO == 1 */
 126
 127
 128 #if (GCD_ALGO == 2)
 129
 130 // Schulmethode:
 131 //   (gcd a b) :==
 132 //   [a:=(abs a), b:=(abs b), while b>0 do (a,b) := (b,(mod a b)), -> a]
 133 // verbessert:
 134 // a=0 -> (abs b)
 135 // b=0 -> (abs a)
 136 // a=1 -> 1
 137 // b=1 -> 1
 138 // a:=(abs a), b:=(abs b)
 139 // Falls a=b: return a; falls a<b: vertausche a und b.
 140 // (*) {Hier a>b>0}
 141 // Falls b=1, return 1. {spart eine Division durch 1}
 142 // Sonst dividieren (divide a b), a:=b, b:=Rest.
 143 //       Falls b=0, return a, sonst goto (*).
 144
 145   const cl_I gcd (const cl_I& a, const cl_I& b)
 146     { if (eq(a,1)) { return 1; } // a=1 -> 1
 147       if (eq(b,1)) { return 1; } // b=1 -> 1
 148       if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
 149       if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
 150       // Beträge nehmen:
 151      {var cl_I abs_a = abs(a);
 152       var cl_I abs_b = abs(b);
 153       var cl_I& a = abs_a;
 154       var cl_I& b = abs_b;
 155       if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums >0
 156         { // bleibt Fixnum, da (gcd a b) <= (min a b)
 157           return L_to_FN(gcd(FN_to_UL(a),FN_to_UL(b)));
 158         }
 159       { var cl_signean vergleich = compare(a,b);
 160         if (vergleich == 0) { return a; } // a=b -> fertig
 161         if (vergleich < 0) { var cl_I tmp = a; a = b; b = a; } // a<b -> a,b vertauschen
 162       }
 163       loop // Hier a > b > 0
 164         { if (eq(b,1)) { return 1; } // b=1 -> Ergebnis 1
 165           { var cl_I_div_t div = cl_divide(a,b); a = b; b = div.remainder; }
 166           if (eq(b,0)) { return a; }
 167         }
 168     }}
 169
 170 #endif /* GCD_ALGO == 2 */
 171
 172
 173 #if (GCD_ALGO == 3)
 174
 175 // Lehmer-Methode:
 176 // vgl. [ D. E. Knuth: The Art of Computer Programming, Vol. 2: Seminumerical
 177 //        Algorithms, Sect. 4.5.2., Algorithm L ]
 178 // und [ Collins, Loos: SAC-2, Algorithms IGCD, DPCC ].
 179 // (gcd a b) :==
 180 // a=0 -> (abs b)
 181 // b=0 -> (abs a)
 182 // a=1 -> 1
 183 // b=1 -> 1
 184 // a:=(abs a), b:=(abs b)
 185 // (*) {Hier a,b>0}
 186 // Falls a=b: return a; falls a<b: vertausche a und b.
 187 // {Hier a>b>0}
 188 // Falls (- (integer-length a) (integer-length b)) >= intDsize/2,
 189 //   lohnt sich eine Division: (a,b) := (b , a mod b). Falls b=0: return a.
 190 // Falls dagegen 0 <= (- (integer-length a) (integer-length b)) < intDsize/2,
 191 //   seien a' die führenden intDsize Bits von a
 192 //   (2^(intDsize-1) <= a' < 2^intDsize) und b' die entsprechenden Bits von b
 193 //   (2^(intDsize/2) <= b' <= a' < 2^intDsize).
 194 //   Rechne den Euklid-Algorithmus mit Beifaktoren für ALLE Zahlen (a,b) aus,
 195 //   die mit a' bzw. b' anfangen; das liefert x1,y1,x2,y2, so daß
 196 //   ggT(a,b) = ggT(x1*a-y1*b,-x2*a+y2*b) und x1*a-y1*b>=0,-x2*a+y2*b>=0.
 197 //   Genauer: Mit offensichtlicher Skalierung betrachten wir
 198 //            a als beliebiges Element des Intervalls [a',a'+1) und
 199 //            b als beliebiges Element des Intervalls [b',b'+1) und
 200 //            führen den Euklid-Algorithmus schrittweise durch:
 201 //            (x1,y1,z1) := (1,0,a'), (x2,y2,z2) := (0,1,b'),
 202 //            Schleife:
 203 //            {Hier x1*a'-y1*b'=z1, x1*a-y1*b in [z1-y1,z1+x1), z1-y1>=0, z1>0,
 204 //             und -x2*a'+y2*b'=z2, -x2*a+y2*b in [z2-x2,z2+y2), z2-x2>=0, z2>0,
 205 //             x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a'.}
 206 //            Falls z1-y1>=z2+y2:
 207 //              (x1,y1,z1) := (x1+x2,y1+y2,z1-z2), goto Schleife.
 208 //            Falls z2-x2>=z1+x1:
 209 //              (x2,y2,z2) := (x2+x1,y2+y1,z2-z1), goto Schleife.
 210 //            Sonst muß man abbrechen.
 211 //            {Zu den Schleifeninvarianten:
 212 //             1. Die Gleichungen x1*a'-y1*b'=z1, -x2*a'+y2*b'=z2,
 213 //                x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a' mit Induktion.
 214 //             2. Die Ungleichungen x1>0, y1>=0, x2>=0, y2>0 mit Induktion.
 215 //             3. Die Ungleichungen z1>=0, z2>=0 nach Fallauswahl.
 216 //             4. Die Ungleichung x1+x2>0 aus x1*z2+x2*z1=b'>0,
 217 //                die Ungleichung y1+y2>0 aus y1*z2+y2*z1=a'>0.
 218 //             5. Die Ungleichung z1>0 wegen Fallauswahl und y1+y2>0,
 219 //                Die Ungleichung z2>0 wegen Fallauswahl und x1+x2>0.
 220 //             6. Die Ungleichungen z1-y1>=0, z2-x2>=0 wegen Fallauswahl.
 221 //             7. Die Ungleichung max(z1,z2) <= a' mit Induktion.
 222 //             8. Die Ungleichung x1+x2 <= x1*z2+x2*z1 = b',
 223 //                die Ungleichung y1+y2 <= y1*z2+y2*z1 = a'.
 224 //             Damit bleiben alle Größen im Intervall [0,beta), kein Überlauf.
 225 //             9. Die Ungleichungen z1+x1<=beta, z2+y2<=beta mit Induktion.
 226 //             10. x1*a-y1*b in (z1-y1,z1+x1) (bzw. [z1,z1+x1) bei y1=0),
 227 //                -x2*a+y2*b in (z2-x2,z2+y2) (bzw. [z2,z2+y2) bei x2=0),
 228 //                da a in a'+[0,1) und b in b'+[0,1).
 229 //                Jedenfalls 0 < x1*a-y1*b < z1+x1 <= x2*z1+x1*z2 = b' falls x2>0,
 230 //                und        0 < -x2*a+y2*b < z2+y2 <= y1*z2+y2*z1 = a' falls y1>0.}
 231 //            Man kann natürlich auch mehrere Subtraktionsschritte auf einmal
 232 //            durchführen:
 233 //            Falls q := floor((z1-y1)/(z2+y2)) > 0 :
 234 //              (x1,y1,z1) := (x1+q*x2,y1+q*y2,z1-q*z2), goto Schleife.
 235 //            Falls q := floor((z2-x2)/(z1+x1)) > 0 :
 236 //              (x2,y2,z2) := (x2+q*x1,y2+q*y1,z2-q*z1), goto Schleife.
 237 //            {Am Schluß gilt -(x1+x2) < z1-z2 < y1+y2 und daher
 238 //             z2-x2 <= b'/(x1+x2) < z1+x1, z1-y1 <= a'/(y1+y2) < z2+y2,
 239 //             und - unter Berücksichtigung von x1*y2-x2*y1=1 -
 240 //             z1-y1 <= b'/(x1+x2) < z2+y2, z2-x2 <= a'/(y1+y2) < z1+x1,
 241 //             also  max(z1-y1,z2-x2) <= min(b'/(x1+x2),a'/(y1+y2))
 242 //                          <= max(b'/(x1+x2),a'/(y1+y2)) < min(z1+x1,z2+y2).}
 243 //   Im Fall y1=x2=0 => x1=y2=1 (der nur bei a'=z1=z2=b' eintreten kann)
 244 //     ersetze (a,b) := (a-b,b). {Beide >0, da a>b>0 war.}
 245 //   Der Fall y1=0,x2>0 => x1=y2=1 => a' = z1 < z2+x2*z1 = b'
 246 //     kann nicht eintreten.
 247 //   Im Fall x2=0,y1>0 => x1=y2=1 ersetze (a,b) := (a-y1*b,b).
 248 //     {Das ist OK, da 0 <= z1-y1 = a'-y1*(b'+1) < a-y1*b < a.}
 249 //   Sonst (y1>0,x2>0) ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
 250 //     {Das ist OK, da 0 <= z1-y1 = x1*a'-y1*(b'+1) < x1*a-y1*b
 251 //                  und 0 <= z2-x2 = -x2*(a'+1)+y2*b' < -x2*a+y2*b
 252 //      und x1*a-y1*b < x1*(a'+1)-y1*b' = z1+x1 <= x2*z1+x1*z2 = b' <= b
 253 //      und -x2*a+y2*b < -x2*a'+y2*(b'+1) = z2+y2 <= y1*z2+y2*z1 = a' <= a.}
 254 // goto (*).
 255
 256   // Define this to 1 in order to use double-word sized a' and b'.
 257   // This gives better x1,y1,x2,y2, because normally the values x1,y1,x2,y2
 258   // have only about intDsize/2 bits and so half of the multiplication work
 259   // is lost. Actually, this flag multiplies the gcd speed by 1.5, not 2.0.
 260   #define DOUBLE_SPEED 1
 261   // Speed increases only for large integers. For small ones, computing
 262   // with double-word sized a' and b' is too costly. The threshold is
 263   // between 12 and 20, around 15.
 264   #define cl_gcd_double_threshold  16
 265
 266   // gcd of two single-word numbers >0
 267   static uintD gcdD (uintD a, uintD b)
 268     { var uintD bit_j = (a | b); // endet mit einer 1 und j Nullen
 269       bit_j = bit_j ^ (bit_j - 1); // Maske = bit(j) | bit(j-1) | ... | bit(0)
 270       if (!((a & bit_j) ==0))
 271         { if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
 272       if (!((b & bit_j) ==0)) goto even_odd;
 273       NOTREACHED;
 274       loop
 275         { odd_odd: // a,b >0, beide ungerade
 276           // Vergleiche a und b:
 277           if (a == b) break; // a=b>0 -> fertig
 278           if (a > b) // a>b ?
 279             { a = a-b;
 280               even_odd: // a,b >0, a gerade, b ungerade
 281               do { a = a>>1; } while ((a & bit_j) ==0);
 282             }
 283             else // a<b
 284             { b = b-a;
 285               odd_even: // a,b >0, a ungerade, b gerade
 286               do { b = b>>1; } while ((b & bit_j) ==0);
 287             }
 288         }
 289       // a=b>0
 290       return a;
 291     }
 292
 293   const cl_I gcd (const cl_I& a, const cl_I& b)
 294     { if (eq(a,1)) { return 1; } // a=1 -> 1
 295       if (eq(b,1)) { return 1; } // b=1 -> 1
 296       if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
 297       if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
 298       if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums /=0
 299         { var sintL a_ = FN_to_L(a);
 300           if (a_ < 0) { a_ = -a_; }
 301           var sintL b_ = FN_to_L(b);
 302           if (b_ < 0) { b_ = -b_; }
 303           return UL_to_I(gcd((uint32)a_,(uint32)b_));
 304         }
 305       CL_ALLOCA_STACK;
 306       var uintD* a_MSDptr;
 307       var uintC a_len;
 308       var uintD* a_LSDptr;
 309       var uintD* b_MSDptr;
 310       var uintC b_len;
 311       var uintD* b_LSDptr;
 312       // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
 313       #define I_abs_to_NUDS(x)  \
 314         I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
 315         if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren:                        */\
 316           { neg_loop_lsp(x##_LSDptr,x##_len); }                                              \
 317         if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen)       */\
 318           { msshrink(x##_MSDptr); x##_len--; }
 319       I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
 320       I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
 321       // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
 322       // beides NUDS, und a_len>0, b_len>0.
 323       // Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
 324       {var uintD* divroomptr; // Platz für Divisionsergebnis
 325        var uintD* c_LSDptr;
 326        var uintD* d_LSDptr;
 327        {var uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
 328         num_stack_alloc(c_len,divroomptr=,c_LSDptr=);
 329         num_stack_alloc(c_len,,d_LSDptr=);
 330         // Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
 331        }
 332        loop
 333          { // Hier a,b>0, beides NUDS.
 334            // Vergleiche a und b:
 335            if (a_len > b_len) goto a_greater_b; // a>b ?
 336            if (a_len == b_len)
 337              { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
 338                if (vergleich > 0) goto a_greater_b; // a>b ?
 339                if (vergleich == 0) break; // a=b ?
 340              }
 341            // a<b -> a,b vertauschen:
 342            swap(uintD*, a_MSDptr,b_MSDptr);
 343            swap(uintC, a_len,b_len);
 344            swap(uintD*, a_LSDptr,b_LSDptr);
 345            a_greater_b:
 346            // Hier a>b>0, beides NUDS.
 347            if (b_len==1) // Beschleunigung eines häufigen Falles
 348              { var uintD b0 = mspref(b_MSDptr,0);
 349                if (b0==1)
 350                   // a>b=1 -> Ergebnis 1.
 351                  { return 1; }
 352                // a>b>1 -> evtl. Division durch b
 353                var uintD a0;
 354                if (a_len==1)
 355                  { a0 = mspref(a_MSDptr,0); }
 356                  else
 357                  { a0 = divu_loop_msp(b0,a_MSDptr,a_len);
 358                    if (a0==0)
 359                      { return UD_to_I(b0); }
 360                  }
 361                return UD_to_I(gcdD(a0,b0));
 362              }
 363            // Entscheidung, ob Division oder Linearkombination:
 364            { var uintD a_msd; // führende intDsize Bits von a
 365              var uintD b_msd; // entsprechende Bits von b
 366              #if DOUBLE_SPEED
 367              var uintD a_nsd; // nächste intDsize Bits von a
 368              var uintD b_nsd; // entsprechende Bits von b
 369              #endif
 370              { var uintC len_diff = a_len-b_len; // Längendifferenz
 371                if (len_diff > 1) goto divide; // >=2 -> Bitlängendifferenz>intDsize -> dividieren
 372                #define bitlendiff_limit  (intDsize/2) // sollte >0,<intDsize sein
 373               {var uintC a_msd_size;
 374                a_msd = mspref(a_MSDptr,0); // führendes Digit von a
 375                integerlengthD(a_msd,a_msd_size=); // dessen Bit-Länge (>0,<=intDsize) berechnen
 376                b_msd = mspref(b_MSDptr,0);
 377                #if HAVE_DD
 378                {var uintDD b_msdd = // 2 führende Digits von b
 379                   (len_diff==0
 380                    ? highlowDD(b_msd, mspref(b_MSDptr,1))
 381                    : (uintDD)b_msd
 382                   );
 383                 // a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
 384                 b_msd = lowD(b_msdd >> a_msd_size);
 385                 if (b_msd < (uintD)bit(intDsize-bitlendiff_limit)) goto divide;
 386                 #if DOUBLE_SPEED
 387                 b_nsd = lowD(highlowDD(lowD(b_msdd), (b_len<=2-len_diff ? 0 : mspref(b_MSDptr,2-len_diff))) >> a_msd_size);
 388                 #endif
 389                }
 390                {var uintDD a_msdd = // 2 führende Digits von a
 391                   highlowDD(a_msd, mspref(a_MSDptr,1));
 392                 a_msd = lowD(a_msdd >> a_msd_size);
 393                 #if DOUBLE_SPEED
 394                 a_nsd = lowD(highlowDD(lowD(a_msdd), (a_len<=2 ? 0 : mspref(a_MSDptr,2))) >> a_msd_size);
 395                 #endif
 396                }
 397                if (a_msd == b_msd) goto subtract;
 398                #else
 399                if (len_diff==0)
 400                  { // a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
 401                    if ((a_msd_size > bitlendiff_limit)
 402                        && (b_msd < (uintD)bit(a_msd_size-bitlendiff_limit))
 403                       )
 404                      goto divide;
 405                    // Entscheidung für Linearkombination ist gefallen.
 406                    // a_msd und b_msd so erweitern, daß a_msd die führenden
 407                    // intDsize Bits von a enthält:
 408                   {var uintC shiftcount = intDsize-a_msd_size; // Shiftcount nach links (>=0, <intDsize)
 409                    if (shiftcount>0)
 410                      { a_msd = a_msd << shiftcount;
 411                        b_msd = b_msd << shiftcount;
 412                        a_msd |= mspref(a_MSDptr,1) >> a_msd_size;
 413                        b_msd |= mspref(b_MSDptr,1) >> a_msd_size;
 414                      }
 415                    if (a_msd == b_msd) goto subtract;
 416                    #if DOUBLE_SPEED
 417                    a_nsd = mspref(a_MSDptr,1);
 418                    b_nsd = mspref(b_MSDptr,1);
 419                    if (shiftcount>0)
 420                      { a_nsd = a_nsd << shiftcount;
 421                        b_nsd = b_nsd << shiftcount;
 422                        if (a_len>2)
 423                          { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
 424                            b_nsd |= mspref(b_MSDptr,2) >> a_msd_size;
 425                      }   }
 426                    #endif
 427                  }}
 428                  else
 429                  // len_diff=1
 430                  { // a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
 431                    if ((a_msd_size >= bitlendiff_limit)
 432                        || (b_msd < (uintD)bit(a_msd_size+intDsize-bitlendiff_limit))
 433                       )
 434                      goto divide;
 435                    // Entscheidung für Linearkombination ist gefallen.
 436                    // a_msd und b_msd so erweitern, daß a_msd die führenden
 437                    // intDsize Bits von a enthält:
 438                    // 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
 439                    a_msd = (a_msd << (intDsize-a_msd_size)) | (mspref(a_MSDptr,1) >> a_msd_size);
 440                    #if DOUBLE_SPEED
 441                    a_nsd = mspref(a_MSDptr,1) << (intDsize-a_msd_size);
 442                    b_nsd = b_msd << (intDsize-a_msd_size);
 443                    a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
 444                    b_nsd |= mspref(b_MSDptr,1) >> a_msd_size;
 445                    #endif
 446                    b_msd = b_msd >> a_msd_size;
 447                  }
 448                #endif
 449                #undef bitlendiff_limit
 450              }}
 451              // Nun ist a_msd = a' > b' = b_msd.
 452              { // Euklid-Algorithmus auf den führenden Digits durchführen:
 453                var partial_gcd_result likobi;
 454                #if DOUBLE_SPEED
 455                if (a_len >= cl_gcd_double_threshold)
 456                  {
 457                    #if HAVE_DD
 458                    partial_gcd(highlowDD(a_msd,a_nsd),highlowDD(b_msd,b_nsd),&likobi); // liefert x1,y1,x2,y2
 459                    #else
 460                    partial_gcd(a_msd,a_nsd,b_msd,b_nsd,&likobi); // liefert x1,y1,x2,y2
 461                    #endif
 462                  }
 463                else
 464                #endif
 465                  { partial_gcd(a_msd,b_msd,&likobi); } // liefert x1,y1,x2,y2, aber nur halb so gut
 466                // Hier y1>0.
 467                if (likobi.x2==0)
 468                  { // Ersetze (a,b) := (a-y1*b,b).
 469                    if (likobi.y1==1) goto subtract; // einfacherer Fall
 470                    // Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
 471                    if (a_len == b_len) { lsprefnext(a_MSDptr) = 0; a_len++; }
 472                    // und y1*b von a subtrahieren:
 473                    mspref(a_MSDptr,0) -= mulusub_loop_lsp(likobi.y1,b_LSDptr,a_LSDptr,b_len);
 474                  }
 475                  else
 476                  { // Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
 477                    // Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
 478                    if (!(a_len==b_len)) { lsprefnext(b_MSDptr) = 0; b_len++; }
 479                    // c := x1*a-y1*b bilden:
 480                    mulu_loop_lsp(likobi.x1,a_LSDptr,c_LSDptr,a_len);
 481                    /* lspref(c_LSDptr,a_len) -= */
 482                      mulusub_loop_lsp(likobi.y1,b_LSDptr,c_LSDptr,a_len);
 483                    // d := -x2*a+y2*b bilden:
 484                    mulu_loop_lsp(likobi.y2,b_LSDptr,d_LSDptr,a_len);
 485                    /* lspref(d_LSDptr,a_len) -= */
 486                      mulusub_loop_lsp(likobi.x2,a_LSDptr,d_LSDptr,a_len);
 487                    // Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
 488                    // lspref(c_LSDptr,a_len) und lspref(d_LSDptr,a_len) =0 sein.
 489                    // a := c und b := d kopieren:
 490                    copy_loop_lsp(c_LSDptr,a_LSDptr,a_len);
 491                    copy_loop_lsp(d_LSDptr,b_LSDptr,a_len);
 492                    // b normalisieren:
 493                    while (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
 494              }   }
 495              if (cl_false)
 496                { subtract: // Ersetze (a,b) := (a-b,b).
 497                  if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
 498                    // Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
 499                    { mspref(a_MSDptr,0) -= 1; }
 500                }
 501              // a normalisieren:
 502              while (mspref(a_MSDptr,0)==0) { msshrink(a_MSDptr); a_len--; }
 503            }
 504            if (cl_false)
 505              { divide: // Ersetze (a,b) := (b , a mod b).
 506               {var uintD* old_a_LSDptr = a_LSDptr;
 507                var DS q;
 508                var DS r;
 509                cl_UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
 510                a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; // a := b
 511                b_len = r.len; if (b_len==0) break; // b=0 -> fertig
 512                b_LSDptr = old_a_LSDptr; // b übernimmt den vorherigen Platz von a
 513                b_MSDptr = copy_loop_lsp(r.LSDptr,b_LSDptr,b_len); // b := r kopieren
 514                goto a_greater_b; // Nun ist a>b>0
 515              }}
 516       }  }
 517       return NUDS_to_I(a_MSDptr,a_len); // NUDS a als Ergebnis
 518       #undef I_abs_to_NUDS
 519     }
 520
 521 #endif /* GCD_ALGO == 3 */
 522
 523 }  // namespace cln