sysdeps/ieee754/dbl-64/halfulp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2005, 2011 Free Software Foundation
   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, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  19  */
  20 /************************************************************************/
  21 /*                                                                      */
  22 /* MODULE_NAME:halfulp.c                                                */
  23 /*                                                                      */
  24 /*  FUNCTIONS:halfulp                                                   */
  25 /*  FILES NEEDED: mydefs.h dla.h endian.h                               */
  26 /*                uroot.c                                               */
  27 /*                                                                      */
  28 /*Routine halfulp(double x, double y) computes x^y where result does    */
  29 /*not need rounding. If the result is closer to 0 than can be           */
  30 /*represented it returns 0.                                             */
  31 /*     In the following cases the function does not compute anything    */
  32 /*and returns a negative number:                                        */
  33 /*1. if the result needs rounding,                                      */
  34 /*2. if y is outside the interval [0,  2^20-1],                         */
  35 /*3. if x can be represented by  x=2**n for some integer n.             */
  36 /************************************************************************/
  37
  38 #include "endian.h"
  39 #include "mydefs.h"
  40 #include <dla.h>
  41 #include "math_private.h"
  42
  43 #ifndef SECTION
  44 # define SECTION
  45 #endif
  46
  47 static const int4 tab54[32] = {
  48    262143, 11585, 1782, 511, 210, 107, 63, 42,
  49        30,    22,   17,  14,  12,  10,  9,  7,
  50         7,     6,    5,   5,   5,   4,  4,  4,
  51         3,     3,    3,   3,   3,   3,  3,  3 };
  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     v.x = y;
  66     if (v.i[LOW_HALF] != 0) return -10.0;
  67     v.x = x;
  68     if (v.i[LOW_HALF] != 0) return -10.0;
  69     if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10;   /* if x =2 ^ n */
  70     k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023;         /* find this n */
  71     z = (double) k;
  72     return (z*y == -1075.0)?0: -10.0;
  73   }
  74                               /* if y > 0  */
  75   v.x = y;
  76     if (v.i[LOW_HALF] != 0) return -10.0;
  77
  78   v.x=x;
  79                               /*  case where x = 2**n for some integer n */
  80   if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
  81     k=(v.i[HIGH_HALF]>>20)-1023;
  82     return (((double) k)*y == -1075.0)?0:-10.0;
  83   }
  84
  85   v.x = y;
  86   k = v.i[HIGH_HALF];
  87   m = k<<12;
  88   l = 0;
  89   while (m)
  90     {m = m<<1; l++; }
  91   n = (k&0x000fffff)|0x00100000;
  92   n = n>>(20-l);                       /*   n is the odd integer of y    */
  93   k = ((k>>20) -1023)-l;               /*   y = n*2**k                   */
  94   if (k>5) return -10.0;
  95   if (k>0) for (;k>0;k--) n *= 2;
  96   if (n > 34) return -10.0;
  97   k = -k;
  98   if (k>5) return -10.0;
  99
 100                             /*   now treat x        */
 101   while (k>0) {
 102     z = __ieee754_sqrt(x);
 103     EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
 104     if (((u-x)+uu) != 0) break;
 105     x = z;
 106     k--;
 107  }
 108   if (k) return -10.0;
 109
 110   /* it is impossible that n == 2,  so the mantissa of x must be short  */
 111
 112   v.x = x;
 113   if (v.i[LOW_HALF]) return -10.0;
 114   k = v.i[HIGH_HALF];
 115   m = k<<12;
 116   l = 0;
 117   while (m) {m = m<<1; l++; }
 118   m = (k&0x000fffff)|0x00100000;
 119   m = m>>(20-l);                       /*   m is the odd integer of x    */
 120
 121             /*   now check whether the length of m**n is at most 54 bits */
 122
 123   if  (m > tab54[n-3]) return -10.0;
 124
 125              /* yes, it is - now compute x**n by simple multiplications  */
 126
 127   u = x;
 128   for (k=1;k<n;k++) u = u*x;
 129   return u;
 130 }