s_sin.c

Glibc 2.3.2源代码(解压后有100多M)

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001 Free Software Foundation
 *
 * 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
/****************************************************************************/
/*                                                                          */
/* MODULE_NAME:usncs.c                                                      */
/*                                                                          */
/* FUNCTIONS: usin                                                          */
/*            ucos                                                          */
/*            slow                                                          */
/*            slow1                                                         */
/*            slow2                                                         */
/*            sloww                                                         */
/*            sloww1                                                        */
/*            sloww2                                                        */
/*            bsloww                                                        */
/*            bsloww1                                                       */
/*            bsloww2                                                       */
/*            cslow2                                                        */
/*            csloww                                                        */
/*            csloww1                                                       */
/*            csloww2                                                       */
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h                     */
/*               branred.c sincos32.c dosincos.c mpa.c                      */
/*               sincos.tbl                                                 */
/*                                                                          */
/* An ultimate sin and  routine. Given an IEEE double machine number x       */
/* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
/* Assumption: Machine arithmetic operations are performed in               */
/* round to nearest mode of IEEE 754 standard.                              */
/*                                                                          */
/****************************************************************************/


#include "endian.h"
#include "mydefs.h"
#include "usncs.h"
#include "MathLib.h"
#include "sincos.tbl"
#include "math_private.h"

static const double
          sn3 = -1.66666666666664880952546298448555E-01,
          sn5 =  8.33333214285722277379541354343671E-03,
          cs2 =  4.99999999999999999999950396842453E-01,
          cs4 = -4.16666666666664434524222570944589E-02,
          cs6 =  1.38888874007937613028114285595617E-03;

