Math functions

Previous

Functions re-implemented in libspot

libspot does not depend on the standard library, so its math. In particular, the SPOT algorithm needs both the exponential and the natural logarithm.

Natural logarithm

Several strategies exist to compute log. Currently we prefer accuracy against speed so we have compared different methods and parameters so as to be as close as possible as the standard library output (on my laptop :grin: but you can also check it with online compiler).

Our simple implementation is actually faster than the standard library while being as accurate as it.

Continued fraction

We currently use the following continued fraction12:

\[ \log\left(1 + z\right) = \dfrac{2 z}{2 + z + \displaystyle\KK_{m = 1}^{\infty}\left(\frac{-m^2 z^2}{(2 m + 1) (2 + z)}\right)} \]

where

\[ \KK_{m = 1}^{\infty}\left(\frac{-m^2 z^2}{(2 m + 1) (2 + z)}\right) = -\cfrac{z^2}{3 (2 + z) - \cfrac{4 z^2}{5 (2 + z) - \cfrac{9 z^2}{7 (2 + z) - \cfrac{16 z^2}{9 (2 + z) - ...}}}} \]

In practice, we truncate this expansion to depth \(d\). For instance, the following pyton script leverages sympy and the recursive behavior of the continued fraction to output a truncated version of the expansion.

import sympy as sp

def K(x, d: int, m: int = 1):
    if m == d:
        return 0

    a = -m * m
    b = 2 * m + 1
    return a * x * x / (b * (2 + x) + K(x, d, m + 1))

def log_cf(x, d: int):
    z = x - 1
    return 2 * z / (2 + z + K(z, d))

x = sp.Symbol("x")
print(log_cf(x, d=11))

Truncation

The following code is used so as to find the "best" depth. We basically retrieve the max mantissa error and the relative error. Indeed we leverage the IEEE754 representation so as to turn the computation \(\log(x), x\in\RR\) into \(\log(x), x\in[1, 2]\).

benchmark/log_cf_accuracy.c 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double const _NAN = 0.0 / 0.0;

double const LOG2 = 0x1.62e42fefa39efp-1;

/*
 * See en.wikipedia.org/wiki/Double-precision_floating-point_format
 */
typedef union {
    double d;
    struct {
        unsigned long long significand : 52;
        unsigned long long exponent : 11;
        unsigned long long sign : 1;
    } bits;
} double_cast;

struct accuracy_errors {
    unsigned long long significand;
    double x_significand;
    double relative;
    double x_relative;
};

void reset_errors(struct accuracy_errors *error) {
    // clean the structure
    error->significand = 0;
    error->relative = 0.0;
    error->x_relative = _NAN;
    error->x_significand = _NAN;
}

double K(double x, unsigned int d, unsigned int m) {
    if (m == d) {
        return 0.;
    }
    double md = m;
    double a = -md * md;
    double b = md + md + 1.;
    return a * x * x / (b * (2. + x) + K(x, d, m + 1));
}

double log_cf(double x, unsigned int d) {
    double z = x - 1.;
    return 2 * z / (2 + z + K(z, d, 1));
}

double _log(double x, unsigned int d) {
    double_cast casted = {.d = x};
    if (casted.bits.exponent == 1022 || casted.bits.exponent == 1023) {
        // 1/2 < x < 1
        // when x is close to 1 (but <1)., it is better
        // to call the function directly instead of setting exponent to 1023
        return log_cf(x, d);
    }
    double exponent = casted.bits.exponent;
    casted.bits.exponent = 1023;
    // when exponent is 1023, 1 <= casted.d <= 2
    // return logd(casted.d, depth) + (exponent - 1023.0) * LOG2;
    return log_cf(casted.d, d) + (exponent - 1023.0) * LOG2;
}

