lib/sinl.c

   1 /* sin (sine) function with 'long double' argument.
   2
   3    Copyright (C) 2003-2024 Free Software Foundation, Inc.
   4
   5    This file is free software: you can redistribute it and/or modify
   6    it under the terms of the GNU Lesser General Public License as
   7    published by the Free Software Foundation, either version 3 of the
   8    License, or (at your option) any later version.
   9
  10    This file is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13    GNU Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public License
  16    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  17
  18 /* s_sinl.c -- long double version of s_sin.c.
  19  * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
  20  */
  21
  22 /*
  23  * ====================================================
  24  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  25  *
  26  * Developed at SunPro, a Sun Microsystems, Inc. business.
  27  * Permission to use, copy, modify, and distribute this
  28  * software is freely granted, provided that this notice
  29  * is preserved.
  30  * ====================================================
  31  */
  32
  33 #include <config.h>
  34
  35 /* Specification.  */
  36 #include <math.h>
  37
  38 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  39
  40 long double
  41 sinl (long double x)
  42 {
  43   return sin (x);
  44 }
  45
  46 #else
  47
  48 /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_sinl.c.  */
  49
  50 /* sinl(x)
  51  * Return sine function of x.
  52  *
  53  * kernel function:
  54  *      __kernel_sinl           ... sine function on [-pi/4,pi/4]
  55  *      __kernel_cosl           ... cosine function on [-pi/4,pi/4]
  56  *      __ieee754_rem_pio2l     ... argument reduction routine
  57  *
  58  * Method.
  59  *      Let S,C and T denote the sin, cos and tan respectively on
  60  *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  61  *      in [-pi/4 , +pi/4], and let n = k mod 4.
  62  *      We have
  63  *
  64  *          n        sin(x)      cos(x)        tan(x)
  65  *     ----------------------------------------------------------
  66  *          0          S           C             T
  67  *          1          C          -S            -1/T
  68  *          2         -S          -C             T
  69  *          3         -C           S            -1/T
  70  *     ----------------------------------------------------------
  71  *
  72  * Special cases:
  73  *      Let trig be any of sin, cos, or tan.
  74  *      trig(+-INF)  is NaN, with signals;
  75  *      trig(NaN)    is that NaN;
  76  *
  77  * Accuracy:
  78  *      TRIG(x) returns trig(x) nearly rounded
  79  */
  80
  81 # include "trigl.h"
  82
  83 long double
  84 sinl (long double x)
  85 {
  86   long double y[2], z = 0.0L;
  87   int n;
  88
  89   /* sinl(NaN) is NaN */
  90   if (isnanl (x))
  91     return x;
  92
  93   /* |x| ~< pi/4 */
  94   if (x >= -0.7853981633974483096156608458198757210492
  95       && x <= 0.7853981633974483096156608458198757210492)
  96     return kernel_sinl (x, z, 0);
  97
  98     /* sinl(Inf) is NaN, sinl(0) is 0 */
  99   else if (x + x == x)
 100     return x - x;               /* NaN */
 101
 102   /* argument reduction needed */
 103   else
 104     {
 105       n = ieee754_rem_pio2l (x, y);
 106       switch (n & 3)
 107         {
 108         case 0:
 109           return kernel_sinl (y[0], y[1], 1);
 110         case 1:
 111           return kernel_cosl (y[0], y[1]);
 112         case 2:
 113           return -kernel_sinl (y[0], y[1], 1);
 114         default:
 115           return -kernel_cosl (y[0], y[1]);
 116         }
 117     }
 118 }
 119
 120 #endif
 121
 122 #if 0
 123 int
 124 main (void)
 125 {
 126   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492));
 127   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *29));
 128   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *2));
 129   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *30));
 130   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *4));
 131   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *32));
 132   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *2/3));
 133   printf ("%.16Lg\n", sinl (0.7853981633974483096156608458198757210492 *4/3));
 134 }
 135 #endif