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

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2019 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 <https://www.gnu.org/licenses/>.
 */
/********************************************************************/
/*                                                                  */
/* MODULE_NAME: dosincos.c                                          */
/*                                                                  */
/*                                                                  */
/* FUNCTIONS:   dubsin                                              */
/*              dubcos                                              */
/*              docos                                               */
/* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h                 */
/*               sincos.tbl                                         */
/*                                                                  */
/* Routines compute sin() and cos() as Double-Length numbers         */
/********************************************************************/


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

#ifndef SECTION
# define SECTION
#endif

extern const union
{
  int4 i[880];
  double x[440];
} __sincostab attribute_hidden;

/***********************************************************************/
/* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
/* as Double-Length number and store it at array v .It computes it by  */
/* arithmetic action on Double-Length numbers                          */
/*(x+dx) between 0 and PI/4                                            */
/***********************************************************************/

void
SECTION
__dubsin (double x, double dx, double v[])
{
  double r, s, c, cc, d, dd, d2, dd2, e, ee,
	 sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
  double p, hx, tx, hy, ty, q;
#endif
  mynumber u;
  int4 k;

  u.x = x + big.x;
  k = u.i[LOW_HALF] << 2;
  x = x - (u.x - big.x);
  d = x + dx;
  dd = (x - d) + dx;
  /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
  MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
  sn = __sincostab.x[k];       /*                                  */
  ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
  cs = __sincostab.x[k + 2];   /*                                  */
  ccs = __sincostab.x[k + 3];  /*                                  */
  /* Taylor series for sin ds=sin(t) */
  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, d, dd, ds, dss, r, s);

  /* Taylor series for cos dc=cos(t) */
  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);

  MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  SUB2 (e, ee, dc, dcc, e, ee, r, s);
  ADD2 (e, ee, sn, ssn, e, ee, r, s);                    /* e+ee=sin(x+dx) */

  v[0] = e;
  v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
/* as Double-Length number and store it in array v .It computes it by */
/* arithmetic action on Double-Length numbers                         */
/*(x+dx) between 0 and PI/4                                           */
/**********************************************************************/

void
SECTION
__dubcos (double x, double dx, double v[])
{
  double r, s, c, cc, d, dd, d2, dd2, e, ee,
	 sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
  double p, hx, tx, hy, ty, q;
#endif
  mynumber u;
  int4 k;
  u.x = x + big.x;
  k = u.i[LOW_HALF] << 2;
  x = x - (u.x - big.x);
  d = x + dx;
  dd = (x - d) + dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
  MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
  sn = __sincostab.x[k];     /*                                  */
  ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
  cs = __sincostab.x[k + 2];   /*                                  */
  ccs = __sincostab.x[k + 3];  /*                                  */
  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, d, dd, ds, dss, r, s);

  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);

  MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);

  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
  MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (ds, dss, d, dd, ds, dss, r, s);
  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
  MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
  MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  ADD2 (e, ee, dc, dcc, e, ee, r, s);
  SUB2 (cs, ccs, e, ee, e, ee, r, s);

  v[0] = e;
  v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