void accuracy(unsigned int d, unsigned int N, struct accuracy_errors *error) {
    double_cast cast_log = {.d = 0.0};
    double_cast cast_log_cf = {.d = 0.0};
    double x;

    // clean the errors
    reset_errors(error);

    // errors
    unsigned long long significand_err = 0;
    double relative_err = 0.0;

    // linspace
    double a = -8.;
    double b = 8.;
    double eps = (b - a) / (double)(N - 1);
    double u = a;

    for (unsigned int i = 0; i < N; i++) {
        x = pow(10., u);
        cast_log.d = log(x);
        cast_log_cf.d = _log(x, d);

        if (cast_log.bits.significand >= cast_log_cf.bits.significand) {
            significand_err =
                cast_log.bits.significand - cast_log_cf.bits.significand;
        } else {
            significand_err =
                cast_log_cf.bits.significand - cast_log.bits.significand;
        }

        // significand error
        if (significand_err > error->significand) {
            error->significand = significand_err;
            error->x_significand = x;
        }

        // relative error
        if (cast_log.d != 0.0) {
            relative_err = (cast_log_cf.d - cast_log.d) / cast_log.d;
            if (relative_err > error->relative) {
                error->relative = relative_err;
                error->x_relative = x;
            }
        }

        u += eps;
    }
}

int main() {
    unsigned int const N = 2000000;
    struct accuracy_errors error;

    printf("depth | max mantissa error | max relative error\n");
    printf("------|--------------------|-------------------\n");
    for (unsigned d = 5; d < 15; d++) {
        accuracy(d, N, &error);
        printf("%-6d|%20llu|%19E\n", d, error.significand, error.relative);
    }
}









The results on my laptop are presented below. Depths 9, 10 or 11 may be good candidates depending on the accuracy need.






depth
max mantissa error
max relative error






5
204513805
3.275499E-08




6
6065190
9.713924E-10




7
179516
2.875076E-11




8
5307
8.500201E-13




9
157
2.514663E-14




10
6
8.008480E-16




11
2
4.022239E-16




12
2
4.022239E-16




13
2
4.022239E-16




14
2
4.022239E-16






Speed


Let us present the speed of our implementation. In the following benchmark, we use an "inline" representation of the continued fraction with depth = 11 (instead of the recursive definitions given above).


We compare computation time of libspot vs the standard library (on 10M runs, \(\scriptsize 10^{-8}<x<10^{8}\)).



benchmark/log_cf_speed.c



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define SEED 0

// we use -O2 optimization flag
#pragma GCC optimize("O2")

double const _NAN = 0.0 / 0.0;
double const DMAX = RAND_MAX;
double const LOG2 = 0x1.62e42fefa39efp-1;

double const CPS = CLOCKS_PER_SEC;

// U(0, 1)
double runif() { return (double)rand() / DMAX; }

/*
 * See en.wikipedia.org/wiki/Double-precision_floating-point_format
 */
typedef union {
    double d;
    struct {
        unsigned long long significand : 52;
        unsigned long long exponent : 11;
        unsigned long long sign : 1;
    } bits;
} double_cast;

double _log_cf_11(double z) {
    double x = z - 1;
    double xx = x + 2;
    double x2 = x * x;

    double xx2 = xx + xx;
    double xx3 = xx + xx2;
    double xx5 = xx3 + xx2;
    double xx7 = xx5 + xx2;
    double xx9 = xx7 + xx2;
    double xx11 = xx9 + xx2;
    double xx13 = xx11 + xx2;
    double xx15 = xx13 + xx2;
    double xx17 = xx15 + xx2;
    double xx19 = xx17 + xx2;
    double xx21 = xx19 + xx2;

    return 2 * x /
           (-x2 / (-4 * x2 /
                       (-9 * x2 /
                            (-16 * x2 /
                                 (-25 * x2 /
                                      (-36 * x2 /
                                           (-49 * x2 /
                                                (-64 * x2 /
                                                     (-81 * x2 /
                                                          (-100 * x2 / xx21 +
                                                           xx19) +
                                                      xx17) +
                                                 xx15) +
                                            xx13) +
                                       xx11) +
                                  xx9) +
                             xx7) +
                        xx5) +
                   xx3) +
            xx);
}

double _log(double x) {
    double_cast casted = {.d = x};
    if (casted.bits.exponent == 1022 || casted.bits.exponent == 1023) {
        // 1/2 < x < 1
        // when x is close to 1 (but <1)., it is better
        // to call the function directly instead of setting exponent to 1023
        return _log_cf_11(x);
    }
    double exponent = casted.bits.exponent;
    casted.bits.exponent = 1023;
    // when exponent is 1023, 1 <= casted.d <= 2
    // return logd(casted.d, depth) + (exponent - 1023.0) * LOG2;
    return _log_cf_11(casted.d) + (exponent - 1023.0) * LOG2;
}

