glibc/sysdeps/ieee754/dbl-64/halfulp.c

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, see <http://www.gnu.org/licenses/>.
 */
/************************************************************************/
/*                                                                      */
/* MODULE_NAME:halfulp.c                                                */
/*                                                                      */
/*  FUNCTIONS:halfulp                                                   */
/*  FILES NEEDED: mydefs.h dla.h endian.h                               */
/*                uroot.c                                               */
/*                                                                      */
/*Routine halfulp(double x, double y) computes x^y where result does    */
/*not need rounding. If the result is closer to 0 than can be           */
/*represented it returns 0.                                             */
/*     In the following cases the function does not compute anything    */
/*and returns a negative number:                                        */
/*1. if the result needs rounding,                                      */
/*2. if y is outside the interval [0,  2^20-1],                         */
/*3. if x can be represented by  x=2**n for some integer n.             */
/************************************************************************/

#include "endian.h"
#include "mydefs.h"
#include <dla.h>
#include <math_private.h>

#ifndef SECTION
# define SECTION
#endif

static const int4 tab54[32] = {
  262143, 11585,  1782, 511, 210, 107, 63, 42,
  30,     22,     17,   14,  12,  10,  9, 7,
  7,      6,      5,    5,   5,   4,   4, 4,
  3,      3,      3,    3,   3,   3,   3, 3
};


double
SECTION
__halfulp (double x, double y)
{
  mynumber v;
  double z, u, uu;
#ifndef DLA_FMS
  double j1, j2, j3, j4, j5;
#endif
  int4 k, l, m, n;
  if (y <= 0)                 /*if power is negative or zero */
    {
      v.x = y;
      if (v.i[LOW_HALF] != 0)
	return -10.0;
      v.x = x;
      if (v.i[LOW_HALF] != 0)
	return -10.0;
      if ((v.i[HIGH_HALF] & 0x000fffff) != 0)
	return -10;                                     /* if x =2 ^ n */
      k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */
      z = (double) k;
      return (z * y == -1075.0) ? 0 : -10.0;
    }
  /* if y > 0  */
  v.x = y;
  if (v.i[LOW_HALF] != 0)
    return -10.0;

  v.x = x;
  /*  case where x = 2**n for some integer n */
  if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0)
    {
      k = (v.i[HIGH_HALF] >> 20) - 1023;
      return (((double) k) * y == -1075.0) ? 0 : -10.0;
    }

  v.x = y;
  k = v.i[HIGH_HALF];
  m = k << 12;
  l = 0;
  while (m)
    {
      m = m << 1; l++;
    }
  n = (k & 0x000fffff) | 0x00100000;
  n = n >> (20 - l);                       /*   n is the odd integer of y    */
  k = ((k >> 20) - 1023) - l;               /*   y = n*2**k                   */
  if (k > 5)
    return -10.0;
  if (k > 0)
    for (; k > 0; k--)
      n *= 2;
  if (n > 34)
    return -10.0;
  k = -k;
  if (k > 5)
    return -10.0;

  /*   now treat x        */
  while (k > 0)
    {
      z = __ieee754_sqrt (x);
      EMULV (z, z, u, uu, j1, j2, j3, j4, j5);
      if (((u - x) + uu) != 0)
	break;
      x = z;
      k--;
    }
  if (k)
    return -10.0;

  /* it is impossible that n == 2,  so the mantissa of x must be short  */

  v.x = x;
  if (v.i[LOW_HALF])
    return -10.0;
  k = v.i[HIGH_HALF];
  m = k << 12;
  l = 0;
  while (m)
    {
      m = m << 1; l++;
    }
  m = (k & 0x000fffff) | 0x00100000;
  m = m >> (20 - l);                       /*   m is the odd integer of x    */

  /*   now check whether the length of m**n is at most 54 bits */

  if (m > tab54[n - 3])
    return -10.0;

  /* yes, it is - now compute x**n by simple multiplications  */

  u = x;
  for (k = 1; k < n; k++)
    u = u * x;
  return u;
}