aboutsummaryrefslogtreecommitdiff
path: root/math/single/s_tanf.c
blob: 0e13d95af845126ad1f2edc106fc7df20c885cb2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
/*
 * s_tanf.c - single precision tangent function
 *
 * Copyright (c) 2009-2018, Arm Limited.
 * SPDX-License-Identifier: MIT
 */

/*
 * Source: my own head, and Remez-generated polynomial approximations.
 */

#include <math.h>
#include "math_private.h"
#include <errno.h>
#include <fenv.h>
#include "rredf.h"

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

float tanf(float x)
{
    int q;

    /*
     * Range-reduce x to the range [-pi/4,pi/4].
     */
    {
        /*
         * I enclose the call to __mathlib_rredf in braces so that
         * the address-taken-ness of qq does not propagate
         * throughout the rest of the function, for what that might
         * be worth.
         */
        int qq;
        x = __mathlib_rredf(x, &qq);
        q = qq;
    }
    if (__builtin_expect(q < 0, 0)) { /* this signals tiny, inf, or NaN */
        unsigned k = fai(x) << 1;
        if (k < 0xFF000000)            /* tiny */
            return FLOAT_CHECKDENORM(x);
        else if (k == 0xFF000000)      /* inf */
            return MATHERR_TANF_INF(x);
        else                           /* NaN */
            return FLOAT_INFNAN(x);
    }

    /*
     * We use a direct polynomial approximation for tan(x) on
     * [-pi/4,pi/4], and then take the negative reciprocal of the
     * result if we're in an interval surrounding an odd rather than
     * even multiple of pi/2.
     *
     * Coefficients generated by the command

./auxiliary/remez.jl --variable=x2 --suffix=f -- '0' '(pi/BigFloat(4))^2' 5 0 'x==0 ? 1/BigFloat(3) : (tan(sqrt(x))-sqrt(x))/sqrt(x^3)' 'sqrt(x^3)'

     */
    {
        float x2 = x*x;
        x += x * (x2 * (
                      3.333294809182307633621540045249152105330074691488121206914336806061620616979305e-01f+x2*(1.334274588580033216191949445078951865160600494428914956688702429547258497367525e-01f+x2*(5.315177279765676178198868818834880279286012428084733419724267810723468887753723e-02f+x2*(2.520300881849204519070372772571624013984546591252791443673871814078418474596388e-02f+x2*(2.051177187082974766686645514206648277055233230110624602600687812103764075834307e-03f+x2*(9.943421494628597182458186353995299429948224864648292162238582752158235742046109e-03f)))))
                      ));
        if (q & 1)
            x = -1.0f/x;

        return x;
    }
}

#ifdef __cplusplus
} /* end of extern "C" */
#endif /* __cplusplus */

/* end of s_tanf.c */