void speed(unsigned int N) {
    unsigned long long reference = 0;
    unsigned long long custom = 0;
    double p = 0.0;
    clock_t start;

    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 20. * (2. * runif() - 1.));
        log(p);
    }
    reference = clock() - start;

    // reference
    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 20. * (2. * runif() - 1.));
        _log(p);
    }
    custom = clock() - start;

    //     printf(" STDLIB: %llu\nLIBSPOT: %llu (x%.1f)\n", reference, custom,
    //            (double)custom / (double)reference);
    printf("LIBSPOT: %.3fs\n STDLIB: %.3fs (x%.1f)\n", (double)custom / CPS,
           (double)reference / CPS, (double)reference / (double)custom);
}

int main() {
    unsigned int const N = 10000000;
    speed(N);
}














When we turn on optimization flags (like -O2) we see that the libspot implementation is faster.


Exponential


Once again we use a continued fraction to compute exponential. 


Continued fraction


We currently use the following continued fraction3:


\[
\exp\left(z\right) = 1 + \cfrac{2 z}{2 - z + 2 \displaystyle\KK\limits_{m = 1}^{\infty}\left(\cfrac{a_m z^2}{1}\right)}
\]


with 
$$
a_m = \dfrac{1}{4 (2 m - 1) (2 m + 1)}
$$
so if we expand the continued fraction
$$
\KK\limits_{m = 1}^{\infty}\left(\cfrac{a_m z^2}{1}\right) = \cfrac{z^2 / 12}{1 + \cfrac{z^2 / 60}{1 + \cfrac{z^2 / 140}{1 + \cfrac{z^2 / 252}{1 + \cfrac{z^2 / 396}{1 + ...}}}}}
$$


import sympy as sp

def K(z2, d: int, m: int = 1):
    if m == d:
        return 0

    mm = 2 * m
    a = 4 * (mm - 1) * (mm + 1)
    return z2 / (a * (1 + K(z2, d, m + 1)))

def exp_cf(z, d: int):
    z2 = z * z
    return 1 + 2 * z / (2 - z + 2 * K(z2, d, 1))

z = sp.Symbol("z")
print(exp_cf(z, d=6))




Truncation


The following benchmark leads us to choose d = 6.



benchmark/exp_cf_accuracy.c



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double const _NAN = 0.0 / 0.0;
double const LOG2 = 0x1.62e42fefa39efp-1;

/*
 * See en.wikipedia.org/wiki/Double-precision_floating-point_format
 */
typedef union {
    double d;
    struct {
        unsigned long long significand : 52;
        unsigned long long exponent : 11;
        unsigned long long sign : 1;
    } bits;
} double_cast;

struct accuracy_errors {
    unsigned long long significand;
    double x_significand;
    double relative;
    double x_relative;
};

void reset_errors(struct accuracy_errors *error) {
    // clean the structure
    error->significand = 0;
    error->relative = 0.0;
    error->x_relative = _NAN;
    error->x_significand = _NAN;
}

double K(double z2, unsigned int d, unsigned int m) {
    if (m == d) {
        return 0.;
    }
    double mm = 2 * m;
    double a = 1. / (4 * (mm - 1) * (mm + 1));
    return z2 * a / (1 + K(z2, d, m + 1));
}

double exp_cf(double z, unsigned int d) {
    double z2 = z * z;
    return 1 + 2 * z / (2 - z + 2 * K(z2, d, 1));
}

double _exp(double x, unsigned int d) {
    if (x < 0) {
        return 1 / _exp(-x, d);
    }
    if (x > LOG2) {
        unsigned long long k = x / LOG2;
        double r = x - LOG2 * (double)k;
        // significand, exponent, sign
        // double_cast z = {.bits = {0, k + 1023, 0}};

        double_cast w = {.d = _exp(r, d)};
        w.bits.exponent += k;
        return w.d;
        // return _exp(r, d) * z.d;
    }

    return exp_cf(x, d);
}

