libgfortran/intrinsics/trigd.c

   1 /* Implementation of the degree trignometric functions COSD, SIND, TAND.
   2    Copyright (C) 2020 Free Software Foundation, Inc.
   3    Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
   4
   5 This file is part of the GNU Fortran runtime library (libgfortran).
   6
   7 Libgfortran is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU General Public
   9 License as published by the Free Software Foundation; either
  10 version 3 of the License, or (at your option) any later version.
  11
  12 Libgfortran is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 Under Section 7 of GPL version 3, you are granted additional
  18 permissions described in the GCC Runtime Library Exception, version
  19 3.1, as published by the Free Software Foundation.
  20
  21 You should have received a copy of the GNU General Public License and
  22 a copy of the GCC Runtime Library Exception along with this program;
  23 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  24 <http://www.gnu.org/licenses/>.  */
  25
  26 #include "libgfortran.h"
  27
  28 #include <math.h>
  29
  30
  31 /*
  32    For real x, let {x}_P or x_P be the closest representible number in the
  33    floating point representation which uses P binary bits of fractional
  34    precision (with IEEE rounding semantics).
  35
  36    Similarly, let f_P(x) be shorthand for {f(x)}_P.
  37
  38    Let ulp_P(x) be the unit of least precision for x: in other words the
  39    maximal value of |a_P - b_P| where a_P <= x <= b_P and a_P != b_P.
  40
  41    Let x  ~= y <-> | x - y | <  ulp_P(x - y).
  42
  43    Let deg(x) be the value of x radians in degrees.
  44
  45    Values for each precision P were selected as follows.
  46
  47
  48    COSD_SMALL = 2**{-N} such that for all x <= COSD_SMALL:
  49
  50      * cos(deg(x)) ~= 1, or equivalently:
  51
  52        |      1 - cos(deg(x))  | < ulp_P(1).
  53
  54    Unfortunately for SIND (and therefore TAND) a similar relation is only
  55    possible for REAL(4) and REAL(8). With REAL(10) and REAL(16), enough
  56    precision is available such that sin_P(x) != x_P for some x less than any
  57    value. (There are values where this equality holds, but the distance has
  58    inflection points.)
  59
  60    For REAL(4) and REAL(8), we can select SIND_SMALL such that:
  61
  62      * sin(deg(x)) ~= deg(x), or equivalently:
  63
  64        | deg(x) - sin(deg(x)) | < ulp_P(deg(x)).
  65
  66  */
  67
  68 /* Build _gfortran_sind_r4, _gfortran_cosd_r4, and _gfortran_tand_r4  */
  69
  70 #define FTYPE       GFC_REAL_4
  71 #define SIND        sind_r4
  72 #define COSD        cosd_r4
  73 #define TAND        tand_r4
  74 #define SUFFIX(x)   x ## f
  75
  76 #define TINY        0x1.p-100f  /* ~= 7.889e-31 */
  77 #define COSD_SMALL  0x1.p-7f    /*  = 7.8125e-3 */
  78 #define SIND_SMALL  0x1.p-5f    /*  = 3.125e-2 */
  79 #define COSD30      8.66025388e-01f
  80
  81 #define PIO180H     1.74560547e-02f     /* high 12 bits.  */
  82 #define PIO180L    -2.76216747e-06f     /* Next 24 bits.  */
  83
  84 #include "trigd_lib.inc"
  85
  86 #undef FTYPE
  87 #undef TINY
  88 #undef COSD_SMALL
  89 #undef SIND_SMALL
  90 #undef COSD30
  91 #undef PIO180H
  92 #undef PIO180L
  93 #undef SIND
  94 #undef COSD
  95 #undef TAND
  96 #undef SUFFIX
  97
  98
  99 /* Build _gfortran_sind_r8, _gfortran_cosd_r8, and _gfortran_tand_r8.  */
 100
 101 #define FTYPE       GFC_REAL_8
 102 #define SIND        sind_r8
 103 #define COSD        cosd_r8
 104 #define TAND        tand_r8
 105 #define SUFFIX(x)   x
 106
 107 #define TINY        0x1.p-1000  /* ~= 9.33e-302 (min exp -1074) */
 108 #define COSD_SMALL  0x1.p-21    /* ~= 4.768e-7 */
 109 #define SIND_SMALL  0x1.p-19    /* ~= 9.537e-7 */
 110 #define COSD30      8.6602540378443860e-01
 111
 112 #define PIO180H     1.7453283071517944e-02      /* high 21 bits.  */
 113 #define PIO180L     9.4484253514332993e-09      /* Next 53 bits.  */
 114
 115 #include "trigd_lib.inc"
 116
 117 #undef FTYPE
 118 #undef TINY
 119 #undef COSD_SMALL
 120 #undef SIND_SMALL
 121 #undef COSD30
 122 #undef PIO180H
 123 #undef PIO180L
 124 #undef SIND
 125 #undef COSD
 126 #undef TAND
 127 #undef SUFFIX
 128
 129
 130 /* Build _gfortran_sind_r10, _gfortran_cosd_r10, and _gfortran_tand_r10.  */
 131
 132 #ifdef HAVE_GFC_REAL_10
 133
 134 #define FTYPE       GFC_REAL_10
 135 #define SIND        sind_r10
 136 #define COSD        cosd_r10
 137 #define TAND        tand_r10
 138 #define SUFFIX(x)   x ## l      /* L */
 139
 140 #define TINY        0x1.p-16400L        /* ~= 1.28e-4937 (min exp -16494) */
 141 #define COSD_SMALL  0x1.p-26L   /* ~= 1.490e-8 */
 142 #undef  SIND_SMALL              /* not precise */
 143 #define COSD30       8.66025403784438646787e-01L
 144
 145 #define PIO180H     1.74532925229868851602e-02L /* high 32 bits */
 146 #define PIO180L    -3.04358939097084072823e-12L /* Next 64 bits */
 147
 148 #include "trigd_lib.inc"
 149 #undef FTYPE
 150 #undef TINY
 151 #undef COSD_SMALL
 152 #undef SIND_SMALL
 153 #undef COSD30
 154 #undef PIO180H
 155 #undef PIO180L
 156 #undef SIND
 157 #undef COSD
 158 #undef TAND
 159 #undef SUFFIX
 160 #endif /* HAVE_GFC_REAL_10 */
 161
 162
 163 /* Build _gfortran_sind_r16, _gfortran_cosd_r16, and _gfortran_tand_r16.  */
 164
 165 #ifdef HAVE_GFC_REAL_16
 166
 167 #define FTYPE       GFC_REAL_16
 168 #define SIND        sind_r16
 169 #define COSD        cosd_r16
 170 #define TAND        tand_r16
 171
 172 #ifdef GFC_REAL_16_IS_FLOAT128  /* libquadmath.  */
 173 #define SUFFIX(x) x ## q
 174 #else
 175 #define SUFFIX(x) x ## l
 176 #endif /* GFC_REAL_16_IS_FLOAT128  */
 177
 178 #define TINY        SUFFIX(0x1.p-16400) /* ~= 1.28e-4937 */
 179 #define COSD_SMALL  SUFFIX(0x1.p-51)    /* ~= 4.441e-16 */
 180 #undef  SIND_SMALL              /* not precise */
 181 #define COSD30      SUFFIX(8.66025403784438646763723170752936183e-01)
 182 #define PIO180H     SUFFIX(1.74532925199433197605003442731685936e-02)
 183 #define PIO180L     SUFFIX(-2.39912634365882824665106671063098954e-17)
 184
 185 #include "trigd_lib.inc"
 186
 187 #undef FTYPE
 188 #undef COSD_SMALL
 189 #undef SIND_SMALL
 190 #undef COSD30
 191 #undef PIO180H
 192 #undef PIO180L
 193 #undef PIO180
 194 #undef D2R
 195 #undef CPYSGN
 196 #undef FABS
 197 #undef FMOD
 198 #undef SIN
 199 #undef COS
 200 #undef TAN
 201 #undef SIND
 202 #undef COSD
 203 #undef TAND
 204 #undef SUFFIX
 205 #endif /* HAVE_GFC_REAL_16 */