Home | History | Annotate | Line # | Download | only in noieee_src 
      1 /*	$NetBSD: n_asincos.c,v 1.11 2024/06/09 13:35:38 riastradh Exp $	*/
      2 /*
      3  * Copyright (c) 1985, 1993
      4  *	The Regents of the University of California.  All rights reserved.
      5  *
      6  * Redistribution and use in source and binary forms, with or without
      7  * modification, are permitted provided that the following conditions
      8  * are met:
      9  * 1. Redistributions of source code must retain the above copyright
     10  *    notice, this list of conditions and the following disclaimer.
     11  * 2. Redistributions in binary form must reproduce the above copyright
     12  *    notice, this list of conditions and the following disclaimer in the
     13  *    documentation and/or other materials provided with the distribution.
     14  * 3. Neither the name of the University nor the names of its contributors
     15  *    may be used to endorse or promote products derived from this software
     16  *    without specific prior written permission.
     17  *
     18  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
     19  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
     22  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     23  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     24  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     25  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     26  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     27  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     28  * SUCH DAMAGE.
     29  */
     30 
     31 #include <sys/cdefs.h>
     32 __RCSID("$NetBSD: n_asincos.c,v 1.11 2024/06/09 13:35:38 riastradh Exp $");
     33 
     34 #ifndef lint
     35 #if 0
     36 static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";
     37 #endif
     38 #endif /* not lint */
     39 
     40 /* ASIN(X)
     41  * RETURNS ARC SINE OF X
     42  * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
     43  * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
     44  *
     45  * Required system supported functions:
     46  *	copysign(x,y)
     47  *	sqrt(x)
     48  *
     49  * Required kernel function:
     50  *	atan2(y,x)
     51  *
     52  * Method :
     53  *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
     54  *		  computed as follows
     55  *			1-x*x                     if x <  0.5,
     56  *			2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
     57  *
     58  * Special cases:
     59  *	if x is NaN, return x itself;
     60  *	if |x|>1, return NaN.
     61  *
     62  * Accuracy:
     63  * 1)  If atan2() uses machine PI, then
     64  *
     65  *	asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
     66  *	and PI is the exact pi rounded to machine precision (see atan2 for
     67  *      details):
     68  *
     69  *	in decimal:
     70  *		pi = 3.141592653589793 23846264338327 .....
     71  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
     72  *    56 bits   PI = 3.141592653589793 227020265 ..... ,
     73  *
     74  *	in hexadecimal:
     75  *		pi = 3.243F6A8885A308D313198A2E....
     76  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
     77  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
     78  *
     79  *	In a test run with more than 200,000 random arguments on a VAX, the
     80  *	maximum observed error in ulps (units in the last place) was
     81  *	2.06 ulps.      (comparing against (PI/pi)*(exact asin(x)));
     82  *
     83  * 2)  If atan2() uses true pi, then
     84  *
     85  *	asin(x) returns the exact asin(x) with error below about 2 ulps.
     86  *
     87  *	In a test run with more than 1,024,000 random arguments on a VAX, the
     88  *	maximum observed error in ulps (units in the last place) was
     89  *      1.99 ulps.
     90  */
     91 
     92 #include "namespace.h"
     93 #include "mathimpl.h"
     94 
     95 __weak_alias(acos, _acos)
     96 __weak_alias(acosf, _asinf)
     97 __weak_alias(asin, _asin)
     98 __weak_alias(asinf, _asinf)
     99 
    100 double
    101 asin(double x)
    102 {
    103 	double s,t,one=1.0;
    104 #if !defined(__vax__)&&!defined(tahoe)
    105 	if(x!=x) return(x);	/* x is NaN */
    106 #endif	/* !defined(__vax__)&&!defined(tahoe) */
    107 	s=copysign(x,one);
    108 	if(s <= 0.5)
    109 	    return(atan2(x,sqrt(one-x*x)));
    110 	else
    111 	    { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
    112 
    113 }
    114 
    115 float
    116 asinf(float x)
    117 {
    118 	return (float)asin(x);
    119 }
    120 
    121 /* ACOS(X)
    122  * RETURNS ARC COS OF X
    123  * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
    124  * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
    125  *
    126  * Required system supported functions:
    127  *	copysign(x,y)
    128  *	sqrt(x)
    129  *
    130  * Required kernel function:
    131  *	atan2(y,x)
    132  *
    133  * Method :
    134  *			      ________
    135  *                           / 1 - x
    136  *	acos(x) = 2*atan2(  / -------- , 1 ) .
    137  *                        \/   1 + x
    138  *
    139  * Special cases:
    140  *	if x is NaN, return x itself;
    141  *	if |x|>1, return NaN.
    142  *
    143  * Accuracy:
    144  * 1)  If atan2() uses machine PI, then
    145  *
    146  *	acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
    147  *	and PI is the exact pi rounded to machine precision (see atan2 for
    148  *      details):
    149  *
    150  *	in decimal:
    151  *		pi = 3.141592653589793 23846264338327 .....
    152  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
    153  *    56 bits   PI = 3.141592653589793 227020265 ..... ,
    154  *
    155  *	in hexadecimal:
    156  *		pi = 3.243F6A8885A308D313198A2E....
    157  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
    158  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
    159  *
    160  *	In a test run with more than 200,000 random arguments on a VAX, the
    161  *	maximum observed error in ulps (units in the last place) was
    162  *	2.07 ulps.      (comparing against (PI/pi)*(exact acos(x)));
    163  *
    164  * 2)  If atan2() uses true pi, then
    165  *
    166  *	acos(x) returns the exact acos(x) with error below about 2 ulps.
    167  *
    168  *	In a test run with more than 1,024,000 random arguments on a VAX, the
    169  *	maximum observed error in ulps (units in the last place) was
    170  *	2.15 ulps.
    171  */
    172 
    173 double
    174 acos(double x)
    175 {
    176 	double t,one=1.0;
    177 #if !defined(__vax__)&&!defined(tahoe)
    178 	if(x!=x) return(x);
    179 #endif	/* !defined(__vax__)&&!defined(tahoe) */
    180 	if( x != -1.0)
    181 	    t=atan2(sqrt((one-x)/(one+x)),one);
    182 	else
    183 	    t=atan2(one,0.0);	/* t = PI/2 */
    184 	return(t+t);
    185 }
    186 
    187 float
    188 acosf(float x)
    189 {
    190 	return (float)acos(x);
    191 }
    192