6262 float64_t MAXfloat64 = np.inf
6363
6464 float64_t NaN = < float64_t> np.nan
65+ float64_t EpsF64 = np.finfo(np.float64).eps
66+
67+ # Consider an operation ill-conditioned if
68+ # it will only have up to 3 significant digits in base 10 remaining.
69+ # https://en.wikipedia.org/wiki/Condition_number
70+ float64_t InvConditionNumber = EpsF64 * 1e3
6571
6672cdef bint is_monotonic_increasing_start_end_bounds(
6773 ndarray[int64_t, ndim= 1 ] start, ndarray[int64_t, ndim= 1 ] end
@@ -532,7 +538,7 @@ cdef void add_skew(float64_t val, int64_t *nobs,
532538 cdef:
533539 float64_t n, delta, delta_n, term1, m3_update, new_m3
534540
535- # This formulas are adapted from
541+ # Formulas adapted from
536542 # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
537543
538544 # Not NaN
@@ -545,7 +551,7 @@ cdef void add_skew(float64_t val, int64_t *nobs,
545551
546552 m3_update = delta_n * fma(term1, n - 2.0 , - 3.0 * m2[0 ])
547553 new_m3 = m3[0 ] + m3_update
548- if fabs(m3_update) + fabs(m3[0 ]) > 1e10 * fabs(new_m3):
554+ if ( fabs(m3_update) + fabs(m3[0 ])) * InvConditionNumber > fabs(new_m3):
549555 # possible catastrophic cancellation
550556 numerically_unstable[0 ] = True
551557
@@ -591,7 +597,7 @@ cdef void remove_skew(float64_t val, int64_t *nobs,
591597 m3_update = delta_n * fma(term1, n + 2.0 , - 3.0 * m2[0 ])
592598 new_m3 = m3[0 ] - m3_update
593599
594- if fabs(m3_update) + fabs(m3[0 ]) > 1e10 * fabs(new_m3):
600+ if ( fabs(m3_update) + fabs(m3[0 ])) * InvConditionNumber > fabs(new_m3):
595601 # possible catastrophic cancellation
596602 numerically_unstable[0 ] = True
597603
0 commit comments