Home | History | Annotate | Line # | Download | only in ld80 
      1 /*-
      2  * Copyright (c) 2017, 2023 Steven G. Kargl
      3  * All rights reserved.
      4  *
      5  * Redistribution and use in source and binary forms, with or without
      6  * modification, are permitted provided that the following conditions
      7  * are met:
      8  * 1. Redistributions of source code must retain the above copyright
      9  *    notice unmodified, this list of conditions, and the following
     10  *    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  *
     15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
     16  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
     17  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
     18  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
     19  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
     20  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     21  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
     22  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
     24  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     25  */
     26 
     27 /*
     28  * See ../src/s_tanpi.c for implementation details.
     29  */
     30 
     31 #include <stdint.h>
     32 
     33 #include "math.h"
     34 #include "math_private.h"
     35 
     36 static const double
     37 pi_hi =  3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
     38 pi_lo = -2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */
     39 
     40 static inline long double
     41 __kernel_tanpil(long double x)
     42 {
     43 	long double hi, lo, t;
     44 
     45 	if (x < 0.25) {
     46 		hi = (float)x;
     47 		lo = x - hi;
     48 		lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
     49 		hi *= pi_hi;
     50 		_2sumF(hi, lo);
     51 		t = __kernel_tanl(hi, lo, -1);
     52 	} else if (x > 0.25) {
     53 		x = 0.5 - x;
     54 		hi = (float)x;
     55 		lo = x - hi;
     56 		lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
     57 		hi *= pi_hi;
     58 		_2sumF(hi, lo);
     59 		t = - __kernel_tanl(hi, lo, 1);
     60 	} else
     61 		t = 1;
     62 
     63 	return (t);
     64 }
     65 
     66 static volatile const double vzero = 0;
     67 
     68 long double
     69 tanpil(long double x)
     70 {
     71 	long double ax, hi, lo, odd, t;
     72 	uint64_t lx;
     73 	uint32_t j0;
     74 	uint16_t hx, ix;
     75 
     76 	EXTRACT_LDBL80_WORDS(hx, lx, x);
     77 	ix = hx & 0x7fff;
     78 	INSERT_LDBL80_WORDS(ax, ix, lx);
     79 
     80 	ENTERI();
     81 
     82 	if (ix < 0x3fff) {			/* |x| < 1 */
     83 		if (ix < 0x3ffe) {		/* |x| < 0.5 */
     84 			if (ix < 0x3fdd) {	/* |x| < 0x1p-34 */
     85 				if (x == 0)
     86 					RETURNI(x);
     87 				INSERT_LDBL80_WORDS(hi, hx,
     88 				    lx & 0xffffffff00000000ull);
     89 				hi *= 0x1p63L;
     90 				lo = x * 0x1p63L - hi;
     91 				t = (pi_lo + pi_hi) * lo + pi_lo * hi +
     92 				    pi_hi * hi;
     93 				RETURNI(t * 0x1p-63L);
     94 			}
     95 			t = __kernel_tanpil(ax);
     96 		} else if (ax == 0.5)
     97 			t = 1 / vzero;
     98 		else
     99 			t = -__kernel_tanpil(1 - ax);
    100 		RETURNI((hx & 0x8000) ? -t : t);
    101 	}
    102 
    103 	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
    104 		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
    105 		odd = (uint64_t)x & 1 ? -1 : 1;
    106 		ax -= x;
    107 		EXTRACT_LDBL80_WORDS(ix, lx, ax);
    108 
    109 		if (ix < 0x3ffe)		/* |x| < 0.5 */
    110 			t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax);
    111 		else if (ax == 0.5L)
    112 			t = odd / vzero;
    113 		else
    114 			t = -__kernel_tanpil(1 - ax);
    115 		RETURNI((hx & 0x8000) ? -t : t);
    116 	}
    117 
    118 	/* x = +-inf or nan. */
    119 	if (ix >= 0x7fff)
    120 		RETURNI(vzero / vzero);
    121 
    122 	/*
    123 	 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
    124 	 * or odd integer to set t = +0 or -0.
    125 	 * For |x| >= 0x1p64, it is always an even integer, so t = 0.
    126 	 */
    127 	t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1));
    128 	RETURNI((hx & 0x8000) ? -t : t);
    129 }
    130