Repository: scipy/scipy
File: scipy/special/gen_harmonic.h (C/C++)

Code before:
            // IEEE: pow(1.0, nan) is 1.0.
            return 1.0;
        }
        return NAN;
    }
    if (std::isinf(a)) {
        if (a > 0) {
            // a = +inf
            return 1.0;
        }
        else {
            // a = -inf
            if (n == 1) {
                // IEEE: pow(1.0, inf) = 1
                return 1.0;
            }
            else {
                return INFINITY;
            }
        }
    }
    if (a == 0) {
        return static_cast<double>(n);
    }
    if (a <= 1) {
        return sum_powers(1, n, a);
    }
    else {
        // If here, we know a is finite and a > 1, and n >= 1.
        //
        // For a > 1, we can use the formula zeta(a, 1) - zeta(a, n + 1),
        // where zeta(a, k) is the Hurwitiz zeta function.
        // But if zeta(a, 1) and zeta(a, n + 1) are close, precision is lost
        // in the subtraction, so we use the explicit sum instead. We consider
        // the values "close" if zeta(a, n + 1)/zeta(a, 1) > 0.5.
        //
        double z1 = xsf::cephes::zeta(a, static_cast<double>(1));
        double znp1 =  xsf::cephes::zeta(a, static_cast<double>(n + 1));
        if (znp1 / z1 <= 0.5) {
            return z1 - znp1;
        }
        else {
            return sum_powers(1, n, a);
        }
    }
}


//
// Computes
//

Review comment (suggestion):
Is there a reason not to do this? Same can be done in the normalized version.
```suggestion
    if (n == 1) {
        // IEEE: pow(1.0, _) is 1.0.
        return 1.0;
    }
    if (std::isnan(a)) {
        return NAN;
    }
    if (std::isinf(a)) {
        if (a > 0) {
            // a = +inf
            return 1.0;
        }
        else {
            // a = -inf
            return INFINITY;
        }
    }
```

Code after:
//
static inline double
gen_harmonic(int64_t n, double a)
{
    if (n < 1) {
        sf_error("_gen_harmonic", SF_ERROR_DOMAIN,
                 "n >= 1 is required, but got n = %" PRId64, n);
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (n == 1) {
        // IEEE: pow(1.0, _) is 1.0.
        return 1.0;
    }
    if (std::isnan(a)) {
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (std::isinf(a)) {
        if (a > 0) {
            // a = +inf
            return 1.0;
        }
        else {
            // a = -inf
            return std::numeric_limits<double>::infinity();
        }
    }
    if (a == 0) {
        return static_cast<double>(n);
    }
    if (a <= 1) {
        return sum_powers(1, n, a);
    }
    else {
        // If here, we know a is finite and a > 1, and n > 1.
        //
        // For a > 1, we can use the formula zeta(a, 1) - zeta(a, n + 1),
        // where zeta(a, k) is the Hurwitiz zeta function.
        // But if zeta(a, 1) and zeta(a, n + 1) are close, precision is lost
        // in the subtraction, so we use the explicit sum instead. We consider
        // the values "close" if zeta(a, n + 1)/zeta(a, 1) > zeta_ratio_threshold.
        //
        double z1 = xsf::cephes::zeta(a, static_cast<double>(1));
        double znp1 =  xsf::cephes::zeta(a, static_cast<double>(n + 1));
        if (znp1 / z1 <= zeta_ratio_threshold) {
            return z1 - znp1;
        }
        else {
            return sum_powers(1, n, a);
        }
    }
}