1240116Smarcel/*- 2240116Smarcel * Copyright (c) 2017-2023 Steven G. Kargl 3240116Smarcel * All rights reserved. 4240116Smarcel * 5240116Smarcel * Redistribution and use in source and binary forms, with or without 6240116Smarcel * modification, are permitted provided that the following conditions 7240116Smarcel * are met: 8240116Smarcel * 1. Redistributions of source code must retain the above copyright 9240116Smarcel * notice unmodified, this list of conditions, and the following 10240116Smarcel * disclaimer. 11240116Smarcel * 2. Redistributions in binary form must reproduce the above copyright 12240116Smarcel * notice, this list of conditions and the following disclaimer in the 13240116Smarcel * documentation and/or other materials provided with the distribution. 14240116Smarcel * 15240116Smarcel * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16240116Smarcel * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17240116Smarcel * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18240116Smarcel * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19240116Smarcel * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20240116Smarcel * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21240116Smarcel * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22240116Smarcel * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23240116Smarcel * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24240116Smarcel * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25240116Smarcel */ 26240116Smarcel 27240116Smarcel/* 28240116Smarcel * See ../src/s_cospi.c for implementation details. 29240116Smarcel */ 30240116Smarcel 31240116Smarcel#include "math.h" 32240116Smarcel#include "math_private.h" 33240116Smarcel 34240116Smarcel/* 35240116Smarcel * pi_hi contains the leading 56 bits of a 169 bit approximation for pi. 36240116Smarcel */ 37240116Smarcelstatic const long double 38240116Smarcelpi_hi = 3.14159265358979322702026593105983920e+00L, 39240116Smarcelpi_lo = 1.14423774522196636802434264184180742e-17L; 40240116Smarcel 41240116Smarcel#include "k_cospil.h" 42240116Smarcel#include "k_sinpil.h" 43240116Smarcel 44260029Sjmmvstatic volatile const double vzero = 0; 45240116Smarcel 46240116Smarcellong double 47240116Smarcelcospil(long double x) 48240116Smarcel{ 49240116Smarcel long double ai, ar, ax, c; 50240116Smarcel 51240116Smarcel ax = fabsl(x); 52240116Smarcel 53240116Smarcel if (ax <= 1) { 54240116Smarcel if (ax < 0.25) { 55240116Smarcel if (ax < 0x1p-60) { 56240116Smarcel if ((int)x == 0) 57240116Smarcel return (1); 58240116Smarcel } 59240116Smarcel return (__kernel_cospil(ax)); 60240116Smarcel } 61240116Smarcel 62240116Smarcel if (ax < 0.5) 63240116Smarcel c = __kernel_sinpil(0.5 - ax); 64240116Smarcel else if (ax < 0.75) { 65240116Smarcel if (ax == 0.5) 66260029Sjmmv return (0); 67240116Smarcel c = -__kernel_sinpil(ax - 0.5); 68240116Smarcel } else 69240116Smarcel c = -__kernel_cospil(1 - ax); 70240116Smarcel return (c); 71240116Smarcel } 72240116Smarcel 73240116Smarcel if (ax < 0x1p112) { 74240116Smarcel /* Split ax = ai + ar with 0 <= ar < 1. */ 75240116Smarcel FFLOORL128(ax, ai, ar); 76240116Smarcel 77240116Smarcel if (ar < 0.5) { 78240116Smarcel if (ar < 0.25) 79240116Smarcel c = ar == 0 ? 1 : __kernel_cospil(ar); 80240116Smarcel else 81240116Smarcel c = __kernel_sinpil(0.5 - ar); 82240116Smarcel } else { 83240116Smarcel if (ar < 0.75) { 84240116Smarcel if (ar == 0.5) 85240116Smarcel return (0); 86240116Smarcel c = -__kernel_sinpil(ar - 0.5); 87240116Smarcel } else 88240116Smarcel c = -__kernel_cospil(1 - ar); 89240116Smarcel } 90240116Smarcel return (fmodl(ai, 2.L) == 0 ? c : -c); 91240116Smarcel } 92240116Smarcel 93240116Smarcel if (isinf(x) || isnan(x)) 94240116Smarcel return (vzero / vzero); 95240116Smarcel 96240116Smarcel /* 97240116Smarcel * For |x| >= 0x1p113, it is always an even integer, so return 1. 98240116Smarcel */ 99240116Smarcel if (ax >= 0x1p113) 100240116Smarcel return (1); 101240116Smarcel /* 102240116Smarcel * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even 103240116Smarcel * or odd integer to return 1 or -1. 104240116Smarcel */ 105240116Smarcel 106240116Smarcel return (fmodl(ax, 2.L) == 0 ? 1 : -1); 107240116Smarcel} 108240116Smarcel