src/integer/gcd/cl_I_gcd_aux2.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 "cln/integer.h"
  13 #include "cl_D.h"
  14
  15 namespace cln {
  16
  17 // Dasselbe wie partial_gcd(z1,z2,erg), nur daß z1 und z2 Doppelworte sind.
  18 // Bevor im Ergebnis erg ein Überlauf eintritt, wird abgebrochen.
  19
  20 #if HAVE_DD
  21
  22 static inline uintDD muluDD_unchecked(uintD q, uintDD a)
  23 {
  24         return muluD(q,lowD(a)) + highlowDD_0(muluD_unchecked(q,highD(a)));
  25 }
  26
  27 // Division: liefert min(floor(x / y), 2^intDsize-1).
  28 // Vorausgesetzt wird, daß  x >= 2 * y > 0.
  29 static uintD floorDD (uintDD x, uintDD y)
  30 {
  31         // vgl. Algorithmus für divu_3232_3232().
  32         var uintD q;
  33         if (y < ((uintDD)1 << intDsize)) {
  34                 if (highD(x) >= y)
  35                         q = ~(uintD)0; // instead of overflow
  36                 else
  37                         divuD(x, (uintD)y, q =, );
  38                 return q;
  39         }
  40         {
  41                 var uintC shift;
  42                 integerlengthD(highD(y), shift=);
  43                 // NB: 0 < shift < intDsize since 2^intDsize <= y < 2^(2*intDsize-1).
  44                 // Determine q := floor((x>>shift) / ((y>>shift)+1)).
  45                 var uintD y_shifted = y >> shift;
  46                 y_shifted += 1;
  47                 if (y_shifted == 0)
  48                         q = highD(x) >> shift;
  49                 else
  50                         divuD(x,y_shifted, q =, );
  51         }
  52         // May need to increment q at most twice.
  53         {
  54                 var uintDD p = muluDD_unchecked(q,y);
  55                 #ifdef DEBUG_GCD
  56                 if (x < p)
  57                         cl_abort();
  58                 #endif
  59                 x -= p;
  60         }
  61         if (x >= y) {
  62                 q += 1;
  63                 x -= y;
  64                 if (x >= y) {
  65                         q += 1;
  66                         #ifdef DEBUG_GCD
  67                         x -= y;
  68                         if (x >= y)
  69                                 cl_abort();
  70                         #endif
  71                 }
  72         }
  73         return q;
  74 }
  75
  76 void partial_gcd (uintDD z1, uintDD z2, partial_gcd_result* erg)
  77 {
  78     var uintD x1 = 1;
  79     var uintD y1 = 0;
  80     var uintD x2 = 0;
  81     var uintD y2 = 1;
  82     for (;;) {
  83         // Hier ist z1-y1>=z2+y2.
  84         // Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
  85         {
  86                 var uintDD zaehler = z1 - (uintDD)y1;
  87                 var uintDD nenner = z2 + (uintDD)y2;
  88                 // z2+y2 <= z1-y1 < beta^2 !
  89                 if (x2 > (~x1) >> 3) // x1 + 8*x2 >= beta ?
  90                         goto do_subtract_1;
  91                 if (y2 > (~y1) >> 3) // y1 + 8*y2 >= beta ?
  92                         goto do_subtract_1;
  93                 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_1.
  94                 if ((zaehler >> 3) < nenner)
  95                         goto do_subtract_1;
  96                 if (1) {
  97                         var uintD q = floorDD(zaehler,nenner);
  98                     repeat_1:
  99                         var uintDD qx = muluD(q,x2);
 100                         if (qx > (uintDD)(~x1)) {
 101                                 // Choose a smaller value for q, to avoid overflow of x1.
 102                                 q = floorD(~x1,x2);
 103                                 goto repeat_1;
 104                         }
 105                         var uintDD qy = muluD(q,y2);
 106                         if (qy > (uintDD)(~y1)) {
 107                                 // Choose a smaller value for q, to avoid overflow of y1.
 108                                 q = floorD(~y1,y2);
 109                                 goto repeat_1;
 110                         }
 111                         x1 += (uintD)qx;
 112                         y1 += (uintD)qy;
 113                         z1 -= muluDD_unchecked(q,z2);
 114                 } else {
 115                     do_subtract_1:
 116                         do {
 117                                 // Checks to avoid overflow.
 118                                 if (x2 > ~x1) goto done;
 119                                 if (y2 > ~y1) goto done;
 120                                 #ifdef DEBUG_GCD
 121                                 if (z1 < z2) cl_abort();
 122                                 #endif
 123                                 // Now really subtract.
 124                                 x1 += x2;
 125                                 y1 += y2;
 126                                 z1 -= z2;
 127                         } while (z1 - (uintDD)y1 >= nenner);
 128                 }
 129         }
 130         if (z2 - (uintDD)x2 <= z1 + (uintDD)(x1-1)) goto done;
 131         // Hier ist z2-x2>=z1+x1.
 132         // Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
 133         {
 134                 var uintDD zaehler = z2 - (uintDD)x2;
 135                 var uintDD nenner = z1 + (uintDD)x1;
 136                 // z1+x1 <= z2-x2 < beta^2 !
 137                 if (x1 > (~x2) >> 3) // x2 + 8*x1 >= beta ?
 138                         goto do_subtract_2;
 139                 if (y1 > (~y2) >> 3) // y2 + 8*y1 >= beta ?
 140                         goto do_subtract_2;
 141                 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_2.
 142                 if ((zaehler >> 3) < nenner)
 143                         goto do_subtract_2;
 144                 if (1) {
 145                         var uintD q = floorDD(zaehler,nenner);
 146                     repeat_2:
 147                         var uintDD qx = muluD(q,x1);
 148                         if (qx > (uintDD)(~x2)) {
 149                                 // Choose a smaller value for q, to avoid overflow of x2.
 150                                 q = floorD(~x2,x1);
 151                                 goto repeat_2;
 152                         }
 153                         var uintDD qy = muluD(q,y1);
 154                         if (qy > (uintDD)(~y2)) {
 155                                 // Choose a smaller value for q, to avoid overflow of y2.
 156                                 q = floorD(~y2,y1);
 157                                 goto repeat_2;
 158                         }
 159                         x2 += (uintD)qx;
 160                         y2 += (uintD)qy;
 161                         z2 -= muluDD_unchecked(q,z1);
 162                 } else {
 163                     do_subtract_2:
 164                         do {
 165                                 // Checks to avoid overflow.
 166                                 if (x1 > ~x2) goto done;
 167                                 if (y1 > ~y2) goto done;
 168                                 #ifdef DEBUG_GCD
 169                                 if (z2 < z1) cl_abort();
 170                                 #endif
 171                                 // Now really subtract.
 172                                 x2 += x1;
 173                                 y2 += y1;
 174                                 z2 -= z1;
 175                         } while (z2 - (uintDD)x2 >= nenner);
 176                 }
 177         }
 178         if (z1 - (uintDD)y1 <= z2 + (uintDD)(y2-1)) goto done;
 179     }
 180     done:
 181         // Keine Subtraktion (ohne Überlauf) mehr möglich.
 182         erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; // Ergebnis
 183 }
 184
 185 #else
 186
 187 // Division: liefert min(floor(xhi|xlo / yhi|ylo), 2^intDsize-1).
 188 // Vorausgesetzt wird, daß  xhi|xlo >= 2 * yhi|ylo > 0.
 189 static uintD floorDD (uintD xhi, uintD xlo, uintD yhi, uintD ylo)
 190 {
 191         // vgl. Algorithmus für divu_3232_3232().
 192         var uintD q;
 193         if (yhi == 0) {
 194                 if (xhi >= ylo)
 195                         q = ~(uintD)0; // instead of overflow
 196                 else
 197                         divuD(xhi,xlo, ylo, q =, );
 198                 return q;
 199         }
 200         {
 201                 var uintC shift;
 202                 integerlengthD(yhi, shift=);
 203                 // NB: 0 < shift < intDsize since 2^intDsize <= y < 2^(2*intDsize-1).
 204                 // Determine q := floor((x>>shift) / ((y>>shift)+1)).
 205                 var uintD y_shifted = (uintD)(yhi << (intDsize-shift)) | (ylo >> shift);
 206                 y_shifted += 1;
 207                 if (y_shifted == 0)
 208                         q = xhi >> shift;
 209                 else
 210                         divuD(xhi >> shift, (uintD)(xhi << (intDsize-shift)) | (xlo >> shift),
 211                               y_shifted,
 212                               q =, );
 213         }
 214         // May need to increment q at most twice.
 215         {
 216                 var uintD phi;
 217                 var uintD plo;
 218                 muluD(q,ylo, phi =, plo =);
 219                 muluD(q,yhi,      , phi +=);
 220                 #ifdef DEBUG_GCD
 221                 if ((xhi < phi) || ((xhi == phi) && (xlo < plo)))
 222                         cl_abort();
 223                 #endif
 224                 xhi = xhi - phi;
 225                 if (xlo < plo)
 226                         xhi -= 1;
 227                 xlo = xlo - plo;
 228         }
 229         if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo))) {
 230                 q += 1;
 231                 xhi = xhi - yhi;
 232                 if (xlo < ylo)
 233                         xhi -= 1;
 234                 xlo = xlo - ylo;
 235                 if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo))) {
 236                         q += 1;
 237                         #ifdef DEBUG_GCD
 238                         xhi = xhi - yhi;
 239                         if (xlo < ylo)
 240                                 xhi -= 1;
 241                         xlo = xlo - ylo;
 242                         if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo)))
 243                                 cl_abort();
 244                         #endif
 245                 }
 246         }
 247         return q;
 248 }
 249
 250 void partial_gcd (uintD z1hi, uintD z1lo, uintD z2hi, uintD z2lo, partial_gcd_result* erg)
 251 {
 252     var uintD x1 = 1;
 253     var uintD y1 = 0;
 254     var uintD x2 = 0;
 255     var uintD y2 = 1;
 256     for (;;) {
 257         // Hier ist z1-y1>=z2+y2.
 258         // Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
 259         {
 260                 var uintD zaehlerhi = z1hi;
 261                 var uintD zaehlerlo = z1lo - y1;
 262                 if (zaehlerlo > z1lo) { zaehlerhi--; }
 263                 var uintD nennerhi = z2hi;
 264                 var uintD nennerlo = z2lo + y2;
 265                 if (nennerlo < z2lo) { nennerhi++; }
 266                 // z2+y2 <= z1-y1 < beta^2 !
 267                 if (x2 > (~x1) >> 3) // x1 + 8*x2 >= beta ?
 268                         goto do_subtract_1;
 269                 if (y2 > (~y1) >> 3) // y1 + 8*y2 >= beta ?
 270                         goto do_subtract_1;
 271                 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_1.
 272                 if ((zaehlerhi >> 3) < nennerhi)
 273                         goto do_subtract_1;
 274                 if ((zaehlerhi >> 3) == nennerhi)
 275                         if (((uintD)(zaehlerhi << (intDsize-3)) | (zaehlerlo >> 3)) < nennerlo)
 276                                 goto do_subtract_1;
 277                 if (1) {
 278                         var uintD q = floorDD(zaehlerhi,zaehlerlo,nennerhi,nennerlo);
 279                     repeat_1:
 280                         var uintD qx;
 281                         var uintD qy;
 282                         {
 283                                 var uintD qxhi;
 284                                 muluD(q,x2, qxhi =, qx =);
 285                                 if ((qxhi > 0) || (qx > ~x1)) {
 286                                         // Choose a smaller value for q, to avoid overflow of x1.
 287                                         q = floorD(~x1,x2);
 288                                         goto repeat_1;
 289                                 }
 290                         }
 291                         {
 292                                 var uintD qyhi;
 293                                 muluD(q,y2, qyhi =, qy =);
 294                                 if ((qyhi > 0) || (qy > ~y1)) {
 295                                         // Choose a smaller value for q, to avoid overflow of y1.
 296                                         q = floorD(~y1,y2);
 297                                         goto repeat_1;
 298                                 }
 299                         }
 300                         x1 += qx;
 301                         y1 += qy;
 302                         {
 303                                 var uintD qzhi;
 304                                 var uintD qzlo;
 305                                 muluD(q,z2lo, qzhi =, qzlo =);
 306                                 muluD(q,z2hi,       , qzhi +=);
 307                                 z1hi -= qzhi;
 308                                 if (z1lo < qzlo)
 309                                         z1hi -= 1;
 310                                 z1lo -= qzlo;
 311                         }
 312                 } else {
 313                     do_subtract_1:
 314                         for (;;) {
 315                                 // Checks to avoid overflow.
 316                                 if (x2 > ~x1) goto done;
 317                                 if (y2 > ~y1) goto done;
 318                                 #ifdef DEBUG_GCD
 319                                 if (z1hi < z2hi) cl_abort();
 320                                 if (z1hi == z2hi) if (z1lo < z2lo) cl_abort();
 321                                 #endif
 322                                 // Now really subtract.
 323                                 x1 += x2;
 324                                 y1 += y2;
 325                                 z1hi -= z2hi;
 326                                 if (z1lo < z2lo)
 327                                         z1hi -= 1;
 328                                 z1lo -= z2lo;
 329                                 var uintD z1dec_hi = z1hi;
 330                                 var uintD z1dec_lo = z1lo - y1;
 331                                 if (z1lo < y1)
 332                                         z1dec_hi -= 1;
 333                                 if (z1dec_hi < nennerhi)
 334                                         break;
 335                                 if (z1dec_hi == nennerhi)
 336                                         if (z1dec_lo < nennerlo)
 337                                                 break;
 338                         }
 339                 }
 340         }
 341         {
 342                 var uintD z1inc_hi = z1hi;
 343                 var uintD z1inc_lo = z1lo + x1-1;
 344                 if (z1inc_lo < z1lo)
 345                         z1inc_hi += 1;
 346                 var uintD z2dec_hi = z2hi;
 347                 var uintD z2dec_lo = z2lo - x2;
 348                 if (z2dec_lo > z2lo)
 349                         z2dec_hi -= 1;
 350                 if (z2dec_hi < z1inc_hi) goto done;
 351                 if (z2dec_hi == z1inc_hi) if (z2dec_lo <= z1inc_lo) goto done;
 352         }
 353         // Hier ist z2-x2>=z1+x1.
 354         // Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
 355         {
 356                 var uintD zaehlerhi = z2hi;
 357                 var uintD zaehlerlo = z2lo - x2;
 358                 if (zaehlerlo > z2lo) { zaehlerhi--; }
 359                 var uintD nennerhi = z1hi;
 360                 var uintD nennerlo = z1lo + x1;
 361                 if (nennerlo < z1lo) { nennerhi++; }
 362                 // z1+x1 <= z2-x2 < beta^2 !
 363                 if (x1 > (~x2) >> 3) // x2 + 8*x1 >= beta ?
 364                         goto do_subtract_2;
 365                 if (y1 > (~y2) >> 3) // y2 + 8*y1 >= beta ?
 366                         goto do_subtract_2;
 367                 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_2.
 368                 if ((zaehlerhi >> 3) < nennerhi)
 369                         goto do_subtract_2;
 370                 if ((zaehlerhi >> 3) == nennerhi)
 371                         if (((uintD)(zaehlerhi << (intDsize-3)) | (zaehlerlo >> 3)) < nennerlo)
 372                                 goto do_subtract_2;
 373                 if (1) {
 374                         var uintD q = floorDD(zaehlerhi,zaehlerlo,nennerhi,nennerlo);
 375                     repeat_2:
 376                         var uintD qx;
 377                         var uintD qy;
 378                         {
 379                                 var uintD qxhi;
 380                                 muluD(q,x1, qxhi =, qx =);
 381                                 if ((qxhi > 0) || (qx > ~x2)) {
 382                                         // Choose a smaller value for q, to avoid overflow of x2.
 383                                         q = floorD(~x2,x1);
 384                                         goto repeat_2;
 385                                 }
 386                         }
 387                         {
 388                                 var uintD qyhi;
 389                                 muluD(q,y1, qyhi =, qy =);
 390                                 if ((qyhi > 0) || (qy > ~y2)) {
 391                                         // Choose a smaller value for q, to avoid overflow of y2.
 392                                         q = floorD(~y2,y1);
 393                                         goto repeat_2;
 394                                 }
 395                         }
 396                         x2 += qx;
 397                         y2 += qy;
 398                         {
 399                                 var uintD qzhi;
 400                                 var uintD qzlo;
 401                                 muluD(q,z1lo, qzhi =, qzlo =);
 402                                 muluD(q,z1hi,       , qzhi +=);
 403                                 z2hi -= qzhi;
 404                                 if (z2lo < qzlo)
 405                                         z2hi -= 1;
 406                                 z2lo -= qzlo;
 407                         }
 408                 } else {
 409                     do_subtract_2:
 410                         for (;;) {
 411                                 // Checks to avoid overflow.
 412                                 if (x1 > ~x2) goto done;
 413                                 if (y1 > ~y2) goto done;
 414                                 #ifdef DEBUG_GCD
 415                                 if (z2hi < z1hi) cl_abort();
 416                                 if (z2hi == z1hi) if (z2lo < z1lo) cl_abort();
 417                                 #endif
 418                                 // Now really subtract.
 419                                 x2 += x1;
 420                                 y2 += y1;
 421                                 z2hi -= z1hi;
 422                                 if (z2lo < z1lo)
 423                                         z2hi -= 1;
 424                                 z2lo -= z1lo;
 425                                 var uintD z2dec_hi = z2hi;
 426                                 var uintD z2dec_lo = z2lo - x2;
 427                                 if (z2lo < x2)
 428                                         z2dec_hi -= 1;
 429                                 if (z2dec_hi < nennerhi)
 430                                         break;
 431                                 if (z2dec_hi == nennerhi)
 432                                         if (z2dec_lo < nennerlo)
 433                                                 break;
 434                         }
 435                 }
 436         }
 437         {
 438                 var uintD z2inc_hi = z2hi;
 439                 var uintD z2inc_lo = z2lo + y2-1;
 440                 if (z2inc_lo < z2lo)
 441                         z2inc_hi += 1;
 442                 var uintD z1dec_hi = z1hi;
 443                 var uintD z1dec_lo = z1lo - y1;
 444                 if (z1dec_lo > z1lo)
 445                         z1dec_hi -= 1;
 446                 if (z1dec_hi < z2inc_hi) goto done;
 447                 if (z1dec_hi == z2inc_hi) if (z1dec_lo <= z2inc_lo) goto done;
 448         }
 449     }
 450     done:
 451         // Keine Subtraktion (ohne Überlauf) mehr möglich.
 452         erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; // Ergebnis
 453 }
 454
 455 #endif
 456
 457 }  // namespace cln