/* as Double-Length number and store it in array v                    */
/**********************************************************************/
void
SECTION
__docos (double x, double dx, double v[])
{
  double y, yy, p, w[2];
  if (x > 0)
    {
      y = x; yy = dx;
    }
  else
    {
      y = -x; yy = -dx;
    }
  if (y < 0.5 * hp0.x)                                 /*  y< PI/4    */
    {
      __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
    }
  else if (y < 1.5 * hp0.x)                        /* y< 3/4 * PI */
    {
      p = hp0.x - y; /* p = PI/2 - y */
      yy = hp1.x - yy;
      y = p + yy;
      yy = (p - y) + yy;
      if (y > 0)
	{
	  __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
	}
      /* cos(x) = sin ( 90 -  x ) */
      else
	{
	  __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
	}
    }
  else   /* y>= 3/4 * PI */
    {
      p = 2.0 * hp0.x - y; /* p = PI- y */
      yy = 2.0 * hp1.x - yy;
      y = p + yy;
      yy = (p - y) + yy;
      __dubcos (y, yy, w);
      v[0] = -w[0];
      v[1] = -w[1];
    }
}
Commit	Line	Data
e4d82761 UD	1	/*
e4d82761 UD	2	* IBM Accurate Mathematical Library
aeb25823	3	* written by International Business Machines Corp.
04277e02	4	* Copyright (C) 2001-2019 Free Software Foundation, Inc.
e4d82761 UD	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
cc7375ce	8	* the Free Software Foundation; either version 2.1 of the License, or
e4d82761	9	* (at your option) any later version.
50944bca	10	*
e4d82761 UD	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
c6c6dd48	14	* GNU Lesser General Public License for more details.
e4d82761 UD	15	*
e4d82761 UD	16	* You should have received a copy of the GNU Lesser General Public License
5a82c748	17	* along with this program; if not, see <https://www.gnu.org/licenses/>.
e4d82761 UD	18	*/
e4d82761 UD	19	/********************************************************************/
50944bca	20	/* */
e4d82761 UD	21	/* MODULE_NAME: dosincos.c */
	22	/* */
	23	/* */
	24	/* FUNCTIONS: dubsin */
	25	/* dubcos */
	26	/* docos */
	27	/* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */
	28	/* sincos.tbl */
	29	/* */
	30	/* Routines compute sin() and cos() as Double-Length numbers */
	31	/********************************************************************/
	32
	33
	34
	35	#include "endian.h"
50944bca	36	#include "mydefs.h"
c8b3296b	37	#include <dla.h>
e4d82761	38	#include "dosincos.h"
1ed0291c	39	#include <math_private.h>
15b3c029	40
31d3cc00 UD	41	#ifndef SECTION
	42	# define SECTION
	43	#endif
	44
af968f62 UD	45	extern const union
	46	{
	47	int4 i[880];
	48	double x[440];
	49	} __sincostab attribute_hidden;
	50
e4d82761 UD	51	/***********************************************************************/
	52	/* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
	53	/* as Double-Length number and store it at array v .It computes it by */
	54	/* arithmetic action on Double-Length numbers */
50944bca	55	/(x+dx) between 0 and PI/4 /
e4d82761 UD	56	/***********************************************************************/
e4d82761 UD	57
31d3cc00 UD	58	void
31d3cc00 UD	59	SECTION
c5d5d574 OB	60	__dubsin (double x, double dx, double v[])
	61	{
	62	double r, s, c, cc, d, dd, d2, dd2, e, ee,
	63	sn, ssn, cs, ccs, ds, dss, dc, dcc;
58985aa9	64	#ifndef DLA_FMS
c5d5d574	65	double p, hx, tx, hy, ty, q;
50944bca	66	#endif
e4d82761 UD	67	mynumber u;
e4d82761 UD	68	int4 k;
50944bca	69
c5d5d574 OB	70	u.x = x + big.x;
	71	k = u.i[LOW_HALF] << 2;
	72	x = x - (u.x - big.x);
	73	d = x + dx;
	74	dd = (x - d) + dx;
	75	/* sin(x+dx)=sin(Xi+t)=sin(Xi)cos(t) + cos(Xi)sin(t) where t ->0 /
	76	MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
	77	sn = __sincostab.x[k]; /* */
	78	ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
	79	cs = __sincostab.x[k + 2]; /* */
	80	ccs = __sincostab.x[k + 3]; /* */
	81	/* Taylor series for sin ds=sin(t) */
	82	MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	83	ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
	84	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	85	ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
	86	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	87	MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	88	ADD2 (ds, dss, d, dd, ds, dss, r, s);
	89
	90	/* Taylor series for cos dc=cos(t) */
	91	MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	92	ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
	93	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	94	ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
	95	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	96	ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
	97	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	98
	99	MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
	100	MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	101	SUB2 (e, ee, dc, dcc, e, ee, r, s);
	102	ADD2 (e, ee, sn, ssn, e, ee, r, s); /* e+ee=sin(x+dx) */
	103
	104	v[0] = e;
	105	v[1] = ee;
e4d82761 UD	106	}
	107	/**********************************************************************/
	108	/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
	109	/* as Double-Length number and store it in array v .It computes it by */
	110	/* arithmetic action on Double-Length numbers */
	111	/(x+dx) between 0 and PI/4 /
	112	/**********************************************************************/
	113
