Home | History | Annotate | Download | only in single 
      1 /*
      2  * s_sincosf.c - single precision sine and cosine functions
      3  *
      4  * Copyright (c) 2009-2018, Arm Limited.
      5  * SPDX-License-Identifier: MIT
      6  */
      7 
      8 /*
      9  * Source: my own head, and Remez-generated polynomial approximations.
     10  */
     11 
     12 #include <fenv.h>
     13 #include <math.h>
     14 #include <errno.h>
     15 #include "rredf.h"
     16 #include "math_private.h"
     17 
     18 #ifdef __cplusplus
     19 extern "C" {
     20 #endif /* __cplusplus */
     21 
     22 #ifndef COSINE
     23 #define FUNCNAME sinf
     24 #define SOFTFP_FUNCNAME __softfp_sinf
     25 #define DO_SIN (!(q & 1))
     26 #define NEGATE_SIN ((q & 2))
     27 #define NEGATE_COS ((q & 2))
     28 #define TRIVIAL_RESULT(x) FLOAT_CHECKDENORM(x)
     29 #define ERR_INF MATHERR_SINF_INF
     30 #else
     31 #define FUNCNAME cosf
     32 #define SOFTFP_FUNCNAME __softfp_cosf
     33 #define DO_SIN (q & 1)
     34 #define NEGATE_SIN (!(q & 2))
     35 #define NEGATE_COS ((q & 2))
     36 #define TRIVIAL_RESULT(x) 1.0f
     37 #define ERR_INF MATHERR_COSF_INF
     38 #endif
     39 
     40 float FUNCNAME(float x)
     41 {
     42     int q;
     43 
     44     /*
     45      * Range-reduce x to the range [-pi/4,pi/4].
     46      */
     47     {
     48         /*
     49          * I enclose the call to __mathlib_rredf in braces so that
     50          * the address-taken-ness of qq does not propagate
     51          * throughout the rest of the function, for what that might
     52          * be worth.
     53          */
     54         int qq;
     55         x = __mathlib_rredf(x, &qq);
     56         q = qq;
     57     }
     58     if (__builtin_expect(q < 0, 0)) { /* this signals tiny, inf, or NaN */
     59         unsigned k = fai(x) << 1;
     60         if (k < 0xFF000000)            /* tiny */
     61             return TRIVIAL_RESULT(x);
     62         else if (k == 0xFF000000)      /* inf */
     63             return ERR_INF(x);
     64         else                           /* NaN */
     65             return FLOAT_INFNAN(x);
     66     }
     67 
     68     /*
     69      * Depending on the quadrant we were in, we may have to compute
     70      * a sine-like function (f(0)=0) or a cosine-like one (f(0)=1),
     71      * and we may have to negate it.
     72      */
     73     if (DO_SIN) {
     74         float x2 = x*x;
     75         /*
     76          * Coefficients generated by the command
     77 
     78 ./auxiliary/remez.jl --variable=x2 --suffix=f -- '0' 'atan(BigFloat(1))^2' 2 0 'x==0 ? -1/BigFloat(6) : (x->(sin(x)-x)/x^3)(sqrt(x))' 'sqrt(x^3)'
     79 
     80          */
     81         x += x * (x2 * (
     82                       -1.666665066929417292436220415142244613956015227491999719404711781344783392564922e-01f+x2*(8.331978663157089651408875887703995477889496917296385733254577121461421466427772e-03f+x2*(-1.949563623766929906511886482584265500187314705496861877317774185883215997931494e-04f))
     83                       ));
     84         if (NEGATE_SIN)
     85             x = -x;
     86         return x;
     87     } else {
     88         float x2 = x*x;
     89         /*
     90          * Coefficients generated by the command
     91 
     92 ./auxiliary/remez.jl --variable=x2 --suffix=f -- '0' 'atan(BigFloat(1))^2' 2 0 'x==0 ? -1/BigFloat(6) : (x->(cos(x)-1)/x^2)(sqrt(x))' 'x'
     93 
     94          */
     95         x = 1.0f + x2*(
     96             -4.999989478137016757327030935768632852012781143541026304540837816323349768666875e-01f+x2*(4.165629457842617238353362092016541041535652603456375154392942188742496860024377e-02f+x2*(-1.35978231111049428748566568960423202948250988565693107969571193763372093404347e-03f))
     97             );
     98         if (NEGATE_COS)
     99             x = -x;
    100         return x;
    101     }
    102 }
    103 
    104 
    105 #ifdef __cplusplus
    106 } /* end of extern "C" */
    107 #endif /* __cplusplus */
    108 
    109 /* end of sincosf.c */
    110