libgcc/config/libbid/div_macros.h

   1 /* Copyright (C) 2007  Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 2, or (at your option) any later
   8 version.
   9
  10 In addition to the permissions in the GNU General Public License, the
  11 Free Software Foundation gives you unlimited permission to link the
  12 compiled version of this file into combinations with other programs,
  13 and to distribute those combinations without any restriction coming
  14 from the use of this file.  (The General Public License restrictions
  15 do apply in other respects; for example, they cover modification of
  16 the file, and distribution when not linked into a combine
  17 executable.)
  18
  19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  21 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  22 for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with GCC; see the file COPYING.  If not, write to the Free
  26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  27 02110-1301, USA.  */
  28
  29 #ifndef _DIV_MACROS_H_
  30 #define _DIV_MACROS_H_
  31
  32 #include "bid_internal.h"
  33
  34 #define FENCE __fence
  35 //#define FENCE
  36
  37 //#define DOUBLE_EXTENDED_ON
  38
  39 #if DOUBLE_EXTENDED_ON
  40
  41
  42 __BID_INLINE__ void
  43 __div_128_by_128 (UINT128 * pCQ, UINT128 * pCR, UINT128 CX, UINT128 CY) {
  44   UINT128 CB, CB2, CB4, CB8, CQB, CA;
  45   int_double d64, dm64, ds;
  46   int_float t64;
  47   double dx, dq, dqh;
  48   long double lq, lx, ly;
  49   UINT64 Rh, R, B2, B4, Ph, Ql, Ql2, carry, Qh;
  50
  51   if (!CY.w[1]) {
  52     pCR->w[1] = 0;
  53
  54     if (!CX.w[1]) {
  55       pCQ->w[0] = CX.w[0] / CY.w[0];
  56       pCQ->w[1] = 0;
  57       pCR->w[1] = 0;
  58       pCR->w[0] = CX.w[0] - pCQ->w[0] * CY.w[0];
  59     } else {
  60
  61       // This path works for CX<2^116 only
  62
  63       // 2^64
  64       d64.i = 0x43f0000000000000;
  65       // 2^64
  66       dm64.i = 0x3bf0000000000000;
  67       // 1.5*2^(-52)
  68       ds.i = 0x3cb8000000000000;
  69       dx = (long double) CX.w[1] * d64.d + (long double) CX.w[0];
  70       dq = dx / (long double) CY.w[0];
  71       dq -= dq * (ds.d);
  72       dqh = dq * dm64.d;
  73       Qh = (UINT64) dqh;
  74       Ql = (UINT64) (dq - ((double) Qh) * d64.d);
  75
  76       Rh = CX.w[0] - Ql * CY.w[0];
  77       Ql2 = Rh / CY.w[0];
  78       pCR->w[0] = Rh - Ql2 * CY.w[0];
  79       __add_carry_out ((pCQ->w[0]), carry, Ql, Ql2);
  80       pCQ->w[1] = Qh + carry;
  81
  82     }
  83     return;
  84   }
  85   // now CY.w[1] > 0
  86
  87   // 2^64
  88   t64.i = 0x5f800000;
  89   lx =
  90     (long double) CX.w[1] * (long double) t64.d + (long double) CX.w[0];
  91   ly =
  92     (long double) CY.w[1] * (long double) t64.d + (long double) CY.w[0];
  93   lq = lx / ly;
  94   pCQ->w[0] = (UINT64) lq;
  95
  96   pCQ->w[1] = 0;
  97
  98   if (!pCQ->w[0]) {
  99     /*if(__unsigned_compare_ge_128(CX,CY))
 100        {
 101        pCQ->w[0] = 1;
 102        __sub_128_128((*pCR), CX, CY);
 103        }
 104        else */
 105     {
 106       pCR->w[1] = CX.w[1];
 107       pCR->w[0] = CX.w[0];
 108     }
 109     return;
 110   }
 111
 112   if (CY.w[1] >= 16 || pCQ->w[0] <= 0x1000000000000000ull) {
 113     pCQ->w[0] = (UINT64) lq - 1;
 114     __mul_64x128_full (Ph, CQB, (pCQ->w[0]), CY);
 115     __sub_128_128 (CA, CX, CQB);
 116     if (__unsigned_compare_ge_128 (CA, CY)) {
 117       __sub_128_128 (CA, CA, CY);
 118       pCQ->w[0]++;
 119       if (__unsigned_compare_ge_128 (CA, CY)) {
 120         __sub_128_128 (CA, CA, CY);
 121         pCQ->w[0]++;
 122       }
 123     }
 124     pCR->w[1] = CA.w[1];
 125     pCR->w[0] = CA.w[0];
 126   } else {
 127     pCQ->w[0] = (UINT64) lq - 6;
 128
 129     __mul_64x128_full (Ph, CQB, (pCQ->w[0]), CY);
 130     __sub_128_128 (CA, CX, CQB);
 131
 132     CB8.w[1] = (CY.w[1] << 3) | (CY.w[0] >> 61);
 133     CB8.w[0] = CY.w[0] << 3;
 134     CB4.w[1] = (CY.w[1] << 2) | (CY.w[0] >> 62);
 135     CB4.w[0] = CY.w[0] << 2;
 136     CB2.w[1] = (CY.w[1] << 1) | (CY.w[0] >> 63);
 137     CB2.w[0] = CY.w[0] << 1;
 138
 139     if (__unsigned_compare_ge_128 (CA, CB8)) {
 140       pCQ->w[0] += 8;
 141       __sub_128_128 (CA, CA, CB8);
 142     }
 143     if (__unsigned_compare_ge_128 (CA, CB4)) {
 144       pCQ->w[0] += 4;
 145       __sub_128_128 (CA, CA, CB4);
 146     }
 147     if (__unsigned_compare_ge_128 (CA, CB2)) {
 148       pCQ->w[0] += 2;
 149       __sub_128_128 (CA, CA, CB2);
 150     }
 151     if (__unsigned_compare_ge_128 (CA, CY)) {
 152       pCQ->w[0] += 1;
 153       __sub_128_128 (CA, CA, CY);
 154     }
 155
 156     pCR->w[1] = CA.w[1];
 157     pCR->w[0] = CA.w[0];
 158   }
 159 }
 160
 161
 162
 163
 164
 165
 166 __BID_INLINE__ void
 167 __div_256_by_128 (UINT128 * pCQ, UINT256 * pCA4, UINT128 CY) {
 168   UINT256 CQ2Y;
 169   UINT128 CQ2, CQ3Y;
 170   UINT64 Q3, carry64;
 171   int_double d64;
 172   long double lx, ly, lq, l64, l128;
 173
 174   // 2^64
 175   d64.i = 0x43f0000000000000ull;
 176   l64 = (long double) d64.d;
 177   // 2^128
 178   l128 = l64 * l64;
 179
 180   lx =
 181     ((long double) (*pCA4).w[3] * l64 +
 182      (long double) (*pCA4).w[2]) * l128 +
 183     (long double) (*pCA4).w[1] * l64 + (long double) (*pCA4).w[0];
 184   ly = (long double) CY.w[1] * l128 + (long double) CY.w[0] * l64;
 185
 186   lq = lx / ly;
 187   CQ2.w[1] = (UINT64) lq;
 188   lq = (lq - CQ2.w[1]) * l64;
 189   CQ2.w[0] = (UINT64) lq;
 190
 191   // CQ2*CY
 192   __mul_128x128_to_256 (CQ2Y, CY, CQ2);
 193
 194   // CQ2Y <= (*pCA4) ?
 195   if (CQ2Y.w[3] < (*pCA4).w[3]
 196       || (CQ2Y.w[3] == (*pCA4).w[3]
 197           && (CQ2Y.w[2] < (*pCA4).w[2]
 198               || (CQ2Y.w[2] == (*pCA4).w[2]
 199                   && (CQ2Y.w[1] < (*pCA4).w[1]
 200                       || (CQ2Y.w[1] == (*pCA4).w[1]
 201                           && (CQ2Y.w[0] <= (*pCA4).w[0]))))))) {
 202
 203     // (*pCA4) -CQ2Y, guaranteed below 5*2^49*CY < 5*2^(49+128)
 204     __sub_borrow_out ((*pCA4).w[0], carry64, (*pCA4).w[0], CQ2Y.w[0]);
 205     __sub_borrow_in_out ((*pCA4).w[1], carry64, (*pCA4).w[1], CQ2Y.w[1],
 206                          carry64);
 207     (*pCA4).w[2] = (*pCA4).w[2] - CQ2Y.w[2] - carry64;
 208
 209     lx = ((long double) (*pCA4).w[2] * l128 +
 210           ((long double) (*pCA4).w[1] * l64 +
 211            (long double) (*pCA4).w[0])) * l64;
 212     lq = lx / ly;
 213     Q3 = (UINT64) lq;
 214
 215     if (Q3) {
 216       Q3--;
 217       __mul_64x128_short (CQ3Y, Q3, CY);
 218       __sub_borrow_out ((*pCA4).w[0], carry64, (*pCA4).w[0], CQ3Y.w[0]);
 219       (*pCA4).w[1] = (*pCA4).w[1] - CQ3Y.w[1] - carry64;
 220
 221       if ((*pCA4).w[1] > CY.w[1]
 222           || ((*pCA4).w[1] == CY.w[1] && (*pCA4).w[0] >= CY.w[0])) {
 223         Q3++;
 224         __sub_borrow_out ((*pCA4).w[0], carry64, (*pCA4).w[0], CY.w[0]);
 225         (*pCA4).w[1] = (*pCA4).w[1] - CY.w[1] - carry64;
 226         if ((*pCA4).w[1] > CY.w[1]
 227             || ((*pCA4).w[1] == CY.w[1] && (*pCA4).w[0] >= CY.w[0])) {
 228           Q3++;
 229           __sub_borrow_out ((*pCA4).w[0], carry64, (*pCA4).w[0],
 230                             CY.w[0]);
 231           (*pCA4).w[1] = (*pCA4).w[1] - CY.w[1] - carry64;
 232         }
 233       }
 234       // add Q3 to Q2
 235       __add_carry_out (CQ2.w[0], carry64, Q3, CQ2.w[0]);
 236       CQ2.w[1] += carry64;
 237     }
 238   } else {
 239     // CQ2Y - (*pCA4), guaranteed below 5*2^(49+128)
 240     __sub_borrow_out ((*pCA4).w[0], carry64, CQ2Y.w[0], (*pCA4).w[0]);
 241     __sub_borrow_in_out ((*pCA4).w[1], carry64, CQ2Y.w[1], (*pCA4).w[1],
 242                          carry64);
 243     (*pCA4).w[2] = CQ2Y.w[2] - (*pCA4).w[2] - carry64;
 244
 245     lx =
 246       ((long double) (*pCA4).w[2] * l128 +
 247        (long double) (*pCA4).w[1] * l64 +
 248        (long double) (*pCA4).w[0]) * l64;
 249     lq = lx / ly;
 250     Q3 = 1 + (UINT64) lq;
 251
 252     __mul_64x128_short (CQ3Y, Q3, CY);
 253     __sub_borrow_out ((*pCA4).w[0], carry64, CQ3Y.w[0], (*pCA4).w[0]);
 254     (*pCA4).w[1] = CQ3Y.w[1] - (*pCA4).w[1] - carry64;
 255
 256     if ((SINT64) (*pCA4).w[1] > (SINT64) CY.w[1]
 257         || ((*pCA4).w[1] == CY.w[1] && (*pCA4).w[0] >= CY.w[0])) {
 258       Q3--;
 259       __sub_borrow_out ((*pCA4).w[0], carry64, (*pCA4).w[0], CY.w[0]);
 260       (*pCA4).w[1] = (*pCA4).w[1] - CY.w[1] - carry64;
 261     } else if ((SINT64) (*pCA4).w[1] < 0) {
 262       Q3++;
 263       __add_carry_out ((*pCA4).w[0], carry64, (*pCA4).w[0], CY.w[0]);
 264       (*pCA4).w[1] = (*pCA4).w[1] + CY.w[1] + carry64;
 265     }
 266     // subtract Q3 from Q2
 267     __sub_borrow_out (CQ2.w[0], carry64, CQ2.w[0], Q3);
 268     CQ2.w[1] -= carry64;
 269   }
 270
 271   // (*pCQ) + CQ2 + carry
 272   __add_carry_out ((*pCQ).w[0], carry64, CQ2.w[0], (*pCQ).w[0]);
 273   (*pCQ).w[1] = (*pCQ).w[1] + CQ2.w[1] + carry64;
 274
 275
 276 }
 277 #else
 278
 279 __BID_INLINE__ void
 280 __div_128_by_128 (UINT128 * pCQ, UINT128 * pCR, UINT128 CX0, UINT128 CY) {
 281   UINT128 CY36, CY51, CQ, A2, CX, CQT;
 282   UINT64 Q;
 283   int_double t64, d49, d60;
 284   double lx, ly, lq;
 285
 286   if (!CX0.w[1] && !CY.w[1]) {
 287     pCQ->w[0] = CX0.w[0] / CY.w[0];
 288     pCQ->w[1] = 0;
 289     pCR->w[1] = pCR->w[0] = 0;
 290     pCR->w[0] = CX0.w[0] - pCQ->w[0] * CY.w[0];
 291     return;
 292   }
 293
 294   CX.w[1] = CX0.w[1];
 295   CX.w[0] = CX0.w[0];
 296
 297   // 2^64
 298   t64.i = 0x43f0000000000000ull;
 299   lx = (double) CX.w[1] * t64.d + (double) CX.w[0];
 300   ly = (double) CY.w[1] * t64.d + (double) CY.w[0];
 301   lq = lx / ly;
 302
 303   CY36.w[1] = CY.w[0] >> (64 - 36);
 304   CY36.w[0] = CY.w[0] << 36;
 305
 306   CQ.w[1] = CQ.w[0] = 0;
 307
 308   // Q >= 2^100 ?
 309   if (!CY.w[1] && !CY36.w[1] && (CX.w[1] >= CY36.w[0])) {
 310     // then Q >= 2^100
 311
 312     // 2^(-60)*CX/CY
 313     d60.i = 0x3c30000000000000ull;
 314     lq *= d60.d;
 315     Q = -4 + (UINT64) lq;
 316
 317     // Q*CY
 318     __mul_64x64_to_128 (A2, Q, CY.w[0]);
 319
 320     // A2 <<= 60
 321     A2.w[1] = (A2.w[1] << 60) | (A2.w[0] >> (64 - 60));
 322     A2.w[0] <<= 60;
 323
 324     __sub_128_128 (CX, CX, A2);
 325
 326     lx = (double) CX.w[1] * t64.d + (double) CX.w[0];
 327     lq = lx / ly;
 328
 329     CQ.w[1] = Q >> (64 - 60);
 330     CQ.w[0] = Q << 60;
 331   }
 332
 333
 334   CY51.w[1] = (CY.w[1] << 51) | (CY.w[0] >> (64 - 51));
 335   CY51.w[0] = CY.w[0] << 51;
 336
 337   if (CY.w[1] < (UINT64) (1 << (64 - 51))
 338       && (__unsigned_compare_gt_128 (CX, CY51))) {
 339     // Q > 2^51
 340
 341     // 2^(-49)*CX/CY
 342     d49.i = 0x3ce0000000000000ull;
 343     lq *= d49.d;
 344
 345     Q = -1 + (UINT64) lq;
 346
 347     // Q*CY
 348     __mul_64x64_to_128 (A2, Q, CY.w[0]);
 349     A2.w[1] += Q * CY.w[1];
 350
 351     // A2 <<= 49
 352     A2.w[1] = (A2.w[1] << 49) | (A2.w[0] >> (64 - 49));
 353     A2.w[0] <<= 49;
 354
 355     __sub_128_128 (CX, CX, A2);
 356
 357     CQT.w[1] = Q >> (64 - 49);
 358     CQT.w[0] = Q << 49;
 359     __add_128_128 (CQ, CQ, CQT);
 360
 361     lx = (double) CX.w[1] * t64.d + (double) CX.w[0];
 362     lq = lx / ly;
 363   }
 364
 365   Q = (UINT64) lq;
 366
 367   __mul_64x64_to_128 (A2, Q, CY.w[0]);
 368   A2.w[1] += Q * CY.w[1];
 369
 370   __sub_128_128 (CX, CX, A2);
 371   if ((SINT64) CX.w[1] < 0) {
 372     Q--;
 373     CX.w[0] += CY.w[0];
 374     if (CX.w[0] < CY.w[0])
 375       CX.w[1]++;
 376     CX.w[1] += CY.w[1];
 377     if ((SINT64) CX.w[1] < 0) {
 378       Q--;
 379       CX.w[0] += CY.w[0];
 380       if (CX.w[0] < CY.w[0])
 381         CX.w[1]++;
 382       CX.w[1] += CY.w[1];
 383     }
 384   } else if (__unsigned_compare_ge_128 (CX, CY)) {
 385     Q++;
 386     __sub_128_128 (CX, CX, CY);
 387   }
 388
 389   __add_128_64 (CQ, CQ, Q);
 390
 391
 392   pCQ->w[1] = CQ.w[1];
 393   pCQ->w[0] = CQ.w[0];
 394   pCR->w[1] = CX.w[1];
 395   pCR->w[0] = CX.w[0];
 396   return;
 397 }
 398
 399
 400 __BID_INLINE__ void
 401 __div_256_by_128 (UINT128 * pCQ, UINT256 * pCA4, UINT128 CY) {
 402   UINT256 CA4, CA2, CY51, CY36;
 403   UINT128 CQ, A2, A2h, CQT;
 404   UINT64 Q, carry64;
 405   int_double t64, d49, d60;
 406   double lx, ly, lq, d128, d192;
 407
 408   // the quotient is assumed to be at most 113 bits,
 409   // as needed by BID128 divide routines
 410
 411   // initial dividend
 412   CA4.w[3] = (*pCA4).w[3];
 413   CA4.w[2] = (*pCA4).w[2];
 414   CA4.w[1] = (*pCA4).w[1];
 415   CA4.w[0] = (*pCA4).w[0];
 416   CQ.w[1] = (*pCQ).w[1];
 417   CQ.w[0] = (*pCQ).w[0];
 418
 419   // 2^64
 420   t64.i = 0x43f0000000000000ull;
 421   d128 = t64.d * t64.d;
 422   d192 = d128 * t64.d;
 423   lx = (double) CA4.w[3] * d192 + ((double) CA4.w[2] * d128 +
 424                                    ((double) CA4.w[1] * t64.d +
 425                                     (double) CA4.w[0]));
 426   ly = (double) CY.w[1] * t64.d + (double) CY.w[0];
 427   lq = lx / ly;
 428
 429   CY36.w[2] = CY.w[1] >> (64 - 36);
 430   CY36.w[1] = (CY.w[1] << 36) | (CY.w[0] >> (64 - 36));
 431   CY36.w[0] = CY.w[0] << 36;
 432
 433   CQ.w[1] = (*pCQ).w[1];
 434   CQ.w[0] = (*pCQ).w[0];
 435
 436   // Q >= 2^100 ?
 437   if (CA4.w[3] > CY36.w[2]
 438       || (CA4.w[3] == CY36.w[2]
 439           && (CA4.w[2] > CY36.w[1]
 440               || (CA4.w[2] == CY36.w[1] && CA4.w[1] >= CY36.w[0])))) {
 441     // 2^(-60)*CA4/CY
 442     d60.i = 0x3c30000000000000ull;
 443     lq *= d60.d;
 444     Q = -4 + (UINT64) lq;
 445
 446     // Q*CY
 447     __mul_64x128_to_192 (CA2, Q, CY);
 448
 449     // CA2 <<= 60
 450     // CA2.w[3] = CA2.w[2] >> (64-60);
 451     CA2.w[2] = (CA2.w[2] << 60) | (CA2.w[1] >> (64 - 60));
 452     CA2.w[1] = (CA2.w[1] << 60) | (CA2.w[0] >> (64 - 60));
 453     CA2.w[0] <<= 60;
 454
 455     // CA4 -= CA2
 456     __sub_borrow_out (CA4.w[0], carry64, CA4.w[0], CA2.w[0]);
 457     __sub_borrow_in_out (CA4.w[1], carry64, CA4.w[1], CA2.w[1],
 458                          carry64);
 459     CA4.w[2] = CA4.w[2] - CA2.w[2] - carry64;
 460
 461     lx = ((double) CA4.w[2] * d128 +
 462           ((double) CA4.w[1] * t64.d + (double) CA4.w[0]));
 463     lq = lx / ly;
 464
 465     CQT.w[1] = Q >> (64 - 60);
 466     CQT.w[0] = Q << 60;
 467     __add_128_128 (CQ, CQ, CQT);
 468   }
 469
 470   CY51.w[2] = CY.w[1] >> (64 - 51);
 471   CY51.w[1] = (CY.w[1] << 51) | (CY.w[0] >> (64 - 51));
 472   CY51.w[0] = CY.w[0] << 51;
 473
 474   if (CA4.w[2] > CY51.w[2] ||
 475       ((CA4.w[2] == CY51.w[2])
 476        && (__unsigned_compare_gt_128 (CA4, CY51)))) {
 477     // Q > 2^51
 478
 479     // 2^(-49)*CA4/CY
 480     d49.i = 0x3ce0000000000000ull;
 481     lq *= d49.d;
 482
 483     Q = -1 + (UINT64) lq;
 484
 485     // Q*CY
 486     __mul_64x64_to_128 (A2, Q, CY.w[0]);
 487     __mul_64x64_to_128 (A2h, Q, CY.w[1]);
 488     A2.w[1] += A2h.w[0];
 489     if (A2.w[1] < A2h.w[0])
 490       A2h.w[1]++;
 491
 492     // A2 <<= 49
 493     CA2.w[2] = (A2h.w[1] << 49) | (A2.w[1] >> (64 - 49));
 494     CA2.w[1] = (A2.w[1] << 49) | (A2.w[0] >> (64 - 49));
 495     CA2.w[0] = A2.w[0] << 49;
 496
 497     __sub_borrow_out (CA4.w[0], carry64, CA4.w[0], CA2.w[0]);
 498     __sub_borrow_in_out (CA4.w[1], carry64, CA4.w[1], CA2.w[1],
 499                          carry64);
 500     CA4.w[2] = CA4.w[2] - CA2.w[2] - carry64;
 501
 502     CQT.w[1] = Q >> (64 - 49);
 503     CQT.w[0] = Q << 49;
 504     __add_128_128 (CQ, CQ, CQT);
 505
 506     lx = ((double) CA4.w[2] * d128 +
 507           ((double) CA4.w[1] * t64.d + (double) CA4.w[0]));
 508     lq = lx / ly;
 509   }
 510
 511   Q = (UINT64) lq;
 512
 513   __mul_64x64_to_128 (A2, Q, CY.w[0]);
 514   A2.w[1] += Q * CY.w[1];
 515
 516   __sub_128_128 (CA4, CA4, A2);
 517   if ((SINT64) CA4.w[1] < 0) {
 518     Q--;
 519     CA4.w[0] += CY.w[0];
 520     if (CA4.w[0] < CY.w[0])
 521       CA4.w[1]++;
 522     CA4.w[1] += CY.w[1];
 523     if ((SINT64) CA4.w[1] < 0) {
 524       Q--;
 525       CA4.w[0] += CY.w[0];
 526       if (CA4.w[0] < CY.w[0])
 527         CA4.w[1]++;
 528       CA4.w[1] += CY.w[1];
 529     }
 530   } else if (__unsigned_compare_ge_128 (CA4, CY)) {
 531     Q++;
 532     __sub_128_128 (CA4, CA4, CY);
 533   }
 534
 535   __add_128_64 (CQ, CQ, Q);
 536
 537   pCQ->w[1] = CQ.w[1];
 538   pCQ->w[0] = CQ.w[0];
 539   pCA4->w[1] = CA4.w[1];
 540   pCA4->w[0] = CA4.w[0];
 541   return;
 542
 543
 544
 545 }
 546
 547 #endif
 548 #endif