void accuracy(unsigned int d, unsigned int N, struct accuracy_errors *error) {
    double_cast cast_exp = {.d = 0.0};
    double_cast cast_exp_cf = {.d = 0.0};
    double x;

    // clean the errors
    reset_errors(error);

    // errors
    unsigned long long significand_err = 0;
    double relative_err = 0.0;

    // linspace
    double a = -10;
    double b = 2;
    double eps = (b - a) / (double)(N - 1);
    double u = a;

    for (unsigned int i = 0; i < N; i++) {
        x = pow(10.0, u);
        cast_exp.d = exp(x);
        cast_exp_cf.d = _exp(x, d);

        if (cast_exp.bits.significand >= cast_exp_cf.bits.significand) {
            significand_err =
                cast_exp.bits.significand - cast_exp_cf.bits.significand;
        } else {
            significand_err =
                cast_exp_cf.bits.significand - cast_exp.bits.significand;
        }

        // significand error
        if (significand_err > error->significand) {
            error->significand = significand_err;
            error->x_significand = x;
        }

        // relative error
        if (cast_exp.d != 0.0) {
            relative_err = (cast_exp_cf.d - cast_exp.d) / cast_exp.d;
            if (relative_err > error->relative) {
                error->relative = relative_err;
                error->x_relative = x;
            }
        }

        u += eps;
    }
}

int main() {
    unsigned int const N = 2000000;
    struct accuracy_errors error;

    printf("depth | max mantissa error | max relative error\n");
    printf("------|--------------------|-------------------\n");
    for (unsigned d = 3; d < 10; d++) {
        accuracy(d, N, &error);
        printf("%-6d|%20llu|%19E\n", d, error.significand, error.relative);
    }
}


















depth
max mantissa error
max relative error






3
6998321454
7.769708E-07




4
13281258
1.014192E-14




5
16126
1.790451E-12




6
81
1.017038E-14




7
90
1.017038E-14




8
90
1.017038E-14




9
90
1.017038E-14






Speed


Like in the previous log benchmark, we use an "inline" representation of the continued fraction with depth = 6. We compare computation time of libspot vs the standard library (on 10M runs, \(\scriptsize 10^{-8}<x<10^{2}\)).



benchmark/exp_cf_speed.c



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define SEED 0

// we use -O2 optimization flag
#pragma GCC optimize("O2")

double const _NAN = 0.0 / 0.0;
double const DMAX = RAND_MAX;
double const LOG2 = 0x1.62e42fefa39efp-1;

double const CPS = CLOCKS_PER_SEC;

// U(0, 1)
double runif() { return (double)rand() / DMAX; }

/*
 * See en.wikipedia.org/wiki/Double-precision_floating-point_format
 */
typedef union {
    double d;
    struct {
        unsigned long long significand : 52;
        unsigned long long exponent : 11;
        unsigned long long sign : 1;
    } bits;
} double_cast;

double _exp_cf_6(double z) {
    double z2 = z * z;

    return 2 * z /
               (2 * z2 /
                    (12 * z2 / (z2 * (z2 + 396) / (6 * (z2 + 154)) + 60) +
                     12) -
                z + 2) +
           1;
}

double _exp(double x) {
    if (x > LOG2) {
        unsigned long long k = x / LOG2;
        double r = x - LOG2 * (double)k;
        double_cast w = {.d = _exp_cf_6(r)};
        w.bits.exponent += k;
        return w.d;
    }

    return _exp_cf_6(x);
}

void speed(unsigned int N) {
    unsigned long long reference = 0;
    unsigned long long custom = 0;
    double p = 0.0;
    clock_t start;

    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 10. * runif() - 8.);
        exp(p);
    }
    reference = clock() - start;

    // reference
    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 10. * runif() - 8.);
        _exp(p);
    }
    custom = clock() - start;

    printf("LIBSPOT: %.3fs\n STDLIB: %.3fs (x%.1f)\n", (double)custom / CPS,
           (double)reference / CPS, (double)reference / (double)custom);
}

int main() {
    unsigned int const N = 10000000;
    speed(N);
}














Once again, libspot implementation is faster when optimization flags are set.


Power


