in C89
Arc tangent functions
An Orthodox Implementation Using Padé approximations.
Slow but Precise, Compact, Portable, and Easy to Understand.

May 20, 2025
Takayuki HOSODA

IEEE 754 Double precision (64-bit)

Listing 1
_atan - arc tangent function in C89 (ANSI-C)
#include <math.h>
#include <float.h>

#ifndef M_SQRT2
#define M_SQRT2  1.41421356237309504880168872421
#endif
#ifndef M_PI_2
#define M_PI_2   1.57079632679489661923132169164
#endif
#ifndef M_PI_4
#define M_PI_4   0.78539816339744830961566084582
#endif
#ifndef DBL_EPSILON
#define DBL_EPSILON    2.2204460492503131E-16
#endif

double _atan(double x);
/* _atan(x)  computes the principal value of the arc tangent of x. */
/* _atan Rev.1.2 (June 3, 2025) (c) 2025 Takayuki HOSODA */
double _atan(double x) {
    double q, s, v, w;
    double offset = 0;
    if (x >  (2 / DBL_EPSILON)) return  M_PI_2;
    if (x < -(2 / DBL_EPSILON)) return -M_PI_2;
    if ((v = fabs(x)) >= 1) {
        if (v > M_SQRT2 + 1) {
            offset = -M_PI_2; v = 1 / v;
        } else {
            offset =  M_PI_4; v = (v - 1) / (1 + v);
        }
    } else if (v >= M_SQRT2 - 1) {
            offset = -M_PI_4; v = (1 - v) / (1 + v);
    }
    s = v * v;
    q =           2401244.6571; w =           2073564.8997;
    q = q * s +  52026974.7295; w = w * s +  32801340.0004;
    q = q * s + 312161850.;     w = w * s + 136306170.;
    q = q * s + 758107350.;     w = w * s + 205633428.;
    q = q * s + 800224425.;     w = w * s + 101846745.;
    q = q * s + 305540235.;     w = w * s;
    if (q != w) v -= w / q * v; /* |err| < 3.2e-17, x in [0, √2-1] */
    s = v + offset;
    return (offset * x >= 0) ? (offset != 0 || x >= 0) ? s : -v : -s;
}


Fig. 1: Relative error of _atan(x), shown as _atan(x) / atan(x) - 1, where atan(x) is a rounded high-precision reference (≈40 digits), not a libm result.


Fig. 2: Theoretical approximation error of _atan(x), shown as _atan(x) / atan(x) - 1, where atan(x) is a rounded high-precision reference (≈40 digits), not a libm result.

X86 Extended-precision (80-bit)

Listing 2
_atanl - arc tangent function in C89 (ANSI-C)
#include <math.h>
#include <float.h>

#define Q_SQRT2  1.41421356237309504880168872421L
#define Q_PI_2   1.57079632679489661923132169164L
#define Q_PI_4   0.78539816339744830961566084582L
#ifndef LDBL_EPSILON
#define LDBL_EPSILON   1.0842021724855044340E-19L
#endif

long double _atanl(long double x);
/* _atanl(x) computes the principal value of the arc tangent of x. */
/* _atanl Rev.1.2 (June 3, 2025) (c) 2025 Takayuki HOSODA */
long double _atanl(long double x) {
    long double q, s, v, w;
    long double offset = 0;

    if (x >  (2 / LDBL_EPSILON)) return  Q_PI_2;
    if (x < -(2 / LDBL_EPSILON)) return -Q_PI_2;
    if ((v = fabs(x)) >= 1) {
        if (v > Q_SQRT2 + 1.0L) {
            offset = -Q_PI_2; v = 1 / v;
        } else {
            offset =  Q_PI_4; v = (v - 1) / (1 + v);
        }
    } else if (v > Q_SQRT2 - 1.0L) {
            offset = -Q_PI_4; v = (1 - v) / (1 + v);
    }
    s = v * v;
    w =               29360128.L; q =            27923620725.L;
    w = w * s +    23930643317.L; q = q * s +   742768311285.L;
    w = w * s +   501165956970.L; q = q * s +  6684914801565.L;
    w = w * s +  3525222266475.L; q = q * s + 28472785265925.L;
    w = w * s + 11362119869100.L; q = q * s + 64710875604375.L;
    w = w * s + 18420316154955.L; q = q * s + 80639706522375.L;
    w = w * s + 14615037006570.L; q = q * s + 51967810869975.L;
    w = w * s +  4512611027925.L; q = q * s + 13537833083775.L;
    w = w * s;
    if (q != w) v -= w / q * v; /* |err| < 1.43e-21, x ∈ [0, √2-1] */
    s = v + offset;
    return (offset * x >= 0) ? (offset != 0 || x >= 0) ? s : -v: -s;
}


Fig. 3: Relative error of _atanl(x), shown as _atanl(x) / atanl(x) - 1, where atanl(x) is a rounded high-precision reference (≈40 digits), not a libm result.

Appendix — Relations and Infinite series of the arctangent functions

\displaystyle\arctan \left({\frac {1}{x}}\right) &=& \left\{ \begin{array}{ll} {\displaystyle\frac {\pi }{2}}-\arctan(x) & {\text{if }} x > 0 \\ -{\displaystyle\frac {\pi }{2}}-\arctan(x) & {\text{if }} x < 0 \end{array}\\ \displaystyle \arctan(u)\pm \arctan(v)&=&\arctan \left({\frac {u\pm v}{1\mp uv}}\right){\pmod {\pi }}\,,\quad uv\neq 1\\ \\ \displaystyle \arctan(z)&=&{\frac {z}{1+z^{2}}}\sum _{n=0}^{\infty }\prod _{k=1}^{n}{\frac {2kz^{2}}{(2k+1)(1+z^{2})}}

Reference

  1. Jean-Michel Muller, "Elementary Functions: Algorithms and Implementation", Birkhauser, 3rd Edition, 2016.
  2. Okumura, Haruhiko; Shudo, Kazuyuki; Sugiura, Masaki; Tsuchimura, Nobuyuki; Tsuru, Kazuo; Hosoda, Takayuki; Matsui, Yoshimitsu; & Mitsunari, Shigeo. (2003). Java ni yoru Arugorizumu Jiten [Encyclopedia of Algorithms with Java]. Gijutsu Hyoron Sha. ISBN4-7741-1729-3 [in Japanese]

External links

www.finetune.co.jp [Mail] © 2000 Takayuki HOSODA.