Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. /* Copyright (C) 2007-2015 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 3, or (at your option) any later
  8. version.
  9.  
  10. GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  11. WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. for more details.
  14.  
  15. Under Section 7 of GPL version 3, you are granted additional
  16. permissions described in the GCC Runtime Library Exception, version
  17. 3.1, as published by the Free Software Foundation.
  18.  
  19. You should have received a copy of the GNU General Public License and
  20. a copy of the GCC Runtime Library Exception along with this program;
  21. see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  22. <http://www.gnu.org/licenses/>.  */
  23.  
  24. #define BID_128RES
  25. #include "bid_internal.h"
  26.  
  27. /*****************************************************************************
  28.  *  BID128 nextup
  29.  ****************************************************************************/
  30.  
  31. #if DECIMAL_CALL_BY_REFERENCE
  32. void
  33. bid128_nextup (UINT128 * pres,
  34.                UINT128 *
  35.                px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  36.   UINT128 x = *px;
  37. #else
  38. UINT128
  39. bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  40.                _EXC_INFO_PARAM) {
  41. #endif
  42.  
  43.   UINT128 res;
  44.   UINT64 x_sign;
  45.   UINT64 x_exp;
  46.   int exp;
  47.   BID_UI64DOUBLE tmp1;
  48.   int x_nr_bits;
  49.   int q1, ind;
  50.   UINT128 C1;                   // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
  51.  
  52.   BID_SWAP128 (x);
  53.   // unpack the argument
  54.   x_sign = x.w[1] & MASK_SIGN;  // 0 for positive, MASK_SIGN for negative
  55.   C1.w[1] = x.w[1] & MASK_COEFF;
  56.   C1.w[0] = x.w[0];
  57.  
  58.   // check for NaN or Infinity
  59.   if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  60.     // x is special
  61.     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
  62.       // if x = NaN, then res = Q (x)
  63.       // check first for non-canonical NaN payload
  64.       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  65.           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
  66.            && (x.w[0] > 0x38c15b09ffffffffull))) {
  67.         x.w[1] = x.w[1] & 0xffffc00000000000ull;
  68.         x.w[0] = 0x0ull;
  69.       }
  70.       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
  71.         // set invalid flag
  72.         *pfpsf |= INVALID_EXCEPTION;
  73.         // return quiet (x)
  74.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
  75.         res.w[0] = x.w[0];
  76.       } else {  // x is QNaN
  77.         // return x
  78.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
  79.         res.w[0] = x.w[0];
  80.       }
  81.     } else {    // x is not NaN, so it must be infinity
  82.       if (!x_sign) {    // x is +inf
  83.         res.w[1] = 0x7800000000000000ull;       // +inf
  84.         res.w[0] = 0x0000000000000000ull;
  85.       } else {  // x is -inf
  86.         res.w[1] = 0xdfffed09bead87c0ull;       // -MAXFP = -999...99 * 10^emax
  87.         res.w[0] = 0x378d8e63ffffffffull;
  88.       }
  89.     }
  90.     BID_RETURN (res);
  91.   }
  92.   // check for non-canonical values (treated as zero)
  93.   if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {      // G0_G1=11
  94.     // non-canonical
  95.     x_exp = (x.w[1] << 2) & MASK_EXP;   // biased and shifted left 49 bits
  96.     C1.w[1] = 0;        // significand high
  97.     C1.w[0] = 0;        // significand low
  98.   } else {      // G0_G1 != 11
  99.     x_exp = x.w[1] & MASK_EXP;  // biased and shifted left 49 bits
  100.     if (C1.w[1] > 0x0001ed09bead87c0ull ||
  101.         (C1.w[1] == 0x0001ed09bead87c0ull
  102.          && C1.w[0] > 0x378d8e63ffffffffull)) {
  103.       // x is non-canonical if coefficient is larger than 10^34 -1
  104.       C1.w[1] = 0;
  105.       C1.w[0] = 0;
  106.     } else {    // canonical
  107.       ;
  108.     }
  109.   }
  110.  
  111.   if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  112.     // x is +/-0
  113.     res.w[1] = 0x0000000000000000ull;   // +1 * 10^emin
  114.     res.w[0] = 0x0000000000000001ull;
  115.   } else {      // x is not special and is not zero
  116.     if (x.w[1] == 0x5fffed09bead87c0ull
  117.         && x.w[0] == 0x378d8e63ffffffffull) {
  118.       // x = +MAXFP = 999...99 * 10^emax
  119.       res.w[1] = 0x7800000000000000ull; // +inf
  120.       res.w[0] = 0x0000000000000000ull;
  121.     } else if (x.w[1] == 0x8000000000000000ull
  122.                && x.w[0] == 0x0000000000000001ull) {
  123.       // x = -MINFP = 1...99 * 10^emin
  124.       res.w[1] = 0x8000000000000000ull; // -0
  125.       res.w[0] = 0x0000000000000000ull;
  126.     } else {    // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  127.       // can add/subtract 1 ulp to the significand
  128.  
  129.       // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
  130.       // q1 = nr. of decimal digits in x
  131.       // determine first the nr. of bits in x
  132.       if (C1.w[1] == 0) {
  133.         if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  134.           // split the 64-bit value in two 32-bit halves to avoid rnd errors
  135.           if (C1.w[0] >= 0x0000000100000000ull) {       // x >= 2^32
  136.             tmp1.d = (double) (C1.w[0] >> 32);  // exact conversion
  137.             x_nr_bits =
  138.               33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
  139.                     0x3ff);
  140.           } else {      // x < 2^32
  141.             tmp1.d = (double) (C1.w[0]);        // exact conversion
  142.             x_nr_bits =
  143.               1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
  144.                    0x3ff);
  145.           }
  146.         } else {        // if x < 2^53
  147.           tmp1.d = (double) C1.w[0];    // exact conversion
  148.           x_nr_bits =
  149.             1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  150.         }
  151.       } else {  // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  152.         tmp1.d = (double) C1.w[1];      // exact conversion
  153.         x_nr_bits =
  154.           65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  155.       }
  156.       q1 = nr_digits[x_nr_bits - 1].digits;
  157.       if (q1 == 0) {
  158.         q1 = nr_digits[x_nr_bits - 1].digits1;
  159.         if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  160.             || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  161.                 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  162.           q1++;
  163.       }
  164.       // if q1 < P34 then pad the significand with zeros
  165.       if (q1 < P34) {
  166.         exp = (x_exp >> 49) - 6176;
  167.         if (exp + 6176 > P34 - q1) {
  168.           ind = P34 - q1;       // 1 <= ind <= P34 - 1
  169.           // pad with P34 - q1 zeros, until exponent = emin
  170.           // C1 = C1 * 10^ind
  171.           if (q1 <= 19) {       // 64-bit C1
  172.             if (ind <= 19) {    // 64-bit 10^ind and 64-bit C1
  173.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  174.             } else {    // 128-bit 10^ind and 64-bit C1
  175.               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  176.             }
  177.           } else {      // C1 is (most likely) 128-bit
  178.             if (ind <= 14) {    // 64-bit 10^ind and 128-bit C1 (most likely)
  179.               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
  180.             } else if (ind <= 19) {     // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
  181.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  182.             } else {    // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
  183.               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  184.             }
  185.           }
  186.           x_exp = x_exp - ((UINT64) ind << 49);
  187.         } else {        // pad with zeros until the exponent reaches emin
  188.           ind = exp + 6176;
  189.           // C1 = C1 * 10^ind
  190.           if (ind <= 19) {      // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
  191.             if (q1 <= 19) {     // 64-bit C1, 64-bit 10^ind
  192.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  193.             } else {    // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
  194.               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
  195.             }
  196.           } else {      // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
  197.             // 64-bit C1, 128-bit 10^ind
  198.             __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  199.           }
  200.           x_exp = EXP_MIN;
  201.         }
  202.       }
  203.       if (!x_sign) {    // x > 0
  204.         // add 1 ulp (add 1 to the significand)
  205.         C1.w[0]++;
  206.         if (C1.w[0] == 0)
  207.           C1.w[1]++;
  208.         if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) {     // if  C1 = 10^34
  209.           C1.w[1] = 0x0000314dc6448d93ull;      // C1 = 10^33
  210.           C1.w[0] = 0x38c15b0a00000000ull;
  211.           x_exp = x_exp + EXP_P1;
  212.         }
  213.       } else {  // x < 0
  214.         // subtract 1 ulp (subtract 1 from the significand)
  215.         C1.w[0]--;
  216.         if (C1.w[0] == 0xffffffffffffffffull)
  217.           C1.w[1]--;
  218.         if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) {       // if  C1 = 10^33 - 1
  219.           C1.w[1] = 0x0001ed09bead87c0ull;      // C1 = 10^34 - 1
  220.           C1.w[0] = 0x378d8e63ffffffffull;
  221.           x_exp = x_exp - EXP_P1;
  222.         }
  223.       }
  224.       // assemble the result
  225.       res.w[1] = x_sign | x_exp | C1.w[1];
  226.       res.w[0] = C1.w[0];
  227.     }   // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  228.   }     // end x is not special and is not zero
  229.   BID_RETURN (res);
  230. }
  231.  
  232. /*****************************************************************************
  233.  *  BID128 nextdown
  234.  ****************************************************************************/
  235.  
  236. #if DECIMAL_CALL_BY_REFERENCE
  237. void
  238. bid128_nextdown (UINT128 * pres,
  239.                  UINT128 *
  240.                  px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  241.   UINT128 x = *px;
  242. #else
  243. UINT128
  244. bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  245.                  _EXC_INFO_PARAM) {
  246. #endif
  247.  
  248.   UINT128 res;
  249.   UINT64 x_sign;
  250.   UINT64 x_exp;
  251.   int exp;
  252.   BID_UI64DOUBLE tmp1;
  253.   int x_nr_bits;
  254.   int q1, ind;
  255.   UINT128 C1;                   // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
  256.  
  257.   BID_SWAP128 (x);
  258.   // unpack the argument
  259.   x_sign = x.w[1] & MASK_SIGN;  // 0 for positive, MASK_SIGN for negative
  260.   C1.w[1] = x.w[1] & MASK_COEFF;
  261.   C1.w[0] = x.w[0];
  262.  
  263.   // check for NaN or Infinity
  264.   if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  265.     // x is special
  266.     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
  267.       // if x = NaN, then res = Q (x)
  268.       // check first for non-canonical NaN payload
  269.       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  270.           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
  271.            && (x.w[0] > 0x38c15b09ffffffffull))) {
  272.         x.w[1] = x.w[1] & 0xffffc00000000000ull;
  273.         x.w[0] = 0x0ull;
  274.       }
  275.       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
  276.         // set invalid flag
  277.         *pfpsf |= INVALID_EXCEPTION;
  278.         // return quiet (x)
  279.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
  280.         res.w[0] = x.w[0];
  281.       } else {  // x is QNaN
  282.         // return x
  283.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
  284.         res.w[0] = x.w[0];
  285.       }
  286.     } else {    // x is not NaN, so it must be infinity
  287.       if (!x_sign) {    // x is +inf
  288.         res.w[1] = 0x5fffed09bead87c0ull;       // +MAXFP = +999...99 * 10^emax
  289.         res.w[0] = 0x378d8e63ffffffffull;
  290.       } else {  // x is -inf
  291.         res.w[1] = 0xf800000000000000ull;       // -inf
  292.         res.w[0] = 0x0000000000000000ull;
  293.       }
  294.     }
  295.     BID_RETURN (res);
  296.   }
  297.   // check for non-canonical values (treated as zero)
  298.   if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {      // G0_G1=11
  299.     // non-canonical
  300.     x_exp = (x.w[1] << 2) & MASK_EXP;   // biased and shifted left 49 bits
  301.     C1.w[1] = 0;        // significand high
  302.     C1.w[0] = 0;        // significand low
  303.   } else {      // G0_G1 != 11
  304.     x_exp = x.w[1] & MASK_EXP;  // biased and shifted left 49 bits
  305.     if (C1.w[1] > 0x0001ed09bead87c0ull ||
  306.         (C1.w[1] == 0x0001ed09bead87c0ull
  307.          && C1.w[0] > 0x378d8e63ffffffffull)) {
  308.       // x is non-canonical if coefficient is larger than 10^34 -1
  309.       C1.w[1] = 0;
  310.       C1.w[0] = 0;
  311.     } else {    // canonical
  312.       ;
  313.     }
  314.   }
  315.  
  316.   if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  317.     // x is +/-0
  318.     res.w[1] = 0x8000000000000000ull;   // -1 * 10^emin
  319.     res.w[0] = 0x0000000000000001ull;
  320.   } else {      // x is not special and is not zero
  321.     if (x.w[1] == 0xdfffed09bead87c0ull
  322.         && x.w[0] == 0x378d8e63ffffffffull) {
  323.       // x = -MAXFP = -999...99 * 10^emax
  324.       res.w[1] = 0xf800000000000000ull; // -inf
  325.       res.w[0] = 0x0000000000000000ull;
  326.     } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) {   // +MINFP
  327.       res.w[1] = 0x0000000000000000ull; // +0
  328.       res.w[0] = 0x0000000000000000ull;
  329.     } else {    // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  330.       // can add/subtract 1 ulp to the significand
  331.  
  332.       // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
  333.       // q1 = nr. of decimal digits in x
  334.       // determine first the nr. of bits in x
  335.       if (C1.w[1] == 0) {
  336.         if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  337.           // split the 64-bit value in two 32-bit halves to avoid rnd errors
  338.           if (C1.w[0] >= 0x0000000100000000ull) {       // x >= 2^32
  339.             tmp1.d = (double) (C1.w[0] >> 32);  // exact conversion
  340.             x_nr_bits =
  341.               33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
  342.                     0x3ff);
  343.           } else {      // x < 2^32
  344.             tmp1.d = (double) (C1.w[0]);        // exact conversion
  345.             x_nr_bits =
  346.               1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
  347.                    0x3ff);
  348.           }
  349.         } else {        // if x < 2^53
  350.           tmp1.d = (double) C1.w[0];    // exact conversion
  351.           x_nr_bits =
  352.             1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  353.         }
  354.       } else {  // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  355.         tmp1.d = (double) C1.w[1];      // exact conversion
  356.         x_nr_bits =
  357.           65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  358.       }
  359.       q1 = nr_digits[x_nr_bits - 1].digits;
  360.       if (q1 == 0) {
  361.         q1 = nr_digits[x_nr_bits - 1].digits1;
  362.         if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  363.             || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  364.                 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  365.           q1++;
  366.       }
  367.       // if q1 < P then pad the significand with zeros
  368.       if (q1 < P34) {
  369.         exp = (x_exp >> 49) - 6176;
  370.         if (exp + 6176 > P34 - q1) {
  371.           ind = P34 - q1;       // 1 <= ind <= P34 - 1
  372.           // pad with P34 - q1 zeros, until exponent = emin
  373.           // C1 = C1 * 10^ind
  374.           if (q1 <= 19) {       // 64-bit C1
  375.             if (ind <= 19) {    // 64-bit 10^ind and 64-bit C1
  376.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  377.             } else {    // 128-bit 10^ind and 64-bit C1
  378.               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  379.             }
  380.           } else {      // C1 is (most likely) 128-bit
  381.             if (ind <= 14) {    // 64-bit 10^ind and 128-bit C1 (most likely)
  382.               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
  383.             } else if (ind <= 19) {     // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
  384.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  385.             } else {    // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
  386.               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  387.             }
  388.           }
  389.           x_exp = x_exp - ((UINT64) ind << 49);
  390.         } else {        // pad with zeros until the exponent reaches emin
  391.           ind = exp + 6176;
  392.           // C1 = C1 * 10^ind
  393.           if (ind <= 19) {      // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
  394.             if (q1 <= 19) {     // 64-bit C1, 64-bit 10^ind
  395.               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
  396.             } else {    // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
  397.               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
  398.             }
  399.           } else {      // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
  400.             // 64-bit C1, 128-bit 10^ind
  401.             __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
  402.           }
  403.           x_exp = EXP_MIN;
  404.         }
  405.       }
  406.       if (x_sign) {     // x < 0
  407.         // add 1 ulp (add 1 to the significand)
  408.         C1.w[0]++;
  409.         if (C1.w[0] == 0)
  410.           C1.w[1]++;
  411.         if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) {     // if  C1 = 10^34
  412.           C1.w[1] = 0x0000314dc6448d93ull;      // C1 = 10^33
  413.           C1.w[0] = 0x38c15b0a00000000ull;
  414.           x_exp = x_exp + EXP_P1;
  415.         }
  416.       } else {  // x > 0
  417.         // subtract 1 ulp (subtract 1 from the significand)
  418.         C1.w[0]--;
  419.         if (C1.w[0] == 0xffffffffffffffffull)
  420.           C1.w[1]--;
  421.         if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) {       // if  C1 = 10^33 - 1
  422.           C1.w[1] = 0x0001ed09bead87c0ull;      // C1 = 10^34 - 1
  423.           C1.w[0] = 0x378d8e63ffffffffull;
  424.           x_exp = x_exp - EXP_P1;
  425.         }
  426.       }
  427.       // assemble the result
  428.       res.w[1] = x_sign | x_exp | C1.w[1];
  429.       res.w[0] = C1.w[0];
  430.     }   // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  431.   }     // end x is not special and is not zero
  432.   BID_RETURN (res);
  433. }
  434.  
  435. /*****************************************************************************
  436.  *  BID128 nextafter
  437.  ****************************************************************************/
  438.  
  439. #if DECIMAL_CALL_BY_REFERENCE
  440. void
  441. bid128_nextafter (UINT128 * pres, UINT128 * px,
  442.                   UINT128 *
  443.                   py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
  444. {
  445.   UINT128 x = *px;
  446.   UINT128 y = *py;
  447.   UINT128 xnswp = *px;
  448.   UINT128 ynswp = *py;
  449. #else
  450. UINT128
  451. bid128_nextafter (UINT128 x,
  452.                   UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  453.                   _EXC_INFO_PARAM) {
  454.   UINT128 xnswp = x;
  455.   UINT128 ynswp = y;
  456. #endif
  457.  
  458.   UINT128 res;
  459.   UINT128 tmp1, tmp2, tmp3;
  460.   FPSC tmp_fpsf = 0;            // dummy fpsf for calls to comparison functions
  461.   int res1, res2;
  462.   UINT64 x_exp;
  463.  
  464.  
  465.   BID_SWAP128 (x);
  466.   BID_SWAP128 (y);
  467.   // check for NaNs
  468.   if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
  469.       || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
  470.     // x is special or y is special
  471.     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
  472.       // if x = NaN, then res = Q (x)
  473.       // check first for non-canonical NaN payload
  474.       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  475.           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
  476.            && (x.w[0] > 0x38c15b09ffffffffull))) {
  477.         x.w[1] = x.w[1] & 0xffffc00000000000ull;
  478.         x.w[0] = 0x0ull;
  479.       }
  480.       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
  481.         // set invalid flag
  482.         *pfpsf |= INVALID_EXCEPTION;
  483.         // return quiet (x)
  484.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
  485.         res.w[0] = x.w[0];
  486.       } else {  // x is QNaN
  487.         // return x
  488.         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
  489.         res.w[0] = x.w[0];
  490.         if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {        // y is SNAN
  491.           // set invalid flag
  492.           *pfpsf |= INVALID_EXCEPTION;
  493.         }
  494.       }
  495.       BID_RETURN (res)
  496.     } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {       // y is NAN
  497.       // if x = NaN, then res = Q (x)
  498.       // check first for non-canonical NaN payload
  499.       if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  500.           (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
  501.            && (y.w[0] > 0x38c15b09ffffffffull))) {
  502.         y.w[1] = y.w[1] & 0xffffc00000000000ull;
  503.         y.w[0] = 0x0ull;
  504.       }
  505.       if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {  // y is SNAN
  506.         // set invalid flag
  507.         *pfpsf |= INVALID_EXCEPTION;
  508.         // return quiet (x)
  509.         res.w[1] = y.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
  510.         res.w[0] = y.w[0];
  511.       } else {  // x is QNaN
  512.         // return x
  513.         res.w[1] = y.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
  514.         res.w[0] = y.w[0];
  515.       }
  516.       BID_RETURN (res)
  517.     } else {    // at least one is infinity
  518.       if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {        // x = inf
  519.         x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
  520.         x.w[0] = 0x0ull;
  521.       }
  522.       if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {        // y = inf
  523.         y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
  524.         y.w[0] = 0x0ull;
  525.       }
  526.     }
  527.   }
  528.   // neither x nor y is NaN
  529.  
  530.   // if not infinity, check for non-canonical values x (treated as zero)
  531.   if ((x.w[1] & MASK_ANY_INF) != MASK_INF) {    // x != inf
  532.     if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {    // G0_G1=11
  533.       // non-canonical
  534.       x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
  535.       x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
  536.       x.w[0] = 0x0ull;
  537.     } else {    // G0_G1 != 11
  538.       x_exp = x.w[1] & MASK_EXP;        // biased and shifted left 49 bits
  539.       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  540.           ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  541.            && x.w[0] > 0x378d8e63ffffffffull)) {
  542.         // x is non-canonical if coefficient is larger than 10^34 -1
  543.         x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
  544.         x.w[0] = 0x0ull;
  545.       } else {  // canonical
  546.         ;
  547.       }
  548.     }
  549.   }
  550.   // no need to check for non-canonical y
  551.  
  552.   // neither x nor y is NaN
  553.   tmp_fpsf = *pfpsf;    // save fpsf
  554. #if DECIMAL_CALL_BY_REFERENCE
  555.   bid128_quiet_equal (&res1, &xnswp,
  556.                       &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  557.                       _EXC_INFO_ARG);
  558.   bid128_quiet_greater (&res2, &xnswp,
  559.                         &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  560.                         _EXC_INFO_ARG);
  561. #else
  562.   res1 =
  563.     bid128_quiet_equal (xnswp,
  564.                         ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  565.                         _EXC_INFO_ARG);
  566.   res2 =
  567.     bid128_quiet_greater (xnswp,
  568.                           ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  569.                           _EXC_INFO_ARG);
  570. #endif
  571.   *pfpsf = tmp_fpsf;    // restore fpsf
  572.  
  573.   if (res1) {   // x = y
  574.     // return x with the sign of y
  575.     res.w[1] =
  576.       (x.w[1] & 0x7fffffffffffffffull) | (y.
  577.                                           w[1] & 0x8000000000000000ull);
  578.     res.w[0] = x.w[0];
  579.   } else if (res2) {    // x > y
  580. #if DECIMAL_CALL_BY_REFERENCE
  581.     bid128_nextdown (&res,
  582.                      &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  583.                      _EXC_INFO_ARG);
  584. #else
  585.     res =
  586.       bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
  587.                        _EXC_INFO_ARG);
  588. #endif
  589.     BID_SWAP128 (res);
  590.   } else {      // x < y
  591. #if DECIMAL_CALL_BY_REFERENCE
  592.     bid128_nextup (&res,
  593.                    &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  594. #else
  595.     res =
  596.       bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  597. #endif
  598.     BID_SWAP128 (res);
  599.   }
  600.   // if the operand x is finite but the result is infinite, signal
  601.   // overflow and inexact
  602.   if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
  603.       && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
  604.     // set the inexact flag
  605.     *pfpsf |= INEXACT_EXCEPTION;
  606.     // set the overflow flag
  607.     *pfpsf |= OVERFLOW_EXCEPTION;
  608.   }
  609.   // if the result is in (-10^emin, 10^emin), and is different from the
  610.   // operand x, signal underflow and inexact
  611.   tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull;
  612.   tmp1.w[LOW_128W] = 0x38c15b0a00000000ull;     // +100...0[34] * 10^emin
  613.   tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
  614.   tmp2.w[LOW_128W] = res.w[0];
  615.   tmp3.w[HIGH_128W] = res.w[1];
  616.   tmp3.w[LOW_128W] = res.w[0];
  617.   tmp_fpsf = *pfpsf;    // save fpsf
  618. #if DECIMAL_CALL_BY_REFERENCE
  619.   bid128_quiet_greater (&res1, &tmp1,
  620.                         &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
  621.                         _EXC_INFO_ARG);
  622.   bid128_quiet_not_equal (&res2, &xnswp,
  623.                           &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
  624.                           _EXC_INFO_ARG);
  625. #else
  626.   res1 =
  627.     bid128_quiet_greater (tmp1,
  628.                           tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
  629.                           _EXC_INFO_ARG);
  630.   res2 =
  631.     bid128_quiet_not_equal (xnswp,
  632.                             tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
  633.                             _EXC_INFO_ARG);
  634. #endif
  635.   *pfpsf = tmp_fpsf;    // restore fpsf
  636.   if (res1 && res2) {
  637.     // set the inexact flag
  638.     *pfpsf |= INEXACT_EXCEPTION;
  639.     // set the underflow flag
  640.     *pfpsf |= UNDERFLOW_EXCEPTION;
  641.   }
  642.   BID_RETURN (res);
  643. }
  644.