Finally, we need to compute \(x^\alpha\) when \(x \ge 0\) and \(\alpha \in\RR\). Currently we basically use our implementation of exp and log as \(x^\alpha = \exp\left(\alpha \log x\right)\).



benchmark/pow_accuracy_speed.c



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define SEED 0

double const _NAN = 0.0 / 0.0;
double const DMAX = RAND_MAX;
double const LOG2 = 0x1.62e42fefa39efp-1;
double const CPS = CLOCKS_PER_SEC;

/*
 * See en.wikipedia.org/wiki/Double-precision_floating-point_format
 */
typedef union {
    double d;
    struct {
        unsigned long long significand : 52;
        unsigned long long exponent : 11;
        unsigned long long sign : 1;
    } bits;
} double_cast;

double runif() { return (double)rand() / DMAX; }

struct accuracy_errors {
    unsigned long long significand;
    double x_significand;
    double relative;
    double x_relative;
};

struct speeds {
    double stdlib;
    double libspot;
};

void reset_errors(struct accuracy_errors *error) {
    // clean the structure
    error->significand = 0;
    error->relative = 0.0;
    error->x_relative = _NAN;
    error->x_significand = _NAN;
}

void reset_speeds(struct speeds *s) {
    s->libspot = 0.0;
    s->stdlib = 0.0;
};

double _log_cf_11(double z) {
    double x = z - 1;
    double xx = x + 2;
    double x2 = x * x;

    double xx2 = xx + xx;
    double xx3 = xx + xx2;
    double xx5 = xx3 + xx2;
    double xx7 = xx5 + xx2;
    double xx9 = xx7 + xx2;
    double xx11 = xx9 + xx2;
    double xx13 = xx11 + xx2;
    double xx15 = xx13 + xx2;
    double xx17 = xx15 + xx2;
    double xx19 = xx17 + xx2;
    double xx21 = xx19 + xx2;

    return 2 * x /
           (-x2 / (-4 * x2 /
                       (-9 * x2 /
                            (-16 * x2 /
                                 (-25 * x2 /
                                      (-36 * x2 /
                                           (-49 * x2 /
                                                (-64 * x2 /
                                                     (-81 * x2 /
                                                          (-100 * x2 / xx21 +
                                                           xx19) +
                                                      xx17) +
                                                 xx15) +
                                            xx13) +
                                       xx11) +
                                  xx9) +
                             xx7) +
                        xx5) +
                   xx3) +
            xx);
}

double _exp_cf_6(double z) {
    double z2 = z * z;

    return 2 * z /
               (2 * z2 /
                    (12 * z2 /
                         (60 * z2 / (140 * z2 / (7 * z2 / 11 + 252) + 140) +
                          60) +
                     12) -
                z + 2) +
           1;
}

double _exp(double x) {
    if (x < 0) {
        return 1.0 / _exp(-x);
    }
    if (x > LOG2) {
        unsigned long long k = x / LOG2;
        double r = x - LOG2 * (double)k;
        double_cast w = {.d = _exp_cf_6(r)};
        w.bits.exponent += k;
        return w.d;
    }

    return _exp_cf_6(x);
}

double _log(double x) {
    double_cast casted = {.d = x};
    if (casted.bits.exponent == 1022 || casted.bits.exponent == 1023) {
        // 1/2 < x < 1
        // when x is close to 1 (but <1)., it is better
        // to call the function directly instead of setting exponent to 1023
        return _log_cf_11(x);
    }
    double exponent = casted.bits.exponent;
    casted.bits.exponent = 1023;
    // when exponent is 1023, 1 <= casted.d <= 2
    // return logd(casted.d, depth) + (exponent - 1023.0) * LOG2;
    return _log_cf_11(casted.d) + (exponent - 1023.0) * LOG2;
}

double _pow(double x, double alpha) { return _exp(alpha * _log(x)); }

double _powbis(double x, double alpha) {
    double_cast casted = {.d = x};
    double p = alpha * (casted.bits.exponent - 1023);
    long long k = p;
    double beta = p - (double)k;
    casted.bits.exponent = 1023;
    double_cast t = {.d = _pow(casted.d, alpha) * _pow(2, beta)};
    t.bits.exponent += k;
    return t.d;
}

