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