31d3cc00 UD	114	void
31d3cc00 UD	115	SECTION
c5d5d574 OB	116	__dubcos (double x, double dx, double v[])
	117	{
	118	double r, s, c, cc, d, dd, d2, dd2, e, ee,
	119	sn, ssn, cs, ccs, ds, dss, dc, dcc;
58985aa9	120	#ifndef DLA_FMS
c5d5d574	121	double p, hx, tx, hy, ty, q;
a1a87169	122	#endif
e4d82761	123	mynumber u;
50944bca	124	int4 k;
c5d5d574 OB	125	u.x = x + big.x;
	126	k = u.i[LOW_HALF] << 2;
	127	x = x - (u.x - big.x);
	128	d = x + dx;
	129	dd = (x - d) + dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
	130	MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
	131	sn = __sincostab.x[k]; /* */
	132	ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
	133	cs = __sincostab.x[k + 2]; /* */
	134	ccs = __sincostab.x[k + 3]; /* */
	135	MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	136	ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
	137	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	138	ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
	139	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	140	MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	141	ADD2 (ds, dss, d, dd, ds, dss, r, s);
	142
	143	MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	144	ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
	145	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	146	ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
	147	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	148	ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
	149	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	150
	151	MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
	152	MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	153
	154	MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	155	ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
	156	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	157	ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
	158	MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	159	MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
	160	ADD2 (ds, dss, d, dd, ds, dss, r, s);
	161	MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	162	ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
	163	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	164	ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
	165	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	166	ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
	167	MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	168	MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
	169	MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
	170	ADD2 (e, ee, dc, dcc, e, ee, r, s);
	171	SUB2 (cs, ccs, e, ee, e, ee, r, s);
	172
	173	v[0] = e;
	174	v[1] = ee;
e4d82761 UD	175	}
	176	/**********************************************************************/
	177	/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
	178	/* as Double-Length number and store it in array v */
	179	/**********************************************************************/
31d3cc00 UD	180	void
31d3cc00 UD	181	SECTION
c5d5d574 OB	182	__docos (double x, double dx, double v[])
	183	{
	184	double y, yy, p, w[2];
	185	if (x > 0)
	186	{
	187	y = x; yy = dx;
	188	}
	189	else
	190	{
	191	y = -x; yy = -dx;
	192	}
	193	if (y < 0.5 * hp0.x) /* y< PI/4 */
	194	{
	195	__dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
	196	}
	197	else if (y < 1.5 * hp0.x) /* y< 3/4 * PI */
	198	{
	199	p = hp0.x - y; /* p = PI/2 - y */
	200	yy = hp1.x - yy;
	201	y = p + yy;
	202	yy = (p - y) + yy;
	203	if (y > 0)
	204	{
	205	__dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
	206	}
	207	/* cos(x) = sin ( 90 - x ) */
	208	else
	209	{
	210	__dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
	211	}
	212	}
	213	else /* y>= 3/4 * PI */
	214	{
	215	p = 2.0 * hp0.x - y; /* p = PI- y */
	216	yy = 2.0 * hp1.x - yy;
	217	y = p + yy;
	218	yy = (p - y) + yy;
	219	__dubcos (y, yy, w);
	220	v[0] = -w[0];
	221	v[1] = -w[1];
	222	}
e4d82761	223	}