void __dubsin(double x, double dx, double w[]);
void __docos(double x, double dx, double w[]);
double __mpsin(double x, double dx);
double __mpcos(double x, double dx);
double __mpsin1(double x);
double __mpcos1(double x);
static double slow(double x);
static double slow1(double x);
static double slow2(double x);
static double sloww(double x, double dx, double orig);
static double sloww1(double x, double dx, double orig);
static double sloww2(double x, double dx, double orig, int n);
static double bsloww(double x, double dx, double orig, int n);
static double bsloww1(double x, double dx, double orig, int n);
static double bsloww2(double x, double dx, double orig, int n);
int __branred(double x, double *a, double *aa);
static double cslow2(double x);
static double csloww(double x, double dx, double orig);
static double csloww1(double x, double dx, double orig);
static double csloww2(double x, double dx, double orig, int n);
/*******************************************************************/
/* An ultimate sin routine. Given an IEEE double machine number x   */
/* it computes the correctly rounded (to nearest) value of sin(x)  */
/*******************************************************************/
double __sin(double x){
	double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
#if 0
	double w[2];
#endif
	mynumber u,v;
	int4 k,m,n;
#if 0
	int4 nn;
#endif

	u.x = x;
	m = u.i[HIGH_HALF];
	k = 0x7fffffff&m;              /* no sign           */
	if (k < 0x3e500000)            /* if x->0 =>sin(x)=x */
	 return x;
 /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
	else  if (k < 0x3fd00000){
	  xx = x*x;
	  /*Taylor series */
	  t = ((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*(xx*x);
	  res = x+t;
	  cor = (x-res)+t;
	  return (res == res + 1.07*cor)? res : slow(x);
	}    /*  else  if (k < 0x3fd00000)    */
/*---------------------------- 0.25<|x|< 0.855469---------------------- */
	else if (k < 0x3feb6000)  {
	  u.x=(m>0)?big.x+x:big.x-x;
	  y=(m>0)?x-(u.x-big.x):x+(u.x-big.x);
	  xx=y*y;
	  s = y + y*xx*(sn3 +xx*sn5);
	  c = xx*(cs2 +xx*(cs4 + xx*cs6));
	  k=u.i[LOW_HALF]<<2;
	  sn=(m>0)?sincos.x[k]:-sincos.x[k];
	  ssn=(m>0)?sincos.x[k+1]:-sincos.x[k+1];
	  cs=sincos.x[k+2];
	  ccs=sincos.x[k+3];
	  cor=(ssn+s*ccs-sn*c)+cs*s;
	  res=sn+cor;
	  cor=(sn-res)+cor;
	  return (res==res+1.025*cor)? res : slow1(x);
	}    /*   else  if (k < 0x3feb6000)    */

/*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
	else if (k <  0x400368fd ) {

	  y = (m>0)? hp0.x-x:hp0.x+x;
	  if (y>=0) {
	    u.x = big.x+y;
	    y = (y-(u.x-big.x))+hp1.x;
	  }
	  else {
	    u.x = big.x-y;
	    y = (-hp1.x) - (y+(u.x-big.x));
	  }
	  xx=y*y;
	  s = y + y*xx*(sn3 +xx*sn5);
	  c = xx*(cs2 +xx*(cs4 + xx*cs6));
	  k=u.i[LOW_HALF]<<2;
	  sn=sincos.x[k];
	  ssn=sincos.x[k+1];
	  cs=sincos.x[k+2];
	  ccs=sincos.x[k+3];
	  cor=(ccs-s*ssn-cs*c)-sn*s;
	  res=cs+cor;
	  cor=(cs-res)+cor;
	  return (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x);
	} /*   else  if (k < 0x400368fd)    */

/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
	else if (k < 0x419921FB ) {
	  t = (x*hpinv.x + toint.x);
	  xn = t - toint.x;
	  v.x = t;
	  y = (x - xn*mp1.x) - xn*mp2.x;
	  n =v.i[LOW_HALF]&3;
	  da = xn*mp3.x;
	  a=y-da;
	  da = (y-a)-da;
	  eps = ABS(x)*1.2e-30;

	  switch (n) { /* quarter of unit circle */
	  case 0:
	  case 2:
	    xx = a*a;
	    if (n) {a=-a;da=-da;}
	    if (xx < 0.01588) {
                      /*Taylor series */
	      t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
	      res = a+t;
	      cor = (a-res)+t;
	      cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
	      return (res == res + cor)? res : sloww(a,da,x);
	    }
	    else  {
	      if (a>0)
		{m=1;t=a;db=da;}
	      else
		{m=0;t=-a;db=-da;}
	      u.x=big.x+t;
	      y=t-(u.x-big.x);
	      xx=y*y;
	      s = y + (db+y*xx*(sn3 +xx*sn5));
	      c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
	      k=u.i[LOW_HALF]<<2;
	      sn=sincos.x[k];
	      ssn=sincos.x[k+1];
	      cs=sincos.x[k+2];
	      ccs=sincos.x[k+3];
	      cor=(ssn+s*ccs-sn*c)+cs*s;
	      res=sn+cor;
	      cor=(sn-res)+cor;
	      cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
	      return (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x);
	    }
	    break;

	  case 1:
	  case 3:
	    if (a<0)
	      {a=-a;da=-da;}
	    u.x=big.x+a;
	    y=a-(u.x-big.x)+da;
	    xx=y*y;
	    k=u.i[LOW_HALF]<<2;
	    sn=sincos.x[k];
	    ssn=sincos.x[k+1];
	    cs=sincos.x[k+2];
	    ccs=sincos.x[k+3];
	    s = y + y*xx*(sn3 +xx*sn5);
	    c = xx*(cs2 +xx*(cs4 + xx*cs6));
	    cor=(ccs-s*ssn-cs*c)-sn*s;
	    res=cs+cor;
	    cor=(cs-res)+cor;
	    cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
	    return (res==res+cor)? ((n&2)?-res:res) : sloww2(a,da,x,n);

	    break;

	  }

	}    /*   else  if (k <  0x419921FB )    */

/*---------------------105414350 <|x|< 281474976710656 --------------------*/
	else if (k < 0x42F00000 ) {
	  t = (x*hpinv.x + toint.x);
	  xn = t - toint.x;
	  v.x = t;
	  xn1 = (xn+8.0e22)-8.0e22;
	  xn2 = xn - xn1;
	  y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x);
	  n =v.i[LOW_HALF]&3;
	  da = xn1*pp3.x;
	  t=y-da;
	  da = (y-t)-da;
	  da = (da - xn2*pp3.x) -xn*pp4.x;
	  a = t+da;
	  da = (t-a)+da;
	  eps = 1.0e-24;

	  switch (n) {
	  case 0:
	  case 2:
	    xx = a*a;
	    if (n) {a=-a;da=-da;}
	    if (xx < 0.01588) {
              /* Taylor series */
	      t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
	      res = a+t;
	      cor = (a-res)+t;
	      cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
	      return (res == res + cor)? res : bsloww(a,da,x,n);
	    }
	    else  {
	      if (a>0) {m=1;t=a;db=da;}
	      else {m=0;t=-a;db=-da;}
	      u.x=big.x+t;
	      y=t-(u.x-big.x);
	      xx=y*y;
	      s = y + (db+y*xx*(sn3 +xx*sn5));
	      c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
	      k=u.i[LOW_HALF]<<2;
	      sn=sincos.x[k];
	      ssn=sincos.x[k+1];
	      cs=sincos.x[k+2];
	      ccs=sincos.x[k+3];
	      cor=(ssn+s*ccs-sn*c)+cs*s;
	      res=sn+cor;
	      cor=(sn-res)+cor;
	      cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
	      return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
		   }
	    break;

	  case 1:
	  case 3:
	    if (a<0)
	      {a=-a;da=-da;}
	    u.x=big.x+a;
	    y=a-(u.x-big.x)+da;
	    xx=y*y;
	    k=u.i[LOW_HALF]<<2;
	    sn=sincos.x[k];
	    ssn=sincos.x[k+1];
	    cs=sincos.x[k+2];
	    ccs=sincos.x[k+3];
	    s = y + y*xx*(sn3 +xx*sn5);
	    c = xx*(cs2 +xx*(cs4 + xx*cs6));
	    cor=(ccs-s*ssn-cs*c)-sn*s;
	    res=cs+cor;
	    cor=(cs-res)+cor;
	    cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
	    return (res==res+cor)? ((n&2)?-res:res) : bsloww2(a,da,x,n);

	    break;

	  }

	}    /*   else  if (k <  0x42F00000 )   */

/* -----------------281474976710656 <|x| <2^1024----------------------------*/
	else if (k < 0x7ff00000) {

	  n = __branred(x,&a,&da);
	  switch (n) {
	  case 0:
	    if (a*a < 0.01588) return bsloww(a,da,x,n);
	    else return bsloww1(a,da,x,n);
	    break;
	  case 2:
	    if (a*a < 0.01588) return bsloww(-a,-da,x,n);
	    else return bsloww1(-a,-da,x,n);
	    break;

	  case 1:
	  case 3:
	    return  bsloww2(a,da,x,n);
	    break;
	  }

	}    /*   else  if (k <  0x7ff00000 )    */

/*--------------------- |x| > 2^1024 ----------------------------------*/
	else return x / x;
	return 0;         /* unreachable */
}


/*******************************************************************/
/* An ultimate cos routine. Given an IEEE double machine number x   */
/* it computes the correctly rounded (to nearest) value of cos(x)  */
/*******************************************************************/

double __cos(double x)
{
  double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
  mynumber u,v;
  int4 k,m,n;

  u.x = x;
  m = u.i[HIGH_HALF];
  k = 0x7fffffff&m;

  if (k < 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */

  else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
    y=ABS(x);
    u.x = big.x+y;
    y = y-(u.x-big.x);
    xx=y*y;
    s = y + y*xx*(sn3 +xx*sn5);
    c = xx*(cs2 +xx*(cs4 + xx*cs6));
    k=u.i[LOW_HALF]<<2;
    sn=sincos.x[k];
    ssn=sincos.x[k+1];
    cs=sincos.x[k+2];
    ccs=sincos.x[k+3];
    cor=(ccs-s*ssn-cs*c)-sn*s;
    res=cs+cor;
    cor=(cs-res)+cor;
    return (res==res+1.020*cor)? res : cslow2(x);

}    /*   else  if (k < 0x3feb6000)    */

  else if (k <  0x400368fd ) {/* 0.855469  <|x|<2.426265  */;
    y=hp0.x-ABS(x);
    a=y+hp1.x;
    da=(y-a)+hp1.x;
    xx=a*a;