src/integer/gcd/cl_I_xgcd.cc

   1 // xgcd().
   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 == 2)
  23
  24 // Schulmethode:
  25 //   (gcd A B) :==
  26 //   [a:=(abs A), b:=(abs B), while b>0 do (a,b) := (b,(mod a b)), -> a]
  27 // verbessert:
  28 // A=1 -> return g=1, (u,v)=(1,0)
  29 // B=1 -> return g=1, (u,v)=(0,1)
  30 // a:=(abs A), ua:=(signum A), va:=0
  31 // b:=(abs B), ub:=0, vb:=(signum B)
  32 // A=0 -> return g=b, (u,v) = (ub,vb)
  33 // B=0 -> return g=a, (u,v) = (ua,va)
  34 // {Stets ua*A+va*B=a, ub*A+vb*B=b, ua*vb-ub*va = +/- 1.}
  35 // Falls a=b: return a,ua,va;
  36 // falls a<b: vertausche a und b, ua und ub, va und vb.
  37 // (*) {Hier a>b>0}
  38 // Falls b=1, return 1,ub,vb. {spart eine Division durch 1}
  39 // Sonst dividieren (divide a b) -> q,r.
  40 //       Falls r=0, return b,ub,vb.
  41 //       a:=b, b := Rest r = a-q*b, (ua,va,ub,vb) := (ub,vb,ua-q*ub,va-q*vb).
  42 //       goto (*).
  43
  44   const cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v)
  45     { if (eq(a,1)) // a=1 -> g=1, (u,v)=(1,0)
  46         { *u = 1; *v = 0; return a; }
  47       if (eq(b,1)) // b=1 -> g=1, (u,v)=(0,1)
  48         { *u = 0; *v = 1; return b; }
  49       // Vorzeichen nehmen:
  50       var cl_I ua = (minusp(a) ? cl_I(-1) : cl_I(1)); // ua := +/- 1
  51       var cl_I va = 0;
  52       var cl_I ub = 0;
  53       var cl_I vb = (minusp(b) ? cl_I(-1) : cl_I(1)); // vb := +/- 1
  54       // Beträge nehmen:
  55      {var cl_I abs_a = abs(a);
  56       var cl_I abs_b = abs(b);
  57       var cl_I& a = abs_a;
  58       var cl_I& b = abs_b;
  59       if (eq(b,0)) // b=0 -> g=a, (u,v) = (ua,va)
  60         { *u = ua; *v = va; return a; }
  61       if (eq(a,0)) // a=0 -> g=b, (u,v) = (ub,vb)
  62         { *u = ub; *v = vb; return b; }
  63       { var cl_signean vergleich = compare(a,b);
  64         if (vergleich == 0) // a=b -> fertig
  65           { *u = ua; *v = va; return a; }
  66         if (vergleich < 0) // a<b -> a,b vertauschen
  67           { swap(cl_I,a,b); swap(cl_I,ua,ub); swap(cl_I,va,vb); }
  68       }
  69       loop // Hier a>b>0
  70         { if (eq(b,1)) // b=1 -> g=b, (u,v) = (ub,vb)
  71             { *u = ub; *v = vb; return b; }
  72           var cl_I_div_t div = cl_divide(a,b); // Division a / b
  73           var cl_I& q = div.quotient;
  74           var cl_I& r = div.remainder;
  75           if (eq(r,0)) // r=0 -> fertig
  76             { *u = ub; *v = vb; return b; }
  77           { var cl_I x = ua-q*ub; ua = ub; ub = x; }
  78           { var cl_I x = va-q*vb; va = vb; vb = x; }
  79           a = b; b = r;
  80         }
  81     }}
  82
  83 #endif /* GCD_ALGO == 2 */
  84
  85
  86 #if (GCD_ALGO == 3)
  87 // (xgcd A B) :==
  88 // wie oben bei (gcd A B).
  89 // Zusätzlich werden Variablen sA,sB,sk,uAa,uBa,uAb,uBb geführt,
  90 // wobei sA,sB,sk Vorzeichen (+/- 1) und uAa,uBa,uAb,uBb Integers >=0 sind mit
  91 //     uAa * sA*A - uBa * sB*B = a,
  92 //   - uAb * sA*A + uBb * sB*B = b,
  93 // ferner  uAa * uBb - uAb * uBa = sk  und daher (Cramersche Regel)
  94 //   uBb * a + uBa * b = sk*sA*A, uAb * a + uAa * b = sk*sB*B.
  95 // Zu Beginn (a,b) := (|A|,|B|), (sA,sB) := ((signum A), (signumB)),
  96 //           (uAa,uBa,uAb,uBb) := (1,0,0,1).
  97 // Beim Ersetzen (a,b) := (a-b,b)
  98 //   ersetzt man (uAa,uBa,uAb,uBb) := (uAa+uAb,uBa+uBb,uAb,uBb).
  99 // Beim Ersetzen (a,b) := (a-y1*b,b)
 100 //   ersetzt man (uAa,uBa,uAb,uBb) := (uAa+y1*uAb,uBa+y1*uBb,uAb,uBb).
 101 // Beim Ersetzen (a,b) := (x1*a-y1*b,-x2*a+y2*b) mit x1*y2-x2*y1=1
 102 //   ersetzt man (uAa,uBa,uAb,uBb) :=
 103 //               (x1*uAa+y1*uAb,x1*uBa+y1*uBb,x2*uAa+y2*uAb,x2*uBa+y2*uBb).
 104 // Beim Ersetzen (a,b) := (b,a)
 105 //   ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa,uBa),
 106 //               sk := -sk, (sA,sB) := (-sA,-sB).
 107 // Beim Ersetzen (a,b) := (b,a-q*b)
 108 //   ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa+q*uAb,uBa+q*uBb),
 109 //               sk := -sk, (sA,sB) := (-sA,-sB).
 110 // Zum Schluß ist a der ggT und a = uAa*sA * A + -uBa*sB * B
 111 // die gewünschte Linearkombination.
 112 // Da stets gilt sk*sA*A = |A|, sk*sB*B = |B|, a>=1, b>=1,
 113 // folgt 0 <= uAa <= |B|, 0 <= uAb <= |B|, 0 <= uBa <= |A|, 0 <= uBb <= |A|.
 114 // Ferner wird sk nie benutzt, braucht also nicht mitgeführt zu werden.
 115
 116   // Define this to 1 in order to use double-word sized a' and b'.
 117   // This gives better x1,y1,x2,y2, because normally the values x1,y1,x2,y2
 118   // have only about intDsize/2 bits and so half of the multiplication work
 119   // is lost. Actually, this flag multiplies the gcd speed by 1.5, not 2.0.
 120   #define DOUBLE_SPEED 1
 121
 122   // Bildet u := u + v, wobei für u genügend Platz sei:
 123   // (Benutzt v.MSDptr nicht.)
 124   static void NUDS_likobi0_NUDS (DS* u, DS* v)
 125     { var uintC u_len = u->len;
 126       var uintC v_len = v->len;
 127       if (u_len >= v_len)
 128         { if (!( addto_loop_lsp(v->LSDptr,u->LSDptr,v_len) ==0))
 129             { if (!( inc_loop_lsp(u->LSDptr lspop v_len,u_len-v_len) ==0))
 130                 { lsprefnext(u->MSDptr) = 1; u->len++; }
 131         }   }
 132         else // u_len <= v_len
 133         { u->MSDptr = copy_loop_lsp(v->LSDptr lspop u_len,u->LSDptr lspop u_len,v_len-u_len);
 134           u->len = v_len;
 135           if (!( addto_loop_lsp(v->LSDptr,u->LSDptr,u_len) ==0))
 136             { if (!( inc_loop_lsp(u->LSDptr lspop u_len,v_len-u_len) ==0))
 137                 { lsprefnext(u->MSDptr) = 1; u->len++; }
 138         }   }
 139     }
 140
 141   // Bildet u := u + q*v, wobei für u genügend Platz sei:
 142   // (Dabei sei nachher u>0.)
 143   static void NUDS_likobi1_NUDS (DS* u, DS* v, uintD q)
 144     { var uintC v_len = v->len;
 145       if (v_len>0) // nur nötig, falls v /=0
 146         { var uintC u_len = u->len;
 147           var uintD carry;
 148           if (u_len <= v_len) // evtl. u vergrößern
 149             { u->MSDptr = clear_loop_lsp(u->MSDptr,v_len-u_len+1);
 150               u->len = u_len = v_len+1;
 151             } // Nun ist u_len > v_len.
 152           carry = muluadd_loop_lsp(q,v->LSDptr,u->LSDptr,v_len);
 153           if (!(carry==0))
 154             { var uintD* ptr = u->LSDptr lspop v_len;
 155               if ((lspref(ptr,0) += carry) < carry)
 156                 { if (!( inc_loop_lsp(ptr lspop 1,u_len-v_len-1) ==0))
 157                     { lsprefnext(u->MSDptr) = 1; u->len++; }
 158             }   }
 159           while (mspref(u->MSDptr,0)==0) { msshrink(u->MSDptr); u->len--; } // normalisieren
 160     }   }
 161
 162   // Bildet (u,v) := (x1*u+y1*v,x2*u+y2*v), wobei für u,v genügend Platz sei:
 163   // (Dabei sei u>0 oder v>0, nachher u>0 und v>0.)
 164   static void NUDS_likobi2_NUDS (DS* u, DS* v, partial_gcd_result* q, uintD* c_LSDptr, uintD* d_LSDptr)
 165     { var uintC u_len = u->len;
 166       var uintC v_len = v->len;
 167       var uintC c_len;
 168       var uintC d_len;
 169       if (u_len >= v_len)
 170         { mulu_loop_lsp(q->x1,u->LSDptr,c_LSDptr,u_len); c_len = u_len+1;
 171           mulu_loop_lsp(q->x2,u->LSDptr,d_LSDptr,u_len); d_len = u_len+1;
 172           if (!(v_len==0))
 173             {{var uintD carry =
 174                 muluadd_loop_lsp(q->y1,v->LSDptr,c_LSDptr,v_len);
 175               if (!(carry==0))
 176                 { var uintD* ptr = c_LSDptr lspop v_len;
 177                   if ((lspref(ptr,0) += carry) < carry)
 178                     { if (!( inc_loop_lsp(ptr lspop 1,u_len-v_len) ==0))
 179                         { lspref(c_LSDptr,c_len) = 1; c_len++; }
 180              }  }   }
 181              {var uintD carry =
 182                 muluadd_loop_lsp(q->y2,v->LSDptr,d_LSDptr,v_len);
 183               if (!(carry==0))
 184                 { var uintD* ptr = d_LSDptr lspop v_len;
 185                   if ((lspref(ptr,0) += carry) < carry)
 186                     { if (!(inc_loop_lsp(ptr lspop 1,u_len-v_len) ==0))
 187                         { lspref(d_LSDptr,d_len) = 1; d_len++; }
 188             }}  }   }
 189         }
 190         else
 191         { mulu_loop_lsp(q->y1,v->LSDptr,c_LSDptr,v_len); c_len = v_len+1;
 192           mulu_loop_lsp(q->y2,v->LSDptr,d_LSDptr,v_len); d_len = v_len+1;
 193           if (!(u_len==0))
 194             {{var uintD carry =
 195                 muluadd_loop_lsp(q->x1,u->LSDptr,c_LSDptr,u_len);
 196               if (!(carry==0))
 197                 { var uintD* ptr = c_LSDptr lspop u_len;
 198                   if ((lspref(ptr,0) += carry) < carry)
 199                     { if (!( inc_loop_lsp(ptr lspop 1,v_len-u_len) ==0))
 200                         { lspref(c_LSDptr,c_len) = 1; c_len++; }
 201              }  }   }
 202              {var uintD carry =
 203                 muluadd_loop_lsp(q->x2,u->LSDptr,d_LSDptr,u_len);
 204               if (!(carry==0))
 205                 { var uintD* ptr = d_LSDptr lspop u_len;
 206                   if ((lspref(ptr,0) += carry) < carry)
 207                     { if (!( inc_loop_lsp(ptr lspop 1,v_len-u_len) ==0))
 208                         { lspref(d_LSDptr,d_len) = 1; d_len++; }
 209             }}  }   }
 210         }
 211       u->MSDptr = copy_loop_lsp(c_LSDptr,u->LSDptr,c_len);
 212       while (mspref(u->MSDptr,0)==0) { msshrink(u->MSDptr); c_len--; }
 213       u->len = c_len;
 214       v->MSDptr = copy_loop_lsp(d_LSDptr,v->LSDptr,d_len);
 215       while (mspref(v->MSDptr,0)==0) { msshrink(v->MSDptr); d_len--; }
 216       v->len = d_len;
 217     }
 218
 219   // Los geht's:
 220   const cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v)
 221     { if (eq(a,1)) // a=1 -> g=1, (u,v)=(1,0)
 222         { *u = 1; *v = 0; return a; }
 223       if (eq(b,1)) // b=1 -> g=1, (u,v)=(0,1)
 224         { *u = 0; *v = 1; return b; }
 225       var sintL sA = (minusp(a) ? ~0 : 0); // Vorzeichen von A
 226       var sintL sB = (minusp(b) ? ~0 : 0); // Vorzeichen von B
 227       CL_ALLOCA_STACK;
 228       var uintD* a_MSDptr;
 229       var uintC a_len;
 230       var uintD* a_LSDptr;
 231       var uintD* b_MSDptr;
 232       var uintC b_len;
 233       var uintD* b_LSDptr;
 234       // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
 235       #define I_abs_to_NUDS(x,zero_statement)  \
 236         I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
 237         if (x##_len == 0) { zero_statement } /* falls =0, fertig      */\
 238         if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren:   */\
 239           { neg_loop_lsp(x##_LSDptr,x##_len); }                         \
 240         if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
 241           { msshrink(x##_MSDptr); x##_len--; }
 242       I_abs_to_NUDS(a, // (abs A) als NUDS erzeugen
 243                     // A=0 -> g=|B|, (u,v) = (0,sB)
 244                     { *u = 0; *v = (sB==0 ? cl_I(1) : cl_I(-1));
 245                       return abs(b);
 246                     });
 247       I_abs_to_NUDS(b, // (abs B) als NUDS erzeugen
 248                     // B=0 -> g=|A|, (u,v) = (sA,0)
 249                     { *u = (sA==0 ? cl_I(1) : cl_I(-1)); *v = 0;
 250                       return abs(a);
 251                     });
 252       // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
 253       // beides NUDS, und a_len>0, b_len>0.
 254       {// Beifaktoren:
 255        var DS uAa;
 256        var DS uBa;
 257        var DS uAb;
 258        var DS uBb;
 259        // Rechenregister:
 260        var uintD* divroomptr; // Platz für Divisionsergebnis
 261        var uintD* c_LSDptr;
 262        var uintD* d_LSDptr;
 263        // Platz für uAa,uBa,uAb,uBb besorgen:
 264        {var uintC u_len = b_len+1;
 265         num_stack_alloc(u_len,,uAa.LSDptr=); uAa.MSDptr = uAa.LSDptr;
 266         num_stack_alloc(u_len,,uAb.LSDptr=); uAb.MSDptr = uAb.LSDptr;
 267        }
 268        {var uintC u_len = a_len+1;
 269         num_stack_alloc(u_len,,uBa.LSDptr=); uBa.MSDptr = uBa.LSDptr;
 270         num_stack_alloc(u_len,,uBb.LSDptr=); uBb.MSDptr = uBb.LSDptr;
 271        }
 272        lsprefnext(uAa.MSDptr) = 1; uAa.len = 1; // uAa := 1
 273        uBa.len = 0; // uBa := 0
 274        uAb.len = 0; // uAb := 0
 275        lsprefnext(uBb.MSDptr) = 1; uBb.len = 1; // uBb := 1
 276        // Jetzt ist uAa = uAa.MSDptr/uAa.len/uAa.LSDptr,
 277        //           uBa = uBa.MSDptr/uBa.len/uBa.LSDptr,
 278        //           uAb = uAb.MSDptr/uAb.len/uAb.LSDptr,
 279        //           uBb = uBb.MSDptr/uBb.len/uBb.LSDptr,
 280        // alles NUDS.
 281        // Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
 282        {var uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
 283         num_stack_alloc(c_len,,c_LSDptr=);
 284         num_stack_alloc(c_len,divroomptr=,d_LSDptr=);
 285         // Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
 286        }
 287        loop
 288          { // Hier a,b>0, beides NUDS.
 289            // Vergleiche a und b:
 290            if (a_len > b_len) goto a_greater_b; // a>b ?
 291            if (a_len == b_len)
 292              { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
 293                if (vergleich > 0) goto a_greater_b; // a>b ?
 294                if (vergleich == 0) break; // a=b ?
 295              }
 296            // a<b -> a,b vertauschen:
 297            swap(uintD*, a_MSDptr,b_MSDptr);
 298            swap(uintC, a_len,b_len);
 299            swap(uintD*, a_LSDptr,b_LSDptr);
 300            a_greater_b_swap:
 301            swap(DS, uAa,uAb); // und uAa und uAb vertauschen
 302            swap(DS, uBa,uBb); // und uBa und uBb vertauschen
 303            sA = ~sA; sB = ~sB; // und sA und sB umdrehen
 304            a_greater_b:
 305            // Hier a>b>0, beides NUDS.
 306            // Entscheidung, ob Division oder Linearkombination:
 307            { var uintD a_msd; // führende intDsize Bits von a
 308              var uintD b_msd; // entsprechende Bits von b
 309              #if DOUBLE_SPEED
 310              var uintD a_nsd; // nächste intDsize Bits von a
 311              var uintD b_nsd; // entsprechende Bits von b
 312              #endif
 313              { var uintC len_diff = a_len-b_len; // Längendifferenz
 314                if (len_diff > 1) goto divide; // >=2 -> Bitlängendifferenz>intDsize -> dividieren
 315                #define bitlendiff_limit  (intDsize/2) // sollte >0,<intDsize sein
 316               {var uintC a_msd_size;
 317                a_msd = mspref(a_MSDptr,0); // führendes Digit von a
 318                integerlengthD(a_msd,a_msd_size=); // dessen Bit-Länge (>0,<=intDsize) berechnen
 319                b_msd = mspref(b_MSDptr,0);
 320                #if HAVE_DD
 321                {var uintDD b_msdd = // 2 führende Digits von b
 322                   (len_diff==0
 323                    ? highlowDD(b_msd, (b_len==1 ? 0 : mspref(b_MSDptr,1)))
 324                    : (uintDD)b_msd
 325                   );
 326                 // a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
 327                 b_msd = lowD(b_msdd >> a_msd_size);
 328                 if (b_msd < (uintD)bit(intDsize-bitlendiff_limit)) goto divide;
 329                 #if DOUBLE_SPEED
 330                 b_nsd = lowD(highlowDD(lowD(b_msdd), (b_len<=2-len_diff ? 0 : mspref(b_MSDptr,2-len_diff))) >> a_msd_size);
 331                 #endif
 332                }
 333                {var uintDD a_msdd = // 2 führende Digits von a
 334                   highlowDD(a_msd, (a_len==1 ? 0 : mspref(a_MSDptr,1)));
 335                 a_msd = lowD(a_msdd >> a_msd_size);
 336                 #if DOUBLE_SPEED
 337                 a_nsd = lowD(highlowDD(lowD(a_msdd), (a_len<=2 ? 0 : mspref(a_MSDptr,2))) >> a_msd_size);
 338                 #endif
 339                }
 340                if (a_msd == b_msd) goto subtract;
 341                #else
 342                if (len_diff==0)
 343                  { // a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
 344                    if ((a_msd_size > bitlendiff_limit)
 345                        && (b_msd < (uintD)bit(a_msd_size-bitlendiff_limit))
 346                       )
 347                      goto divide;
 348                    // Entscheidung für Linearkombination ist gefallen.
 349                    // a_msd und b_msd so erweitern, daß a_msd die führenden
 350                    // intDsize Bits von a enthält:
 351                   {var uintC shiftcount = intDsize-a_msd_size; // Shiftcount nach links (>=0, <intDsize)
 352                    if (shiftcount>0)
 353                      { a_msd = a_msd << shiftcount;
 354                        b_msd = b_msd << shiftcount;
 355                        if (a_len>1)
 356                          { a_msd |= mspref(a_MSDptr,1) >> a_msd_size;
 357                            b_msd |= mspref(b_MSDptr,1) >> a_msd_size;
 358                      }   }
 359                    if (a_msd == b_msd) goto subtract;
 360                    #if DOUBLE_SPEED
 361                    if (a_len>1)
 362                      { a_nsd = mspref(a_MSDptr,1);
 363                        b_nsd = mspref(b_MSDptr,1);
 364                        if (shiftcount>0)
 365                          { a_nsd = a_nsd << shiftcount;
 366                            b_nsd = b_nsd << shiftcount;
 367                            if (a_len>2)
 368                              { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
 369                                b_nsd |= mspref(b_MSDptr,2) >> a_msd_size;
 370                      }   }   }
 371                    else
 372                      { a_nsd = 0; b_nsd = 0; }
 373                    #endif
 374                  }}
 375                  else
 376                  // len_diff=1
 377                  { // a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
 378                    if ((a_msd_size >= bitlendiff_limit)
 379                        || (b_msd < (uintD)bit(a_msd_size+intDsize-bitlendiff_limit))
 380                       )
 381                      goto divide;
 382                    // Entscheidung für Linearkombination ist gefallen.
 383                    // a_msd und b_msd so erweitern, daß a_msd die führenden
 384                    // intDsize Bits von a enthält:
 385                    // 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
 386                    a_msd = (a_msd << (intDsize-a_msd_size)) | (mspref(a_MSDptr,1) >> a_msd_size);
 387                    #if DOUBLE_SPEED
 388                    a_nsd = mspref(a_MSDptr,1) << (intDsize-a_msd_size);
 389                    b_nsd = b_msd << (intDsize-a_msd_size);
 390                    if (a_len>2)
 391                      { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
 392                        b_nsd |= mspref(b_MSDptr,1) >> a_msd_size;
 393                      }
 394                    #endif
 395                    b_msd = b_msd >> a_msd_size;
 396                  }
 397                #endif
 398                #undef bitlendiff_limit
 399              }}
 400              // Nun ist a_msd = a' > b' = b_msd.
 401              { // Euklid-Algorithmus auf den führenden Digits durchführen:
 402                var partial_gcd_result likobi;
 403                #if DOUBLE_SPEED
 404                #if HAVE_DD
 405                partial_gcd(highlowDD(a_msd,a_nsd),highlowDD(b_msd,b_nsd),&likobi); // liefert x1,y1,x2,y2
 406                #else
 407                partial_gcd(a_msd,a_nsd,b_msd,b_nsd,&likobi); // liefert x1,y1,x2,y2
 408                #endif
 409                #else
 410                partial_gcd(a_msd,b_msd,&likobi); // liefert x1,y1,x2,y2, aber nur halb so gut
 411                #endif
 412                // Hier y1>0.
 413                if (likobi.x2==0)
 414                  { // Ersetze (a,b) := (a-y1*b,b).
 415                    if (likobi.y1==1) goto subtract; // einfacherer Fall
 416                    // Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
 417                    if (a_len == b_len) { lsprefnext(a_MSDptr) = 0; a_len++; }
 418                    // und y1*b von a subtrahieren:
 419                    mspref(a_MSDptr,0) -= mulusub_loop_lsp(likobi.y1,b_LSDptr,a_LSDptr,b_len);
 420                    NUDS_likobi1_NUDS(&uAa,&uAb,likobi.y1); // uAa := uAa + y1 * uAb
 421                    NUDS_likobi1_NUDS(&uBa,&uBb,likobi.y1); // uBa := uBa + y1 * uBb
 422                  }
 423                  else
 424                  { // Ersetze (uAa,uAb) := (x1*uAa+y1*uAb,x2*uAa+y2*uAb) :
 425                    NUDS_likobi2_NUDS(&uAa,&uAb,&likobi,c_LSDptr,d_LSDptr);
 426                    // Ersetze (uBa,uBb) := (x1*uBa+y1*uBb,x2*uBa+y2*uBb) :
 427                    NUDS_likobi2_NUDS(&uBa,&uBb,&likobi,c_LSDptr,d_LSDptr);
 428                    // Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
 429                    // Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
 430                    if (!(a_len==b_len)) { lsprefnext(b_MSDptr) = 0; b_len++; }
 431                    // c := x1*a-y1*b bilden:
 432                    mulu_loop_lsp(likobi.x1,a_LSDptr,c_LSDptr,a_len);
 433                    /* lspref(c_LSDptr,a_len) -= */
 434                      mulusub_loop_lsp(likobi.y1,b_LSDptr,c_LSDptr,a_len);
 435                    // d := -x2*a+y2*b bilden:
 436                    mulu_loop_lsp(likobi.y2,b_LSDptr,d_LSDptr,a_len);
 437                    /* lspref(d_LSDptr,a_len) -= */
 438                      mulusub_loop_lsp(likobi.x2,a_LSDptr,d_LSDptr,a_len);
 439                    // Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
 440                    // lspref(c_LSDptr,a_len) und lspref(d_LSDptr,a_len) =0 sein.
 441                    // a := c und b := d kopieren:
 442                    copy_loop_lsp(c_LSDptr,a_LSDptr,a_len);
 443                    copy_loop_lsp(d_LSDptr,b_LSDptr,a_len);
 444                    // b normalisieren:
 445                    while (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
 446              }   }
 447              if (cl_false)
 448                { subtract: // Ersetze (a,b) := (a-b,b).
 449                  NUDS_likobi0_NUDS(&uAa,&uAb); // uAa := uAa + uAb
 450                  NUDS_likobi0_NUDS(&uBa,&uBb); // uBa := uBa + uBb
 451                  if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
 452                    // Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
 453                    { mspref(a_MSDptr,0) -= 1; }
 454                }
 455              // a normalisieren:
 456              while (mspref(a_MSDptr,0)==0) { msshrink(a_MSDptr); a_len--; }
 457            }
 458            if (cl_false)
 459              { divide: // Ersetze (a,b) := (b , a mod b).
 460               {var uintD* old_a_LSDptr = a_LSDptr;
 461                var DS q;
 462                var DS r;
 463                cl_UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
 464                a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; // a := b
 465                b_len = r.len; if (b_len==0) goto return_a_coeffsb; // b=0 -> fertig
 466                b_LSDptr = old_a_LSDptr; // b übernimmt den vorherigen Platz von a
 467                b_MSDptr = copy_loop_lsp(r.LSDptr,b_LSDptr,b_len); // b := r kopieren
 468                // (uAa,uAb) := (uAb,uAa+q*uAb) :
 469                if (!(uAb.len==0))
 470                  { cl_UDS_mul(q.LSDptr,q.len,uAb.LSDptr,uAb.len,c_LSDptr); // q * uAb
 471                    var DS c;
 472                    c.LSDptr = c_LSDptr; c.len = q.len + uAb.len;
 473                    if (lspref(c_LSDptr,c.len-1)==0) { c.len--; } // normalisieren
 474                    NUDS_likobi0_NUDS(&uAa,&c); // zu uAa addieren
 475                  } // noch uAa,uAb vertauschen (später)
 476                // (uBa,uBb) := (uBb,uBa+q*uBb) :
 477                if (!(uBb.len==0))
 478                  { cl_UDS_mul(q.LSDptr,q.len,uBb.LSDptr,uBb.len,c_LSDptr); // q * uBb
 479                    var DS c;
 480                    c.LSDptr = c_LSDptr; c.len = q.len + uBb.len;
 481                    if (lspref(c_LSDptr,c.len-1)==0) { c.len--; } // normalisieren
 482                    NUDS_likobi0_NUDS(&uBa,&c); // zu uBa addieren
 483                  } // noch uBa,uBb vertauschen (später)
 484                goto a_greater_b_swap; // Nun ist a>b>0
 485              }}
 486          }
 487        // Nun ist a = b. Wähle diejenige der beiden Linearkombinationen
 488        //   a =  uAa*sA * A + -uBa*sB * B
 489        //   b = -uAb*sA * A +  uBb*sB * B
 490        // die die betragsmäßig kleinsten Koeffizienten hat.
 491        // Teste auf uBa < uBb. (Das kann auftreten, z.B. bei
 492        // A=560014183, B=312839871 wird a=b=1, uAa < uAb, uBa < uBb.)
 493        // Falls uBa = uBb, teste auf uAa < uAb. (Das kann auftreten, z.B. bei
 494        // A=2, B=3 wird a=b=1, uAa < uAb, uBa = uBb.)
 495        if (uBb.len > uBa.len) goto return_a_coeffsa;
 496        if (uBb.len < uBa.len) goto return_a_coeffsb;
 497        // (uBb.len == uBa.len)
 498        { var cl_signean vergleich = compare_loop_msp(uBb.MSDptr,uBa.MSDptr,uBb.len);
 499          if (vergleich > 0) goto return_a_coeffsa;
 500          if (vergleich < 0) goto return_a_coeffsb;
 501        }
 502        if (uAb.len > uAa.len) goto return_a_coeffsa;
 503        if (uAb.len < uAa.len) goto return_a_coeffsb;
 504        // (uAb.len == uAa.len)
 505        if (compare_loop_msp(uAb.MSDptr,uAa.MSDptr,uAb.len) > 0)
 506          return_a_coeffsa:
 507          { // uAa mit Vorfaktor sA versehen:
 508            lsprefnext(uAa.MSDptr) = 0; uAa.len++;
 509            if (!(sA==0)) { neg_loop_lsp(uAa.LSDptr,uAa.len); }
 510            // uBa mit Vorfaktor -sB versehen:
 511            lsprefnext(uBa.MSDptr) = 0; uBa.len++;
 512            if (sB==0) { neg_loop_lsp(uBa.LSDptr,uBa.len); }
 513            *u = DS_to_I(uAa.MSDptr,uAa.len); // DS uAa als Vorfaktor von A
 514            *v = DS_to_I(uBa.MSDptr,uBa.len); // DS uBa als Vorfaktor von B
 515          }
 516          else
 517          return_a_coeffsb:
 518          { // uAb mit Vorfaktor -sA versehen:
 519            lsprefnext(uAb.MSDptr) = 0; uAb.len++;
 520            if (sA==0) { neg_loop_lsp(uAb.LSDptr,uAb.len); }
 521            // uBb mit Vorfaktor sB versehen:
 522            lsprefnext(uBb.MSDptr) = 0; uBb.len++;
 523            if (!(sB==0)) { neg_loop_lsp(uBb.LSDptr,uBb.len); }
 524            *u = DS_to_I(uAb.MSDptr,uAb.len); // DS uAb als Vorfaktor von A
 525            *v = DS_to_I(uBb.MSDptr,uBb.len); // DS uBb als Vorfaktor von B
 526          }
 527       }
 528       return NUDS_to_I(a_MSDptr,a_len); // NUDS a als ggT
 529       #undef I_abs_to_NUDS
 530     }
 531
 532 #endif /* GCD_ALGO == 3 */
 533
 534 }  // namespace cln