void accuracy(double alpha, unsigned int N, struct accuracy_errors *error) {
    double_cast cast_pow = {.d = 0.0};
    double_cast cast_pow_cf = {.d = 0.0};
    double x;

    // clean the errors
    reset_errors(error);

    // errors
    unsigned long long significand_err = 0;
    double relative_err = 0.0;

    // linspace
    double a = -1.;
    double b = 1.;
    double eps = (b - a) / (double)(N - 1);
    double u = a;

    for (unsigned int i = 0; i < N; i++) {
        x = pow(10., u);
        cast_pow.d = pow(x, alpha);
        cast_pow_cf.d = _pow(x, alpha);

        if (cast_pow.bits.significand >= cast_pow_cf.bits.significand) {
            significand_err =
                cast_pow.bits.significand - cast_pow_cf.bits.significand;
        } else {
            significand_err =
                cast_pow_cf.bits.significand - cast_pow.bits.significand;
        }

        // significand error
        if (significand_err > error->significand) {
            error->significand = significand_err;
            error->x_significand = x;
        }

        // relative error
        if (cast_pow.d != 0.0) {
            relative_err = (cast_pow_cf.d - cast_pow.d) / cast_pow.d;
            if (relative_err > error->relative) {
                error->relative = relative_err;
                error->x_relative = x;
            }
        }

        u += eps;
    }
}

void speed(double alpha, unsigned int N, struct speeds *s) {
    unsigned long long reference = 0;
    unsigned long long custom = 0;
    double p = 0.0;
    clock_t start;

    reset_speeds(s);

    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 10. * runif() - 8.);
        pow(p, alpha);
    }
    reference = clock() - start;

    // reference
    srand(SEED);
    start = clock();
    for (unsigned int k = 0; k < N; k++) {
        p = pow(10., 10. * runif() - 8.);
        _pow(p, alpha);
    }
    custom = clock() - start;

    s->libspot = (double)custom / CPS;
    s->stdlib = (double)reference / CPS;
}

int main() {
    unsigned int const N = 1000000;
    struct accuracy_errors error;
    struct speeds s;

    srand(SEED);

    double alpha[] = {
        pow(10., -3. + runif()), pow(10., -2. + runif()),
        pow(10., -1. + runif()), pow(10., runif()),
        pow(10., 1. + runif()),
    };

    const size_t A = sizeof(alpha) / sizeof(double);

    printf("ACCURACY ===\n");
    printf("  alpha  | max mantissa error | max relative error\n");
    printf("---------|--------------------|-------------------\n");

    for (size_t i = 0; i < A; i++) {
        accuracy(alpha[i], N, &error);
        printf("%-9.5f|%20llu|%19E\n", alpha[i], error.significand,
               error.relative);
    }

    printf("\nSPEED ===\n");
    printf("  alpha  | libspot |  stdlib | stdlib/libspot \n");
    printf("---------|---------|---------|----------------\n");
    for (size_t i = 0; i < A; i++) {
        speed(alpha[i], N, &s);
        printf("%-9.5f|%9.4f|%9.4f|%16.2f\n", alpha[i], s.libspot, s.stdlib,
               s.stdlib / s.libspot);
    }
}














This current choice leads to worse performances than the standard library. This function is then likely to be improved in the future.








    

  1. 

    Formula (4.5) p.111
    

    

    Khovanskii, A. N. (1963). The application of continued fractions and their generalizations to problems in approximation theory (p. 144). Groningen: Noordhoff.
    

    

    ↩
    

    2. 

  2. 

    Formula (11.2.3) p.196
    

    

    Cuyt, A. A., Petersen, V., Verdonk, B., Waadeland, H., & Jones, W. B. (2008). Handbook of continued fractions for special functions. Springer Science & Business Media.
    

    

    ↩
    

    3. 

  3. 

    Formula (11.1.2) p.194
    

    

    Cuyt, A. A., Petersen, V., Verdonk, B., Waadeland, H., & Jones, W. B. (2008). Handbook of continued fractions for special functions. Springer Science & Business Media.
    

    

    ↩
    

    4. 




  



                                

    
    
        
            
            
            
        
        Initial fit