gdtoa/_ldbl_util.c

   1 /*
   2  * Copyright (c) 2004, 2008 Apple Inc. All rights reserved.
   3  *
   4  * @APPLE_LICENSE_HEADER_START@
   5  *
   6  * This file contains Original Code and/or Modifications of Original Code
   7  * as defined in and that are subject to the Apple Public Source License
   8  * Version 2.0 (the 'License'). You may not use this file except in
   9  * compliance with the License. Please obtain a copy of the License at
  10  * http:www.opensource.apple.com/apsl/ and read it before using this
  11  * file.
  12  *
  13  * The Original Code and all software distributed under the License are
  14  * distributed on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER
  15  * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
  16  * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
  17  * FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR NON-INFRINGEMENT.
  18  * Please see the License for the specific language governing rights and
  19  * limitations under the License.
  20  *
  21  * @APPLE_LICENSE_HEADER_END@
  22  */
  23
  24 #include <sys/cdefs.h>
  25 #include <stdint.h>
  26 #include <strings.h>
  27 #include <float.h>
  28 #include <math.h>
  29
  30 #include "fpmath.h"
  31
  32 #define BITS64          64
  33 #define DBL_BIAS        (DBL_MAX_EXP - 1)
  34 #define DBL_SUBEXP      (-DBL_BIAS + 1)
  35 #define LL_BITS         (8 * sizeof(int64_t))
  36 #define LL_HIGHBIT      (1LL << 63)
  37
  38 __private_extern__ int
  39 _ldbl2array32dd(union IEEEl2bits u, uint32_t *a)
  40 {
  41         int bit, shift, highbit, dexp;
  42         uint64_t a64[2];
  43         int64_t t64;
  44         int extrabit = 0;
  45
  46         if(u.d[0] == 0.0) {
  47             a[0] = a[1] = a[2] = a[3] = 0;
  48             return 0;
  49         }
  50
  51         bzero(a64, sizeof(a64));
  52
  53         switch (__fpclassifyd(u.d[0])) {
  54         case FP_NORMAL:
  55                 /*
  56                  * special case: if the head double only has the high (hidden)
  57                  * bit set, and the tail double is non-zero and is opposite
  58                  * in sign, then we increment extrabit to keep 106 bit
  59                  * precision in the results.
  60                  */
  61                 if(u.bits.manh == 0 && u.d[1] != 0 && u.bits.sign != u.bits.sign2)
  62                     extrabit++;
  63                 a64[1] = (1LL << (LDBL_MANT_DIG - BITS64 - 1 + extrabit));
  64                 a64[1] |= ((uint64_t)u.bits.manh >> (BITS64 - LDBL_MANL_SIZE - extrabit));
  65                 a64[0] = ((uint64_t)u.bits.manh << (LDBL_MANL_SIZE + extrabit));
  66                 break;
  67         case FP_SUBNORMAL:
  68                 a64[1] |= ((uint64_t)u.bits.manh >> (BITS64 - LDBL_MANL_SIZE));
  69                 a64[0] = ((uint64_t)u.bits.manh << LDBL_MANL_SIZE);
  70                 /* the tail double will be zero, so we are done */
  71                 goto done;
  72         default:
  73                 goto done;
  74         }
  75
  76         dexp = (int)u.bits.exp - (int)u.bits.exp2;
  77         /*
  78          * if the tail double is so small to not fit in LDBL_MANT_DIG bits,
  79          * then just skip it.
  80          */
  81         if (dexp >= LDBL_MANT_DIG + extrabit) {
  82                 reshift:
  83                 if (extrabit) {
  84                         bit = a64[1] & 1;
  85                         a64[1] >>= 1;
  86                         a64[0] >>= 1;
  87                         a64[0] |= ((uint64_t)bit) << (BITS64 - 1);
  88                         extrabit = 0;
  89                 }
  90                 goto done;
  91         }
  92
  93         switch (__fpclassifyd(u.d[1])) {
  94         case FP_NORMAL:
  95                 bit = LDBL_MANT_DIG - dexp - 1 + extrabit;
  96                 t64 = (1LL << bit);
  97                 break;
  98         case FP_SUBNORMAL:
  99                 bit = LDBL_MANT_DIG - (int)u.bits.exp + extrabit;
 100                 t64 = 0;
 101                 break;
 102         default:
 103                 /* should never get here */
 104                 goto reshift;
 105         }
 106         shift = LDBL_MANL_SIZE - bit - 1;
 107         if (shift >= 0)
 108                 t64 |= (u.bits.manl >> shift);
 109         else
 110                 t64 |= (u.bits.manl << (-shift));
 111         highbit = ((a64[0] & LL_HIGHBIT) != 0);
 112         if (u.bits.sign == u.bits.sign2) {
 113                 a64[0] += t64;
 114                 if (highbit && !(a64[0] & LL_HIGHBIT))  /* carry */
 115                         a64[1]++;
 116         } else {
 117                 a64[0] -= t64;
 118                 if (!highbit && (a64[0] & LL_HIGHBIT))  /* borrow */
 119                         a64[1]--;
 120         }
 121
 122         done:
 123         a[0] = (uint32_t)a64[0];
 124         a[1] = (uint32_t)(a64[0] >> 32);
 125         a[2] = (uint32_t)a64[1];
 126         a[3] = (uint32_t)(a64[1] >> 32);
 127         return extrabit;
 128 }