sysdeps/ieee754/dbl-64/halfulp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2015 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /************************************************************************/
  20 /*                                                                      */
  21 /* MODULE_NAME:halfulp.c                                                */
  22 /*                                                                      */
  23 /*  FUNCTIONS:halfulp                                                   */
  24 /*  FILES NEEDED: mydefs.h dla.h endian.h                               */
  25 /*                uroot.c                                               */
  26 /*                                                                      */
  27 /*Routine halfulp(double x, double y) computes x^y where result does    */
  28 /*not need rounding. If the result is closer to 0 than can be           */
  29 /*represented it returns 0.                                             */
  30 /*     In the following cases the function does not compute anything    */
  31 /*and returns a negative number:                                        */
  32 /*1. if the result needs rounding,                                      */
  33 /*2. if y is outside the interval [0,  2^20-1],                         */
  34 /*3. if x can be represented by  x=2**n for some integer n.             */
  35 /************************************************************************/
  36
  37 #include "endian.h"
  38 #include "mydefs.h"
  39 #include <dla.h>
  40 #include <math_private.h>
  41
  42 #ifndef SECTION
  43 # define SECTION
  44 #endif
  45
  46 static const int4 tab54[32] = {
  47   262143, 11585,  1782, 511, 210, 107, 63, 42,
  48   30,     22,     17,   14,  12,  10,  9, 7,
  49   7,      6,      5,    5,   5,   4,   4, 4,
  50   3,      3,      3,    3,   3,   3,   3, 3
  51 };
  52
  53
  54 double
  55 SECTION
  56 __halfulp (double x, double y)
  57 {
  58   mynumber v;
  59   double z, u, uu;
  60 #ifndef DLA_FMS
  61   double j1, j2, j3, j4, j5;
  62 #endif
  63   int4 k, l, m, n;
  64   if (y <= 0)                 /*if power is negative or zero */
  65     {
  66       v.x = y;
  67       if (v.i[LOW_HALF] != 0)
  68         return -10.0;
  69       v.x = x;
  70       if (v.i[LOW_HALF] != 0)
  71         return -10.0;
  72       if ((v.i[HIGH_HALF] & 0x000fffff) != 0)
  73         return -10;                                     /* if x =2 ^ n */
  74       k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */
  75       z = (double) k;
  76       return (z * y == -1075.0) ? 0 : -10.0;
  77     }
  78   /* if y > 0  */
  79   v.x = y;
  80   if (v.i[LOW_HALF] != 0)
  81     return -10.0;
  82
  83   v.x = x;
  84   /*  case where x = 2**n for some integer n */
  85   if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0)
  86     {
  87       k = (v.i[HIGH_HALF] >> 20) - 1023;
  88       return (((double) k) * y == -1075.0) ? 0 : -10.0;
  89     }
  90
  91   v.x = y;
  92   k = v.i[HIGH_HALF];
  93   m = k << 12;
  94   l = 0;
  95   while (m)
  96     {
  97       m = m << 1; l++;
  98     }
  99   n = (k & 0x000fffff) | 0x00100000;
 100   n = n >> (20 - l);                       /*   n is the odd integer of y    */
 101   k = ((k >> 20) - 1023) - l;               /*   y = n*2**k                   */
 102   if (k > 5)
 103     return -10.0;
 104   if (k > 0)
 105     for (; k > 0; k--)
 106       n *= 2;
 107   if (n > 34)
 108     return -10.0;
 109   k = -k;
 110   if (k > 5)
 111     return -10.0;
 112
 113   /*   now treat x        */
 114   while (k > 0)
 115     {
 116       z = __ieee754_sqrt (x);
 117       EMULV (z, z, u, uu, j1, j2, j3, j4, j5);
 118       if (((u - x) + uu) != 0)
 119         break;
 120       x = z;
 121       k--;
 122     }
 123   if (k)
 124     return -10.0;
 125
 126   /* it is impossible that n == 2,  so the mantissa of x must be short  */
 127
 128   v.x = x;
 129   if (v.i[LOW_HALF])
 130     return -10.0;
 131   k = v.i[HIGH_HALF];
 132   m = k << 12;
 133   l = 0;
 134   while (m)
 135     {
 136       m = m << 1; l++;
 137     }
 138   m = (k & 0x000fffff) | 0x00100000;
 139   m = m >> (20 - l);                       /*   m is the odd integer of x    */
 140
 141   /*   now check whether the length of m**n is at most 54 bits */
 142
 143   if (m > tab54[n - 3])
 144     return -10.0;
 145
 146   /* yes, it is - now compute x**n by simple multiplications  */
 147
 148   u = x;
 149   for (k = 1; k < n; k++)
 150     u = u * x;
 151   return u;
 152 }