Consolidate multiple precision sin/cos functions

[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / sincos32.c

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2013 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: sincos32.c                                     */
/*                                                              */
/*  FUNCTIONS: ss32                                             */
/*             cc32                                             */
/*             c32                                              */
/*             sin32                                            */
/*             cos32                                            */
/*             mpsin                                            */
/*             mpcos                                            */
/*             mpranred                                         */
/*             mpsin1                                           */
/*             mpcos1                                           */
/*                                                              */
/* FILES NEEDED: endian.h mpa.h sincos32.h                      */
/*               mpa.c                                          */
/*                                                              */
/* Multi Precision sin() and cos() function with p=32  for sin()*/
/* cos() arcsin() and arccos() routines                         */
/* In addition mpranred() routine  performs range  reduction of */
/* a double number x into multi precision number   y,           */
/* such that y=x-n*pi/2, abs(y)<pi/4,  n=0,+-1,+-2,....         */
/****************************************************************/
#include "endian.h"
#include "mpa.h"
#include "sincos32.h"
#include <math_private.h>

#ifndef SECTION
# define SECTION
#endif

/* Compute Multi-Precision sin() function for given p.  Receive Multi Precision
   number x and result stored at y.  */
static void
SECTION
ss32 (mp_no *x, mp_no *y, int p)
{
  int i;
  double a;
  mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
  for (i = 1; i <= p; i++)
    mpk.d[i] = 0;

  __sqr (x, &x2, p);
  __cpy (&oofac27, &gor, p);
  __cpy (&gor, &sum, p);
  for (a = 27.0; a > 1.0; a -= 2.0)
    {
      mpk.d[1] = a * (a - 1.0);
      __mul (&gor, &mpk, &mpt1, p);
      __cpy (&mpt1, &gor, p);
      __mul (&x2, &sum, &mpt1, p);
      __sub (&gor, &mpt1, &sum, p);
    }
  __mul (x, &sum, y, p);
}

/* Compute Multi-Precision cos() function for given p. Receive Multi Precision
   number x and result stored at y.  */
static void
SECTION
cc32 (mp_no *x, mp_no *y, int p)
{
  int i;
  double a;
  mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
  for (i = 1; i <= p; i++)
    mpk.d[i] = 0;

  __sqr (x, &x2, p);
  mpk.d[1] = 27.0;
  __mul (&oofac27, &mpk, &gor, p);
  __cpy (&gor, &sum, p);
  for (a = 26.0; a > 2.0; a -= 2.0)
    {
      mpk.d[1] = a * (a - 1.0);
      __mul (&gor, &mpk, &mpt1, p);
      __cpy (&mpt1, &gor, p);
      __mul (&x2, &sum, &mpt1, p);
      __sub (&gor, &mpt1, &sum, p);
    }
  __mul (&x2, &sum, y, p);
}

/* Compute both sin(x), cos(x) as Multi precision numbers.  */
void
SECTION
__c32 (mp_no *x, mp_no *y, mp_no *z, int p)
{
  mp_no u, t, t1, t2, c, s;
  int i;
  __cpy (x, &u, p);
  u.e = u.e - 1;
  cc32 (&u, &c, p);
  ss32 (&u, &s, p);
  for (i = 0; i < 24; i++)
    {
      __mul (&c, &s, &t, p);
      __sub (&s, &t, &t1, p);
      __add (&t1, &t1, &s, p);
      __sub (&mptwo, &c, &t1, p);
      __mul (&t1, &c, &t2, p);
      __add (&t2, &t2, &c, p);
    }
  __sub (&mpone, &c, y, p);
  __cpy (&s, z, p);
}

/* Receive double x and two double results of sin(x) and return result which is
   more accurate, computing sin(x) with multi precision routine c32.  */
double
SECTION
__sin32 (double x, double res, double res1)
{
  int p;
  mp_no a, b, c;
  p = 32;
  __dbl_mp (res, &a, p);
  __dbl_mp (0.5 * (res1 - res), &b, p);
  __add (&a, &b, &c, p);
  if (x > 0.8)
    {
      __sub (&hp, &c, &a, p);
      __c32 (&a, &b, &c, p);
    }
  else
    __c32 (&c, &a, &b, p);	/* b=sin(0.5*(res+res1))  */
  __dbl_mp (x, &c, p);		/* c = x  */
  __sub (&b, &c, &a, p);
  /* if a > 0 return min (res, res1), otherwise return max (res, res1).  */
  if (a.d[0] > 0)
    return (res < res1) ? res : res1;
  else
    return (res > res1) ? res : res1;
}

/* Receive double x and two double results of cos(x) and return result which is
   more accurate, computing cos(x) with multi precision routine c32.  */
double
SECTION
__cos32 (double x, double res, double res1)
{
  int p;
  mp_no a, b, c;
  p = 32;
  __dbl_mp (res, &a, p);
  __dbl_mp (0.5 * (res1 - res), &b, p);
  __add (&a, &b, &c, p);
  if (x > 2.4)
    {
      __sub (&pi, &c, &a, p);
      __c32 (&a, &b, &c, p);
      b.d[0] = -b.d[0];
    }
  else if (x > 0.8)
    {
      __sub (&hp, &c, &a, p);
      __c32 (&a, &c, &b, p);
    }
  else
    __c32 (&c, &b, &a, p);	/* b=cos(0.5*(res+res1))  */
  __dbl_mp (x, &c, p);		/* c = x                  */
  __sub (&b, &c, &a, p);
  /* if a > 0 return max (res, res1), otherwise return min (res, res1).  */
  if (a.d[0] > 0)
    return (res > res1) ? res : res1;
  else
    return (res < res1) ? res : res1;
}

/* Compute sin() of double-length number (X + DX) as Multi Precision number and
   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
   original input and DX is ignored.  */
double
SECTION
__mpsin (double x, double dx, bool reduce_range)
{
  double y;
  mp_no a, b, c, s;
  int n;
  int p = 32;

  if (reduce_range)
    {
      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
      __c32 (&a, &c, &s, p);
    }
  else
    {
      n = -1;
      __dbl_mp (x, &b, p);
      __dbl_mp (dx, &c, p);
      __add (&b, &c, &a, p);
      if (x > 0.8)
        {
          __sub (&hp, &a, &b, p);
          __c32 (&b, &s, &c, p);
        }
      else
        __c32 (&a, &c, &s, p);	/* b = sin(x+dx)  */
    }

  /* Convert result based on which quarter of unit circle y is in.  */
  switch (n)
    {
    case 1:
      __mp_dbl (&c, &y, p);
      break;

    case 3:
      __mp_dbl (&c, &y, p);
      y = -y;
      break;

    case 2:
      __mp_dbl (&s, &y, p);
      y = -y;
      break;

    /* Quadrant not set, so the result must be sin (X + DX), which is also in
       S.  */
    case 0:
    default:
      __mp_dbl (&s, &y, p);
    }
  return y;
}

/* Compute cos() of double-length number (X + DX) as Multi Precision number and
   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
   original input and DX is ignored.  */
double
SECTION
__mpcos (double x, double dx, bool reduce_range)
{
  double y;
  mp_no a, b, c, s;
  int n;
  int p = 32;

  if (reduce_range)
    {
      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
      __c32 (&a, &c, &s, p);
    }
  else
    {
      n = -1;
      __dbl_mp (x, &b, p);
      __dbl_mp (dx, &c, p);
      __add (&b, &c, &a, p);
      if (x > 0.8)
        {
          __sub (&hp, &a, &b, p);
          __c32 (&b, &s, &c, p);
        }
      else
        __c32 (&a, &c, &s, p);	/* a = cos(x+dx)     */
    }

  /* Convert result based on which quarter of unit circle y is in.  */
  switch (n)
    {
    case 1:
      __mp_dbl (&s, &y, p);
      y = -y;
      break;

    case 3:
      __mp_dbl (&s, &y, p);
      break;

    case 2:
      __mp_dbl (&c, &y, p);
      y = -y;
      break;

    /* Quadrant not set, so the result must be cos (X + DX), which is also
       stored in C.  */
    case 0:
    default:
      __mp_dbl (&c, &y, p);
    }
  return y;
}

/* Perform range reduction of a double number x into multi precision number y,
   such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ...
   Return int which indicates in which quarter of circle x is.  */
int
SECTION
__mpranred (double x, mp_no *y, int p)
{
  number v;
  double t, xn;
  int i, k, n;
  mp_no a, b, c;

  if (ABS (x) < 2.8e14)
    {
      t = (x * hpinv.d + toint.d);
      xn = t - toint.d;
      v.d = t;
      n = v.i[LOW_HALF] & 3;
      __dbl_mp (xn, &a, p);
      __mul (&a, &hp, &b, p);
      __dbl_mp (x, &c, p);
      __sub (&c, &b, y, p);
      return n;
    }
  else
    {
      /* If x is very big more precision required.  */
      __dbl_mp (x, &a, p);
      a.d[0] = 1.0;
      k = a.e - 5;
      if (k < 0)
	k = 0;
      b.e = -k;
      b.d[0] = 1.0;
      for (i = 0; i < p; i++)
	b.d[i + 1] = toverp[i + k];
      __mul (&a, &b, &c, p);
      t = c.d[c.e];
      for (i = 1; i <= p - c.e; i++)
	c.d[i] = c.d[i + c.e];
      for (i = p + 1 - c.e; i <= p; i++)
	c.d[i] = 0;
      c.e = 0;
      if (c.d[1] >= HALFRAD)
	{
	  t += 1.0;
	  __sub (&c, &mpone, &b, p);
	  __mul (&b, &hp, y, p);
	}
      else
	__mul (&c, &hp, y, p);
      n = (int) t;
      if (x < 0)
	{
	  y->d[0] = -y->d[0];
	  n = -n;
	}
      return (n & 3);
    }
}