Skip to content

Commit

Permalink
Move iv_ratio.h into special/ directory.
Browse files Browse the repository at this point in the history
  • Loading branch information
fancidev committed May 7, 2024
1 parent 2debbfe commit 3639f63
Show file tree
Hide file tree
Showing 6 changed files with 29 additions and 25 deletions.
9 changes: 0 additions & 9 deletions scipy/special/_add_newdocs.py
Original file line number Diff line number Diff line change
Expand Up @@ -12845,12 +12845,3 @@ def add_newdoc(name, doc):
scalar or ndarray
""")

add_newdoc("_iv_ratio",
r"""
_iv_ratio(v, x)
Internal function, do not use.
Return `iv(v, x) / iv(v-1, x)` for `v >= 1` and `x >= 0`.
""")
2 changes: 1 addition & 1 deletion scipy/special/_generate_pyx.py
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@
'_spherical_kn', '_spherical_kn_d', 'airy', 'airye', 'bei', 'beip', 'ber', 'berp',
'binom', 'exp1', 'expi', 'expit', 'exprel', 'gamma', 'gammaln', 'hankel1',
'hankel1e', 'hankel2', 'hankel2e', 'hyp2f1', 'it2i0k0', 'it2j0y0', 'it2struve0',
'itairy', 'iti0k0', 'itj0y0', 'itmodstruve0', 'itstruve0', 'iv', 'ive', 'jv',
'itairy', 'iti0k0', 'itj0y0', 'itmodstruve0', 'itstruve0', 'iv', '_iv_ratio', 'ive', 'jv',
'jve', 'kei', 'keip', 'kelvin', 'ker', 'kerp', 'kv', 'kve', 'log_expit',
'log_wright_bessel', 'loggamma', 'logit', 'mathieu_a', 'mathieu_b', 'mathieu_cem',
'mathieu_modcem1', 'mathieu_modcem2', 'mathieu_modsem1', 'mathieu_modsem2',
Expand Down
8 changes: 8 additions & 0 deletions scipy/special/_special_ufuncs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "special/fresnel.h"
#include "special/gamma.h"
#include "special/hyp2f1.h"
#include "special/iv_ratio.h"
#include "special/kelvin.h"
#include "special/lambertw.h"
#include "special/legendre.h"
Expand Down Expand Up @@ -139,6 +140,7 @@ extern const char *hankel2_doc;
extern const char *hankel2e_doc;
extern const char *hyp2f1_doc;
extern const char *iv_doc;
extern const char *iv_ratio_doc;
extern const char *ive_doc;
extern const char *jv_doc;
extern const char *jve_doc;
Expand Down Expand Up @@ -413,6 +415,12 @@ PyMODINIT_FUNC PyInit__special_ufuncs() {
);
PyModule_AddObjectRef(_special_ufuncs, "iv", iv);

PyObject *iv_ratio = SpecFun_NewUFunc(
{static_cast<func_dd_d_t>(special::iv_ratio)},
"_iv_ratio", iv_ratio_doc
);
PyModule_AddObjectRef(_special_ufuncs, "_iv_ratio", iv_ratio);

PyObject *ive = SpecFun_NewUFunc(
{static_cast<func_ff_f_t>(special::cyl_bessel_ie), static_cast<func_dd_d_t>(special::cyl_bessel_ie),
static_cast<func_fF_F_t>(special::cyl_bessel_ie), static_cast<func_dD_D_t>(special::cyl_bessel_ie)},
Expand Down
8 changes: 8 additions & 0 deletions scipy/special/_special_ufuncs_docs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1760,6 +1760,14 @@ const char *iv_doc = R"(
)";

const char *iv_ratio_doc = R"(
_iv_ratio(v, x, out=None)
Internal function, do not use.
Return `iv(v, x) / iv(v-1, x)` for `v >= 1` and `x >= 0`.
)";

const char *ive_doc = R"(
ive(v, z, out=None)
Expand Down
5 changes: 0 additions & 5 deletions scipy/special/functions.json
Original file line number Diff line number Diff line change
Expand Up @@ -1302,10 +1302,5 @@
"hypergeom_skewness_float": "fff->f",
"hypergeom_skewness_double": "ddd->d"
}
},
"_iv_ratio": {
"_iv_ratio.h++": {
"iv_ratio": "dd->d"
}
}
}
22 changes: 12 additions & 10 deletions scipy/special/_iv_ratio.h → scipy/special/special/iv_ratio.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@

#pragma once

#include "special/tools.h"
#include "special/error.h"
#include <cstdint> // for uint64_t
#include <cmath> // for frexp, ldexp, fmax, fma
#include <utility> // for pair
#include "config.h"
#include "tools.h"
#include "error.h"

namespace special {

/* Generates the "tail" of Perron's continued fraction for `iv(v,x)/iv(v-1,x)`
* for v >= 1 and x >= 0.
Expand Down Expand Up @@ -58,12 +58,12 @@ inline double iv_ratio(double v, double x) {
return std::numeric_limits<double>::quiet_NaN();
}
if (v < 1 || x < 0) {
special::set_error("iv_ratio", SF_ERROR_DOMAIN, NULL);
set_error("iv_ratio", SF_ERROR_DOMAIN, NULL);
return std::numeric_limits<double>::signaling_NaN();
}
if (std::isinf(v) && std::isinf(x)) {
// There is not a unique limit as both v and x tends to infinity.
special::set_error("iv_ratio", SF_ERROR_DOMAIN, NULL);
set_error("iv_ratio", SF_ERROR_DOMAIN, NULL);
return std::numeric_limits<double>::signaling_NaN();
}
if (std::isinf(v)) {
Expand All @@ -81,16 +81,18 @@ inline double iv_ratio(double v, double x) {
double xc = x * c;

IvRatioCFTailGenerator cf(vc, xc, c);
auto [fc, terms] = special::detail::series_eval_kahan(
special::detail::continued_fraction_series(cf),
auto [fc, terms] = detail::series_eval_kahan(
detail::continued_fraction_series(cf),
std::numeric_limits<double>::epsilon() * 0.5,
1000,
2*vc);

if (terms == 0) { // failed to converge; should not happen
special::set_error("iv_ratio", SF_ERROR_NO_RESULT, NULL);
set_error("iv_ratio", SF_ERROR_NO_RESULT, NULL);
return std::numeric_limits<double>::signaling_NaN();
}

return xc / (xc + fc);
}

} // namespace special

0 comments on commit 3639f63

Please sign in to comment.