Revision history for Stats-LikeR
0.31 2026-08
[read_table]
- `na.strings` / `na_values`:
- `read_table` can now be told which field texts mean "missing". It is one option
under three names — R's `na.strings`, pandas' `na_values`, and `write_table`'s
own `undef.val` — so it can be spelled whichever way the surrounding code
already does. Each takes a string or an array reference, each means the same
thing, and passing more than one is an error the way `sep` and `delim` together
are:
read_table('cohort.csv', 'na.strings' => 'NA');
read_table('cohort.csv', na_values => ['NA', 'N/A', 'NULL', '-']);
read_table('cohort.csv', 'undef.val' => 'NA');
- An empty field has always read as `undef`. Nothing else did, and there was no
way to ask for it, so a file could not be read back the way it had been
written: `write_table` renders an `undef` cell with `undef.val`, which 0.302
listed as a known limitation whose closing sentence proposed exactly this
option. The round trip now closes, and `undef.val` is accepted as the third
name so that both halves of it can use one spelling:
write_table($rows, 'out.csv', 'undef.val' => 'NA');
my $back = read_table('out.csv', 'undef.val' => 'NA'); # undef again
- `write_table`'s `undef.val` is a single token, being what it writes;
`read_table`'s takes a list, being every token it should recognise.
- The consequence had been quiet rather than loud. An unmapped `NA` is an
ordinary string, and Perl numifies a string to `0`, so it did not stop `mean`
or `sd` — it pulled them toward zero, warning once per cell on `STDERR` and
dying only under `warnings FATAL => 'all'`. On the `titanic.csv` in the
repository root, whose 2207 rows carry 3690 literal `NA` cells across six
columns:
- | column | `NA` cells | mean, unmapped | mean, mapped |
|---|---|---|---|
| `age` | 2 | 30.416855 | 30.444444 |
| `fare` | 916 | 19.540347 | 33.404760 |
| `sibsp` | 900 | 0.295877 | 0.499617 |
| `parch` | 900 | 0.228364 | 0.385616 |
- The mapped column is `mean(d[[col]], na.rm=TRUE)` in R 4.6.1 to every digit
shown, for all four.
- Three choices are worth stating, because each follows R where pandas differs,
and each is asserted in the test rather than left to be discovered:
- **It is off by default.** R's `read.table` defaults to `na.strings = "NA"`
and pandas recognises a 19-token set, but `read_table` recognises nothing
beyond the empty field unless asked. A literal `NA` is real data in every
script written against an earlier release, and `t/read_table.t` pins that;
changing the default is a separate decision from adding the option.
- **The list replaces the set rather than extending one.**
`'na.strings' => 'baz'` maps `baz` and leaves `NA` and `NaN` alone. pandas
adds `baz` to its defaults and maps all three.
- **The match is on the exact field text**, case-sensitively and with no
whitespace stripped. `' NA'` is not `'NA'`, and `'-999.000'` is not
`'-999.0'` — pandas numifies both sides and matches them; R compares the
string and does not, and neither does this.
- The mapping happens where the empty-field rule already did, so it reaches every
`output.type`, the `.xlsx` reader, and a value a `filter` writes back. The
header is never mapped, so a column may still be named `NA`. A `filter` runs
after the mapping and so sees `undef` rather than the token, which makes
`sub { defined $_ }` the way to keep the rows that have a value. For `hoh`, a
mapped row-name cell is a missing row name and is refused exactly as an empty
one is.
- Everything above is checked in the new `t/read_table.na_strings.t` (95 tests),
whose cases are taken from the reference suites rather than invented here: R's
`tests/reg-tests-1a.R` — the `read.table(na.strings="foo")` case at 1911-1918
and PR#6781 at 2899-2901 — and `src/library/utils/man/read.table.Rd`, plus
pandas 2.2.3's `pandas/tests/io/parser/test_na_values.py` (`test_string_nas`,
`test_detect_string_na`, `test_non_string_na_values` for gh-3611,
`test_default_na_values`, `test_custom_na_values`, `test_na_trailing_columns`,
`test_na_values_scalar` for gh-12224, and `test_na_values_dict_aliasing`). The
expected values are frozen literals with their provenance in the file header,
so the test needs no R and no Python to run. Three divergences from pandas are
asserted deliberately, the three above; a fourth is that a row short of the
header is still an alignment error naming the row rather than the NA padding
pandas does.
- The full suite is 130 files and 26,444 tests, and `./test.all.perls.pl` passes
on all five local perls — `5.10.1`, `5.12.5` (long double), `5.42.3`, `5.44.0`
and `5.44.0-quadmath` — with no warnings on any of them.
[Distribution functions]
- Eight new functions — `qnorm`, `pt`, `qt`, `pchisq`, `qchisq`, `pf`, `qf` and
`pbinom` — with R's names, R's positional-or-named parameters and R's
`lower.tail` / `log.p` flags. Before this, `pnorm` and `dnorm` were the whole
family, so a Wald bound or a likelihood-ratio test could not be finished from
outside the module: the critical value simply was not expressible. See
[Distribution functions](#distribution-functions) for the details; two are
worth pulling out here.
- No tail is formed as `1 -` the other one. That subtraction discards
everything below machine epsilon, which is the range a p-value is interesting
in, so each function is routed to the parameterisation that computes the
requested side directly — `pt(q)` is `pt_upper(-q)` by symmetry, and `pchisq`
lower goes through a new `igam()` (the lower incomplete gamma that `igamc()`
had always formed internally and thrown away):
pchisq(1e-30, 1); # 7.978845608028654e-16, where 1 - igamc gives 0
- `qf` disagrees with R in the far tail of *F*, and is right. Arbitrated by
`t/distributions.mpmath.py` at `mp.dps = 60`, bisecting the defining equation
rather than calling a library inverse:
- | | R 4.6.1 | this module |
|---|---|---|
| `qf(2^-20, 1, 10)` | 9.9e-4 out | 1.0e-15 out |
| `qf(2^-20, 1, 1)` | 5.2e-5 out | 4.1e-15 out |
| `qf(2^-16, 1, 10)` | 2.7e-6 out | 6.9e-16 out |
- R's own `pf` settles it: `pf` of this module's answer returns the target
probability to 1.6e-15, and `pf` of R's own `qf` answer returns it to 4.9e-4.
It is R's inverse in that tail, not its CDF. Thirteen such rows are asserted
twice in `t/distributions.R.scipy.t` — right against the 60-digit value, and
still disagreeing with R — so if a future R fixes `qf`, or this module
regresses onto R's answer, the row fails rather than quietly agreeing. Three
rows go the other way and are recorded the same way rather than tolerated.
- `log => 1` on a `p*` function is `log()` of the probability already computed,
not a log carried through the series, so a tail that rounded to `1` logs to `0`
where R reports a tiny negative number. Four of SciPy's six `chi2` log
references land in that case and are pinned as exactly `0` with the true value
named. `ncp` and `qbinom` are not implemented.
- `t/distributions.R.scipy.t` is 2041 tests: a frozen table of R 4.6.1 values at
dyadic inputs (generated by the committed `t/distributions.R.scipy.R`), SciPy
1.18.0's own mpmath cases from `TestChi2` and `TestStudentT`, and R's own
round-trip identities from `tests/d-p-q-r-tests.R`. Every tolerance is the
worst error actually measured, printed by the test's own diagnostic.
- Ten figures now sit with the documentation — one per function, `dnorm` and
`pnorm` in their own sections and the eight new ones under
[What each one does, in one picture](#what-each-one-does-in-one-picture). Each
shades the region its function integrates and writes the integral it evaluates,
so the `p*`/`q*` pair reads as one identity answered from opposite ends: given
the boundary, return the area; given the area, return the boundary. `pbinom` is
drawn as bars with a summation sign rather than an integral, because a binomial
is discrete and drawing an integral over a histogram would misstate what it
computes. `distribution.plots.pl` regenerates all ten; like the other plot
scripts it is author-only and `[PruneFiles]` keeps it out of a release.
- Two portability problems surfaced while pinning this and are worth recording,
because neither is about the C code. Perl 5.10.1's string-to-number conversion
does not parse a 17-digit literal exactly — it reads `0.9999847412109375` one
ulp out, and reflecting a probability onto `1 - p` turns that ulp into 7.3e-12
relative — and it reads `9.332636185032188789901e-302` as
`8.999999999999996e-302`, a 3.7% error. So every argument and every
exactly-representable expected value in `t/distributions.R.scipy.t` is written
as a dyadic ratio, `M/2^K`, which rebuilds by integer division on any perl. And
where R's frozen value is `0` because a double underflowed, a long-double or
`__float128` build computes the real number instead — `pbinom(0, 1000, 0.75)`
really is `0.25^1000` = 8.7e-603 — so the test accepts that as agreement with a
reference that could not go lower, rather than failing the wider build for being
right. The same applies to `log => 1` near a probability of 1, where the
resolvable accuracy is `NV_EPSILON / |log p|` and the tolerance is scaled off
measured machine epsilon accordingly.
[Every confidence limit now comes from the same critical value]
- Adding `qnorm` made a discrepancy visible that had always been there. Nine
places in `LikeR.xs` built their own `z` by calling `inverse_normal_cdf()` —
Moro's rational approximation, of which the file's own comment said:
*"good to about 1e-9, which is fine for a confidence limit"* — while `qnorm`
went through the Newton-polished `normal_quantile_hp()`. Which function you asked
therefore decided how many digits you got, and a Wald bound written by hand
from `qnorm` differed from the one `glm` reported in the tenth digit.
- It is one function now. `std_qnorm()` seeds with Moro, polishes with Newton
against the `erfc`-based normal CDF, and reflects a probability above `0.5` onto
`1 - p` before inverting; `glm`, `cor_test`, `prop_test`, `epi_2x2`,
`cmh_test`, `roc`, `survfit`, `coxph` and `wilcox_test` all call it, `qnorm`
*is* it, and two near-duplicates are gone: `wilcox_qnorm()` in the XS, which was
the same Newton loop without the reflection, and the pure-Perl `_qnorm()`
(Acklam's approximation, also `~1e-9`) that `cohen_d` used for its interval.
- Measured against `mpmath` at `mp.dps = 60`, bisecting `erfc(-z/sqrt(2))/2 = p`
rather than calling a library inverse, worst relative error over conf.level
0.8 to 0.9999:
- | | worst relative error | |
|---|---|---|
| 0.302 and earlier, at these nine sites | 1.8e-9 | 8.0e6 ulp |
| R 4.6.1 `qnorm` | 4.8e-16 | 2.1 ulp |
| SciPy 1.18.0 `norm.ppf` | 1.5e-16 | 0.7 ulp |
| 0.303 | 8.2e-17 | 0.4 ulp |
- Seven orders of magnitude, and the result is closer to the truth than either
reference. Nothing in the existing suite failed: every one of these functions is
cross-validated against R at tolerances that had room for `1e-9` in them,
because they had to. What changed is that they no longer need it —
`t/distributions.R.scipy.t`'s `glm` agreement check went from `5e-9`, with its
comment explaining that the gap was the design, to bit-identical.
- The bug this exposed: `0.95` is a `double`:
- With the critical value down to 0.4 ulp, the next-largest error in a confidence
limit turned out to be the confidence level. Eleven functions declared their
default as
NV conf_level = 0.95;
- and `0.95` in C is a literal of type `double`, so on a long-double or
`__float128` perl what lands in the `NV` is `0.95` rounded to 53 bits — 2.2e-17
from the `NV` nearest `0.95`. `p = 1 - (1 - conf_level)/2` inherits all of it,
the normal density at the 95% point is `0.0584`, and the critical value comes
out 3.8e-16 wrong: about 1700 ulp of a long double, and the single largest error
in a 95% interval on those builds. `fisher_test` and `binom_test` had dodged
this for releases by parsing the literal through perl's own number parser; that
idiom is now the `NV_CONF_95` macro and all thirteen sites use it.
- This was invisible while the critical value itself was only good to `1e-9`.
It is the reason `t/qnorm.crit.R.scipy.t` freezes `p` as well as `conf.level`:
`1 - (1 - 0.999)/2` has to round, and it rounds differently at each NV width, so
a shared expected value is only meaningful if every build is asked the same
question.
- `t/qnorm.crit.R.scipy.t` is 166 tests — the value itself against R 4.6.1,
SciPy 1.18.0 and the 60-digit arbiter, and then the same value recovered out of
seven functions' reported intervals at six confidence levels each. It also pins
what `inverse_normal_cdf()` alone returns and asserts it is nowhere near the
tolerance, so the file cannot pass by accident if the routing is ever undone.
[Result keys are dotted everywhere]
- Every key of every returned hash now uses R's dotted spelling. Previously the
module answered to both conventions and it was impossible to guess which:
`p_value` from nine functions but `p.value` from two, `conf_int` from four but
`conf.int` from three, `prop_test` returning `conf.int` and `p_value` in the
same hash, and `shapiro_test` returning both spellings of its own p-value.
- | was | is |
|---|---|
| `p_value` | `p.value` |
| `conf_int`, `conf_level` | `conf.int`, `conf.level` |
| `null_value`, `null_value_name`, `statistic_name` | `null.value`, `null.value.name`, `statistic.name` |
| `estimate_x`, `estimate_y` | `estimate.x`, `estimate.y` |
| `group_stats`, `std_err` | `group.stats`, `std.err` |
| `odds_ratio`, `risk_ratio`, `risk_diff` (+ `_ci`) | `odds.ratio`, `risk.ratio`, `risk.diff` (+ `.ci`) |
| `risk_exposed`, `risk_unexposed` | `risk.exposed`, `risk.unexposed` |
| `n_pos`, `n_neg`, `auc_se`, `auc_ci` | `n.pos`, `n.neg`, `auc.se`, `auc.ci` |
| `n_active`, `n_inactive`, `n_top`, `active_count`, `enrichment_factor`, `rie_min`, `rie_max` | the same with dots |
| `n_risk`, `n_event`, `n_censor` | `n.risk`, `n.event`, `n.censor` |
| `exp_coef`, `lr_stat`, `lr_p_value`, `lr_df`, `loglik_null` | `exp.coef`, `lr.stat`, `lr.p.value`, `lr.df`, `loglik.null` |
| `has_na`, `old_coords` | `has.na`, `old.coords` |
| `p_adjust` (`dunn_test`'s column) | `p.adjust` |
- **This breaks code that reads the old names, and it is the only change here
that does.** `$r->{p_value}` is now `$r->{'p.value'}` — note the quotes, which a
dotted key requires. `table_one`'s output column is renamed with everything
else.
- Argument names are untouched. They are a different surface, where the
underscore spelling exists on purpose as a synonym for the dotted one, so
`conf_level => 0.99` and `'conf.level' => 0.99` both still work everywhere they
did, and `wilcox_test(..., conf_int => 1)` still asks for an interval. The
duplicate stores are gone: `shapiro_test` and `kruskal_test` each wrote their
p-value twice under the two spellings and now write it once.
[Decorative comment rules removed]
- 1,495 lines that were nothing but a comment marker and a run of dashes, and 380
comments whose text trailed off into one, are gone from `LikeR.xs`,
`lib/Stats/LikeR.pm` and `t/`. `/* ===== WILCOXON EXACT NULL DISTRIBUTIONS =====`
is now `/* WILCOXON EXACT NULL DISTRIBUTIONS`, which is what the house style in
`CLAUDE.md` asks for anyway. Nothing after `__DATA__` was touched, no POD
region, and no line holding a `/*` or `*/` was deleted, so no comment can have
been left unterminated. The suite is unchanged at 28,485 tests either side of it.
0.302 2026-08-22 CDT
- `sum`, `min`, `max`, `mean`, `sd`, `var`, `quantile`, `cor` and `cov` are
between 2.5 and 10.7 times faster, and four of the nine are now faster than R.
`merge` is between 1.9 and 5.3 times faster and now beats R's `merge` at every
size measured, and `uniq` is about four times faster at every size measured. Nothing any of them returns changes, except where it was wrong:
five bugs turned up while the measurements were being taken. One had been
hanging the interpreter, one had been segfaulting, and three had been quietly
returning the wrong frame for any input built on a tied array.
- The distribution had also been building without an optimizer, which is where a
third of this came from and none of the credit does.
- At `n = 1,000,000`, one column of `rnorm` values, best of seven, measured by
`scale.pl` against `scale.R` and `scale.py` on the same machine:
- | | 0.301 | 0.302 | R 4.6.1 | NumPy 2.5.2 |
|---|---|---|---|---|
| `sum` | 4.07 ms | 1.37 ms | 0.795 ms | 0.202 ms |
| `min` | 4.29 ms | 1.51 ms | 1.06 ms | 0.160 ms |
| `max` | 4.25 ms | 1.50 ms | 1.06 ms | 0.160 ms |
| `mean` | 4.07 ms | 1.37 ms | 1.59 ms | 0.204 ms |
| `sd` | 8.25 ms | 2.39 ms | 2.80 ms | 0.960 ms |
| `var` | 8.13 ms | 2.40 ms | 2.79 ms | 0.962 ms |
| `quantile` | 160 ms | 14.9 ms | 20.6 ms | 15.0 ms |
| `cor` | 20.5 ms | 8.24 ms | 6.43 ms | 4.67 ms |
| `cov` | 27.8 ms | 8.75 ms | 4.82 ms | 4.79 ms |
- The ratios are better below a million, where the column still fits in cache. At
`n = 100,000`, `sum` and `mean` are 5.1 times faster than 0.301, `sd` and `var`
5.2, `quantile` 8.0, and `sum` overtakes R.
- What is left between this and NumPy is not a constant factor waiting to be
found. A perl array is an array of `SV*`; NumPy holds a contiguous `double[]`
and runs SIMD over it. 1.2 ns per element is about four cycles — load the
pointer, load the flags, branch, load the `NVX`, add — and that is near the
floor for walking a million SVs. Closing it would need a packed-buffer type
holding `NV*` directly, which is a different module rather than a tuning pass.
[The module was being built without an optimizer]
- `compile.sh` ran `perl Makefile.PL OPTIMIZE='-Wall'`. `OPTIMIZE` *replaces*
perl's own `$Config{optimize}` rather than adding to it, so that was not "the
usual flags plus a warning switch" — it was the whole optimize line, and every
build it made, including the one `make install` then installed, was compiled at
`-O0`. `test.all.perls.pl` defaulted to the same string, so the entire support
matrix was built and timed that way too.
- Both say `-O2 -Wall` now. On the same `LikeR.c`, summing `1e6` NVs takes 4.10 ms
at `-O0` against 2.63 ms at `-O2`, and `cor` 25.5 ms against 12.9 ms. Read the
rest of this section as what was left after the compiler was allowed to do its
job.
- `t/build.optimize.t` reads both scripts and fails if either hands `Makefile.PL`
an `OPTIMIZE` with no `-O` in it. There is no way to ask a loaded `.so` what it
was compiled at, so the check has to be on the scripts; neither ships, so the
file skips outside the repository.
[Reading a column: one `av_fetch` per element]
- All nine functions read their input by calling `av_fetch` once per element. That
call is out of line and repeats, for every element, a bounds check, a
negative-index fixup and a magic dispatch that a straight walk of `AvARRAY`
already knows the answers to. On a million-element column it was about half of
what `sum` cost and more than half of `cor`.
- Walking the block directly is only sound while nothing can move it, and `SvNV`
can move it: on an SV carrying get magic it runs `mg_get`, on a reference it may
run an overloaded `0+`, and on a string that is not a number it raises a warning
that a `$SIG{__WARN__}` handler catches. Each of those is perl code, and perl
code may push to the very array being walked, which reallocs the block and
leaves the walk holding freed memory.
- So the fast path is decided per element, not per array. `sv_plain_nv` takes an
SV only when it already holds a number and carries neither magic nor a
reference, which makes `SvNV` a register read that calls nothing; everything
else — a hole, a string, a reference, a tied element — goes back to `av_fetch`,
which is both correct and where the block pointer gets re-derived. A tied array
is refused whole, before the scan starts.
- The scans live in their own functions and are marked `LIKER_NOINLINE`, which is
not decoration. Inlined into the caller, a scan shares a stack frame with the
`av_fetch` loop that finishes whatever it would not take; the accumulators then
have to stay addressable across that call, so the scan spills them to memory
every iteration instead of keeping them in registers. Summing `1e6` NVs: 2.77
ns/element for the `av_fetch` loop being replaced, 2.25 with the scan inlined
into it, 1.29 with the scan compiled on its own. There is no portable spelling
of `noinline`, so it is annotated where the compiler has one and left empty
where it does not, which gives back the inlined speed and never a wrong answer.
- `min` and `max` get a scan each rather than sharing one that tracks both ends.
The running extreme is a loop-carried dependency, each compare-and-move waiting
on the previous element's result, so carrying the end the caller will not return
costs a second chain for nothing: 1.79 ms against 1.35 ms on `1e6` NVs.
[`var`, `sd` and `cov`: Welford to two passes]
- Welford's recurrence updates the running mean with `mean += delta / count`. That
divide sits in a loop-carried dependency chain, so each element waits out the
previous one's latency — 5.34 ns/element against 2.70 for two passes.
- Two passes is also what R computes. `src/library/stats/src/cov.c` takes a first
mean, refines it with a second pass (`tmp = tmp + sum(x - tmp)/n`, the `MEAN`
macro) and only then sums the squared deviations about the refined mean.
Expanding that third pass about the *unrefined* mean leaves
`sum((x-m)^2) - (sum(x-m))^2/n`, so carrying the sum of the deviations alongside
the sum of their squares gives R's answer in one pass fewer. That is the
Chan-Golub-LeVeque correction, so it is no less accurate than Welford on an
ill-centred column.
- `mean` is deliberately not changed: it remains the single-pass `sum/n`, which is
not what R's `mean` computes and never was.
- Reading the arguments twice is visible only to a tied array, whose `FETCH` now
runs once per pass.
[`quantile` orders only the statistics it reads]
- `quantile` sorted the whole column and then indexed into it. Type 7 reads at
most two order statistics per probability — `x[j]` at `j = floor((n-1)p)`, and
`x[j+1]` where the interpolation weight is non-zero — which for the five default
probs is ten values out of a million.
- `nv_select_multi` places just those. It takes the middle index first and
recurses into the two ranges that index separates, costing `O(n log k)` rather
than the `O(n k)` a left-to-right sweep of selects would: each level of the
recursion touches each element at most once, and the index list halves at every
level. Past `k = n/2` distinct indices a full sort is the cheaper way to get the
same array, and that is what still happens there.
- This is also what R does — `quantile.default` hands its indices to
`sort(partial=)`, and NumPy reaches `np.partition` for the same reason. On `1e6`
random doubles the ten order statistics of the five default probs cost 19.2 ms
this way against 69.1 ms for the full introsort.
- The probabilities are parsed before the column is read now, rather than after it
is sorted, so an invalid `probs` vector no longer costs a pass over the data
first.
[`cor` and `cov`: sorts that inline their comparison]
- Spearman ranks its input before correlating it, and `rank_data` did that through
`qsort` with an out-of-line comparator; Kendall's tau-b sorted three times the
same way. That comparator call cannot be inlined, and for a sort of a few
thousand elements it is most of what the sort costs — the module's own `nv_sort`
takes 210us where `qsort` takes 355us on 5,000 NVs, which is why it exists.
- `LIKER_DEFINE_SORT(T, PFX, LESS)` generates that introsort for a given struct
and ordering: median of three left at the midpoint so both partition scans have
a sentinel, recursion into the shorter side and a loop on the longer so the
stack stays `O(log n)` deep, a heapsort fallback past `2*floor(log2 n)` levels
so a median-of-three killer cannot reach `O(n^2)`, and one insertion sort over
the nearly-ordered remainder. It is `kw_sort`, the Kruskal-Wallis sort, made
generic. `rank_data` and Kendall's `(x, y)` sort are generated from it; Kendall's
third sort, a plain `NV` column, calls `nv_sort` directly.
- glibc's `qsort` is also a mergesort that allocates a scratch buffer the size of
the array it is given. Peak RSS for a Spearman correlation of two `2e6` columns
falls from 455 MB to 424 MB — exactly the 32 MB it was allocating to sort `2e6`
sixteen-byte records.
- At `n = 200,000`: `cor(..., 'spearman')` 66.9 ms to 27.2 ms, `cov(...,
'spearman')` 65.9 ms to 26.7 ms, `cor(..., 'kendall')` 70.5 ms to 44.7 ms. Small
repeated calls gain the same factor: ten thousand `cor(..., 'spearman')` calls
over twelve elements, the shape `agg` and `group_by` produce, 6.4 ms to 3.1 ms.
[`cor` read an array of arrays once per column]
- `cor(\@matrix)` extracted one column at a time, and for each column it walked
every row — `av_fetch` on the outer array, `SvRV` to reach the row, then
`av_fetch` for the cell. The outer array was therefore walked `ncols` times
over, and every one of the `ncols * nrows` cells paid a fetch on it that the
column before had already paid.
- The pre-validation pass that already checks every row is an array reference now
keeps the row `AV*`s it resolves instead of discarding them, and the extraction
transposes in a single pass over the rows through `av_extract_or_nan` — the same
reader the vector branch uses, so a short row, a hole, an `undef` and a
non-numeric cell all behave here exactly as they do there. Outer fetches fall
from `ncols * nrows` to `nrows`, and each row's cells are read in the order they
are stored rather than one per pass.
- At 20,000 rows: `cor` of a 40-column matrix 40.6 ms to 26.7 ms, and its
cross-correlation against an 8-column matrix 28.7 ms to 12.4 ms. The transpose
keeps one write stream open per column, which is what a wide matrix pays for
this, and 2,000 x 300 still comes out ahead — 144.6 ms to 134.8 ms.
- A row that was itself a tied array used to read as all-`NA`, because `av_fetch`
on one hands back a deferred `PVLV` whose value only arrives when `mg_get` runs:
`cor(\@rows)` croaked "standard deviation is 0 in x column 0" on data that was
perfectly well defined. `av_extract_or_nan` runs the get magic, so those rows
now give the same answer as the untied equivalent. A cell that is a tied scalar
was already read correctly and still is.
- The zero-variance check that used to be folded into the extraction is now
`nv_all_equal`, one function shared with the vector branch, which returns at the
first pair of values that differ instead of sweeping the column. On the vector
path that replaces two full passes over `x` and `y`, which is where the `cor`
row of the table above moves from 10.5 ms to 8.24 ms — the two passes were a
fifth of what `cor(\@x, \@y)` cost at `n = 1,000,000`.
- `cov` does not gain from any of this: its Pearson path neither sorts nor scans
for zero variance, and 8.75 ms is the same figure it read before. What changed
there is housekeeping — `Newx` in place of a bare `safemalloc`, and the method
string resolved once rather than by three more `strcmp` calls.
[`min` over more than 65,535 arguments never returned]
- `min` counted its arguments in an `unsigned short int`. `items` is the whole
flattened argument list, so `min(@x)` on an array of more than 65,535 scalars
wrapped the counter and the loop never reached its bound:
my @a = (1 .. 70_000);
min(@a); # 0.301: hangs, forever
- It is a `Stack_off_t` now, which is what `max` and `sum` already used.
`power_t_test` carried the same counter; it takes named pairs, so nothing could
reach 65,535 of them, and it is widened anyway.
- A regression here is a hang rather than a wrong answer, so `t/hot_path.t` runs
the 70,000-argument calls under `alarm` where the platform has one. Without that
guard a smoker would sit there instead of failing.
[A tied array read as all-`undef`]
- `sum(\@tied)` croaked `undefined value at array ref index 0`, and so did every
other function here. `av_fetch` on a tied array does not return the element: it
returns a mortal `PVLV` that only acquires the element's value once `mg_get` has
run on it. The `SvOK` test came first, saw an empty `PVLV`, and reported every
element of every tied array as undefined. An ordinary array holding a tied
scalar failed the same way.
- The slow path runs `SvGETMAGIC` before it looks now. This predates 0.302 and is
not fallout from the scan — but the scan is what made it worth finding, because
`sv_plain_nv` refuses every magical SV and so sends all of them down that one
path.
- That slow path then had a bug of its own, and it was worse. `av_slow_at` held
the `SV **` that `av_fetch` returned — a pointer *into* `AvARRAY` — ran the get
magic through it, and read it again afterwards. `mg_get` runs perl, and the
element in `t/hot_path.t` pushes 200 values onto the array being walked, which
reallocates that block: both reads after the magic were reads of freed memory.
It usually returned the right answer because the freed block usually still held
the old pointer, so the test that exists for precisely this case passed on four
of the five perls and, on the fifth, only failed when the suite was run four
perls at a time and something else had claimed the block —
`sum` returned 1085 instead of 10, once. Under `valgrind` it is an "invalid
read of size 8" on every perl, every time.
- It loads the element into a local before running the magic now. The SV itself
does not move; only the array of pointers to it does, so a copy of the pointer
stays good where a second read through the slot does not. The same shape — take
`SV **` from `av_fetch`, run magic, dereference again — is still present in the
`median` and `moment` readers and in three places in `transpose`, and is the
same one-line fix in each; those are untouched here.
- `uniq` was missed by this pass and fixed later in the same release; see
"`uniq`: the same keys, without a perl hash" below.
[`cor` returned `NaN` on an off-centre column]
- `pearson_corr` computed the correlation from raw cross-products,
`(n*sxy - sx*sy) / sqrt((n*sx2 - sx^2) * (n*sy2 - sy^2))`. Each of those three
terms subtracts two nearly equal large numbers, so it loses a digit of the
answer for every digit by which the mean exceeds the spread. On a column with
mean `1e9` and a spread of `0.05` it does not merely lose precision:
`n*sx2 - sx*sx` comes out negative, the square root is `NaN`, and `cor` returns
`NaN` for an ordinary pair of vectors. At mean `1e6` it was wrong in the fifth
decimal place — `0.74945606763640826` against R's `0.74943997817595964`.
- It centres first now, as R's `cov.c` does, with the same compensated two-pass
correction `var` uses. R reports `0.063454463733801467` for the `1e9` case, and
so does this.
- `cov`'s Pearson branch was Welford's covariance recurrence, two divides deep in
the dependency chain; it is `nv_cov2` now, the two-pass form, shared with the
Spearman branch.
[A leak on `cor`'s matrix error path]
- `cor(\@x, \@y)` with a `y` matrix whose rows are empty croaked "y matrix has
zero columns" without freeing the `x` columns it had already extracted. On a
40-column matrix that is 40 buffers and their table per call: 20,000 such croaks
grew RSS by 44 MB. The branch's allocations are now released through one macro
that every croak path calls, and the row tables and the scratch row are freed as
soon as the last value has been read out of them rather than at the end.
- The `nrows < 2` guard was missing the same frees, but it is unreachable — a
one-row matrix has a constant column, so the zero-variance croak above always
gets there first. It frees anyway.
- Seven `croak` formats in `cor` passed a `size_t` to `%lu` with no cast, which on
Windows is a 64-bit argument read through a 32-bit conversion. They use
`%" UVuf "` with a `(UV)` cast now, as the rest of the file does, and `cov`'s
one — cast to `unsigned long`, so well defined but still narrowing on Win64 —
went with them. Every message is unchanged wherever `unsigned long` is 64 bits,
which is everywhere the suite runs.
[Unchanged: `cor` and `cov` still do not see an overloaded object as a number]
- `cor` and `cov` decide what is numeric with `looks_like_number`, which is false
for any reference, an overloaded `0+` included, so a column of overloaded
objects reads as all-`NA`: `cov` returns `NaN` and `cor` reaches its
zero-variance croak. `sum`, `min` and `max` do honour the overload, through
`SvNV`. The two families genuinely disagree. That predates 0.302 and is a
different question from how a column is walked, so `t/hot_path.t` pins the
current behaviour rather than changing it quietly.
[`merge`: a hash join that stopped building a perl hash]
- At `n = 300,000` a `merge($df, $df, how => 'inner', on => 'id')` over a
six-column HoA took 575 ms. It takes 109 ms. Best of ten, one process per
measurement, pinned to one CPU, against `scale.R` and `scale.py` on the same
machine:
- | `n` | 0.301 | 0.302 | R 4.6.1 | pandas 2.x |
|---|---|---|---|---|
| 1,000 | 0.471 ms | 0.245 ms | 0.46 ms | 0.68 ms |
| 10,000 | 5.45 ms | 2.60 ms | 3.35 ms | 0.88 ms |
| 100,000 | 138 ms | 33.3 ms | 80.0 ms | 2.95 ms |
| 300,000 | 575 ms | 109 ms | 155 ms | 9.83 ms |
- Under callgrind the 30,000-row join fell from 358M instructions to 170M. Three
things were paying for that.
- The right-hand index was a perl `HV` keyed by the join key: an `HE`, a
shared `HEK` and an `IV` SV per distinct key, so three allocations per right row
whenever the key is unique, which is exactly what a join on an id column is.
It is the same open-addressed table over a key arena that `drop_duplicates`
already interns into (`dd_ctx`) — no SV, no HEK, one `memcpy` of each distinct
key — so the two functions now share their definition of "the same key" as well
as their code for deciding it. `mg_key` writes straight into that arena instead
of making four `sv_catpvn` calls per key column into a scratch SV, and a
one-column join writes the cell's bytes bare: with a single field there is
nothing for a length prefix to disambiguate, and on an integer id the prefix and
its two separators were about as many bytes again to hash and to compare.
- It emitted each output row as it found it, so the row count was not known
until the join had finished and every output column grew by `av_push`, being
reallocated its way up to the answer. The probe now records `(left row, right
row)` pairs into a flat list and builds nothing; only when that list is complete
— so the count is exact — is each column allocated once at its final size and
filled straight into `AvARRAY`, the way `filter` already builds its columns.
The list costs two `SSize_t` per output *row* against roughly one SV per output
*cell*, so it is a small fraction of the result on any join wide enough to care,
a cross join included. `mg_emit` is gone; `mg_column` and `mg_build` replace it,
and the cross join goes through the same two stages as every other `how`.
- `mg_cell` called `av_fetch` once per cell, which was 7% of the call: a join
reads every key cell of both frames and every cell of every column it keeps.
It reads the block directly now, through the same `av_at` the other frame
functions use.
- What is left is the result itself. A 300,000-row join of two six-column frames
is 3.3 million output SVs, and at 109 ms that is 33 ns each — allocate, copy a
cell into it, store it. pandas is not doing the same work: its columns are
contiguous typed buffers and a join is a `take` over them.
[`uniq`: the same keys, without a perl hash]
- `uniq` on a million-element column of `rnorm` values took 0.80 s, against R's
`unique()` at 0.019 and pandas' `pd.unique()` at 0.029. It takes 0.16 s. Same
method as the tables above — one process per measurement, pinned to one CPU,
the fastest of seven, `plot.scaling.pl` against `scale.R` and `scale.py` on the
same machine — except that the `0.301` column is the 0.301 *code* built at `-O2`,
not the released build, so this is the rewrite on its own and not the optimizer
again:
- | `n` | 0.301 code | 0.302 | R 4.6.1 | pandas 2.2.3 |
|---|---|---|---|---|
| 1,000 | 0.294 ms | 0.069 ms | 0.009 ms | 0.020 ms |
| 10,000 | 3.19 ms | 0.825 ms | 0.092 ms | 0.184 ms |
| 100,000 | 41.2 ms | 9.71 ms | 1.16 ms | 2.30 ms |
| 300,000 | 174 ms | 39.2 ms | 3.71 ms | 5.59 ms |
| 1,000,000 | 801 ms | 160 ms | 19.4 ms | 28.6 ms |
- The slope was never the problem. Fitted over the ladder by `scaling.slopes.tsv`,
all three are linear — 1.12 for `uniq`, 1.11 for R, 1.04 for pandas, and 1.14
for `uniq` before the rewrite — so this was a constant factor from the start,
and none of it was the hashing.
- `SvPV` on an `NV` is a `%.15g` that has to be redone every time. Walking a
million `NV` SVs costs 2 ms; walking them and calling `SvPV` on each costs
268 ms. That is ten times R's whole runtime for the same call, and it is the
floor the old code could not get under however few distinct values there were:
a column of a million elements holding ten distinct values still took 0.226 s.
Worse, `SvPV` leaves the rendered buffer on the caller's own SV without ever
reusing it (`sv_2pv_flags` re-renders an `NV` on the next pass regardless), so
asking a large numeric column for its distinct values grew that column by tens
of megabytes, for good. `nk_num_pv` — already written for `drop_duplicates`,
about four times faster and only taken where the answer is provably the same —
renders into the XSUB's own stack buffer and touches nothing.
- The `seen` hash was a perl `HV`. An `HE` and a copied `HEK` per distinct
key is ~72 MB of scattered small allocations at a million of them, and the
table rehashes at every doubling on the way there. Going from ten distinct
values to a million added 0.47 s of pure insert cost. It is `dd_ctx` now, the
same open-addressed slot array over a key arena that `drop_duplicates` and
`merge` intern into, presized from the element count.
- It hashed each key twice, once for `hv_exists` and again for `hv_store`.
Collapsing that into one `hv_fetch(..., 1)` is the obvious repair and is the
wrong one: the lvalue fetch mints an SV per key, and on the million distinct
doubles it measured 0.919 s against the pair's 0.756 — while also reading
`AvARRAY` directly, so the comparison flatters it. That is why the pair had
survived. One hash of one key into one open-addressed probe is what the arena
gives instead.
- What the rewrite does *not* do is compare doubles by value. That is how R and
pandas get the factor of six to eight they still have, and it is a different
answer:
`0.1 + 0.2` and `0.3` are two doubles that print the same, so they are one
value to `uniq` and two to `unique()`. `uniq` is documented to compare the way
`eq` and `List::Util::uniq` do, so it renders every element and compares the
text; that rendering pass is essentially all of the distance that is left.
Scalar context builds no result list at all and takes about a fifth off again.
- `uniq(\@tied)` croaked `undefined value at array ref index 0`. It was reading
elements through `av_fetch` and testing `SvOK` on the `PVLV` that comes back,
which is the bug "A tied array read as all-`undef`" describes above. `uniq` was
simply not one of the functions that pass reached. It reads through
`av_slow_at` now, like the rest.
[`filter`: measured, and left alone]
- `filter` is not faster, and the reason is worth writing down so it is not
looked for again. At `n = 300,000` on a five-column HoA, `col('x') > 0` keeping
half the rows costs 22.5 ms, of which the predicate is 1.3 ms. The other
21.2 ms is allocating and filling the 750,000 SVs of the result, which the
documented contract requires: an HoA input, or any `hoa` output, builds fresh
arrays and fresh cell values. That is 19 ns per numeric cell and 44 ns per
string cell, and the difference between the two is one `malloc` for the string's
buffer.
- Copy-on-write is the obvious way out of that `malloc` and is not available.
`newSVsv` is `newSVsv_flags(sv, SV_GMAGIC|SV_NOSTEAL)`, which does not pass
`SV_COW_SHARED_HASH_KEYS`, so its string copies are never copy-on-write; and
`SV_DO_COW_SVSETSV`, the flag pair a plain `my $b = $a` gets, is defined as `0`
unless `PERL_CORE` is set — perl's own `sv.h` says "the core is safe for this
COW optimisation, XS code on CPAN may not be". Reaching past that by hand was
measured anyway and is slower in both directions: `filter` 23.0 ms against 28.3
ms, `merge` 103 ms against 126 ms. For a string as short as a data frame's
usually are, the `malloc` costs less than what COW puts in its place — an
out-of-line `sv_setsv`, the `CowREFCNT` increment, and the read-only flip on the
source. The finding is recorded at `flt_cell_copy` rather than only here.
- The one lever that does move it is the copy itself. Sharing the surviving cells
instead — one refcount bump, which is what `drop_duplicates` does for exactly
this reason — was prototyped and measured at 13.6 ms against R's 14.4 ms. It is
not done, because `filter` documents the opposite and someone may be relying on
it; changing that is a decision about the interface, not a tuning pass.
- Two things did change. `flt_num` takes a bare `IOK`/`NOK` cell without calling
`SvGETMAGIC` or `looks_like_number` — neither can say anything about an SV that
already holds a number, and `looks_like_number` is out of line — which is most
of what the predicate pass costs on a numeric column. And the tied-frame bugs
below.
[A tied frame read as all-`undef`, and one that segfaulted]
- `filter`, `merge` and `drop_duplicates` all walk a column or a row through
`AvARRAY`, which is why they are fast and which is wrong for a tied array: its
elements do not exist until `FETCH` has run, and `AvARRAY` on one is not the
block they live in — on an array that has never held a real element it is a null
pointer. `av_fetch` is the way in, and what it returns for a tied element is a
mortal `PVLV` that only acquires the value once `mg_get` has run on it, so an
`SvOK` or `SvROK` test placed before that says "undef" or "not a reference" for
every cell and every row of the frame.
- This is the same pair of mistakes `sum(\@tied)` made, above, and it predates
0.302 in all three functions. Between them they were wrong in four distinct
ways:
tie my @x, 'TiedArray'; @x = (1, -2, 3);
tie my @y, 'TiedArray'; @y = (10, 20, 30);
filter({ x => \@x, y => \@y }, col('x') > 0);
# 0.301: { x => [], y => [] } -- an empty frame, whatever the predicate
# 0.302: { x => [1, 3], y => [10, 30] }
merge({ id => \@x }, { id => [1, 3], w => ['a','b'] }, how => 'inner', on => 'id');
# 0.301: { id => [], w => [] } -- no key matched, because every key read undef
# 0.302: { id => [1, 3], w => ['a', 'b'] }
drop_duplicates({ k => \@x, v => \@y });
# 0.301: { k => [undef], v => [undef] } -- every row identical, so one survived
# 0.302: { k => [1, -2, 3], v => [10, 20, 30] }
drop_duplicates([ map { { k => $_ } } 1 .. 3 ]); # with the AoH itself tied
# 0.301: segmentation fault
- The crash was `_aoh_key_union`, which indexed `AvARRAY` with no bounds check and
no test for a tied array at all; on a tied AoH that is a null dereference on the
first row. An outer join was the quietest of the four — with no key matching, it
returned the two frames as disjoint halves, well-formed and entirely wrong.
- There is one reader now. `av_at` returns an element as it is — block read where
that is sound, `av_fetch` where the array is tied, and the pointer derived per
element rather than hoisted, because copying a cell can run perl and perl can
push to the array being walked. `av_ref_at` adds the `mg_get` and the
"is it a reference to the right thing" test for a row; `av_row_keep` is how a
row gets into a result, handing back the caller's own reference for a plain
array and a fresh reference to the same row for a tied one, since the `PVLV` is
mortal and stays bound to the tie. `dd_cell` runs the get magic before it looks
at the cell. None of it is measurable: `drop_duplicates` is within 2% of 0.301
on all three shapes at 30,000 and 300,000 rows, and `merge`'s numbers above are
with it.
- One consequence is visible and is documented under `drop_duplicates`: a tied HoA
column's cells cannot be *shared* into the result, because the tie has no cell
SV to share, so they are copied. It is the only place the "what survives is
shared" rule cannot hold, and copying is the only reading of it that can be
true.
- A tied *frame* hash — the outer hash of a HoA or HoH, rather than a column or a
row inside it — is still not supported by any of the three. All of them read
its shape and its columns through `HeVAL(hv_iternext(...))`, which for a tied
hash is not the value. That is a larger change than this one and is not made
here; `t/tied.frames.t` says so where it stops.
[Tests]
- Three new test files, and the suite is 129 files and 26,331 tests.
- `t/var_sd_cov.R.t` (162 tests) cross-validates `sum`, `min`, `max`, `mean`,
`var`, `sd`, `cov` and `cor` against R 4.6.1 over ten columns chosen to separate
the algorithms rather than to be representative: an arithmetic ladder, a
constant column, `n = 2`, mixed signs, and the column with mean `1e9` and spread
`0.05` that breaks a naive variance and did break the old `cor`. The expected
values are frozen literals with their provenance in the file header; the
generator, `t/var_sd_cov.R.R`, is committed beside it, and the test itself never
calls R.
- Every value in the corpus is a dyadic rational, so the vectors are the same
numbers at every NV width, and the corpus is rebuilt in perl rather than pasted
in — the `n`, `sum`, `min` and `max` columns are what pin the perl construction
to R's. Those frozen numbers are printed `%.40g`, the double's exact decimal
expansion, because 17 significant digits round-trip a double back to a double
and no further: a long-double perl reads `100000000004.83398` as a different
number from the `100000000004.833984375` R had, which was enough to fail an
exact comparison on two of the five perls. The tolerance on the exact columns is
two ulps of whatever NV the running perl was built with, found by bisection at
run time rather than hardcoded, because perl's string-to-NV conversion is not
correctly rounded on a long-double build — `perl-5.12.5` reads that literal one
ulp short.
- Every case also runs through a tied array, which cannot take the scan and has to
reach the same answer through `av_fetch`. That is the assertion that the fast
path and the slow path are the same function.
- `t/hot_path.t` (135 tests) covers how the input is read rather than what comes
out: the 70,000-argument counter; every kind of element the scan must refuse —
IV, UV above `IV_MAX`, NV, string, `0+` overload, tied element, and one magical
element in the middle of an otherwise plain column; holes and `undef` in all
three of their documented behaviours; and an element whose get magic pushes 200
values onto the array it lives in while that array is being walked. `cor`, `cov`
and `var` are checked for invariance under a location shift of `1e3`, `1e6` and
`1e9`, which is the property the old correlation failed, and `quantile`'s
partial sort is checked against a full sort at every rank on a column with ties.
- `t/tied.frames.t` (81 tests) covers `filter`, `merge` and `drop_duplicates`
over frames built on tied arrays and tied row hashes: every output shape, every
`how`, every `keep`, single and composite keys, an empty tied frame of each
shape, and the sharing rule on both sides of its one exception. Every case runs
the same call over tied input and over an identical plain copy and requires the
two answers to be equal — a fixed expected value would pin the plain answer as
well, which `t/filter.t`, `t/merge.t` and `t/drop_duplicates.t` already do, and
what has to be pinned here is that the two routes agree. `t/merge.t`'s existing
reference join — plain Perl, run over all six input/output shape combinations —
is what checks that the rewritten join still means what it meant.
- `./test.all.perls.pl` passes on all five local perls — `5.10.1`, `5.12.5`
(long double), `5.42.3`, `5.44.0` and `5.44.0-quadmath` — with no warnings on
any of them, and `LikeR.c` compiles clean under `-std=c99 -O2 -Wall` against
every one of their `CORE` directories.
0.301 2026-08-21 CDT
- there are numerous additions of `restrict` keywords, which may or may not improve speed
[kruskal_test]
- `kruskal_test` cross-validated against R 4.6.1's `stats::kruskal.test()` and
SciPy 1.18.0's `TestKruskal`, driven by those suites' own cases rather than by
cases invented here. Six bugs are fixed: five in how the arguments and the data
are read before any ranking happens, and one in the chi-squared tail, which
reaches every function that uses it. The test now needs a third of the memory
and runs in under half the time, and it returns the same answer twice in a row,
which it did not before.
- Everything below is checked in the new `t/kruskal_test.R.scipy.t` (870 tests),
whose expected values are frozen literals with their provenance recorded in the
file header; it needs no R and no Python to run. The generator that produced
them, `t/kruskal_test.R.scipy.R`, is committed beside it. There had been no
`t/kruskal_test.t` at all — the only coverage was six assertions in `t/01.t` on
the single Hollander & Wolfe example, and none of the six bugs would have shown
up in it. The full suite is 125 files and 25,951 tests, and
`./test.all.perls.pl` passes on all five local perls — `5.10.1`, `5.12.5`
(long double), `5.42.3`, `5.44.0` and `5.44.0-quadmath` — with no warnings on
any of them.
- NaN was ranked instead of dropped:
- `looks_like_number` is true for `NaN`, so a `NaN` went into the ranking. R
treats `NaN` as `NA` and `complete.cases(x, g)` removes it before `rank()` ever
sees it:
- | data | was | R 4.6.1 |
|---|---|---|
| `c(1,2,3,4,5,6)` with one `NaN`, n = 7 | `H = 4.5` | `H = 3.8571428571428577` |
| `1:24` with one `NaN` | `H = 17.28` | `H = 16.5` |
- It also handed `cmp_nv3` a comparison that is never true for any pair
involving the `NaN`, which leaves `qsort` without the strict weak ordering it is
entitled to — the same defect `wilcox_test` had fixed in 0.298, where the
comment describing it is still in the file. `+Inf` and `-Inf` are neither `NA`
nor `NaN` to R and a rank test has no trouble with them, so they are still kept
and ranked; that is now pinned rather than incidental.
- The bad value propagated: `table_one`'s `_t1_cont_p` hands its groups straight
to `kruskal_test`, so a single `NaN` among nine observations was reported at
`p = 0.0273` where dropping it gives `0.0439`.
- A group with no data inflated the degrees of freedom:
- The hash-of-arrays form counted every key in `k`, including a key whose array
was empty or whose every element had been dropped. On
`{a => [1,1,1], b => [2,2,2], c => []}` that gave `df = 2` and
`p = 0.0820849986238988` for what is a two-group problem with
`df = 1, p = 0.025347318677468304`. Such a group was already skipped when
forming the statistic and when building `group.stats`; only `df` still counted
it.
- R refuses this case outright — `all groups must contain data` — and so does
`kruskal_test` now, because the alternative is to test the groups that do have
data under a `df` that counts one that does not. R's order of checks came with
it: it filters each group, refuses an empty one, and only then counts what is
left, so `{a => [], b => []}` is `all groups must contain data` and not
`not enough observations`. SciPy takes the other side of this and returns `NaN`
with a `SmallSampleWarning`; the divergence is recorded in the test file rather
than papered over. The `x`/`g` form cannot reach any of this — it mints a group
id the first time an observation survives the filter — so nothing changes there.
- Group labels were truncated at a NUL and lost their UTF-8 flag:
- The `x`/`g` path read the label with `SvPV_nolen` and then took `strlen` of it.
Perl strings are counted, not NUL-terminated, so `"a\0X"` and `"a\0Y"`
collapsed into one group: `kruskal_test([1..6], ["a\0X","a\0X","a\0Y","a\0Y","b","b"])`
came back as two groups with `df = 1, H = 2.4` instead of three with
`df = 2, H = 4.571428571428573`. Both paths also copied the label's bytes while
dropping perl's UTF-8 flag when storing it into `group.stats`, so a label
outside latin-1 came back as mojibake and the two input paths disagreed with
each other about labels inside it. The length now travels with the string and
carries the flag in its sign, which is `hv_store`'s own convention, so a label
comes back `eq` to what went in. Dropping the `strlen` also drops a pass over
every label.
- A trailing named argument read past the argument stack:
- The named-argument loop took `ST(arg_idx + 1)` without checking that there was
one, so an odd argument list read one slot past the top of the stack — and what
it found there changed which branch ran: `kruskal_test(\%h, 'x')` came back
complaining that `'h'` cannot be mixed with `'x'`/`'g'`, because `x_sv` had been
assigned whatever was past the end. `binom_test`, `chisq_test`, `fisher_test`,
`wilcox_test`, `var_test` and `prcomp` all guard this; `kruskal_test` was the
one that did not. It now croaks `odd number of named arguments`.
- An infinite chi-squared statistic gave no p-value:
- `get_p_value` short-circuits a statistic at or below zero and otherwise goes to
`igamc`. `+Inf` is neither, so it reached the continued fraction, where the
first `1/d` is `1/Inf = 0` and then `del = 0 * Inf` is `NaN` — an overwhelmingly
significant result reported as no result at all. R's
`pchisq(Inf, df, lower.tail = FALSE)` is `0`, and so is this now. `NaN` in gives
`NaN` out, as R does, rather than running the continued fraction to its full
10,000-iteration safety bound first.
- This is reachable from `kruskal_test`. When a sample has no variation at all
the tie correction is `(n^3 - n)/(n^3 - n)`, and once `n^3` is past `2^53` the
subtraction of `n` is lost from one side or the other, so an inexact zero is
divided by an exact zero. R has the same problem and returns `+Inf`, `-Inf` or
`NaN` depending on which way `n` rounded: `NaN` at `n = 250000`, `-Inf` at
`300000`, `NaN` at `400000`, `+Inf` at `500000`, `NaN` at `750000`, `+Inf` at
`1000000`, `NaN` at `1500000` and `+Inf` at `2000000`. `kruskal_test` now agrees
with it on all eight. Below that the correction is exact and both give `NaN`,
which the corpus pins. `get_p_value` is shared, so `chisq_test`, `prop_test`,
`mcnemar_test`, `friedman_test`, `cmh_test`, `logrank_test` and `coxph` get the
same fix.
- Three times less memory, twice the speed:
- The ranking no longer goes through `RankInfo` and `rank_and_count_ties`.
`kruskal_test` wants per-group rank sums, not the ranks themselves, so it sorts
a 16-byte `(value, group)` pair and adds each tie block's averaged rank straight
into the group sums, instead of storing an `NV` rank per observation for a
second pass to read.
- It also no longer calls `qsort`. glibc's `qsort` is a mergesort that allocates a
scratch buffer the size of the whole array — measured on glibc 2.39 as a
`VmHWM` of 116 MB going to 230 MB across one sort of a 114 MB array — which was
half of the function's peak memory, and its comparison goes through a function
pointer that cannot be inlined. In its place is a median-of-three introsort that
recurses on the smaller partition and loops on the larger, so the stack stays
`O(log n)`, with a heapsort fallback past a depth of `2*floor(log2(n))` so an
adversarial input cannot drive it to `O(n^2)`, and insertion sort for short
runs. Sorting the same five-million-element array takes 0.375s against `qsort`'s
1.01s and allocates nothing. The third change is the group-label array on the
`x`/`g` path, which was sized at one pointer per *observation* to hold one per
*group* — 40 MB at `n = 5e6` to hold three pointers — and now grows on demand.
- At `n = 5,000,000` over three groups, measured as `VmHWM` either side of the
call:
- | | 0.3 | 0.301 |
|---|---|---|
| peak memory | 228 MB (47.8 B/obs) | 76 MB (15.9 B/obs) |
| `kruskal_test(\@x, \@g)` | 1.20 s | 0.557 s |
| `kruskal_test(\%h)` | 1.11 s | 0.467 s |
- Sorted, reversed, all-equal, organ-pipe and median-of-three-killer inputs all
stay under 0.32s at `n = 2e6`, which is what the depth limit is there for. The
sort is checked against an independent pure-Perl implementation of the whole
test over 748 structured cases — those shapes at every n either side of the
insertion-sort threshold — and 49,712 random ones.
- The same input now gives the same answer:
- `H` moved by up to `1.2e-14` between runs on identical data. Nothing was random:
the sum of `R_i^2 / n_i` walked the groups by group id, and on the
hash-of-arrays path an id is minted in `hv_iternext` order, which is perl's
per-process hash order. Equal values were also left in whatever relative order
the sort happened to leave them, which came from the same place.
- The sort now orders by value and then by group, which makes it a total order,
and the `k` terms of the sum are ordered before they are added — smallest first,
which is the better-conditioned direction as well as a canonical one. `k` is the
number of groups, not the number of observations, so it costs nothing next to
the ranking. `H` is now bit-identical to R on all 37 corpus cases in all four
call forms, and stays so across 60 runs under `PERL_PERTURB_KEYS=1`.
[Documentation]
- `kruskal_test` gains two sections: what happens to non-numeric, undefined,
`NaN` and infinite elements and to a group left with no data, and what the
returned fields are — `statistic`, `parameter`, `method` and the p-value under
both `p.value` and `p.value` from R's `htest`, plus the `size` and `mean`
sub-hashes of `group.stats`, which are computed over the same observations the
statistic used.
0.3 2026-08-16 CDT
[shapiro_test]
- `shapiro_test` rebuilt against R 4.6.1's `src/library/stats/src/swilk.c` — AS R94,
Royston (1995) — driven by R's and SciPy's own test suites rather than by cases
invented here. Four bugs are fixed, one of them a case R's regression suite
tests for by name, and the statistic is now more accurate than R's own on a
sample whose values dwarf its spread.
- Everything below is checked in the new `t/shapiro_test.R.scipy.t` (146 tests),
whose expected values are frozen literals with their provenance recorded in the
file header; it needs no R and no Python to run. The generators that produced
them, `t/shapiro_test.R.scipy.R` and `t/shapiro_test.R.scipy.py`, are committed
beside it. The full suite is 124 files and 25,081 tests, and
`./test.all.perls.pl` passes on all five local perls — `5.10.1`, `5.12.5`
(long double), `5.42.3`, `5.44.0` and `5.44.0-quadmath` — with no warnings on
any of them.
- The p-value could come back negative:
- R's `tests/reg-tests-1b.R` contains exactly one `shapiro.test` assertion, and it
is this:
stopifnot(shapiro.test(c(0,0,1))$p.value >= 0)
- `shapiro_test([0,0,1])` returned `-4.6648135328131477e-15`. At `n = 3` the
p-value is `6/pi * (asin(sqrt(W)) - asin(sqrt(3/4)))` and `W` has an exact floor
of `3/4` that `c(0,0,1)` sits on, so the subtraction lands on zero from
whichever side the constants round to; R clamps the result at 0 and this module
did not. It is the same defect SciPy fixed as gh-18322. The clamp is in, and
`asin(sqrt(3/4)) = pi/3` is now carried to NV width rather than R's 15 digits,
so the same case comes out at `+4.2e-16` before clamping instead of below zero.
- W and the p-value were only good to nine digits:
- The expected normal order statistics that AS R94 weights the sample with came
from `inverse_normal_cdf()`, which is Moro's approximation and good to about
`1e-9`. They go straight into `W`, so nine digits there is nine digits in the
answer — where R reports sixteen. Against the values SciPy pins in
`TestShapiro`, every one of them annotated upstream as *"reference values
generated using R shapiro.test"*:
- | SciPy case | W was | W is | R 4.6.1 |
|---|---|---|---|
| `test_basic` x1 | `0.900472879324135` | `0.900472879317561` | `0.90047287931756` |
| `test_basic` x2 | `0.959026945965277` | `0.959026946032345` | `0.95902694603234` |
| `test_basic2` x4 | `0.834666275331324` | `0.834666275318169` | `0.83466627531817` |
- and the p-values with them:
- | SciPy case | p was | p is | R 4.6.1 |
|---|---|---|---|
| `test_basic` x1 | `0.0420895752342124` | `0.0420895752222577` | `0.04208957522226` |
| `test_basic` x2 | `0.524597929157127` | `0.524597930470668` | `0.5245979304707` |
| `test_basic2` x4 | `0.000913490482316994` | `0.000913490481812984` | `0.000913490481813` |
- Moro's value is still the starting point, but Newton against the `erfc`-based
normal CDF finishes it. The loop stops as soon as another pass could not move
the answer, so a `double` build pays for one refinement and only the wider NVs
pay for a second. Across a 180-sample sweep over normal, uniform, exponential,
log-normal, Cauchy, tied, tiny-scale and grid data at every n from 3 to 5000,
the worst remaining disagreement with R is `2.0e-15` in `W` and `2.6e-12` in the
p-value — the latter is not sloppier arithmetic but the same last bit amplified,
since the p-value is a function of `log(1 - W)` over a sigma of about `0.6`.
- 1 - W was formed by subtracting from 1:
- The p-value depends on `log(1 - W)`, and `W` runs to within `1e-5` of 1 on a
large normal sample, so computing `W = b^2/ssq` and then `1 - W` throws away
exactly the digits the p-value is made of. R does not do this — its `swilk.c`
forms `w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx)` directly and says so
in a comment — and now neither does this module.
- More accurate than R when the values dwarf their own spread:
- R's `swilk.c` divides the sample by its range but never centres it, so `1e9 +
noise` loses most of its significant digits before `W` is ever formed. SciPy
filed the same complaint from the other end as gh-14462 and works around it by
subtracting the median; this module now does that too, which costs one
subtraction per value.
- Measured against a 60-digit `mpmath` evaluation of AS R94 on the identical
doubles, over `1e6 + noise` and `1e9 + noise` at every n from 3 to 5000:
- | | worst relative error in W | worst in the p-value |
|---|---|---|
| R 4.6.1 | `1.9e-8` | `2.6e-7` |
| `shapiro_test` | `1.0e-15` | `1.4e-13` |
- On well-conditioned samples the two still agree to the last few ulp, so this is
a divergence only where R has already lost the digits. `t/shapiro_test.R.scipy.t`
asserts the invariance rather than R's number there, and records why at that
section.
- Faster as well:
- The sort now goes through the module's own introsort rather than `qsort()`,
whose comparator the compiler cannot inline; the order statistics are generated
for half the sample and mirrored, since the weights are antisymmetric; and ten
`pow()` calls became Horner evaluations. `pow()` is a `__float128` call on a
quadmath perl.
- | n | was | is |
|---|---|---|
| 10 | 0.83 µs | 0.78 µs |
| 100 | 4.78 µs | 5.05 µs |
| 1000 | 79.3 µs | 49.8 µs |
| 5000 | 578 µs | 459 µs |
- `n = 100` is the one size that got slower: at that length the accurate
quantiles are most of the work and there is not enough sorting to pay for them.
That trade was taken deliberately.
- One documented value was wrong:
- The hash printed under `shapiro_test` in this README, in `read.me.pod` and in
the module's own POD showed `statistic 0.960870680168535` and `p.value
0.589650577093106` — the pre-fix numbers, and for `[1..19]` while the example
above them calls `shapiro_test([1..5])`. It now shows what `[1..5]` actually
returns, `0.986762155447719` and `0.96717393596804`, which is R's
`shapiro.test(1:5)` to the last digit R prints.
[quantile]
- Interpolation ran between order statistics that were equal:
- R's type 7 interpolates only when the index falls strictly between two order
statistics *that differ* — `index > lo & x[hi] != qs` in `quantile.default`.
This module always evaluated `(1 - g) * x[j] + g * x[j+1]`, which does not
return `v` when both sides are `v`. On a two-valued sample at `n = 999` it
reported `0.99999999999994` for `1`, and on the 602 identical values R's
PR#16672 was filed about it failed to return that value at every prob — which
is the monotonicity failure the PR is about. Both now match R exactly.
- Probabilities a hair outside [0, 1] were refused:
- A probability arrived at by arithmetic rather than written down can land just
outside the interval. R allows `100 * .Machine$double.eps` of overshoot and
clamps to the endpoint — its PR#17891, `quantile(0:1, 1+1e-14) == 1` — where
this module raised an error. It now clamps within the same allowance and still
errors on anything further out. R's constant is used rather than `NV_EPSILON`
on purpose: it is part of what the function *accepts*, so a long-double or
`__float128` build must not reject a `probs` vector R takes.
- Faster:
- The sort was `qsort()` with a function-pointer comparator; it is now the same
introsort `shapiro_test` uses. Ordering 5000 NVs costs about 61,000
comparisons, and paying for an indirect call on every one of them is most of
what a sort of that size costs.
- | n | was | is |
|---|---|---|
| 100 | 3.2 µs | 2.7 µs |
| 1000 | 52.5 µs | 22.8 µs |
| 10,000 | 1.07 ms | 0.72 ms |
| 100,000 | 13.7 ms | 8.8 ms |
- Both fixes and the sort are covered by the new `t/quantile.R.t` (197 tests, or
205 under `EXTENDED_TESTING`), built on the two assertions R's own suite makes
about `quantile` — that `quantile(x, ((1:n)-1)/(n-1))` recovers `sort(x)`, and
that it equals the type-7 interpolation computed by hand off the sorted sample
— run over seven input shapes chosen to break a quicksort (sorted, reversed,
organ pipe, two-valued, tie ladder, sawtooth) at every n either side of the
insertion-sort threshold and the recursion depth limit, plus PR#16672 and
PR#17891 verbatim and 79 frozen R value tables. Its generator,
`t/quantile.R.R`, is committed beside it.
[Documentation]
- Illustrations for three more functions, drawn by `t.test.plots.pl` and
`skew.kurtosis.plots.pl`, both committed:
- **`t_test`** gains six: what the estimate, the standard error and the null
distribution are and which area of it the p-value is; how `conf.int` is the
estimate plus or minus a t quantile and how `conf.level` sets that quantile;
the three `alternative`s side by side with the region each counts and the
interval that goes with it; `p.value` as a function of `mu`, crossing
`1 - conf.level` exactly at the two bounds of `conf.int`; paired,
`var_equal` and Welch on the same data, with the Welch degrees of freedom as
the two spreads separate; and two distributions separating with the interval
retreating from `mu` as the p-value falls.
- **`skew`** gains a left-tailed, a symmetric and a right-tailed sample against
the same `N(0, 1)` curve, with the mean and median of each, which is what the
sign of the statistic is reporting.
- **`kurtosis`** gains a flat-shouldered, a normal and a heavy-tailed sample
with the tails behind each drawn out, since it is the tails and not the peak
that the statistic is measuring.
0.298 2026-08-12 CDT
[wilcox_test]
- A rewrite of `wilcox_test` against R 4.6.1, driven by R's and SciPy's own test
suites rather than by cases invented here. It brings the function up to the
exact conditional inference R gained in 4.6.0, fixes six bugs — two of which
returned confidently wrong p-values on the *default* code path — and adds the
Hodges-Lehmann estimate and confidence interval, `digits.rank`, and the
Edgeworth series.
- Everything below is checked in the new `t/wilcox_test.R.scipy.t` (3,242 tests),
whose expected values are frozen literals with their provenance recorded in the
file header; it needs no R and no Python to run. The full suite is 120 files and
23,149 tests, and `./test.all.perls.pl` passes on all five local perls —
`5.10.1`, `5.12.5` (long double), `5.42.3`, `5.44.0` and `5.44.0-quadmath` —
with no warnings on any of them.
- Exact p-values are now computed when there are ties:
- R 4.6.0 added exact (conditional) inference in the presence of ties, via Torsten
Hothorn's implementation of the Streitberg-Röhmel shift algorithm; R's
`doc/NEWS.Rd` announces it and `tests/reg-tests-1d.R` records the consequence at
its degenerate one-sample cases: *"For R >= 4.6.0 warnings for exact with ties
are gone."* Before that, ties ruled out an exact p-value and both R and this
module fell back to the normal approximation with a warning.
- `wilcox_test` now does what R does. When ties are present the null distribution
is the conditional one given the observed ranks, and the same holds for zero
differences in the signed-rank test. The warnings are gone with them.
- This changes published answers on tied data, including R's own documented
examples:
- | case | was | is (R 4.6.1) |
|---|---|---|
| `?wilcox.test` man-page data, `wilcox_test(\@x, \@y)` | `0.13291945818531886` | `0.12990538872891813` |
| the `airquality` Ozone example (`W = 127.5`) | `1.2080783e-04` | `6.1087351888e-05` |
| `wilcox_test([1,2,2,3], [4,5,5,6], exact => 1)` | `0.02842953599879653` + a warning | `0.028571428571428571` |
| `wilcox_test([1,1])` | `0.34577858615116` | `0.5` |
| `wilcox_test([4,3,2], [3,2,1], paired => 1)` | `0.14891467317876567` | `0.25` |
- Two further consequences are worth knowing about. V itself changes when zero
differences are present, because the exact test ranks `|x - mu|` over every
observation and only afterwards drops the ranks belonging to the zeroes, where
the approximation drops the zeroes first and ranks what is left:
`wilcox_test([-1, 0, 1])` gives `V = 2.5` exactly and `V = 1.5` with
`exact => 0`. R's two branches differ in exactly the same way. And degenerate
inputs that used to be fatal now return a result, as they must for
`tests/reg-tests-1d.R` line 332 to pass: `wilcox_test([0])` gives `V = 0`,
`p = 1`, and so does `wilcox_test([0,0,0,0,0])`, which SciPy pins as
`test_all_zeros_exact`.
- If you need the old numbers, `exact => 0` still asks for the approximation and
is unchanged.
- The exact upper tail was returning zero, on the default path:
- `p_greater` was computed as `1 - CDF(q - 1)`. That subtraction cancels away every
significant digit once the true p falls below `NV_EPSILON`, and then returns a
flat `0`. It did not take a contrived input to reach: two perfectly separated
samples of 30 apiece are inside the automatic exact branch, no `exact => 1`
required.
- | m = n | was | is | R 4.6.1 |
|---|---|---|---|
| 20 | `7.2544192875e-12` | `7.2544445519e-12` | `7.2544445519e-12` |
| 25 | `7.8825834748e-15` | `7.9107286024e-15` | `7.9107286024e-15` |
| 30 | `0` | `8.4556169461e-18` | `8.4556169461e-18` |
| 49 | `0` | `3.9250145965e-29` | `3.9250145965e-29` |
- Both tails are now summed directly. That alone is not enough for the rank-sum
table, whose Gaussian-binomial recurrence is built with subtractions, so far up
the support a count of `1` is the difference of numbers around `C(m+n, n)` and
has already been rounded into noise. The table is folded about its centre before
summing, so only well-conditioned entries are ever touched — the same thing R's
`pwilcox()` does when it folds `q` about `m*n/2` and flips `lower_tail`.
- The signed-rank tail was accurate to `n = 49` by luck (`1 - 2^-49` is exactly
representable) and reached `0` from about `n = 53`; forcing
`exact => 1` on `n = 120` returned `0` where R gives `1.5046327690525337e-36`,
and now returns it too.
- `int m * n` overflowed, and said the samples were identical:
- `exact_pwilcox` took `int m, int n` and computed `int max_u = m * n`. For two
separated samples of 50,000 that wraps negative, every statistic looks out of
range, and the function returns `1.0`:
wilcox_test([1 .. 50000], [50001 .. 100000], exact => 1); # p = 1
- Signed overflow is also undefined behaviour, so a different optimiser was
entitled to do something else entirely. Sizes and indices in the exact
distributions are `size_t` now, the multiplications are checked for wrap before
they happen, and a table that would need more than 16 million cells is refused
outright with a message naming `exact => 0` rather than attempted.
- NaN was ranked instead of dropped:
- `NaN` is `NA` to R, and R drops it. `looks_like_number` accepts it, `d == 0.0`
is false for it, so it went into the rank buffer — and `cmp_nv3` returns `0` for
every comparison involving it, which leaves `qsort` without the strict weak
ordering the C standard entitles it to.
- The visible symptom is R's own regression case, `tests/reg-tests-1d.R` line
3546, which asserts that a paired test is unaffected by pairs whose difference
is `Inf - Inf`:
- | | was | is (and R) |
|---|---|---|
| `1:5` vs `4*(0:4)` | `V = 1`, `p = 0.125` | `V = 1`, `p = 0.125` |
| the same with `+Inf` appended to both | `V = 1`, `p = 0.0625` | `V = 1`, `p = 0.125` |
| the same with `-Inf` and `+Inf` on both | `V = 2`, `p = 0.046875` | `V = 1`, `p = 0.125` |
- `NaN` — in either sample, and however it arises — is now dropped with the other
missing values. `±Inf` is not missing and is kept, since a rank test has no
trouble with it; SciPy's `test_gh_11355b` pins five cases of that and they all
agree.
- An empty `y` ran a different test:
- `wilcox_test([1,2,3], [])` fell through to the one-sample branch and returned a
signed-rank result, silently answering a question nobody asked. It croaks now,
with R's message.
- `mu` was likewise unvalidated: `mu => Inf` or `mu => NaN` turned every
difference into a non-number and produced a confident answer from the wreckage.
Both croak now, as they do in R.
- A dying `$SIG{__WARN__}` handler leaked the rank buffer:
- The warnings in `wilcox_test` were emitted while the `RankInfo` and difference
buffers were held as raw pointers. A `__WARN__` handler that dies — or `warnings
FATAL` at the call site — longjmps straight past the `Safefree`. Under valgrind,
500 iterations of the ties path with such a handler lost 95,616 bytes in 498
blocks. Every allocation now goes through `Newx` plus `SAVEFREEPV`, the idiom
`chisq_test` in the same file already used, so it is released by the save stack
however the call unwinds. The same 500 iterations now report `definitely lost: 0
bytes`, as does a sweep over every croak path and every branch of the function.
- New: `conf.int`, and a Hodges-Lehmann estimate:
- R has returned a distribution-free confidence interval and a point estimate
since PR#1150 in 2001, and `tests/reg-tests-1a.R` has guarded them ever since
with Hollander & Wolfe's published numbers. `wilcox_test` now computes both, by
all four of R's routes — the exact interval from the order statistics of the
Walsh averages or the pairwise differences, the exact interval conditional on
the observed ranks when there are ties, and the asymptotic interval from a root
search:
my $r = wilcox_test(\@y, \@x, paired => 1, conf_int => 1);
# $r->{estimate} == -0.46
# $r->{conf.int} == [-0.786, -0.010]
# $r->{conf.level} == 0.9609375
- Those are Hollander & Wolfe (1999) 2nd ed., pp. 40 and 53, to the digit. So are
the two-sample values from pp. 111 and 126: estimate `-0.305`, interval
`(-0.76, 0.15)`.
- The level a rank test can actually deliver is a step function of the data, so
`conf.level` reports what was achieved rather than echoing what was asked for —
`0.9609375` above, not `0.95`. `conf.level`, `tol.root` and R's alpha-doubling
search for a level the data can support (with its *requested conf.level not
achievable* warning) all behave as R's do.
- New: `digits.rank`, `edgeworth`, and more of R's result fields:
- `digits.rank` rounds each value to a given number of significant digits before
ranking, so that ties are decided on the rounded values. R's man page recommends
it because tie detection is an exact `==` on floating point, and its own worked
example shows `(4:2)/10` against `(3:1)/10` — three differences that ought to be
`0.1` and are three different doubles. Ported from R's `fprec()`, half-to-even
rounding included.
- `edgeworth => 1, 2, 3` adds up to three Edgeworth correction terms to the normal
approximation, the refinement R 4.6.0 reaches through its integer `correct`. It
is ignored on the exact path, and — as in R — ignored when there are ties, or
when the signed-rank test dropped a zero, because the series is derived for
untied ranks.
- The result hash gains `statistic.name` (`"W"` or `"V"`, as R prints), plus
`null.value` and `null.value.name`, and `estimate` / `conf.int` / `conf.level`
when an interval was asked for.
- Three deliberate differences from R:
- Each is asserted in the test file, so that changing one later is a choice rather
than a drift.
- 1. `correct` is a boolean here. R 4.6.0 turned its `correct` into an integer
`0:3`, in which numeric `0` still applies the continuity correction and only
`FALSE` removes it — so in R, `correct = 0` and `correct = FALSE` are
different tests. Keeping that would mean `correct => 0` no longer meaning
"off", which is what it means for every other flag in this module. `correct`
stays a boolean, and R's `correct = k` is `correct => 1, edgeworth => k`.
2. A zero variance is reported, not propagated. With `exact => 0` and every
observation tied there is nothing to divide by. R divides anyway and returns
`NaN`; this warns and returns `p = 1`. The default path no longer reaches it
at all, since the exact test handles all-tied data.
3. An all-tied interval does not raise. R's one-sample code warns and hands
back a `NaN` interval at level `0`; its two-sample code warns and then dies
inside `uniroot` with *missing value where TRUE/FALSE needed*. We give the
one-sample answer in both places.
- There is one place where this module is simply more accurate than R. R's exact
p-values on tied data come from a density it normalises entry by entry;
`wilcox_test` sums the integer permutation counts and divides once. For the
worst case in the corpus — an 11-against-12 tied rank sum whose p-value is
exactly `4/676039` — this returns the correctly rounded double and R is
`1.2e-11` high. Checked against exact rational arithmetic, and recorded in the
test file rather than papered over.
- Testing:
- `t/wilcox_test.R.scipy.t` takes its cases from the references' own suites:
- R's `tests/reg-tests-1a.R` (the PR#1150 Hollander & Wolfe intervals),
`reg-tests-1b.R` (the Wolfgang Huber `wilcox.test(1, 2:60)` case, and the
check that the asymptotic estimate does not move with `alternative`),
`reg-tests-1d.R` (the six degenerate one-sample calls and the `±Inf`
identities), and the man-page examples whose printed output is pinned in
`tests/Examples/stats-Ex.Rout.save`.
- SciPy 1.17.1's `TestMannWhitneyU`, whose header reads *"All magic numbers are
from R wilcox.test"* — `cases_basic`, `cases_continuity`, `cases_9184`,
`cases_2118`, `test_tie_correct`, `test_exact_U_equals_mean`,
`test_gh_11355b` and the 30-against-20 asymptotic cases — and
`TestWilcoxon`'s `test_accuracy_wilcoxon`, `test_wilcoxon_tie`,
`test_onesided`, `test_exact_pval`, `test_exact_p_1`, `test_all_zeros_exact`
and `test_symmetry_gh19872_gh20752`.
- A 663-case sweep generated by `t/wilcox_test.R.scipy.R`, committed next to the
test, crossing four data shapes against every alternative, `exact` state,
`correct` state, `mu` and `conf.int` setting.
- Beyond the file, 960 further randomised calls were compared against R 4.6.1 and
agree everywhere except the three divergences above.
- One lesson from getting that to pass on every NV width is worth recording: the
corpus data has to be exactly representable. Whether two values tie decides
which branch runs, and `1.6 - 2 - 0.5` does not land on the same value in a
`double`, an x87 `long double` and a `__float128`. A corpus of one-decimal
values passed on the default perl and failed on `perl-5.12.5` and quadmath with
a *different statistic*, not merely a different last digit. Every generated
value is now a whole number of quarters or of 1024ths. For the same reason the
asymptotic interval, which is only ever pinned down to `tol.root`, is generated
at `tol.root = 1e-12` rather than freezing wherever Brent's method happened to
stop on one machine.
- A compiler-warning audit of `LikeR.xs` for `-Wint-conversion`, `-Wimplicit-int`,
`-Wreturn-mismatch` and `-Wdeclaration-missing-parameter-type`, and a pass
tightening integer types that can only hold a count or a flag. No behaviour
changed: the full suite (116 files, 18,546 tests) passes, and every function
touched was diffed call-for-call against a build of the previous release, with
`dnorm`, `pnorm`, one- and two-sample `ks_test`, `fisher_test`, `auc` and the set
operations re-checked against R 4.6.1 and found bit-identical.
- All four of those warnings were already clean, and stay clean on a `double`, a
`long double` and a `__float128` build. Two of them cannot be tested with the
GCC most systems still default to: `-Wreturn-mismatch` and
`-Wdeclaration-missing-parameter-type` are GCC 14 additions — where they are
errors rather than warnings — and GCC 13 rejects both as unrecognized options,
so a check that appears to pass on 13 has really only skipped them.
[Dead code removed]
- Turning the audit up to `-Wextra` found two branches that could never run, both
of them a test for negativity on a value whose type is unsigned:
- 1. `r_pow_di` takes `unsigned int n`, so its `if (n < 0) return 1.0 / r_pow_di(x, -n);`
was unreachable — a leftover of R's `R_pow_di`, which takes a signed exponent.
All three callers (in `K2x`, for the exact one-sample Kolmogorov-Smirnov
distribution) pass a non-negative exponent, so the unsigned parameter is the
correct one and the reciprocal branch simply goes.
2. `hoa2aoh` casts `HvUSEDKEYS` to `U32` and then clamps with `if (ncols < 0) ncols = 0;`.
[Types narrowed to what they can actually hold]
- Eighteen `int`s that only ever hold 0 or 1 became `bool`, a convention the file
already followed in some 219 other places; each was confirmed by reading every
call site rather than by name. The flag parameters of `ft_pnhyper`, `K2l`,
`c_dnorm`, `c_pnorm`, `c_pnorm_both`, `set_multiplicity` and `roc_split`, the
`is_cat` field of `AnFac`, the `lower_pos` and `frac_low` locals of `auc`,
`auroc`, `roc` and `bedroc`, and the return types of `mg_key` and
`psmirnov_exact_test`. Several of these were already being handed a `bool` by
their callers — `dnorm`'s and `pnorm`'s `log` and `lower` options, for instance —
so only the helper signatures were behind. `c_pnorm_both`'s loop counter became
`unsigned int`.
- Two that look like flags and are not: `c_pnorm_both`'s `i_tail` is three-valued,
and `set_multiplicity`'s `gimme` carries a Perl `G_*` context value. Both stay
`int`.
- Also six coefficient tables in `c_pnorm_both` written `const static double`,
which puts the storage class after the qualifier and draws
`-Wold-style-declaration`; they are now `static const double`.
[runif argument validation, and every warning names its function]
- `runif` accepts its arguments either positionally or by name, and decided which
was which by asking whether the current argument was a string *and* whether
another argument followed it. A key at the end of the list therefore failed the
second half of that test and fell through to the positional branch, where it was
read as a number: `runif(5, 'min')` took `SvNV("min")`, which is 0, silently set
`min = 0`, and returned five values. The only sign anything was wrong was perl's
own `Argument "min" isn't numeric`, which does not say which function provoked
it. `runif(5, bogus => 1)` went the same way, taking `bogus` as `min` and `1` as
`max`. Every sibling that parses named arguments — `rbinom`, `binom_test`,
`fisher_test`, `dnorm`, `pnorm` — rejects both of those.
- `runif` now does too. A string argument is treated as a key when it is not a
number, which is decidable from the key alone, so a dangling or misspelled key
is an error instead of a silent coercion; a numeric string is still positional,
so `runif("9")` is unchanged. Named values are checked for numerichood before
use, which is what keeps perl's unattributed warning from being the diagnostic.
- `n` is also range-checked now. It was read straight through `SvUV()`, so
`runif(-1)` wrapped to 2**64-1, `av_extend()` read that back as a negative
`SSize_t`, and perl died with `panic: av_extend_guts() negative count (-2)` --
which names neither the function nor the argument at fault. A negative or
over-large `n` now croaks and says so. Non-integer `n` still truncates toward
zero, as R's `runif()` does, and `runif(0)` still returns an empty list.
- Separately, three warnings did not name the function emitting them, unlike every
other warning in the file: one in `ks_test` (the 1-sided exact 1-sample case
falling back to asymptotic) and two in `wilcox_test`'s signed-rank branch (exact
p-value abandoned for ties, and for zeroes). All three now carry the prefix
their siblings already had. The one warning left deliberately bare is the
`warn("%s", m)` in the uninitialized-value catcher, which re-emits somebody
else's warning verbatim and must not add to it.
[Argument-stack indices are now Stack_off_t]
- `-Wextra` reported 58 `-Wsign-compare` warnings, and 33 of them were one idiom:
an index declared `size_t`, `unsigned`, `unsigned int` or `unsigned short int`
and then compared against `items`. `items` is neither of those — XSUB.h's
`dITEMS` declares it `Stack_off_t items = (Stack_off_t)(SP - MARK)`, a *signed*
type, because it is a stack-pointer difference. Every one of those comparisons
was converting the signed side to unsigned.
- The indices are now `Stack_off_t` themselves, which is the type they are
compared against: 25 declarations across 23 functions — `binom_test`,
`ks_test`, `wilcox_test`, `write_table`, `max`, `runif`, `quantile`, `mean`,
`mode`, `sum`, `sd`, `uniq`, `var`, `t_test`, `median`, `matrix`, `fisher_test`,
`power_t_test`, `var_test`, `dnorm`, `value_counts`, `prcomp` and `pnorm`.
That is a retype, not a cast: writing
`(size_t)items` at each comparison would silence the warning just as well, but
it would be wrong the day `Stack_off_t` widens, which is exactly what it exists
to allow. `t_test`'s index was `unsigned short int`, which drew no warning at
all — integer promotion made the comparison signed — and was the same latent
mistake regardless.
- `Stack_off_t` arrived in perl 5.39.2 and this distribution supports 5.010, so
the preamble now carries a shim typedef guarded on `PERL_STACK_OFFSET_DEFINED`,
the macro perl.h defines next to the typedef. On 5.10.1 and 5.12.5 neither the
macro nor the type exists and the shim supplies `I32`, which is what the stack
offset was on every perl before that.
- The 33 warnings are gone, 25 remain, and no warning category increased —
verified by compiling the before and after trees and diffing the warning sets.
The remaining 25 are unrelated signedness pairs (`size_t` against `ssize_t`,
`IV` against `size_t`, `STRLEN` against `ssize_t`) and are left alone. The full
suite passes on perl 5.10.1 and 5.12.5, the two builds that depend on the shim,
as well as on 5.42.3, 5.44.0 and 5.44.0-quadmath; and 94 calls covering all 23
retyped functions — positional and named forms, bare lists against arrayrefs,
`write_table`'s emitted bytes, and the odd-argument and unknown-argument croaks
that this index arithmetic drives — produce identical output before and after.
[NV was being computed at double precision on wide builds]
- Every libm call in `LikeR.xs` was written bare — `sqrt(x)`, `log(x)`,
`lgamma(x)` — and C has no type-generic `<math.h>`. Those functions take a
`double`, so on a perl built with `-Duselongdouble` or `-Dusequadmath` every one
of them converted the `NV` down to 53 bits of mantissa, computed there, and
converted the result back. Nothing warned and nothing failed to compile; the
answers were simply less accurate than the perl running them. On perl-5.12.5
(`long double`), `sd(1..5)` returned exactly the double-rounded `sqrt(2.5)`,
9.5e-17 away from the value perl's own `sqrt` gives.
- All 412 of those calls now go through `nv_*` macros that paste on the suffix for
the width `NV` actually is: none for `double`, `l` for `long double`, `q` for
`__float128`. The 80 `isnan`/`isinf`/`isfinite` calls became `nv_isnan`,
`nv_isinf` and `nv_isfinite`, which classify by comparing against `NV_MAX` — the
largest finite `NV` — rather than calling libm at all. The C99 macros could not
be kept: where a platform does not provide the type-generic versions,
`isfinite()` is a plain `double` function, and narrowing a large-but-finite long
double into it reports the value as infinite rather than merely rounding it.
Perl's own `Perl_isnan`/`Perl_isinf`/`Perl_isfinite` were used up to 0.298 and
could not be kept either: on every perl before 5.22 those route through a
`Perl_fp_class()` block in `perl.h` that has never compiled — the macro is
written with an empty parameter list and compares against `FP_CLASS_*` names no
`<ieeefp.h>` defines. That block is dead code wherever Configure finds
`isinf()`, so it is invisible on Linux and glibc, and live on illumos/Solaris,
where it broke the 0.298 build outright.
- The long-double row is conditional. The `l` variants are C99 but some libms —
the thinner BSD ones especially — do not ship the whole set, so `Makefile.PL`
link-tests all twenty as a unit and defines `LIKER_HAVE_LONG_DOUBLE_MATH` only
if every one resolves; otherwise the build falls back to the `double` functions,
which is exactly what it did before and so cannot regress. `__float128` needs no
probe: `<quadmath.h>` and `-lquadmath` come with the quadmath perl itself, and
the built object was checked with `nm` — it references `lgammaq`, `expq`,
`sqrtq` and no double-width libm symbol at all.
- Accuracy on the long-double build, measured against values that are exact in
binary or known in closed form: `sd(1..5)` is now bit-identical to perl's
`sqrt(2.5)`, and `fisher_test([[3,1],[1,3]])` moves from 1.5e-16 to 6.4e-18
relative error against the exact 17/35. The remaining 6.4e-18 is an accuracy
floor in that function's own summation, not a width problem — the `__float128`
build lands on the same figure.
- This costs time where the wide math is software-emulated: the suite takes 352s
on the quadmath perl, against 67s when it was quietly running on hardware
doubles. The other four perls are unaffected.
[The build ran itself twice, and clobbered its own Makefile doing it]
- `make` had to be run twice or the `.so` came out stamped with the wrong version
and refused to load. The cause: ExtUtils::MakeMaker scans the directory for
`*.PL` files to run during the build, and `dev.Makefile.PL` — a local
convenience wrapper, not part of the distribution — looks like one. It was being
run mid-build as `perl dev.Makefile.PL dev.Makefile`, and since it calls
`WriteMakefile()` it overwrote the real `Makefile` with its own: no `DEFINE`, no
probed C99 flag, and a different `VERSION`. The second `make` then rebuilt from
that. `PL_FILES => {}` turns the scan off; nothing here is generated by a `.PL`
file.
- The version half was a stale literal: the checked-in `Makefile.PL` pinned
`VERSION => "0.28"` while `lib/Stats/LikeR.pm` had moved to 0.298, and
`XSLoader::load()` passes `$VERSION` to a `.so` compiled with `-DXS_VERSION`
from that literal. It now reads `VERSION_FROM => lib/Stats/LikeR.pm`. One `make`
after `perl Makefile.PL` is enough again, and the non-quadmath builds are about
a third faster for not doing the work twice.
[Portability: Solaris, the BSDs, and vendor compilers]
- The C99 flag is now probed instead of guessed. `Makefile.PL` was selecting
`-std=gnu99` on any compiler whose name matched `/\b(?:g?cc|clang)\b/`, and
`$Config{cc}` is plain `cc` for Oracle Studio on Solaris and for aCC on HP-UX —
both of which reject that flag outright, so the build failed there before it
compiled a line. Each candidate is now trial-compiled and the first that works
wins: `-std=gnu99`/`-std=c99` for gcc and clang, `-xc99=all` for Studio,
`-qlanglvl=extc99` for AIX `xlc`, `-AC99` for HP-UX, and nothing at all for a
compiler already in C99 mode. MSVC is skipped outright, since it warns rather
than errors on switches it does not know and would make the probe settle on a
no-op.
- Two things that would have failed to compile off Linux are gone. `<strings.h>`
and its `strcasecmp` — POSIX-only, absent on MSVC — are replaced by a small
`str_ieq_ascii()`, which also drops the locale dependency: `tolower()` under a
Turkish locale maps `I` outside ASCII, which should never decide whether
`"TRUE"` matches `"true"`. And bare C99 `restrict`, used on 151 pointers here,
now has an `#ifdef` mapping it to `__restrict` on MSVC and `__restrict__` on
older gcc, and defining it away where no spelling exists, rather than losing the
annotation.
- `LikeR.xs` also compiles clean under strict `-std=c99` with no GNU extensions,
which is the closest available local proxy for a vendor compiler.
[Dead code: sample()'s private PRNG]
- A splitmix64 generator sat at the top of the file under a comment promising a
PRNG stream separate from `Drand01()`, seeded lazily from `/dev/urandom` with a
`time()^PID` fallback. None of it was true: no seeding code was ever written, no
caller ever existed, and its state started at a fixed 0, so had anything called
it the "random" sample would have been the same sequence in every process.
`sample()` draws from `Drand01()` and always did, which is the behaviour that is
wanted — `srand($seed)` governs it the way `set.seed()` governs R. The generator
and its comment are removed.
[Tests]
- Two files, 273 assertions, and both were checked against a deliberately broken
build rather than merely observed to pass.
- `t/nv_width.t` fails if the math width ever comes undone. Its sharp assertion
needs no tolerance at all: `sd(1..5)` must be the identical NV to perl's
`sqrt(2.5)`, which holds on any width and breaks the moment a `double` gets in
the way. It is width-adaptive rather than skipped on a `double` perl, computing
the NV epsilon of the running build instead of assuming one.
- `t/scale.keywords.t` covers `scale()`'s string options — `"mean"`, `"sd"`,
`"none"`, `"true"`, `"false"`, `""` and their case variants — which had no
coverage at all: `t/01.t` passes only the numeric forms. Expected values come
from R 4.6.1 `base::scale()` at `options(digits=17)` and are frozen in the file,
so it needs no R at run time. Deleting the case fold from `str_ieq_ascii()`
fails 11 of its assertions; usefully, all 11 are the "off" spellings, because an
unmatched string falls through to `SvTRUE` and still means "compute it", so
`"MEAN"` would keep working while `"NONE"` flipped. That is recorded in the file
so the section is not trusted for more than it proves.
- The suite is 118 files and 18,819 tests, passing on perl 5.10.1, 5.12.5, 5.42.3
(threaded), 5.44.0 and 5.44.0-quadmath, with no compiler warnings on any of
them.
0.297 2026-08-10 CDT
- https://www.cpantesters.org/cpan/report/260534ea-9474-11f1-8ca2-bfb68deea6df bug fix
0.296 2026-08-09 CDT
- fixed CPAN bug: https://www.cpantesters.org/cpan/report/fcf32c68-75a5-1014-bc87-8fe0d10910fe
- write_table.announce.t ran its child perl through -e, which cannot carry double quotes or shell metacharacters on Windows; the child program now goes in a file
- chisq_test now matches R 4.6.1 bit-for-bit on the statistic across 170 randomized cross-check cases, and the full suite (116 files, 18,546 tests) passes.
- Bugs found and fixed in LikeR.xs
- 1. A 1×k or k×1 table returned df = 0, p = 1 — no test at all. R collapses a single-row/column matrix to a vector and runs goodness-of-fit (if (min(dim(x)) == 1L) x <- as.vector(x)); now so does this. [[10,20,30]] went from X²=0, df=0, p=1 to X²=10, df=2, p=0.006738.
2. Yates' label was attached even when the correction was zero. R only says "with Yates' continuity correction" when min(0.5, |O−E|) > 0. A table sitting exactly on its expectation, and every zero-margin table, were mislabelled.
3. Yates was computed per cell instead of as R's single whole-table min(0.5, abs(x-E)) — equal in theory on a 2×2, not always in the last bits.
4. No input validation. Negatives, infinities, NaN, strings and undef were silently coerced to 0 and produced garbage or NaN; all-zero data returned NaN; a single element returned df = 0. All now croak with R's wording. Ragged array rows and 2D hash rows with mismatched column keys were silently zero-filled — now fatal.
5. Uniform expectation used n/k instead of R's n * (1/k), and sums were accumulated in a plain NV where R uses a long double. Together these put the statistic 1–2 ulp off R on most inputs; both fixed (ct_acc_t).
6. Hash input was read in Perl's randomized key order, so which row a malformed hash got blamed on was a coin toss. Rows and columns are now sorted, as fisher_test already does.
7. Segfault on sparse arrays (av_fetch returns NULL for a hole) — this one I introduced during the rewrite and caught before finishing; guarded by ct_av_get
- Three of those cross-checks compared the statistic to R's printed value relatively, and on the tables in question R's value is not a statistic. Where a 2×2 has all four |O−E| equal, Yates' min(0.5, |O−E|) cancels every corrected residual, so the exact statistic is 0 and the exact p is 1; what R prints there — 1.4515367733818938e-24 for [[1573,3],[4,0]], 2.9347503914472165e-32 for [[1,2],[3,4]], 7.1842689582627857e-32 for [[1.5,2.5],[3.5,4.5]] — is the leftover of forming E in floating point, the four |O−E| differing in their last bits so that the minimum comes out a hair below the rest. Its size is a property of the NV rather than of the test: a double build reproduces R's digits, and a __float128 build cancels the whole way to 0. Comparing that relatively can only pass on the width R happened to use, and it failed with rel diff = 1 on the quadmath perl and on 5.12.5. Those three cases in t/chisq_test.R.scipy.t now check the statistic against 0 and the p-value against 1 with absolute tolerances of 1e-20 and 1e-11, R's numbers staying in the file as provenance. LikeR.xs is unchanged — the wide-NV answer was the more accurate one. The suite passes on perl 5.10.1, 5.12.5, 5.42.3-thr, 5.44.0 and 5.44.0-quadmath.
0.295 2026-08-08 CDT
- bug fix https://www.cpantesters.org/cpan/report/0f13fed6-92f5-11f1-b043-dc326e8775ea
- Removed `restrict` where it made no difference, or was potentially dangerous
[drop_duplicates, merge, value_counts]
- These three decide what counts as the same row, the same join key, or the same
value by a cell's Perl stringification, and on numeric columns that one
conversion was most of the work they did.
- `sv_2pv_flags()` renders an NV with `snprintf("%.*g", NV_DIG, x)`, about 140 ns
a cell, and — unlike the IV case, where `SvPOK_or_cached_IV` lets the `SvPV`
macros hand back the string perl cached on the SV — it never reuses that PV, so
every pass over a column of doubles paid the conversion again. It does leave the
buffer behind, which is why keying a frame used to grow the caller's own numeric
columns by about 64 bytes a cell, permanently: reading a frame ought to be a
read.
- `nk_num_pv()` now renders bare integers and bare doubles into the caller's own
scratch buffer instead, and leaves the SV untouched. Its double path is `%.15g`
about four times faster than the C library's, and taken only where the answer is
provably the same: the magnitude is scaled into `[1e14, 1e15)` in `long double` —
64 mantissa bits against the double's 53 — which bounds the scaled value's error
under 2e-4, so a fractional part further than 2e-3 from one half rounds exactly
as the true value would. About one cell in 300 lands nearer than that and goes
back through `SvPV`, as do zero, the non-finite values, `use locale`, an x87
control word left at double precision, and any build whose NV is not an IEEE
double. It agreed with the C library's own `%.15g` over 90 million random bit
patterns; `t/drop_duplicates.t`, `t/merge.t` and `t/value_counts.t` now group
tens of thousands of doubles both ways and require the two answers to match.
- Two further changes in `drop_duplicates` alone:
- Its interning table started at 64 slots and doubled, so a pass over 10,000
distinct rows rehashed nine times, each one a scattered walk over a table too
big for L2. The row count is known before the pass starts and bounds the group
count, so it is now used as the hint — capped, so a large frame of few distinct
rows does not pay for a slot per row.
- An HoA result copied every surviving cell, while AoA and AoH already shared the
whole surviving row. It now shares the cells too. **This is a behaviour
change.** The frame, and an HoA's column arrays, are still new, so they can be
reshaped without touching the input; but assigning *through* a survivor —
`$out->{col}[0] = ...` — now writes to the input's cell, exactly as
`$out->[0]{col} = ...` always did for AoA and AoH. Clone the result if you need
full independence.
- Measured on the 10,000-row frame `benchmark.pl` uses (five columns: two doubles,
one integer, two strings), on one machine, with only these paths toggled. Time is
the median of 25 calls in one process; RAM is `benchmark.pl`'s own figure, the
`VmRSS` delta of a forked child running the call once, median of nine. The string
row is there to show where the win is not: it is confined to numeric cells.
- | Call | Time before | Time after | RAM before | RAM after |
|---|---|---|---|---|
| `drop_duplicates($hoa)` | 5.45 ms | 1.88 ms (2.9x) | 5.36 MB | 1.45 MB (3.7x) |
| `merge`, inner join on an integer key | 7.78 ms | 5.62 ms (1.4x) | 7.41 MB | 6.41 MB (1.2x) |
| `merge`, inner join on a double key | 10.28 ms | 5.91 ms (1.7x) | 7.12 MB | 6.42 MB (1.1x) |
| `value_counts` on a double column | 2.98 ms | 1.66 ms (1.8x) | 1.98 MB | 1.55 MB (1.3x) |
| `value_counts` on a string column | 0.215 ms | 0.220 ms | 0.69 MB | 0.71 MB |
- `group_by` and `pivot_table` were left alone: `group_by` hands the cell SV
straight to `hv_fetch_ent`, so perl does the stringification internally and
reaching it means byte-level `hv_*` calls and a change to how UTF-8 keys are
handled, and `pivot_table` is pure Perl.
0.294 2026-08-07 CDT
- bug fixes: https://www.cpantesters.org/cpan/report/368ca238-73ee-1014-a03f-97f1b88bf904
- `binom_test` was cross-validated against R 4.6.1 `stats::binom.test` and SciPy
1.17.1 `scipy.stats.binomtest` using their own test suites rather than cases
invented here: SciPy's `TestBinomTest`, R's `binom.test(c(800,10))` from
`tests/reg-tests-2.R`, the `?binom.test` example, and an R-generated corpus of
383 p-values and 1560 Clopper-Pearson bounds. They are in
`t/binom_test.R.scipy.t`. Two fixes came out of it, both in the incomplete beta
that every tail and confidence bound goes through:
- Its continued fraction stopped after a flat 500 terms, but it needs about
0.25 sqrt(a+b) of them once the shape parameters are large, so it was quietly
cut short at big `n`: `binom_test(10079990, 21000000, p => 0.48)` returned
0.996781946606 where R and SciPy both give 0.9966892187965, i.e. wrong in the
fourth decimal of a printed p-value. The cap now scales with sqrt(a+b), and
the front factor moved off differenced `lgamma` onto the same saddle-point
form `dbinom` already used here. Agreement with R over these cases went from
9.3e-5 to 3.3e-13 relative.
- The Clopper-Pearson bounds are found by bisection, which stopped at an
absolute width of 1e-15, so a bound far below 1 came back with only four
correct digits: `binom_test(1, 1000000000, alternative => 'greater',
conf.level => 0.999)` gave 1.00053299e-12 against R's 1.00050033e-12. The
stopping rule is now relative to where the bracket sits, and such bounds now
hold about 1e-15.
- Both fixes also help `t_test`, `var_test` and `cor_test`, which use the same
function. One limit remains, pinned by the tests rather than left to chance: the
upper bound for a handful of successes in a billion trials still carries about
1e-9 of relative error, because the complement branch of the incomplete beta
cannot resolve a tiny `x` past the spacing of `1-x`.
0.293 2026-08-06 CDT
- Fixed quadmath error https://www.cpantesters.org/cpan/report/83bcd9a2-9123-11f1-aac1-f3cd035a6881
- `fisher_test` was cross-validated against R 4.6.1 `stats::fisher.test` and SciPy
1.17.1 `scipy.stats.fisher_exact` using their own test suites rather than cases
invented here: SciPy's 84-case R-generated corpus
(`scipy/stats/tests/data/fisher_exact_results_from_r.py`, four numbers per case
over two confidence levels and all three alternatives), its
`TestFisherExact`, R's regression suite (`tests/reg-tests-1{a,b,d,e}.R`:
PR#644, PR#1662, PR#4688, PR#10558, PR#18336, PR#17671 and the "exact
fisher.test" entry) and the `?fisher.test` examples. They are in
`t/fisher_test.R.scipy.t`. Three fixes came out of it:
- The 2x2 hypergeometric density was built by differencing `lgamma`, which
costs the back half of a large table's p-value: at a margin of 8.4e7,
`lgamma` is about 1.4e9, where a double's spacing is 2.4e-7, and exponentiating
that turns into a relative error of the same size. SciPy's gh-3014 case came
out right to only seven digits. The density is now assembled from Loader's
saddle-point binomial, which is how R's own `dhyper` avoids this and which
`binom_test` already had in the file; its three terms stay O(1) whatever the
margins are. Worst-case agreement with R over the 84-case corpus went from
2.1e-12 to 5.2e-14 relative, and gh-3014 from 2.2e-07 to 1.5e-16.
- The R x C enumeration charged only its leaves against its safety cap, so a
table wide enough to spend the time in the interior of the tree neither
finished nor stopped: R's PR#4688 table (4x3, N = 16442), whose whole point
upstream is that `fisher.test` must fail rather than return `p = Inf`, ran for
over five minutes here without doing either. Every node is now counted, and
that table is declined in about a second.
- The R x C enumeration now bounds each subtree before walking it. `lgamma(x+1)`
is convex and `a!b! <= (a+b)!`, which together bracket the probability of every
completion of a partial table; when the whole subtree falls inside the tail its
mass is added in closed form (`N'! / (prod R_i! prod C_j!)`, from counting the
remaining observations into rows two ways), and when it falls outside the
subtree is dropped. The margins are also transposed and sorted first, so the
fattest row and column are the ones the enumeration gets for free. R's Job
Satisfaction 4x4 example went from 7.5s to 0.3s and PR#644's 19x2 from 1.0s to
under 0.05s, and the 6x6 table of PR#18336 -- which segfaulted R before 4.2.0
and which R 4.6.1 still declines with `hash key 5e+09 > INT_MAX` -- is now
computable at 0.6322160531, agreeing with R's own 2e6-replicate
`simulate.p.value` fallback to within its sampling error.
- Two behaviours that the two references disagree about are now pinned by tests
rather than left to chance: a table with an empty row or column returns R's
`p = 1` with an odds ratio of 0 and a CI of (0, Inf), not SciPy's NaN odds ratio;
and a table with a single row or column is rejected as R rejects it, rather than
returning SciPy's `p = 1`.
0.292 2026-08-05 CDT
- fixed long-double bug https://www.cpantesters.org/cpan/report/506975f6-906a-11f1-8f30-a201c4f2440e
- `power_t_test` was cross-validated against R 4.6.1 `power.t.test` and against
`scipy.stats.nct` driven by `scipy.optimize.brentq`, over a grid of 288 cases
covering all five solved-for parameters, all three types, both alternatives and
`strict`. Three fixes came out of it:
- The Simpson sum behind the noncentral *t* CDF put a fixed 30000-step grid on
`u = w/(1+w)`, and the chi density it integrates defeats that at both ends.
The density carries `w**(df-1)`, so unless `df` is a whole number some
derivative of it is infinite at `w = 0` and Simpson's error bound does not
hold: two good digits at `df = 1.2` with `sig_level = 1e-4`, five at
`df = 1.2`, nine at `df = 1.8`. Substituting `w = z**m`, with `m` chosen so
that `m*df - 1 >= 3`, restores the bounded derivatives and brings all of
those to machine precision. It also puts the origin's contribution at zero,
which subsumes a separate bug: the sum had been dropping its `u = 0` endpoint
term, worth 7e-7 of absolute power at `df == 1`. `nu` is now also floored the
way R floors it, per sample rather than in total.
- The same density has standard deviation `1/sqrt(2*df)` and so narrows without
bound, while the grid did not. Past `df` of about 1e7 the steps went clean
over the peak: `power_t_test(n => 4e7, delta => 0)` returned 0.138 where the
answer can only be `sig_level/2`, and a large-cohort `n` solved 9% low. Above
`df` of 1e3 the steps now go on `w` across +/- 12 standard deviations of the
mode, with the chi normalisation taken from Stirling's series to keep the peak
height from cancelling away; and above 4e5, where those log terms cancel too
hard for any grid to help, the Abramowitz & Stegun 26.7.10 asymptotic form
takes over -- the same formula, at the same cut-off, that R's `pnt.c` uses.
That is also 25 times quicker than integrating.
- The power was formed as `1 - P(T <= t)`, which loses most of its digits to
cancellation when the power is small. It is now integrated as the upper tail
directly.
- The four inverse solvers were plain bisection stopped at the bracket width,
which capped `n`, `delta`, `sd` and `sig_level` at R's own four or five
significant figures. They now use regula falsi with the Illinois correction
against a relative tolerance, so they match machine-precision `brentq` roots
to ~1e-13 in fewer evaluations than the bisection took. The `tol` default
moved from `1.22e-4` to `1e-12` to match.
- Nothing checked that the bracket held a root, so an unreachable target came
back as a bracket endpoint wearing the requested power: solving for `sd` with
`power => 0.01` returned `delta * 1e7`, and with a negative `delta` returned a
negative standard deviation. Unreachable targets now croak and name the range
searched. `sig_level` and `power` outside `[0, 1]`, an `n` below 2, a negative
`sd`, and an unrecognised `type` or `alternative` are rejected as well --
`type => 'twosample'` used to be read silently as `'two.sample'`.
- New test file `t/power_t_test.R.scipy.t` carries the cross-validated grid.
0.291 2026-08-04 CDT
- POD formatting improvements
[`lm`, `glm`]
- Formula parsing and data reading are now shared between `lm` and `glm` too, so the
two agree on what a formula means and on what a row is called. `lm` had the better
parser and `glm` the better row naming; each now has both.
- `lm` now names rows the way `glm` does — from a `row.names`, `_row`,
`rownames` or `.rownames` column when the data has one, and 1-based integers
otherwise. `lm` previously always used integers, so `fitted.values` and
`residuals` came back keyed `1..n` for data whose rows had names, and did not
match what `glm` or `predict` returned for the same data; the `predict`
documentation already described the shared behaviour. A row-name column is a label
rather than a measurement, so `y ~ .` now excludes it in both.
- Design-matrix construction is now shared between `lm` and `glm`, and decides a
categorical column's coding term by term using R's margin rule: the reference
level is dropped when the term with that column removed is itself in the model.
Three bugs fall out of that, all confirmed against R 4.6.1 and statsmodels
0.14.6.
- Bug fixes:
- Four in `glm`, from the parser it now shares with `lm`. Three of them ended the
same way: a term that names no column evaluates to `NaN` for every row, every row
is dropped as incomplete, and the fit dies with `0 degrees of freedom (too many
NAs or parameters > observations)` — never mentioning the formula.
- **`glm` truncated a formula at 511 characters.** It copied the formula into a
fixed `char[512]`, so a model with enough predictors to overrun that lost the
tail. The buffer now grows with the formula.
- **`glm` did not understand `.`.** It parsed the formula before reading the data,
so there were no column names to expand `.` into and the term stayed a literal
`.`. Formula splitting now happens first and term expansion after the data is
read, so `y ~ .` works in both.
- **`glm` did not understand `+ 0` or a leading `0 +`.** Only `- 1` suppressed the
intercept; the other two spellings R accepts left a term named `0`. All three now
work in both, as do `+ 1` and a leading `1 +`.
- **`glm` read the `-1` inside `I(...)` as intercept suppression.** It searched the
whole right-hand side for the substring, so `y ~ I(x-1)` silently became
`y ~ I(x) - 1`: a different model, fitted without complaint. The scan now steps
over `I(...)`, leaving the term alone. `I()` still supports only `^power`, so
that formula is an error in both rather than a wrong answer in one.
- And the three that fall out of the shared design matrix:
- **A categorical column in a model with no intercept lost a level.** With no
intercept there is no baseline for a reference level to be measured against, so
R codes the factor in full — one column per group, each coefficient that
group's own mean. Both functions dropped the reference level anyway, so
`len ~ supp - 1` fitted `len ~ suppVC - 1`: a model forcing every observation
at the reference level to a fitted value of 0. On R's `ToothGrowth` that meant
a residual sum of squares of 16056 against R's 3247, and an R² of 0.35 against
0.87. Where two categorical main effects appear with no intercept only the
first is coded in full, as in R, since coding both would be rank deficient.
- **An interaction involving a categorical column could not be built.** The
interaction was looked up as a single column literally named `dose:supp`;
finding none, it evaluated to `NaN` for every row, every row was dropped as
incomplete, and the fit died with `0 degrees of freedom (too many NAs or
parameters > observations)`. Interactions now expand to the product of their
components' indicator columns, so `len ~ dose * supp` gives `dose`, `suppVC`
and `dose:suppVC`. `predict` already understood such coefficient names; now
they can be produced.
- **`a*b*c` expanded only its first `*`.** Crossing is associative, so
`y ~ a * b * c` now yields every non-empty subset (`a`, `b`, `c`, `a:b`, `a:c`,
`b:c`, `a:b:c`), ordered by degree as R's `terms()` orders them. Previously the
chunk was split once, producing the unusable terms `b*c` and `a:b*c`, and the
fit died the same way as above. Crossing more than 16 columns now croaks rather
than expanding to 2^n terms.
- **`predict` scored reference-level rows as if the term were absent.** It
registered factor dummies from `levels[1..]` only, on the assumption that a
reference level never has a coefficient — true for a factor coded by contrasts,
but not for one coded in full. Every row at the reference level of a
no-intercept model therefore came back 0.
- **`glm` halved its IRLS step whenever the deviance rose, costing iterations and
accuracy in the standard errors.** R truncates a step only when the deviance
comes out non-finite; a deviance that merely increases is not divergence. The
standard IRLS start puts `mu` at `y + 0.1`, essentially on the data, so the
initial deviance is near zero and the first real step almost always raises it —
on the nine-point poisson fit in `t/glm.t`, from 0.016 to 1.54. That was read as
divergence and the step was halved ten times over, turning R's four iterations
into seven.
- The extra iterations reached the same coefficients, so the symptom appeared
only in the standard errors. They are built from the information matrix of the
*penultimate* iterate — in R because `summary.glm` inverts the QR that
`glm.fit` kept from its last weighted least squares call, and here because the
IRLS sweep leaves that inverse in place — so stopping on a different iteration
than R means reporting a different matrix. Poisson standard errors were 5e-8 to
2e-5 away from R's while the coefficients agreed to twelve digits; they now
agree to about 1e-14. Binomial standard errors were up to 6e-7 out and now
agree to 2e-14, except on a near-separable fit, where the `varmu` floor of
1e-10 (a guard against dividing by an underflowed variance) accounts for the
remaining difference — 1.9e-9 on `am ~ wt * hp`, where three of 32 fitted
probabilities are within 1e-12 of 0 or 1 and R itself warns. Gaussian fits are
unaffected: their weights are all 1, so the matrix is `X'X` either way.
- The same condition had its `isfinite` test on the accepting side, so a
genuinely divergent step producing a non-finite deviance was kept rather than
truncated. That is now the one case that does trigger halving.
- **The negative-binomial theta alternation stopped early and started from the
wrong place.** `MASS::glm.nb` does not simply maximise over theta; it alternates
between an IRLS fit at the current theta and a fresh ML estimate of theta at the
current fitted means, and which fit it lands on depends on the schedule. Four
details of that schedule were wrong here, and all four are now reproduced:
- The alternation stopped on a relative test of the log-likelihood alone,
`|dll| < 1e-7 * (|ll| + 0.1)`. `glm.nb` requires
`(|dLm| / d1 + |dtheta|) < 1e-8` with `d1 = sqrt(2 * max(1, df.residual))`
taken from its Poisson pass — theta itself has to have settled, not just the
log-likelihood. The old test was satisfied roughly 2e-5 of log-likelihood
early, which left theta 8e-7 out and dragged the coefficients 8e-6 with it.
- The first pass now runs as a genuine **Poisson** fit, as `glm.nb`'s does,
rather than a negative-binomial fit at a large stand-in theta. That pass
supplies both the first theta and the `d1` above.
- Later passes are **warm started** from the previous pass's means
(`etastart = log(mu)`), so they converge to the fit `glm.nb` reaches rather
than to the same optimum approached from a cold start.
- Theta is re-estimated at the means each pass **started** from, not the ones it
produced: `glm.nb` calls `theta.ml(Y, mu)` and only then reassigns
`mu <- fit$fitted.values`. `theta.ml` itself now also uses MASS's own stopping
rule, an absolute Newton-step tolerance of `.Machine$double.eps^0.25`.
- Across eighteen fits spanning dispersion from theta 0.41 to theta 69000, theta
now agrees with `glm.nb` to 3.4e-9, coefficients to 5.8e-9, standard errors to
8.4e-10 and deviance to 1.3e-9 — previously 8e-7, 8e-6, 3e-6 and 6e-7. The one
exception is genuinely near-Poisson data, where theta is not identified at all
(its own standard error exceeds the estimate, and `glm.nb`'s `theta.ml` reports
"iteration limit reached"); theta there agrees only to about 4e-6 relative,
while the coefficients still agree to 1.6e-10.
- A separate consequence: a negative-binomial fit with theta supplied was
starting from the Poisson `mustart` of `y + 0.1`, where R's
`negative.binomial()$initialize` sets `y + (y == 0)/6`. Different starting
values walk different iterates, and since the standard errors come from the
penultimate one, that showed as standard errors 6e-7 from R's while the
coefficients agreed to 1e-9. Such fits now match R to 2e-15.
- Note on that comparison: standard errors for a negative-binomial fit hold the
dispersion at 1, which is what `glm.nb` and `summary.negbin` do. R's
`summary.glm`, handed a `negative.binomial` family directly, instead *estimates*
the dispersion and prints standard errors scaled by its square root — 1.0839 on
one of the test data sets, so about 4% larger. Compare against
`summary(fit, dispersion = 1)` to see the values this module reports.
0.29 2026-08-03 CDT
[t_test]
- `t_test` was cross-checked against R's `stats::t.test` and `scipy.stats` case by
case, including the cases their own suites pin: R's regression tests
(`reg-tests-1a.R`, "t.test with one group of size one") and scipy's
`TestTTest_1samp`, `TestTTest_ind.test_special_cases`, `test_ttest_rel_ci_1d`,
`test_1samp_ci_1d` and `test_pvalue_ci`. On 2000 randomised comparisons against
R — all four modes, all three alternatives, random `mu` and `conf.level`, sample
sizes 2 to 40 and data scales spanning 1e-4 to 1e4 — the statistic and the
degrees of freedom agree to 2e-11 and the p-value to 3e-9, holding to eight
digits even where the p-value is subnormal (5e-310). What the comparison did
turn up was seven ways a call could come back wrong rather than loud, all of
them now fixed and covered by `t/t_test.t`.
- **`undef` was coerced to 0 instead of being dropped.** This is the one worth
re-running results over. `t_test` did not filter missing values, so a column
with gaps in it was tested with every gap counted as a zero: R gives
`t.test(c(1,2,NA,4,5))` a `t` of 3.286 on 3 degrees of freedom, and `t_test`
answered 2.588 on 4. No error, no warning, and an answer close enough to the
real one to look right. `undef` and `NaN` are now dropped the way R drops `NA`,
per-vector for a one-sample or unpaired test, and on complete cases when
`paired` so a half-missing pair goes whole rather than contributing a
difference against zero.
- **A `y` of fewer than two observations returned a silent `NaN`.** `var_y`
divided by `ny - 1`, so `t_test(\@x, [$one_value])` propagated `0/0` into the
statistic, the p-value and both interval bounds without raising. The two
thresholds R uses are now both in place: a Welch test needs a variance from
each side and refuses without one, while a pooled test tolerates a side of one
observation, since that side contributes no sum of squares. That second case is
what R's own regression suite pins — `t.test(y=x[1], x=x[-1], var.equal=TRUE)`
is a well-defined test with 8 degrees of freedom, and `t_test` now answers it
instead of returning `NaN` in one direction and croaking in the other. An empty
`y` is caught by the same check.
- **`alternative` was never validated.** The p-value helper fell through to
two-sided for any string it did not recognise, so a typo — `'gerater'` — ran a
different test than the caller asked for and reported nothing. It is now
checked the way R's `match.arg` checks it. `scipy`'s `"two-sided"` spelling is
unambiguous, so it is accepted rather than rejected.
- **A one-sided interval was wrong when `conf.level < 0.5`.** That case needs a
negative t quantile, and `qt_tail` searched upward from zero only, so it
returned roughly zero and collapsed the bound onto `mu`: R puts the upper bound
of `t.test(1:10, mu=5, conf.level=0.3, alternative="less")` at 4.9797 where
`t_test` reported 5.0000000036. `qt_tail` now reduces by symmetry first, so the
root it brackets is always positive.
- **`qt_tail` silently saturated at 1e6.** Past that its doubling loop gave up
and returned the ceiling, so `conf.level` of 0.99999999 and 0.9999999999 came
back with the *identical* interval, ±1048576, against R's ±6.4e7 and ±6.4e9.
The ceiling is gone; the loop now runs until `t * t` would overflow.
- **Interval accuracy no longer depends on the data's scale.** `qt_tail`
bisected to an absolute 1e-8 on the quantile, which is 1e-8 × `std.err` on the
interval — fine for data around 1, an error of 2 units for data around 1e9. It
now bisects to adjacent doubles. Worst interval error across the 2000
randomised cases went from 2.1e-4 to 5.3e-11 relative. At extreme
`conf.level` this makes `t_test` the more accurate of the two: `t.test` asks
for `qt(1 - alpha/2, df)`, and representing a 5e-9 tail as the double
`1 - 5e-9` costs eight significant figures of it, so R's own interval for
`conf.level=0.99999999` is off by 0.7 in the eighth digit. Working in the upper
tail throughout agrees with R's `qt(alpha/2, df, lower.tail=FALSE)` to 15
digits.
- **"Essentially constant" was an absolute test.** Only an exactly-zero variance
was rejected, so a spread below what a double can resolve at the data's own
magnitude was reported as a finding: four values around 1e10 differing by 1e-5
gave `t` = 4e15 and a p-value of 3e-47. The comparison is now relative, as R's
is. The exactly-zero case, where R returns `NaN`, raises here instead.
- **A defined non-array `y` was ignored.** `t_test(\@x, y => 5)` quietly ran a
one-sample test. It now raises. An explicit `undef` still means absent, as R's
`y = NULL` does.
- `qt_tail` is shared with `power_t_test`, which gains the same precision; it is
only ever called there with a tail below 0.5, so nothing about its behaviour
changes. `t_test` remains allocation-free — the missing-value filtering happens
inside the same single Welford pass that was already there, and the result hash
is built after the last error check rather than before the first.
[write_table: `tex.longtable.head`]
- A `longtable` freezes only the header sitting inside `\endfirsthead` /
`\endhead`, and `tex.longtable` never wrote those blocks — its header was an
ordinary first body row, leaving the frozen one to be hand-written by the
caller. That header then had no link to `col.names`: reorder the columns and
the labels at the top of every page keep the old order while the data below
them moves, and the generated header appears again as a duplicate first row.
- `tex.longtable.head` generates the repeat machinery from the table's own header
record, so it cannot drift. A true-but-numeric value emits
`\endfirsthead`/`\endhead`/`\endfoot` with no continuation caption; any other
true value is the caption used on pages after the first, written verbatim.
Implies `tex.longtable`. `tex.longtable` on its own is unchanged.
- The wrapper keeps one static token, the `\hline` closing its `\caption` line,
because a leading `\hline` in an `\input`ed file is a `Misplaced \noalign`
error — TeX has already begun the row by the time it expands the `\input`.
[skew, kurtosis]
- Two new XS functions describing the shape of a sample beyond its spread:
`skew` for the third central moment and `kurtosis` for the fourth. Both take
arguments the way `sd` and `var` do — numbers, array references or a mixture,
flattened into one sample — and both also accept `x => \@data` and
`type => 1|2|3`.
- `type` selects among the three sample conventions, which disagree noticeably on
small samples. The default is `type => 2`: `G1` and `G2`, the estimators
unbiased for a normal sample, as reported by SAS, SPSS, Stata, Excel's `SKEW()`
and `KURT()`, and `scipy.stats` with `bias => FALSE`. `type => 1` is the plain
moment ratio (`moments::skewness`) and `type => 3` is `b1`/`b2`
(`e1071::skewness`'s own default). All three, for both functions, agree with R
to about 1e-15. `kurtosis` returns *excess* kurtosis — 3 is already subtracted,
so a normal sample sits near 0.
- One pass, no allocation: the third and fourth central moments accumulate
through Welford's recurrence extended to higher moments (Terriberry) rather
than the textbook expansion in raw moments. That expansion is not usable on
real data — for a column of values around 1e7, a lab value in the wrong units
or a timestamp, `sum(x**3)/n` is about 1e21 while the third central moment is
single digits, so every significant figure cancels away.
- A constant sample croaks rather than returning a silent `NaN` from `0/0`, and
a `type` whose denominator the sample is too small for (`type => 2` needs
`n >= 3` for `skew` and `n >= 4` for `kurtosis`) says which.
- Both read tied arrays. `av_fetch` on a tied array returns a deferred `PVLV`
rather than the value, and `SvOK` on one of those is false until its
get-magic has run, so without an `SvGETMAGIC` every element of a tied array
looks undefined.
[median]
- `median` now reads tied arrays too. It already had a separate `av_fetch` path
for them — a tied array keeps nothing in `AvARRAY`, so the fast path would read
off a null pointer — but that path was missing the `SvGETMAGIC` described
above, so it rejected every tied array as undefined instead of computing the
answer. `mean`, `sd`, `var`, `sum`, `min` and `max` still reject tied arrays.
They have no `AvARRAY` fast path to guard, so they croak rather than crash, and
the same one-line fix would make each of them work.
[oneway_test]
- `oneway_test` was cross-checked case by case against R's `stats::oneway.test`
(both branches), R's `anova(aov())` for the `Sum Sq` / `Mean Sq` columns,
`statsmodels.stats.oneway.anova_oneway(use_var="unequal")` and
`scipy.stats.f_oneway`. The 37 data sets are R's own built-ins — `chickwts`,
`InsectSprays`, `PlantGrowth`, `iris`, `ToothGrowth`, `mtcars`, `warpbreaks`,
`sleep`, `airquality`, `CO2`, `esoph`, `OrchardSprays`, `faithful`, `quakes` —
plus hand-built numerical edge cases. The statistic and both degrees of freedom
already matched R everywhere; what the comparison turned up was four ways a call
could come back wrong rather than loud, all now fixed and covered by
`t/oneway_test.R.scipy.t`. Statistic, degrees of freedom and p-value now agree
with R to 1.3e-12 relative error across all 37, and on 2000 randomised
comparisons against R — both branches, 2 to 8 groups, sizes 2 to 40,
deliberately heteroscedastic, data scales 1e-4 to 1e4 — the statistic and the
degrees of freedom agree to 1e-12 and the p-value to 8e-11, the worst of those
being a p-value of 2.4e-66.
- **Every p-value below about 1e-16 was returned as a flat 0.** `Pr(>F)` was
built as `1 - pf(F, df1, df2)`, and 1 minus something that close to 1 has no
bits left to carry the answer: `faithful` split at `waiting > 70` should give
`1.2099104551915e-76` under Welch and `5.50783574504386e-103` pooled, and
`oneway_test` reported `0` for both. Anything from about 1e-9 downward was
losing relative precision the same way, quietly — `ToothGrowth` by dose came
back as `9.99200722162641e-16` against R's `9.53272701169993e-16`, off by 4.6%
with nothing to indicate it. The p-value is now evaluated in the upper tail
directly, via the beta symmetry `1 - I_x(a, b) = I_{1-x}(b, a)`, so no
subtraction from 1 happens at any point, and the range down to the smallest
representable double is reported at full precision.
- **An `F` of `Inf` produced a p-value of `NaN` instead of 0.** When every group
is constant but their means differ, the within-group sum of squares is 0 and
`F` is legitimately infinite; R reports `p = 0`. `pf` formed
`df1*f/(df1*f + df2)`, which is `Inf/Inf` — a `NaN` that propagated straight
into `Pr(>F)`. `Inf` and `NaN` are now handled explicitly, matching R's
`p = 0` and `p = NaN` respectively.
- **A `NaN` Welch denominator df was reported as 1e300.** A group with zero
variance gets an infinite Welch weight, which makes R's `tmp` term `NaN` and
its denominator df `NaN` with it. `oneway_test` had a `(tmp > 0.0)` guard that
a `NaN` fails, so it substituted a magic `1e300` — a number that reads as a
real, very large degrees of freedom and would be believed as one. The guard is
gone; `Residuals`/`Df` and `Residuals`/`Mean Sq` are `NaN` there, as in R.
- **`formula` mode read `undef` and non-numeric response cells as 0.0.** The
hash and array-of-arrays shapes were already fixed to die on these (and pinned
by `t/oneway_test.bugs.t`), but the formula path has its own fill loop and was
missed, so
oneway_test({ y => [1, 2, 3, undef, 5, 6], lab => [qw(a a a b b b)] },
formula => 'y ~ lab');
- silently tested group `b` as `(0, 5, 6)` — a mean of 3.67 instead of 5.5 — and
returned an `F` of 0.735 with no complaint. All three input shapes now enforce
the documented contract identically.
- Two places where `oneway_test` is the more accurate side and the reference is
not, now documented rather than treated as disagreements: the sums of squares
are accumulated two-pass, so on two groups near 1e8 `Residuals`/`Sum Sq` is
exactly `10` where R's QR-based `anova(aov())` gives `10.0000000521067`; and
where the exact between-group sum of squares is 0, `oneway_test` returns 0
rather than R's 1e-30-scale residue.
0.281 2026-08-03 CDT
- `median` (LikeR.xs) — the same answers in about an eighth of the time. On the
`benchmark.pl` case (10,000 normals in one array ref) a call went from 0.83 ms
to 0.097 ms, which puts it ahead of the two implementations it was behind:
`numpy.median` at 0.105 ms and R's `median` at 0.196 ms, measured on the same
machine. Small samples — a per-group median under `agg` or `group_by`, which is
where most calls to it come from — went from 647 ns to 251 ns.
- A median is the middle one or two values, so most of the sort the function used
to do was wasted work: `qsort` orders all n elements at a cost of n log n
comparisons, every one an indirect call through a function pointer the compiler
cannot see into, to answer a question that depends on one or two of them. Those
values are now selected instead, and the sample is walked once rather than
twice.
- The selection is introselect, the same shape numpy's `partition` uses:
quickselect with a median-of-three pivot, an insertion sort once a range is
small, and a heapsort fallback past a depth limit, so an input crafted to
defeat the pivot choice degrades to O(n log n) rather than O(n²). The awkward
data people actually have comes out faster than random data rather than
slower — sorted, reversed, all-equal and organ-pipe samples of 100,000 values
each take about 0.25 ms against 0.90 ms for random ones. For an even count the
lower of the middle pair is the largest value left below the upper one, which
a scan of that side finds without a second selection.
- The counting pass is gone. It walked every element through `av_fetch` before
any arithmetic, only to size the buffer; the array lengths give the same
count, exactly, because an undef anywhere still dies.
- The pass that remains reads cells through `AvARRAY` instead of `av_fetch`,
with tied arrays kept on the `av_fetch` path, since only it sees their values.
- A sample of 256 values or fewer is copied to the C stack rather than the heap,
so the common small call no longer pays for a malloc and free at all: 10,000
of them now grow RSS by nothing. Larger samples still copy n values, which is
what leaves the caller's array in its original order — the selection reorders
whatever it works on, and `t/median.t` checks that the input comes back
untouched.
- Error messages that carry an index or a count were unreadable on older perls,
in nineteen places across LikeR.xs, and now are not. `croak` runs perl's own
formatter rather than the C library's, and that formatter does not understand
C99's `z` length modifier: it printed the conversion literally, so
`median(1, undef)` on perl 5.10 or 5.12 said `undefined value at argument index
%zu` instead of naming the argument that was undefined — the one thing the
message existed to say. `min`, `max`, `mode`, `sum`, `sd`, `var`, `median`,
`mcnemar_test`, `friedman_test`, `hoa2hoh` and `oneway_test` were all affected.
They now use `UVuf`, as the rest of the file already did. The `snprintf` calls
elsewhere in the file are unaffected and unchanged: those do go to the C
library, where `%zu` means what it says.
- New `t/croak.messages.t` covers every one of those messages: each is triggered,
checked for the number it should name, and then swept for any conversion left
unexpanded, which is what will catch the next one written with `%zu`. Run
against the code as it stood before this change, thirty-two of its assertions
fail on perl 5.10.
- New `t/median.t`: every length from 1 to 25 and from 254 to 258 — either side of
the insertion-sort cutoff and of the point where the buffer moves off the stack
— across sorted, reversed, all-equal, two-valued, duplicate-heavy, organ-pipe
and median-of-three-killer samples, each checked against a plain Perl sort,
together with the error messages and `Test::LeakTrace` over the stack, heap,
mixed-argument and croak paths.
- fix for threaded Perls https://www.cpantesters.org/cpan/report/2dbacf8f-7138-1014-a1ab-f0f91cf3b922
0.28 2026-08-02 CDT
- `p_adjust` (LikeR.xs) now takes a data frame as well as a flat list of
p-values, and hands the corrected values back in the shape they arrived in. An
AoA, AoH, HoA or HoH goes in and a new frame of the same kind comes out, with
the same rows, columns and row labels; the input is left alone. Everything the
flat form did is unchanged — an arrayref of p-values still returns a list, in
order, with the same numbers.
- `columns => 'p_value'` (or an arrayref of names, or 0-based positions for an
AoA) says which columns hold p-values, and copies the rest of the frame
through untouched, so a results table with a `gene` column no longer has to
be taken apart and put back together around the call. Without `columns`
every cell is treated as a p-value, which is right for a frame that is
nothing but p-values; a label column in one dies with a message naming the
offending value and pointing at `columns`, rather than correcting a string
coerced to zero.
- All the p-values in the frame are corrected as one family, whichever shape
they came in, so the family size is the number of p-value cells.
- The method still reads positionally and may now also be given as
`method => ...`. `none`, which the function has always accepted, is now
documented along with the rest.
- Cells are visited in a fixed order — by row and then column name, or column
name and then row for a HoA — so tied p-values break the same way on every
run instead of following hash iteration order.
- `drop_duplicates`, `filter`, `t_test`, `vals`: speed/RAM improvements
- Incompatible: the `'?'` / `'h'` argument added in 0.27 is gone (lib/Stats/LikeR.pm). `agg('h')`, `read_table('?')` and the fifty-odd other pure-Perl functions that took it no longer print help and die — they treat the string as data, the way the XS functions always have. `h('agg')`, `h(*agg)` and `h(\&agg)` are unchanged and remain the way to ask, for every function in the distribution.
- It was a help route that only half the module had, so what a lone `'h'` meant depended on whether the callee happened to be written in XS or in Perl, and a column, file or option value really named `'h'` needed `$Stats::LikeR::HELP = 0` to get through. That variable is gone too; nothing reads its arguments for a help flag any more.
- `bedroc` still prints its own short XS usage summary for `bedroc('h' | 'H' | '?')`, which is hand-written and predates all of this.
- `merge` (LikeR.xs) — same joins, a third of the time and a fifth of the memory. Nothing about the result changes: every join type, shape combination and edge case produces exactly what it did before, and `t/merge.t` now checks all six input/output paths against a plain-Perl reference join over a randomized corpus.
- The old implementation transposed both frames into arrays of row hashes, joined those, and transposed the result back. A 10,000-row HoA joined to itself therefore built 20,000 throwaway row hashes and copied every cell three times before returning. It now reads each frame where it lies — a HoA column by column, an AoH/HoH row by row — and writes the result straight into the shape being returned, so the only cells copied are the ones the caller keeps.
- The right frame's index is a hash of row numbers chained through a flat array, rather than an array-ref of index scalars per distinct key, and one reused buffer builds every join key instead of one scalar per row.
- Column names are resolved to their column (HoA) or interned once as shared hash keys (AoH/HoH) before the join starts, so the per-row work is a lookup rather than a lookup and a rehash.
- Measured on the `benchmark.pl` case (two 10,000-row frames, six columns, inner join on `id`): 0.052 s and 41.4 MB before, 0.017 s and 7.3 MB after. An outer join of the same frames went from 0.113 s to 0.008 s.
- `write_table` (LikeR.xs) — two changes, one of them incompatible.
- Every format now prints the coloured `wrote <file>` confirmation line, not just LaTeX and `.xlsx`. Delimited output (csv/tsv) was silent before. The line is identical in all cases: the file name in black on cyan, with the SGR codes inline so there is still no `Term::ANSIColor` dependency. Nothing is announced when nothing is written.
- **Incompatible:** `row.names` now defaults to **off** in every format. It previously defaulted **on** everywhere, following R's `write.table`, which meant a call that said nothing about row names got a label column and a leading empty header cell (`,gene,n`) it had not asked for. Pass `row.names => 1` for the old behaviour; `row.names => 'col'` is unchanged.
- New `h2aoh` and `aoh2h` (lib/Stats/LikeR.pm), which add the flat hash to the shapes the conversion family understands. A plain hash is a two-column table folded shut, and until now nothing would unfold it: `value_counts` hands one back, and no frame function would take it.
- `h2aoh(\%h, var_name => .., value_name => ..)` unfolds a flat hash into a two-column AoH, one row per pair, under column names the caller picks. `sort => 'key' | 'value' | 'none'` fixes the row order, which hash iteration otherwise leaves to chance; `'value'` is biggest-first for numbers, so `value_counts` output comes out the way pandas' `Series.value_counts()` orders it.
- `aoh2h` folds a two-column AoH back down, with `duplicates => 'die' | 'first' | 'last'` deciding what a repeated key means. The two are exact inverses under their defaults.
- The column options are named `var_name` / `value_name` after `melt`, which emits the same two columns. R spells this pair `tibble::enframe()` / `deframe()`; pandas spells it `pd.Series(d).reset_index()` and `Series.to_dict()`.
0.27 2026-07-26 CDT
- New `h` function: `h('agg')`, `h(*agg)` or `h(\&agg)` prints that function's section of this document and returns, in the spirit of R's `?function`. `h()` lists every documented function. It covers the XS functions as well as the Perl ones, because it looks the name up in the module's POD instead of reading an argument list — see [Getting help](#getting-help).
- The pure Perl functions also accept `'?'` or `'h'` in place of their arguments, which prints the same text and then dies. `$Stats::LikeR::HELP = 0` switches that off for code that has to pass a column or file really named `'h'`.
- `qcut`'s hand-written usage message was replaced by its section of this document; `qcut('h')` and `qcut('?')` still die, but `qcut('H')` no longer means help.
- speed improvements in calculation of Kendall tau and p-value. Improvement of writing xlsx files that won't show in time, but pure waste was removed.
- Addition of `auc`, `auroc`, `cmh_test`, `epi_2x2`, `roc` functions
- `prcomp` now accepts AoH input
- glm extended (LikeR.xs)
- family => 'poisson' (log link) and family => 'negbin' — negative-binomial θ estimated by ML via a MASS::glm.nb-style outer loop, or fixed with theta =>. Matched R to ~1e-8 (coefs, deviance, null-dev, AIC, SE, θ); exact Poisson limit when data aren't over-dispersed.
- Every non-gaussian family now returns exp (odds/rate/incidence-rate ratios + conf.low/conf.high), link-scale conf.int, conf.level, and theta (negbin). Count families report z-statistics. OR/CI matched R's confint.default exactly.
- New XS tests (all matched R exactly)
- prop_test — 1/2/k-sample proportions (Yates, Wilson & Wald-diff CIs)
- mcnemar_test — matrix or paired vectors; continuity correction; exact => 1 binomial
- friedman_test — repeated-measures rank test, tie-corrected
- dunn_test — post-Kruskal pairwise, 7 adjustment methods
- New Perl functions (lib/Stats/LikeR.pm, matched base-R references)
- Effect sizes: cohen_d (+Hedges g, CI), smd, cramers_v (+Bergsma bias-corrected), eta_squared (η²/partial/ω²)
- vif, hosmer_lemeshow (matches hoslem.test)
- age_standardize — direct standardization + Fay–Feuer gamma CI (matches epitools::ageadjust.direct)
0.26 2026-07-20 CDT
- https://www.cpantesters.org/cpan/report/fc7d01a0-83f4-11f1-b543-8a9ac547de9a
Fixed a long-double issue
0.25 2026-07-19 CDT
- https://www.cpantesters.org/cpan/report/3376f80e-83bf-11f1-a5f3-44496e8775ea
- Fixed a use-after-free in `fisher_test` on the hash (HoH) input path: the "row is missing column key" error freed its scratch arrays and then read the key strings back out of them to build the croak message. This was harmless on glibc but crashed (`SIGBUS`) under stricter allocators such as FreeBSD's, failing `t/fisher_test.t` on CPAN smokers. The key pointers are now captured before the arrays are freed.
0.24 2026-07-19 CDT
- `interpolate`'s numeric core moved from pure Perl to XS (`_interp_column_xs`): ~5× faster for `linear` on large columns, ~11× for `pchip`, and ~50× for the spline methods whose dense solve dominates. Results are unchanged (bit-for-bit versus the former Perl kernels).
- `Ronly` now accepts one or more array references (like `Lonly`), returning the values found only in the last reference; the two-argument form is unchanged, and `Ronly(@refs)` equals `Lonly(reverse @refs)`.
- `interpolate` gains full `pandas.DataFrame.interpolate` method parity: `nearest`, `zero`, `slinear`, `pad`/`ffill`, `bfill`/`backfill`, `quadratic`, `cubic`, `cubicspline`, `pchip`, `akima`, `barycentric`, `krogh`, `polynomial`, `spline`, and `index`/`values`/`time`, plus an `x` argument for custom abscissae and an `order` argument. Matched to pandas/scipy within 1e-6.
- `t/transpose.t` no longer loads `Devel::Confess` in its leak tests: its `$SIG{__DIE__}` stack-trace objects landed in `$@` and were reported as leaks by `Test::LeakTrace` on the croak paths under older perls (e.g. 5.12.3). The die-path leak checks now also clear `$@` so the exception object cannot be miscounted.
- `cfilter` simplification, use of `qr///` filtering on columns
- `summary` output now looks more like `view`, and accepts HoH
- `fisher_test` can compute larger tables than just 2x2
- `read_table` reads xlsx files significantly faster and with less RAM.
- Addition of `bfill`, `drop_duplicates`, `ffill`, `melt`, and `pivot_table`
- Original `Lonly` code removed, as it was a special case of `get_unique`, and `get_unique` was re-named to `Lonly`.
- Removal of `Devel::Confess` from testing and dependencies.
0.23 2026-07-10 CDT
- `rename_cols` takes HoH as input
- `write_table` prints row names as first column; writes longtable with comments
- `assign` gains `map_cell { ... }` for in-place per-cell column edits
[assign]
- `assign` now accepts a third kind of column value, `map_cell { ... }`, for editing an existing column in place — no "copy, substitute, return" boilerplate and no dependence on `s///r` (unavailable on the older perls this module supports).
- Inside a `map_cell` block, `$_` is the **named column's current cell** (not the whole row), the block's return value is **ignored**, and the modified `$_` is stored back: `assign($df, 'Res.' => map_cell { s/^[A-Z]:// })`.
- The row is still available as `$_[0]` (sibling columns), the index as `$_[1]`, and the row key as `$_[2]` (HoH only).
- **Undef/missing cells pass through untouched** (undef in → undef out): the block is skipped for them, so `s///` never warns on an uninitialized value.
- Supported on all three shapes; for HoA the target column must already exist. A plain `sub { ... }` is unchanged, so `map_cell` is purely additive.
- `map_cell` is exported alongside `assign`.
- Tests: `assign.t` (AoH + HoA) and `assign.HoH.t` gained `map_cell` coverage — in-place `s///`, `$_[0]`/`$_[1]`/`$_[2]` context, new-column-from-undef, the missing-HoA-column death path, and `no_leaks_ok` guards. Verified building and passing the full suite on perl 5.10.1, 5.12.5 (long-double), and 5.42.2.
[`group_by`]
- Fixed group_by to honor all filter hashrefs (option 1)
- Root cause: the XS captured only ST(3), so every filter hashref after the first was silently dropped — including the README's documented multi-hashref form.
- Change (LikeR.xs):
- Removed the single-ST(3) capture and the filter_hv PREINIT var.
- Added a FOR_EACH_FILTER(body) macro that walks the arg stack from ST(3) to ST(items-1), iterating every { column => sub } pair and ANDing them together. It iterates the stack directly rather than heap-collecting the hashrefs, so a croaking filter sub still can't leak anything (verified). Non-hashref args are skipped.
- Rewrote the filter loop in all three branches (AoH / HoA / HoH) to use the macro, keeping each branch's own value-fetch logic.
- One build wrinkle worth noting: xsubpp parses every non-# line in the inter-XSUB region as a candidate function signature, so a /* ... */ comment there breaks the build (it tried to parse column => sub / (ST(3)..) as a signature). I moved the macro's documentation into the XSUB body (real C) and left the macro comment-free, matching the existing EVAL_FILTER style.
- Tests (t/group_by.HoH.filter.t, 17 assertions):
- HoH single-column filter, AND filter (both the one-hashref and separate-hashref forms now give identical results), no-match → empty hash, missing/undef target excluded despite passing the filter, and no_leaks_ok
- Mentioning a non-existent column is now fatal.
0.22 2026-07-07 CDT
- returned `Devel::Confess` to required dependencies to fix for CPAN testers.
0.21 2026-07-07 CDT
- Better warning message for undefined data for `aoh2hoh`, `assign`, `dropna`
- addition of `agg`, `concat`, `drop_cols`, `rank`, `rename_cols`, `select_cols` functions
- Improving Kwalitee (sic): added `[PodWeaver]` to dist.ini; as well as `Changes` file
[`assign`]
- `assign` now accepts two kinds of column value, so a function that already returns a whole column (like `rank`) drops in without wrapping.
- **Per-row coderef** (unchanged): called once per row, `$_` is the row, and the single scalar it returns is the cell. A single arrayref return is still stored *as the cell*, so arrayref-valued columns keep working.
- **Whole-column coderef** (new): if the coderef returns a *list* of more than one value, that whole list becomes the column, laid down positionally. This is what makes `'ΔG rank' => sub { rank( vals($df, 'dG_kcal_mol') ) }` work directly — no `[ ... ]` needed.
- **Arrayref value** (new): a ready-made column, e.g. `col => [ rank(...) ]`, copied into the frame.
- The coderef is probed once (row 0 for AoH/HoH, the first synthesized view for HoA) to decide per-row vs whole-column, so per-row code is never run twice on row 0. Every column value is length-checked against the row count and a mismatch dies. HoH is now a supported, documented shape alongside AoH and HoA; whole-column and arrayref values align to sorted key order.
- Tests: `assign.t` (AoH + HoA) and `assign_HoH.t` were expanded to cover every shape × value-kind combination — per-row scalar, whole-column list, arrayref value, single-arrayref-as-cell, `rank()` integration, chaining, `$_[1]` index, `$_[2]` row key (HoH), overwrite, ragged HoA columns, empty frames, length-mismatch and bad-value / odd-arg / non-hash-row death paths, and `no_leaks_ok` guards on the new whole-column and arrayref paths.
[`read_table`]
- Fixed handling of commented-out header lines and made filter columns
referenceable by the name as it appears in the file.
- **Commented-out header recovery.** `_parse_csv_file` treats a line whose
comment marker is followed by whitespace (e.g. `# PDB<TAB>score`) as a
comment and drops it, so a header written that way never reached the
callback and the first *data* row was silently mistaken for the header.
`read_table` now recovers it: the first physical line, if it is
`marker + whitespace` and splits into two or more fields, is held as a
candidate header and confirmed only when its field count matches the first
data row. If the counts disagree the candidate was an ordinary leading
comment and is discarded, so a prose comment that happens to contain the
separator (e.g. `# note, see README`) is never mistaken for a header. A
marker hugging its text (`#id,val`) is delivered by the parser and
un-commented in the callback as before. The marker and any following
whitespace are stripped, so `# PDB` is stored as the clean name `PDB`.
- **Filter columns may be named as written in the file.** Filter keys are
matched against the header by exact name first, then retried with the
leading comment marker (and surrounding whitespace) stripped, so a
commented-header column resolves whether it is referenced as `# PDB` or by
its clean name `PDB`:
read_table(
'regression_rank.tabular.tsv',
filter => { '# PDB' => sub { $_ == 2 } },
);
- **Clearer "column not found" error.** The failure now names the file and
lists the actual header instead of printing it to STDOUT (a library
shouldn't print):
read_table: Filter column 'nope' not found in the header of FILE;
header is: 'PDB', 'score'
0.20 2026-07-05 CDT
- addition of `ncol`, `nrow`, and `pnorm` functions
- `filter` can filter by row names with `$_[1]`
- `view` now accepts array of arrays in addition to AoH, HoA, and HoH
[csort]
- Two behavioural changes, both contained to the `csort` XSUB (the `cs_*` helpers are untouched).
- Row names survive a Hash-of-Hashes sort. Sorting a HoH previously discarded the outer keys. Now each row is folded into a *fresh* row hash (a private container over aliased, read-only cells) that carries its outer key under a `row.name` column, so the name flows into whichever shape you request:
my $hoh = { alpha => { id => 1 }, beta => { id => 2 } };
csort($hoh, 'id'); # AoH: each row gains a row.name field
csort($hoh, 'id', 'hoa'); # HoA: an aligned row.name column
- The column name defaults to `row.name` and can be overridden with an optional 4th argument (mirroring `hoa2hoh`'s named-key style): `csort($df, 'id', 'aoh', 'sample')`.
- The outer key is authoritative — it wins over any pre-existing same-named field in the row.
- Once present, the column is sortable like any other: `csort($hoh, 'row.name')`.
- Because rows are now *copied* rather than shared, the caller's HoH is never mutated by the injection. (Minor behaviour change: output rows are no longer the same refs as the source rows.)
- Clearer usage message. The signature is now `csort(...)`, so xsubpp no longer emits the misleading auto-generated `Usage: Stats::LikeR::csort(data, by, output=&PL_sv_undef)`. Argument count is checked by hand, and the croak now shows both real calling forms:
Usage: csort($df, 'column.name', 'HoA')
or csort($df, sub { $b->{'No.'} <=> $a->{'No.'} }, 'hoa')
(optional 4th arg names the row-name column when sorting a HoH; default 'row.name')
- `data`/`by`/`output` are read as `ST(0..2)`; `output` still defaults to matching the input shape.
- Tightened validation messages. The `$data` croak now reads `hash-ref (HoA or HoH)`, and the `$by` croak includes a concrete example: `a column name (e.g. 'No.') or a comparator code-ref using $a and $b, e.g. sub { $b->{'No.'} <=> $a->{'No.'} }`. Existing HoA croaks (`unequal lengths`, `not found`, `not an array-ref`) are unchanged.
- When sorting, undefined values in the sorting column are placed at the bottom
[cor]
- Fixed an unsigned-integer underflow in `kendall_tau_b` and added a regression test.
- Bug:
- In `kendall_tau_b`, concordant/discordant counts `C` and `D` are declared `size_t` (unsigned). The numerator was computed as:
return (NV)(C - D) / denom;
- The subtraction `C - D` happens in unsigned arithmetic *before* the cast to `NV`. When discordant pairs dominate (`D > C`), the result wraps to a huge positive value instead of going negative.
- For the arrays:
dG_kcal_mol: -7.765, -9.328, -10.326, -9.038, -9.608, -9.779, -9.975, -6.906
anomaly_rank: 154, 155, 161, 188, 76, 172, 173, 69
- there are `C = 9` concordant and `D = 19` discordant pairs (no ties). `9 - 19` wraps to `18446744073709551607`, so the function returned ~`6.6e17` instead of the correct `-10/28 = -0.3571428571`.
- Fix:
- Cast each operand to `NV` before subtracting, so the arithmetic is signed:
return ((NV)C - (NV)D) / denom;
- Only that one line changed. The denominator sums (`C + D + tie_x`, `C + D + tie_y`) are non-negative, so they were left as-is.
- Regression test — `cor.t`:
- Kendall on the offending arrays pinned to `-0.3571428571`.
- Explicit `[-1, 1]` range guard (the real backstop — the pre-fix value `~6.6e17` blows past the bound regardless of exact magnitude), plus a negative-sign assertion.
- Pearson (`-0.4889102301`), Spearman (`-0.4761904762`), and default-method coverage of the three `compute_cor` branches.
- Kendall boundary cases: perfectly concordant (`+1`), perfectly discordant (`-1`), self-correlation (`+1`), and a tie case exercising `tie_x` in the denominator.
- `no_leaks_ok` per method (guarded with `unless $INC{'Devel/Cover.pm'}`).
- Croak paths: length mismatch, unknown method, zero-variance input.
[XS refactor]
- Consolidate helper functions to reduce binary size, find bugs, and back the changes with tests. Every change was validated by translating the XS (`ExtUtils::ParseXS`) and compiling the result
with the module's own `ccflags`.
- Outcome:
- **Net change to the source:** ~154 fewer lines; helper-function count down by 4 (7 removed, 3 added).
- **Genuine bugs fixed:** two instances of the same latent defect (see below). The rest of the work was behavior-preserving consolidation.
- Function consolidation:
- | Change | Before | After |
|---|---|---|
| Three-way `NV` comparator | `compare_rank`, `cmp_rank_item`, `cmp_rank_info`, `compare_NVs` | single `cmp_nv3` (reads the leading `NV` member, valid for `RankInfo`/`RankItem`/raw `NV`) |
| Average-rank routine | `compute_ranks` + `compare_index` restoration sort | existing `rank_data` (scatters ranks into `out[idx]`, no second sort) |
| String comparator | `cmp_string_wt`, `lm_str_qsort` (byte-identical) | single `cmp_string_wt` |
| Multiplicity filter & set difference | `intersection` + `get_unique` (~90% shared); `Lonly`/`Ronly` duplicated bodies; a separate `set_difference()` | one shared `set_multiplicity()` with an "all vs. one" mode flag and a `from_last` flag: `intersection` (all), `Lonly` (one, first array), `Ronly` (one, last array) |
- All merges were confirmed behavior-preserving: the collapsed comparators are
equivalent on ordinary values, `NaN`, and infinities, and `compute_ranks` and
`rank_data` produce identical average ranks.
- Bugs:
- Two comparators stabilized their sort by returning `a->idx - b->idx` directly,
where the index field is an unsigned `size_t`. The subtraction wraps and is then
truncated to `int`, which is implementation-defined and gives the wrong sign
once a difference exceeds `INT_MAX`.
- `compare_index` — removed entirely (the routine that used it, `compute_ranks`, was replaced by `rank_data`).
- `cmp_pval` — the tie-break comparator in the p-adjust path. **Missed in the initial review; found later** via a `-Wconversion` compile of the earlier source. Fixed to compare with the `(a > b) - (a < b)` idiom.
- Caveat on severity: on every mainstream ABI (LP64, LLP64, ILP32), the
low-word truncation happens to reproduce the correct sign for any array smaller
than ~2^31 elements, so this never produces a wrong result at realistic sizes.
It is a portability/UB issue, not a runtime failure, which is why no functional
test detects it (see "Testing", below).
- `LikeR.xs` — consolidated helpers; `compare_index` removed; `cmp_pval` fixed.
[`view`]
- non-ASCII characters now print
[`write_table`]
- new option to output to LaTeX table
0.19 2026-07-01 CDT
- numerous `SSize_t var1 = av_len(var) + 1` are changed to `size_t var1 = av_len(var) + 1` as `size_t`; as the result cannot be negative, in order to expand numerical range
- Addition of `hoa2hoh`, `binom_test`, `chunk`, `get_union`, `get_unique`, `Lonly`, `Ronly`, `qcut`, and 3 tukey functions
- Better warnings when non-array references are given to `intersection`
- `view` now breaks columns into chunks for very wide data sets, more closely matching R's behavior
0.18 2026-06-28 CDT
- `restrict` keyword added to numerous places within `intersection` to decrease CPU time
- fix to dist.ini for dependencies
- fixed POD rendering
0.17 2026-06-23 CDT (approx)
- addition of `assign`, which adds new columns based on calculations from other columns
- addition of `hoa2aoh`, transforming hash of arrays to array of hashes
- addition of `predict`, using results from `aov`, `glm`, and `lm`
- addition of `aoh2hoh` transforming array of hash into hash of hashes, `intersection`, `uniq`, and `vals`
[`aov`]
- Bug fixes:
- **`size_t` underflow on empty arrays.** Three loops were bounded by `av_len(...)`
compared against an unsigned counter; `av_len` returns `-1` for an empty array,
which turned `k <= len` into a `SIZE_MAX` loop. The `stack()` value loop, the `.`
column-expansion loop, and the `group.stats` column loop now use a signed
`SSize_t` bound.
- **HoH row count.** Row count for hash-of-hashes input was taken from the return
value of `hv_iterinit`; it now uses `HvUSEDKEYS(hv)` with a separate
`hv_iterinit`, matching `predict`.
- **Buffer overflow in interaction parsing.** `strcpy(right, colon + 1)` into a
fixed `char right[256]` is now `snprintf(right, sizeof(right), ...)`.
- Performance / memory:
- **Removed the per-row `row_x` scratch allocation.** Design rows are built
directly into `X_mat[valid_n]`; `valid_n` simply does not advance on a rejected
row. Interaction columns read their operands from the same in-progress row, so
the logic is unchanged.
- **`row_names` is no longer dead.** Surviving row names are transferred (pointer
move, no copy) into `surv_names` to key `fitted.values`; rejected rows are freed
in place.
- **Dropped a `restrict` UB.** `orig_data_sv` aliases `data_sv`; the `restrict`
qualifier was removed.
- New, `predict`-compatible output keys:
- **`coefficients`** — OLS estimates recovered by back-substitution on the R factor
left in `X_mat` against Q'y in `Y` (no re-derivation). Keys are the expanded term
names (`Intercept`, continuous names, `base.level` dummies, and `a:b` interaction
products). Aliased columns are reported as `NaN`, which `predict` drops.
- **`fitted.values`** — `Xb` over the non-aliased columns, keyed by surviving row
name. Computed from a snapshot of the design (`Dsav`) taken before the QR
overwrites `X_mat`. Costs one transient copy of the design matrix; negligible for
typical ANOVA where the column count is small.
- **`xlevels`** — sorted level list per factor, index 0 = reference, aligned with
the contrast coding used to build the dummies.
- **`family`** — `"gaussian"`.
- Cleanup-path correctness:
- `xlevels_hv`, `Dsav`, and `surv_names` are freed on both the "0 degrees of
freedom" croak and the normal exit. The interaction-main-effects croak in
PHASE 3 also frees `xlevels_hv`.
- Known limitations (unchanged):
- The intercept-stripping string surgery (`-1`, `+0`, `+1`, ...) operates on the
whole RHS and can still mangle `I(x-1)`-style transforms; treat `I()` with
arithmetic constants carefully.
- Top-level keys `coefficients` / `fitted.values` / `xlevels` / `family` /
`group.stats` share the return hash with the ANOVA rows; a predictor literally
named one of those would collide.
[`predict`]
- New: factor-bearing interaction terms:
- Previously, interaction coefficients such as `GroupB:Sexmale` or `GroupB:x` fell
through to the continuous `evaluate_term` path and died on a nonexistent column.
They are now handled directly:
- **`dummy_hv`** stores each dummy's factor base index (an `IV`) instead of
`&PL_sv_yes`, so a dummy name maps back to its `(base, level)` in O(1)
(`level == name + strlen(base)`). `hv_exists` lookups are unaffected.
- During coefficient caching, any `:` term with at least one factor-dummy component
is routed to a separate list (`icopy` / `ibeta`); pure-continuous interactions
(e.g. `x:z`) stay on the existing `evaluate_term` path, so prior behavior is
preserved.
- Each routed term is parsed once into flat component arrays. Factor components
store a base index and level pointer; continuous components store the term string
and get the same up-front column-existence validation as main terms.
- Per row, each factor's raw level is read once into `raw_lv[]` and reused by both
main effects and interactions (no duplicate `get_data_string_alloc`). An
interaction's value is the product of its components: a factor component
contributes `1.0` iff the row's level matches the dummy's level (reference levels
give `0`), continuous components go through `evaluate_term`.
- This covers factor×factor, factor×continuous, continuous×continuous, and n-way
combinations.
- Other:
- HoH row count uses `HvUSEDKEYS` (already present).
- The unseen-factor-level croak now frees every level string already read for the
current row, not just the current one.
[Tests]
- **`aov.t`** — one-way ANOVA against hand-computed values (Df / Sum Sq / Mean Sq /
F / decomposition); identical results across HoA / HoH / AoH / stacked input;
simple regression; `.` expansion; intercept removal (`-1`); two-way with
interaction (Type I SS on a balanced design); NaN listwise deletion; all croak
paths; leak checks.
- **`predict.t`** — `predict(training) == fitted.values` round-trips for one-way,
regression, factor×factor, factor×continuous, and continuous×continuous models;
explicit predicted values; agreement across HoA / AoH / HoH / flat newdata;
no-newdata path; binomial `link` vs `response`; gaussian identity link; all croak
paths; leak checks.
- Leak tests use `no_leaks_ok` guarded by `unless $INC{'Devel/Cover.pm'}` and skipped
when `Test::LeakTrace` is absent.
- Assumptions worth confirming:
- The NaN-deletion test relies on `evaluate_term` returning `NaN` for a non-finite
response value (an `Inf - Inf` NaN is fed in deterministically).
- The continuous×continuous round-trip relies on `evaluate_term("x:z")` yielding
`x * z` — the same assumption the pre-existing `predict` continuous-interaction
path already made. If that path was untested, this round-trip now exercises it.
[`view`]
- now returns colored output; fixed bug with incorrect widths; undefined values show as `undef` rather than `NA`, as in Data::Printer
[`csort`]
- now accepts Hash of Hashes; addition of `restrict` which should decrease calculation time
[filter]
- **Added hash-of-hashes (HoH) input.** In addition to AoH and HoA, `filter` now accepts an HoH (`{ key => { col => val, ... }, ... }`); each inner hash is one row, and matching keys are preserved by default (HoH -> HoH).
- **Added `output.type`.** `filter($df, $pred, 'output.type' => 'aoh'|'hoa')` selects the returned shape (aliases `out` / `output_type`; a bare positional type also works). When omitted, the input shape is preserved. `hoh` is not a selectable output, since it would require choosing a key column.
- **`col()` reworked, not removed.** Both predicate forms are kept: `col('age') >= 18` still works and is the concise/composable option, while a coderef covers everything else. Internally `col()` is now **pure Perl** — an overloaded class that builds a per-row closure — and `filter` unwraps that closure so `col()` and a coderef share one evaluation path. The previous standalone XS predicate evaluator (`filt_eval`/`filt_ctx`) is gone; delete it if your tree still has it. One consequence: a `col()` comparison now costs the same per row as the equivalent coderef (a Perl call), rather than being evaluated in C.
- **Unchanged guarantees:** the input frame is never modified; `undef` (and, for numeric ops, non-numeric) cells never match a `col()` comparison; AoH/HoH rows are shared rather than copied where possible; keep-all/keep-none shapes are well defined per output type; Perl 5.10 compatibility is retained. A latent `SvTRUE(POPs)` double-evaluation in the per-row call helper (which crashed on perls where `SvTRUE` is a multi-eval macro) was fixed along the way.
[read_table]
- Added an opt-in `auto.row.names` argument so `read_table` can read the file R
produces by default from `write.table(x, sep="\t")`.
- The problem:
- R's `write.table` defaults to `row.names=TRUE, col.names=TRUE`, which writes the
row-names column in every data row but emits no header label for it. So a
frame with N columns comes out as N header fields over N+1 data fields — e.g.
`mtcars` gives 11 headers but 12-field rows. By default `read_table` (correctly)
rejects that as ragged:
Alignment error on mtcars.tsv data row 1 (12 fields vs 11 headers).
- The change:
- `auto.row.names` turns on R's own `read.table` rule: **when, and only when, the
header is exactly one field short of the data rows, treat the first field of
each row as an (unlabelled) row-names column.**
# default: the leading column is named 'row_name'
my $df = read_table('mtcars.tsv', 'auto.row.names' => 1);
# or give it a name
my $df = read_table('mtcars.tsv', 'auto.row.names' => 'model');
- The synthesized column behaves like any other first column: it appears in `aoh`
and `hoa` output, and for `hoh` it becomes the default key (so rows are keyed by
the model name). This also lines up with the existing handling of R's
`col.names=NA` output (a blank leading header), which still produces a
`row_name` column with no flag needed.
- What did not change:
- The strict alignment check is still the default. Without `auto.row.names` the
lopsided file still croaks, and even with it, a row that is off by anything
other than exactly one field still croaks — so the corruption guard only relaxes
for the one case R itself treats specially.
- Tested in `t/read_table.2.t` (16 assertions, Perl 5.10.1 and 5.38): aoh / hoa /
hoh output, custom column name, the already-aligned file (flag is a no-op), the
`col.names=NA` path, and the strict / ragged croak paths.
- additional bugfix:
# This is a comment
id,name,val
1,Alice,10.5
2,Bob,
3,Charlie,15.2
- would not be read correctly using `read_table`, but now is read correctly
[value_counts]
- now accepts array of hashes
0.16 2026-06-17 CDT
- changes to dist.ini, the minimum Perl version disappeared when I fixed other problems
- clarifications between run time and test dependencies
- addition of `csort` function to sort AoH and HoA
- addition of `aoh2hoa` to translate array of hashes into a hash of arrays
- fix of long double functions: https://www.cpantesters.org/cpan/report/5d5d9836-6a5f-11f1-aadb-63fd6d8775ea
[`glm`]
- output residual keys now use names, not integers
[`lm`]
[Bug fixes]
- Memory leak on the zero-degrees-of-freedom error path. When
`valid_n <= p`, the cleanup freed the `valid_row_names` *array* but not the
per-row name strings it held (those had been transferred out of `row_names`,
whose own array was already freed). The strings leaked on every such error.
Added the per-entry `Safefree` loop before freeing the array, matching the
normal path.
- HoH input validated only the first row. Only the first hash value was
checked to be a `HASHREF`; subsequent values were `SvRV`'d unconditionally, so
a malformed row (`{ a => {...}, b => 5 }`) dereferenced a non-reference. Every
row is now validated, with the partial allocations cleaned up before the
`croak`, mirroring the existing AoH path.
- `isspace` on a possibly-signed `char`. `isspace(*src)` is undefined for
byte values ≥ 0x80 on platforms where `char` is signed. Cast to
`(unsigned char)` before the call.
[Speed / RAM improvements]
- Formula buffer is now heap-allocated to fit. `char f_cpy[512]` silently
truncated any longer formula. Replaced with a buffer sized to
`strlen(formula) + 1`, so there is no fixed limit and no truncation.
- `.`-expansion buffer is now a growable heap buffer. `char rhs_expanded[2048]`
silently dropped expanded terms once full. It is now a buffer that doubles on
demand. Appends also went from `strcat` (which rescans from the start every
time — O(n²) over many columns) to an O(1) amortised append that tracks the
write position.
- No more per-row scratch allocation in matrix construction. The original
`safemalloc`'d a `row_x` buffer, filled it, copied it into `X`, and freed it
*for every row* — `n` allocations plus `n*p` copies. Each candidate row is now
written straight into `X` at its prospective commit slot; a row that fails
listwise deletion is simply overwritten by the next candidate. This removes the
`n` allocate/free cycles and the copy loop entirely.
- Categorical levels sorted with `qsort`. The level list used an O(n²) bubble
sort; replaced with `qsort` (relevant only for high-cardinality factors).
- Unused tail of `X` reclaimed after listwise deletion. `X` is allocated for
all `n` rows up front (`valid_n` is unknown until rows are scanned). When rows
are dropped, `X` is now `Renew`ed down to `valid_n * p`, returning the unused
tail to the allocator before the OLS phase.
- Minor robustness. The argument-parsing index was widened from
`unsigned short` to `I32` to match `items`, and the HoH row count now uses
`HvUSEDKEYS` rather than relying on `hv_iterinit`'s return value.
[Known limitations (left unchanged)]
- A multi-way term such as `a*b*c` is split only on the first `*`, so it yields
`a`, `b*c`, and `a:b*c` rather than a full three-way expansion. Deeper
interactions silently fail (the unparsable term evaluates to `NaN` and the
rows are dropped). This matches the documented two-way `*` support.
- HoA input takes the row count from the first column; columns shorter than
that simply contribute dropped rows rather than raising an error.
[`oneway_test`]
- Bug fixes:
- Memory leaks on error paths. Nearly every `croak` after an allocation
leaked memory. `croak` does a `longjmp`, so anything allocated but not yet
freed is lost. Affected paths:
- AoA and hash first-pass errors leaked `sizes` and any `gnames[]` entries
allocated so far.
- Formula-mode "not found as an array ref" errors leaked `lhs` and `rhs`.
- All post-allocation errors now route through a single `fail:` label that frees
every pointer unconditionally. Pointers are initialised to `NULL` and `gnames`
is zero-allocated with `Newxz`, so the cleanup is always safe to run.
- Undefined and non-numeric cells silently coerced to `0.0`. The original
second pass used `(svp && *svp) ? SvNV(*svp) : 0.0`, meaning an `undef` or
non-numeric cell was quietly treated as zero, silently corrupting the
F-statistic. Each cell is now validated with `SvOK` and `looks_like_number`;
the call dies naming the group and observation index, consistent with the rest
of `Stats::LikeR` (`mean`, `sum`, `cor`, etc.).
- Unsigned wraparound on empty array input. `k = (size_t)av_len(in_av) + 1`
cast to `size_t` *before* adding, so an empty array (`av_len` returns `-1`)
produced `SIZE_MAX` rather than `0`. Changed to
`k = (size_t)(av_len(in_av) + 1)` so the `+1` is done in signed arithmetic
before the cast.
- Unreliable group count from `hv_iterinit`. `hv_iterinit` returns the
number of buckets in use rather than the number of keys for tied hashes.
Replaced with `HvUSEDKEYS`, which always returns the correct key count.
- Improvements:
- `var.equal` accepted as an alias for `var_equal`. R users write
`var.equal`; the argument parser now accepts both spellings.
- Perl memory API used throughout. `safemalloc` and manual `memcpy` replaced
with `Newx`, `Newxz`, `savepv`, and `savepvn`. `savepvn` additionally
preserves embedded NUL bytes in group key strings, which the previous
`strlen`-based copies silently truncated.
- Known limitations (not changed):
- A factor column named `Residuals` or `group.stats` in a formula call will
collide with reserved top-level keys in the result hash.
- Group names containing an embedded NUL are stored correctly but are still
truncated at `strlen` when written into the output hash keys.
[`view`]
- default view shifted to 80 characters to match Linux window length
- New features:
- **`rows` is accepted as a synonym for `n`** (the number of rows shown).
Passing both `n` and `rows` is an error.
- **Unknown arguments are now rejected.** `view` validates its argument names
against the documented set (`n`, `rows`, `na`, `max_width`, `ellipsis`,
`gap`, `cols`, `columns`, `to`, `return_only`, `row.names`, `row_names`) and
dies listing any it does not recognise, so a misspelt option (e.g. `widht`)
is caught instead of silently ignored.
- **`n` / `rows` is validated.** It must be a non-negative integer; `undef` or
a non-numeric value now dies with a clear message instead of producing
warnings and being treated as `0`.
- **flat/simple hashes are accepted as input**
- Bug fixes:
- **`n => 0` now still prints the column header.** Column names were collected
only from the rows being shown, so requesting zero rows produced an empty
header line. At least one row is now scanned (when data exists) so the
header always lists the columns.
- **An empty hash (`{}`) no longer dies.** It was rejected as
*"neither ARRAY nor HASH"*; it is now shown as an empty table
(`0 rows x 0 cols`), matching the handling of an empty array.
- **The `row_names` alias now drives the Hash-of-Hashes label header.** The
header for the row-label column consulted only `row.names`, so
`row_names => 'id'` displayed `row_name` instead of `id`. Both spellings are
now honoured consistently.
- **Malformed nested values degrade gracefully.** A Hash-of-Arrays column or
Hash-of-Hashes row whose value is not actually an array/hash reference now
renders as empty cells rather than throwing a dereference error.
- Performance:
- Column gathering no longer sorts once per scanned row. Unique column names
are collected across the scanned rows and sorted a single time (same output
order), and the ellipsis length is computed once rather than per cell.
- Tests:
- `t/view.t` is self-contained (the `view` implementation is inlined; it loads
no other files) and covers the new argument handling, the bug fixes above,
and the existing AoH / HoA / HoH behaviour, alignment, truncation, and
output-path handling.
[`wilcox_test`]
- Corrected four bugs in the `wilcox_test` XSUB plus a portability fix in its exact signed-rank helper. Behaviour on valid input is unchanged: the R-agreement cases (unpaired `W = 58`, `p = 0.13292`; paired one-sided `V = 40`, `p = 0.019531`; separated exact `W = 0`, `p = 0.028571`) all still match R's `wilcox.test`.
- Bug fixes:
- **Invalid `alternative` is now rejected.** Any value other than `less` or `greater` previously fell through to the two-sided branch and returned a two-sided result mislabelled with the bad string, so a typo like `alternative => "twosided"` silently "worked". It now croaks unless `alternative` is one of `two.sided`, `less`, `greater`.
- **Zero/negative variance is guarded.** When every observation is tied the approximation's variance collapses to 0 and the old code divided by `sqrt(0)`: `wilcox_test([5,5,5], [5,5,5])` returned `p = 0` (a "significant" difference between identical samples). It now warns and returns `p = 1`.
- **Two-sided continuity correction at `z = 0`.** R uses `sign(z) * 0.5`, so the correction is `0` when the statistic sits exactly on its mean; the old code used `-0.5`. Example: `wilcox_test([1,4], [2,3], exact => 0)` changed from `p = 0.698535` to `p = 1` (matches R).
- **`exp` no longer shadows libm.** The local `exp` accumulator (mean of the statistic) shadowed the C library `exp()`; renamed to `mean_w` (two-sample) and `mean_v` (signed-rank). No active miscompute, removed as a latent hazard.
- Cosmetic:
- Collapsed a no-op ternary that assigned the same signed-rank exact method string on both branches; the `method` field is now simply `Wilcoxon signed rank exact test`.
- Portability (exact signed-rank helper):
- **`exact_psignrank` no longer calls `powl()`.** The `2^n` normaliser is now built by exact repeated doubling, which has no long-double libm dependency. This fixes an `Undefined symbol "powl"` load failure reported by a CPAN smoker (FreeBSD, perl 5.20, `nvtype=double`) whose libm lacks the long-double math functions; the symbol resolved on glibc, which is why local builds passed. `long double` accumulation in the DP is retained — only the `powl` call was at fault.
- **`int` → `size_t`** for `n`, `max_v`, and the DP loop counters, which also removes a `size_t`-to-`int` narrowing at the call site. The `floor()` result (`k`) stays signed so its negative-`q` sentinel still fires, and is cast to `size_t` only after the `k < 0` check.
- Tests:
- Added `t/wilcox_test.t` (flat, no subtests): R-agreement cases, option handling (`paired`, `correct`, `exact`, `mu`, named/positional `x`/`y`, NA dropping), regressions for all four bug fixes, argument-error and `alternative`-validation checks, output shape, and `no_leaks_ok` coverage of the two-sample, exact, and paired allocation paths.
0.15 2026-06-11 CDT
- `view` function added, similar to R's `head`
- `read_table`:
filter => {
'Testosterone, total (nmol/L)' => sub { defined $_ },
}
- was broken by the change in undefined variables in 0.14, but is back to being `undef`
- `col2col` improvement in sectioning in README
- Numerous changes to prevent quadmath/long double CPAN test failures
- Minimum Scalar::Util version in dist.ini is now 1.22, see https://www.cpantesters.org/cpan/report/6b682236-6567-11f1-a3bc-a055f9c4ba34
- `Digest::SHA` removed as a dependency
[`read_table`]
- Bug fixes:
- **A comment-prefixed header is now read correctly.** `read_table` strips a
leading comment marker from the header line (so a file may begin with
`#id,val`), but that strip was dead code: the XS parser skipped *every* line
beginning with the comment string before the callback ever saw it, so a
commented header was silently dropped and the first data row was mistaken for
the header. The parser now delivers the first content line even when it
begins with the comment marker, and only skips comment lines after the header
has been seen.
- **Carriage returns inside quoted fields are preserved.** The parser stripped
`\r` unconditionally, so a quoted value such as `"x\ry"` lost its carriage
return and would not survive a `write_table` -> `read_table` round-trip. `\r`
is now stripped only as part of a trailing CRLF line ending and as a stray CR
*outside* quotes; inside quotes it is literal data.
- **Duplicate column names no longer corrupt `hoa` output.** With
`output.type => 'hoa'`, a repeated column name pushed the same cell once per
occurrence, so the affected columns came out longer than the others and the
arrays no longer lined up by row. Columns are now keyed by unique header name
(first-seen order preserved, later values win, one warning emitted).
- **A defined non-CODE callback is now an error.** Passing a defined argument
that was not a CODE reference silently fell through to slurp mode and ignored
the argument; it now croaks
(*"callback must be a CODE reference"*).
- **An undefined/empty `hoh` row-name now dies instead of keying on `""`.**
With `output.type => 'hoh'`, a row whose row-name column was empty/undef was
stored under the `''` key and raised *"uninitialized value"* warnings. It now
dies, naming the column and the offending data row.
- **A numeric filter key past the last column now dies.** A 1-based numeric
filter key greater than the column count was accepted, then silently extended
every row through the `$_` write-back. It is now rejected up front with a
message naming the column count.
- **`sep` and `delim` together now die.** Supplying both silently preferred
`delim`; passing both is now an explicit error (`delim` remains an alias for
`sep` when used alone).
- **The library no longer prints to STDOUT.** The unknown-argument path used
`say` to dump the offending names to STDOUT before dying; the names are now
carried in the `die` message itself.
- Better diagnostics:
- Alignment errors now report **which data row** is ragged
(*"Alignment error on FILE data row N (X fields vs Y headers)"*), instead of
only the field/header counts.
- Memory-leak fixes (exception paths):
- The parser allocated its working buffers (`current_row`, `field`, and — in
slurp mode — `data`) in the XS `INIT:` block, i.e. *before* any validation, and
freed them only by falling off the end of the function. Any non-local exit
therefore leaked:
- the open-failure `croak` leaked the row buffer and field (and the slurp
accumulator);
- far more commonly, a `die` thrown **inside the row callback** — which
`read_table` does routinely on alignment errors, bad row names, and filter
exceptions — unwound straight out of the XS frame and leaked the field, the
current row, the line buffer, the slurp accumulator, *and the open file
handle*.
- Allocations now happen in `CODE:` after every croak-able check, and every
long-lived resource (the file handle via `SAVEDESTRUCTOR_X`, the buffers via
`SAVEFREESV`) is tied to the save stack, which an exception unwinds. Measured
with `Test::LeakTrace`: a `die` mid-file went from 5 leaked SVs to 0, and an
open failure from 2 to 0. This is the likely source of the constant-size leaks
seen in CPAN-tester reports for the exception-path tests.
- Performance:
- **~2.5x faster parsing** (57 -> 145 MB/s on a 100k-row quoted file). The core
loop appended one character at a time with `sv_catpvn(field, &ch, 1)`; it now
scans runs of ordinary bytes with `memchr` / a bounded scan and appends each
run in a single `sv_catpvn`, copying field contents in bulk rather than byte
by byte.
- Internal / non-behavioral:
- XS declarations moved from `INIT:` to `PREINIT:`; allocations deferred into
`CODE:` (see the leak fixes above).
- The filter loop now aliases the row hash with `local *_ = \%line_hash`
instead of copying it with `local %_ = %line_hash`. This removes a full
per-row hash copy for every filtered row and fixes a latent staleness bug:
after a filter mutated `$_` and the change was written back, `%_` still
reflected the pre-mutation copy, so a subsequent filter in the same row saw
stale values. With aliasing, `%_` *is* the row, so write-backs are always
visible.
- Known limitation (not changed):
- **`undef.val` does not round-trip back to `undef`.** `write_table` renders an
`undef` cell as an empty field by default, and `read_table` maps an empty
field back to `undef`, so the *default* round-trip is clean. But if a file is
written with a token such as `'undef.val' => 'NA'`, `read_table` has no
inverse option and reads `NA` back as the string `'NA'`. `read_table` also
cannot distinguish a deliberately quoted empty string (`""`) from a missing
value -- both become `undef`. Adding an `na.strings`-style option to
`read_table` (mapping configurable tokens and/or empty fields to `undef`)
would close this gap.
[`write_table`]
- Behavior change:
- **`undef` cells now write as an empty field, not an empty string.** A missing
or `undef` value renders as nothing between separators (`a,,c`) rather than a
quoted empty string (`a,'',c` / `a,"",c`). Supplying `'undef.val' => 'NA'`
(or any other token) still overrides this, exactly as before. This is the
only change that can alter the bytes of an existing output file; if you relied
on the previous default, pass `'undef.val' => ''` to keep an explicit empty
field, or your chosen placeholder.
- Bug fixes:
- **Wide-character / UTF-8 column names and row keys now round-trip.**
Previously, cells were looked up with the raw bytes of the column name
(`hv_fetch(..., SvPV_nolen(name), strlen(name), ...)`), which fails to match a
UTF-8-flagged hash key: the column header printed correctly but every cell
under it came back empty. All lookups now fetch by SV (`hv_fetch_ent`), header
lists are gathered and sorted as SVs (`sortsv` + `sv_cmp`, preserving the
flag) instead of being round-tripped through `char *`, and the `row.names`
column is matched with `sv_eq` rather than `strcmp`. Embedded NUL bytes in
keys are handled correctly as a side effect.
- **`col.names => []` no longer loops forever.** An empty `col.names` array made
`av_len()` return `-1`, which — compared against an unsigned `size_t` loop
index — wrapped to `SIZE_MAX` and ran effectively without end. This was fixed
for flat hashes previously; it was still present for hash-of-hashes,
hash-of-arrays, and array-of-hashes, plus both `row.names` header-filtering
loops. All such loops now use a signed index.
- **Tables wider than 65,535 columns no longer hang.** One header loop used an
`unsigned short` index that silently wrapped past 65,535 and never terminated.
It now uses `size_t` like the rest of the code.
- **Flat-hash cells holding a reference now croak.** Every other input shape
rejects a nested reference with
*"Cannot write nested reference types to table"*; a flat hash instead
stringified it (e.g. `ARRAY(0x55...)`) into the file. It now croaks
consistently.
- **`'undef.val' => undef` is handled cleanly.** It previously called
`SvPV_nolen` on `undef`, raising an *"uninitialized value"* warning and
yielding an empty string by accident. It is now treated explicitly as an empty
field, with no warning.
- Memory-leak fixes (exception paths):
- The row-key list gathered for hash-of-hashes input was leaked when the output
file could not be opened.
- The *"Could not get headers"* croak on hash-of-arrays input leaked both the
already-open filehandle and the headers array.
- Internal / non-behavioral:
- Numeric row labels are now formatted into a reused stack buffer instead of a
per-row `savepv()` / `safefree()` allocation (no functional change; removes a
cast-away-`const` and one allocation per row).
- Several signed/unsigned index types were made consistent (`SSize_t` vs
`size_t`) to match `av_len()` and silence the conditions behind the loop bugs
above.
- Tests:
- `t/write_table.t` expanded from 17 to 69 assertions. New coverage targets each
fix above: the empty-field default and `undef.val => undef` (no warning),
`col.names => []` termination across all four input shapes, the
>65,535-column header loop (gated behind `EXTENDED_TESTING=1`), in-sequence
numeric row labels, nested-reference rejection, CSV quoting corners
(carriage return, separators inside column names, multi-character separators),
empty input writing no file, and UTF-8 column names and row keys. Two leak
assertions cover the exception paths above.
0.14 2026-06-08 CDT
- `filter` function added for rows
- `read_table` reads undefined values to `undef` instead of `NA`, which makes calculations easier
- `write_table` writes undef by default as an empty string `''`
- `hoh2hoa` transforms a hash of hashes into an hash of arrays
- `quantile` uses `NV` instead of `double` to allow for high-precision 128-bit floats to be used on quadmath machines when available: https://www.cpantesters.org/cpan/report/296f4868-631f-11f1-abba-ff15558d240b
- Numerous switches from `double` to `NV` for local precision, like above
- numerous changes to `col2col` for ease of use and working with datasets with numerous undefined values
- dist.ini now links to math library when compiling: https://www.cpantesters.org/cpan/report/785e26d8-6397-11f1-89c0-dc066e8775ea
- `fisher_test` now should be complete, errors with confidence intervals fixed
0.13 2026-06-07 CDT
- `read_table`: speed improvements; commented headers are now allowed
- `write_table`: fix for
Attempt to free temp prematurely: SV 0x56417a2ae610 at t/write_table.t line 182.
main::wrote_ok(",age\x{a}Alice,30\x{a}Bob,25\x{a}", "row.names => 'name' uses that column as labels", HASH(0x56417a272250), "row.names", "name") called at t/write_table.t line 203
Attempt to free unreferenced scalar: SV 0x56417a2ae610 at t/write_table.t line 183.
main::wrote_ok(",age\x{a}Alice,30\x{a}Bob,25\x{a}", "row.names => 'name' uses that column as labels", HASH(0x56417a272250), "row.names", "name") called at t/write_table.t line 203
- `write_table` gives better warnings for incorrect types of data given
- Numerous changes to dist.ini to improve CPAN testing, especially for Win32
0.12 2026-06-08 CDT
- `add_data` can also take hash of arrays, and various mixes of data types
- `ljoin`: Addition of `restrict` keywords in many places; should improve CPU performance
- Better POD formatting, correction of output hash for README's `add_data`
- `chisq_test` can now accept hash of hashes as input
- new `transpose` function for switching 2D hash keys and 2D array indices, and `col2col` for comparing columns against columns
- removed unused function from C helpers
- `value_counts`: addition of restrict keywords in preinit, should improve CPU performance
- MANIFEST.skip changed to MANIFEST.SKIP to improve CPAN testing
- using `is_deeply` for tests of `transpose`, which may or may not work with CPAN testers (experimental)
- Added function name to warnings, so I actually know which function is producing the error
- `write_table` can also take `file` and `data` as args, in addition to positions
- fixed `write_table` as it could hang if given empty `col.names` or `row.names`
- Added `__EXTENSIONS__` to source XS file for better CPAN testing
0.11 2026-06-03 CDT
- better POD formatting for tables
- addition of MANIFEST.skip to get better testing results on CPAN
- `glm`: bugfix for when there is no intercept in the formula, new test cases in t/glm.t
- `write_table` now accepts simple hashes as input, in addition to hash of arrays, hash of hashes, and arrays of hashes
- Better documentation for t-test
0.10 2026-06-01 CDT (approx)
- changes to compilation for CPAN, trying to get this work on Windows
- Addition of `prcomp` and `value_counts`
- `matrix` will work without key names, just like in R. Testing for `matrix` has improved.
0.09 2026-06-01 CDT (approx)
- context changes in XS `dTHX`, `pTHX_`, and `aTHX_` to get better CPAN testing results
- `restrict` keywords added to `lm` to increase speed
0.08 2026-05-26 CDT
- Speed improvement in `summary` of hashes.
- Addition of `add_data`, `dnorm`, `group_by`, `ljoin`, and `mode` functions
- Chi-squared function no longer has Perl wrapper, and all code is in XS, which should result in a minor speed increase with 1 less function call.
- Compiler changes for GNU source and inclusion of `strings.h`, to ensure more CPAN testing works better.
- `read_table` now returns hash-of-hash in {row}{column}
0.07 2026-05-24 CDT
- Addition of `summary` function.
- Formulas can now be omitted from `aov`, resulting in a stacked calculation as R would think.
- Addition of `oneway_test` for multi-group comparisons that does not assume normality like `aov` does.
- `read_table` and `write_table` now automatically set separators for `.csv` files as `,` and `.tsv` files as `"\t"`, respectively, so these values no longer need to be specified separately from the file name.
0.06 2026-05-19 CDT
- Changed compiler options so that Solaris will work
- signed integers changed to unsigned in `glm`
- Added restrict keywords to `power_t_test`, and made `int` to `unsigned int`
0.05 2026-05-08 CDT
- Leak testing for `sample`
- removal of Data::Printer dependency for easier CPAN testing
- switched several `unsigned int` variable to `I32` so that clang doesn't complain
- added restrict keyword for `sample`
0.04 2026-5-17 CDT
- addition of `sample` function
- GNU source, to maximize compatibility and ease installation
- removal of JSON dependency to ease installation
0.03 2026-5-13 CDT
- Compatibility back to Perl 5.10
0.02 2026-5-7 CDT
- back-compatible to Perl 5.10, instead of original 5.40, ensuring more people can use it
- added var_test
- mean, min, sum, median, var, and max die with undefined values, and print the offending indices
- "group.stats" added to aov, for TukeyHSD in the future
- "cor" dies when given data with standard deviation of 0
- `write_table` now has `undef.val` option, which shows how undefined values are printed to tables, which is `NA` by default.