18e77cc64SDavid Schultz /*- 25e53a4f9SPedro F. Giffuni * SPDX-License-Identifier: BSD-2-Clause-FreeBSD 35e53a4f9SPedro F. Giffuni * 48e77cc64SDavid Schultz * Copyright (c) 2007 Steven G. Kargl 58e77cc64SDavid Schultz * All rights reserved. 68e77cc64SDavid Schultz * 78e77cc64SDavid Schultz * Redistribution and use in source and binary forms, with or without 88e77cc64SDavid Schultz * modification, are permitted provided that the following conditions 98e77cc64SDavid Schultz * are met: 108e77cc64SDavid Schultz * 1. Redistributions of source code must retain the above copyright 118e77cc64SDavid Schultz * notice unmodified, this list of conditions, and the following 128e77cc64SDavid Schultz * disclaimer. 138e77cc64SDavid Schultz * 2. Redistributions in binary form must reproduce the above copyright 148e77cc64SDavid Schultz * notice, this list of conditions and the following disclaimer in the 158e77cc64SDavid Schultz * documentation and/or other materials provided with the distribution. 168e77cc64SDavid Schultz * 178e77cc64SDavid Schultz * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 188e77cc64SDavid Schultz * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 198e77cc64SDavid Schultz * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 208e77cc64SDavid Schultz * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 218e77cc64SDavid Schultz * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 228e77cc64SDavid Schultz * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 238e77cc64SDavid Schultz * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 248e77cc64SDavid Schultz * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 258e77cc64SDavid Schultz * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 268e77cc64SDavid Schultz * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 278e77cc64SDavid Schultz */ 288e77cc64SDavid Schultz 298e77cc64SDavid Schultz #include <sys/cdefs.h> 308e77cc64SDavid Schultz __FBSDID("$FreeBSD$"); 318e77cc64SDavid Schultz 328e77cc64SDavid Schultz /* 338e77cc64SDavid Schultz * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows 348e77cc64SDavid Schultz * an accuracy of <= 0.7412 ULP. 358e77cc64SDavid Schultz */ 368e77cc64SDavid Schultz 378e77cc64SDavid Schultz #include <float.h> 38340076f0SSteve Kargl #ifdef __i386__ 39340076f0SSteve Kargl #include <ieeefp.h> 40340076f0SSteve Kargl #endif 418e77cc64SDavid Schultz 42046e2d5dSSteve Kargl #include "fpmath.h" 438e77cc64SDavid Schultz #include "math.h" 448e77cc64SDavid Schultz #include "math_private.h" 458e77cc64SDavid Schultz #if LDBL_MANT_DIG == 64 46c273267eSSteve Kargl #include "../ld80/e_rem_pio2l.h" 47046e2d5dSSteve Kargl static const union IEEEl2bits 48046e2d5dSSteve Kargl pio4u = LD80C(0xc90fdaa22168c235, -00001, 7.85398163397448309628e-01L); 49046e2d5dSSteve Kargl #define pio4 (pio4u.e) 508e77cc64SDavid Schultz #elif LDBL_MANT_DIG == 113 51c273267eSSteve Kargl #include "../ld128/e_rem_pio2l.h" 52046e2d5dSSteve Kargl long double pio4 = 7.85398163397448309615660845819875721e-1L; 538e77cc64SDavid Schultz #else 548e77cc64SDavid Schultz #error "Unsupported long double format" 558e77cc64SDavid Schultz #endif 568e77cc64SDavid Schultz 578e77cc64SDavid Schultz long double 588e77cc64SDavid Schultz cosl(long double x) 598e77cc64SDavid Schultz { 608e77cc64SDavid Schultz union IEEEl2bits z; 61c273267eSSteve Kargl int e0; 62c273267eSSteve Kargl long double y[2]; 638e77cc64SDavid Schultz long double hi, lo; 648e77cc64SDavid Schultz 658e77cc64SDavid Schultz z.e = x; 668e77cc64SDavid Schultz z.bits.sign = 0; 678e77cc64SDavid Schultz 688e77cc64SDavid Schultz /* If x = +-0 or x is a subnormal number, then cos(x) = 1 */ 698e77cc64SDavid Schultz if (z.bits.exp == 0) 708e77cc64SDavid Schultz return (1.0); 718e77cc64SDavid Schultz 728e77cc64SDavid Schultz /* If x = NaN or Inf, then cos(x) = NaN. */ 738e77cc64SDavid Schultz if (z.bits.exp == 32767) 748e77cc64SDavid Schultz return ((x - x) / (x - x)); 758e77cc64SDavid Schultz 76340076f0SSteve Kargl ENTERI(); 77340076f0SSteve Kargl 788e77cc64SDavid Schultz /* Optimize the case where x is already within range. */ 79046e2d5dSSteve Kargl if (z.e < pio4) 80340076f0SSteve Kargl RETURNI(__kernel_cosl(z.e, 0)); 818e77cc64SDavid Schultz 82c273267eSSteve Kargl e0 = __ieee754_rem_pio2l(x, y); 83c273267eSSteve Kargl hi = y[0]; 84c273267eSSteve Kargl lo = y[1]; 858e77cc64SDavid Schultz 868e77cc64SDavid Schultz switch (e0 & 3) { 878e77cc64SDavid Schultz case 0: 888e77cc64SDavid Schultz hi = __kernel_cosl(hi, lo); 898e77cc64SDavid Schultz break; 908e77cc64SDavid Schultz case 1: 918e77cc64SDavid Schultz hi = - __kernel_sinl(hi, lo, 1); 928e77cc64SDavid Schultz break; 938e77cc64SDavid Schultz case 2: 948e77cc64SDavid Schultz hi = - __kernel_cosl(hi, lo); 958e77cc64SDavid Schultz break; 968e77cc64SDavid Schultz case 3: 978e77cc64SDavid Schultz hi = __kernel_sinl(hi, lo, 1); 988e77cc64SDavid Schultz break; 998e77cc64SDavid Schultz } 1008e77cc64SDavid Schultz 101340076f0SSteve Kargl RETURNI(hi); 1028e77cc64SDavid Schultz } 103