/ 
Stats-LikeR-0.316

Changes for version 0.316 - 2026-09-11

  • read_table on .xlsx: the worksheet parser moved from perl to XS
    • Reading an .xlsx was slow out of all proportion to the format. On a 21,845 x 50 workbook (7.8 MB on disk, 36 MB of worksheet XML, 1,013,220 cells, an 8 MB shared-string table of 117,870 entries) `read_table` took 2.47 s and peaked at 263 MB. The same table written out as a 34.7 MB CSV and read back through `_parse_csv_file` took 0.139 s and 156 MB -- from a larger file, into the same array of hashes. Almost the whole gap was the perl side, not the format.
    • Where it went: the nested regexes in `_parse_xlsx_sheet` were 1.68 s of the 2.47 s, or 1.66 us per cell, against 0.099 s for a bare `$ws =~ m{<c\b}g` count of the same string -- so the scan was 6% of it and perl's per-op overhead the rest. Decompressing the worksheet was 0.25 s, the shared strings 0.21 s, and read_table's own per-row closure 0.32 s.
    • The parse is now `xlsx_ws_scan()` in `LikeR.xs`, and `read_table` offers it the same `csv_plan` fast path a delimited file gets, so an `aoh` or `hoa` is assembled in C from the cells the parser already has. The perl callback still reads the header, and still takes every row when a `filter` or a `hoh` shape needs per-row perl, exactly as it does for a CSV.
    • Rows are emitted as they are parsed. They have to be padded to the widest row in the sheet, which the perl parser could not know until it had read the last one, so it buffered the whole sheet as an array of arrays first
      • 76 MB on that workbook, before a single row had been handed on. A first pass now reads nothing but the cell references and the width is known before any row is built. Both passes are the same function so that they cannot see different cells.
    • Shared-string cells are copy-on-write copies of the table entry rather than fresh buffers, which is worth 29 MB there. That needs asking for: `sv.h` defines `SV_DO_COW_SVSETSV` -- what `sv_setsv()` and `newSVsv()` pass -- to the real flags only under `PERL_CORE`, and to 0 otherwise, so an XS copy is an ordinary copy unless it names `SV_COW_SHARED_HASH_KEYS` itself. Nothing here writes through `SvPVX` of a cell, so it is safe to. Perl 5.20 is where copy-on-write started covering plain strings, so on 5.10 and 5.12 this is an ordinary copy as before.
    • `_unzip_member` reads onto the end of its buffer instead of concatenating through a second scalar: `IO::Uncompress::Base::read()` makes its truncating `substr()` a no-op when the offset is already the buffer's length. 0.171 s against 0.243 s on the 36 MB part, and no realloc slack left behind. (`BlockSize => 1<<20` is not the way to do this: 2.08 s.)
    • `xl/sharedStrings.xml` is parsed in XS too, and `_xlsx_col_idx` is now the `xlsx_ref_col()` that places the cells, exposed to perl so that `t/xlsx_col_idx.t` still exercises the letter arithmetic on its own.
    • Together: 2.47 s -> 0.48 s and 263 MB -> 221 MB on that workbook, and 0.083 s -> 0.031 s on a 2,207-row one that uses inline strings and no shared-string table. What is left of the peak is the table the caller asked for, plus the 36 MB worksheet part and the 33 MB shared-string table, which are both live for the length of the parse.
  • What the new parser answers differently
    • Two things, both on input the format does not allow. A cell reference past `XFD` -- the last of the 16,384 columns ECMA-376 gives a worksheet -- is now read as no reference at all, which puts the cell in the next column, the same place a cell with no `r=` goes. Every row is padded to the widest column the sheet mentions, so the reference is what decides what a row costs in memory, and a bad one should cost the 16,384 the format allows rather than the twelve million `ZZZZZ` asks for or the three hundred million of `ZZZZZZ`. The perl parser placed the cell wherever the arithmetic landed.
    • A numeric character reference above `&#x7FFFFFFF;` is now left in the text instead of being decoded. The perl version handed the number to `chr()` unguarded, and that was not one answer but two: `&#999999999999;` came back as thirteen bytes of perl's extended UTF-8 on an ivsize=8 build, and died outright on `5.44.0-i686` with "Use of code point 0xFFFFFFFF is not allowed". XML 1.0 does not allow a character reference above `#x10FFFF` at all, so nothing legal is lost, and the answer is now the same on every perl in the matrix.
    • Everything else is unchanged, and was checked rather than assumed: the two parsers were run side by side over `Affinity Dataset(main).xlsx` and `titanic.xlsx` in all three `output.type` shapes with and without a `filter` and `na.strings`, over fifteen hand-built worksheets covering the awkward layouts, and over 400 generated strings of entity soup. Every `Data::Dumper` of the result is byte-identical, warnings and die messages included.
  • How the new parser was checked
    • A tokenizer is only as good as the malformed input it survives, so it was fuzzed rather than reasoned about: every one of the 362 prefixes of a worksheet part, 20,000 byte-level mutations of it, and 12,000 generated workbooks each read twice -- once through the XS fast path and once through the perl callback -- so that the two could be required to agree cell for cell. They do, on all of them. Three defects in the new code turned up that way and are fixed:
    • `av_store()` already releases whatever was in the slot, so the parser's own `SvREFCNT_dec` for a repeated `r=` in one row was a second free of the same SV. Only malformed input has two cells claiming one column, which is why nothing but a fuzz would have found it.
    • The two passes could see different cells. The second skipped each cell body to `</c>`; the first had no reason to and did not. On a cell whose `</c>` is missing that search runs on to the next cell's, so the pass that skipped saw fewer cells, and rows were then built to a width measured from a different reading of the same file. Both passes now skip -- the scan is the same scan or it is not the same file.
    • The alignment croak that followed segfaulted. `S_fast_row()` frees the row buffer itself before it croaks, because it unwinds past `_parse_csv_file`'s local, and the new parser's save-stack destructor then freed it again. It now hands the buffer over and takes it back only if the call returns.
    • A fourth was found by the new test on `perl-5.10.1` and `perl-5.12.5`, and by nothing else: a cell is placed by its column reference, so a row with a gap leaves the slots between untouched, and perl 5.10 and 5.12 fill a newly extended array with `&PL_sv_undef` where 5.14 and later zero it. `SvCUR()` on the immortal undef dereferences a NULL `SvANY`. Only the first row after a fresh buffer can meet one -- every later row is cleared to NULLs -- and `t/read_table.xlsx.t`'s sparse cells are not in its first data row, so it passed on both perls while the new file segfaulted. The hole test now accepts either value.
  • read_table on .xlsx: a third off the time, straight to Compress::Raw::Zlib
    • With the worksheet parser in XS, decompression is what a read of an .xlsx spends its time on: 0.239 s of the 0.413 s the 21,845 x 50 workbook takes, against 0.015 s to parse its 117,870 shared strings. `_unzip_member_fast()` reads the archive's central directory itself and inflates one member with `Compress::Raw::Zlib`. The 36 MB worksheet part goes from 0.174 s to 0.067 s and the 8 MB shared-string part from 0.052 s to 0.023 s. Reading the whole workbook: 0.481 s to 0.323 s, with peak RSS unchanged at 179 MB.
    • Most of that gap is one thing, and it is not the inflate. `Unzip.pm`'s `ckParams()` sets `crc32 => 1` unconditionally ("unzip always needs crc32"), so every byte goes through `Compress::Raw::Zlib::crc32()` -- 0.076 s on that part, more than the 0.052 s the inflate itself costs -- and the comparison against the stored CRC happens only under `Strict`, which defaults to 0. The check is paid for and never made. This path does not compute it either, which is the same answer for the same money. Asking `Inflate` for `-CRC32 => 1` costs the identical 0.076 s, being the same per-chunk call, so turning it on here would only be worth anything with a croak on mismatch behind it -- and that is a decision about what `read_table` should do with a damaged workbook, not a speed one.
    • `-Bufsize => 1<<16`: 0.052 s against 0.063 s at the 4 KB default on that part, and flat from there (1<<18, 1<<20 and 1<<22 all measure 0.052 s).
    • What it declines, handing the archive back to `IO::Uncompress::Unzip` unchanged: zip64 (whether by the end-of-central-directory sentinels or a member's own), split archives, encrypted members, a compression method that is neither stored nor deflate, a directory that does not check out against its own signatures and lengths, a local header that is not one, and an inflate that does not end where the directory says it should. Declining is not an error -- it reads the central directory and one deflate stream and leaves everything past that to the module that has been ported to all of it.
    • Two things worth knowing about the shape of an archive. `IO::Compress::Zip` with `Append => 1` does not extend one: it writes a second complete archive after the first, and the final end-of-central-directory record describes only that second one, with offsets counted from where it begins. It reads as a single archive only to a reader that scans local headers forwards, which is what `IO::Uncompress::Unzip` does -- so `t/read_table.xlsx.t` and `t/read_table.xlsx.parser.t`, whose fixtures are built that way, go on exercising the fallback. Real workbooks do not have this: Excel, LibreOffice, openpyxl and this module's own `write_table` all take the fast path, which `t/unzip_member.t` pins by reading one that `write_table` wrote. Second, `IO::Uncompress` defaults to `Transparent => 1`, so a file that is not an archive at all comes back as its own raw bytes rather than failing. That has always been true of `_unzip_member`; the fast path declines such a file and changes none of it.
    • `Compress::Raw::Zlib` joins the prerequisites. It is core since 5.8.8 and `IO::Uncompress::Unzip` was already loading it, but it is dual-life and it is now named directly.
    • It also reads members the fallback cannot. `IO::Uncompress::Unzip` 2.020 (perl 5.10.1) and 2.024 (5.12.5) refuse a stored member written in streaming mode -- "Header Error: Streamed Stored content not supported" -- and fail to find any member past the first in a streamed multi-member archive at all, both of which this path reads on every perl in the matrix. That is why the new tests assert an accepted member against the bytes that went into the fixture rather than against what the fallback makes of it, and cross-check the fallback only where it can answer.
    • `t/unzip_member.t` is new, 82 tests: deflated, stored, empty and streamed members, and one whose sizes are in its local header rather than a trailing descriptor; a name another name begins with and one another ends with; an archive comment, and a comment holding the end-of-central-directory signature as a decoy for the backwards scan; an absent member (answered "handled, no such member", so that the sharedStrings.xml plenty of workbooks lack does not pay for a second scan); zip64, a concatenated archive, a file that is not an archive, one that does not exist, an empty one and a truncated one -- each asserted to be declined and then to come back from the fallback exactly as it did before, undef included; and read_table over two workbooks whose directories are intact, one from `write_table` and one with a shared-string table.
  • A malformed .xlsx could ask for an unbounded amount of memory
    • `xlsx_ref_col()` refuses a column reference past XFD, the last of the 16,384 columns ECMA-376 allows, and its comment says why: every row in a sheet is padded to the widest column the sheet mentions, so "a bad one should cost the 16,384 the format allows rather than the twelve million `ZZZZZ` would ask for". It did not. A reference it refuses falls back to "the next column", and that counter had no ceiling at all, so `ZZZZZ1` repeated 20,000 times in one row got 20,000 columns -- the exact outcome the cap was written to prevent, reached through the cap's own fallback.
    • What it cost: a 54 KB workbook of 20,000 such cells, plus 200 ordinary one-cell rows, came back as 264 MB of empty strings in 0.41 s; 60,000 cells made it 804 MB and 1.35 s, and nothing bounded it but the size of the input. The amplification is about 5,000x, so a 16 MB file of that shape would ask for tens of gigabytes. A cell with no `r=` at all reaches the same counter, so it did not even need a malformed reference.
    • `xlsx_ws_scan()` now clamps the column index to XLSX_MAX_COL on that path too. Past the ceiling the cells pile up in the last column, last one winning, which is what a repeated `r=` in one row already did. The same line runs on both passes, so they go on agreeing about the width. The two cases above now stop at 16,384 columns and 234 MB, and stop growing: 60,000 cells cost what 20,000 do.
    • `t/read_table.xlsx.parser.t` covers both routes to the counter -- 20,000 unreadable references and 20,000 cells with no `r=` -- asserted on `_parse_xlsx_sheet_xs` directly, because `read_table` would fold the unnamed columns into one key and hide the width. Both fail without the clamp.
  • read_table on .xlsx: 42 MB less peak RSS
    • `_unzip_member()` returned the decompressed part as a string, which costs a full copy of it: perl cannot hand back a lexical's pad slot, so `return $content` copies. On the 36 MB worksheet part of the 21,845 x 50 workbook that was 34 MB of peak RSS for nothing -- 82.4 MB against 48.5 MB for a scalar reference, and 81 MB still resident after the string went out of scope against 13 MB, the rest being heap the allocator never gave back.
    • It returns a reference now and its four callers dereference. `$$ws` on an argument list pushes the SV itself, so the part still reaches `_parse_xlsx_sheet_xs` and `_xlsx_sst_xs` without a copy. Reading the whole workbook went from 220.2 MB and 0.487 s to 177.9 MB and 0.475 s.
  • The `autodie` dependency is gone
    • `autodie` was a prerequisite for five calls in one file: three `close`s and two `open`s in `lib/Stats/LikeR.pm`, all of them reading this module's own POD or peeking at the first line of a CSV. Nothing in `LikeR.xs` used it -- the XS side opens through `PerlIO_open` and croaks for itself -- and no test loaded it. `open`/`close` now check for themselves through `_open_read()` and `_close()`, and the pragma and the prereq are both dropped.
    • Those two helpers raise the failure at the same points, with the message `autodie::exception` 2.37 would have built: `_format_open` (through `_FORMAT_OPEN` and `_format_open_with_mode`) for the open, `_format_close` for the close, each followed by `add_file_and_line`, which is why `_io_die()` takes the file and line one frame up rather than its own. The text is byte-identical, checked against a script running under `autodie`. The one visible difference is that `$@` now holds a plain string instead of an `autodie::exception` object, which nothing here ever inspected.
    • Two error paths the pragma had been masking are now visible, and both were dead code: `_pod_open`'s `open ... or return undef` and `read_table`'s `open ... or die "read_table: can't open $file: $!"` could never run, because `autodie`'s `open` died first. They are removed rather than revived, so a failed open stays fatal and keeps the wording it had.
  • Prerequisites: one added, three moved to the test phase
    • `IO::Uncompress::Unzip` is now declared. `_unzip_member` has required it for every `.xlsx` read since 0.24 and it was never in the list; it is core since 5.9.4, so the 5.010 floor already guaranteed it, but it is dual-life and a used module belongs in the metadata. `Cwd` stays for the same reason -- `provenance_path()` in `LikeR.xs` reaches it through `load_module()`, which croaks rather than degrades if the require fails, so the tex and xlsx writers genuinely need it.
    • `Test::Exception`, `Test::LeakTrace` and `Test::More` move from `[Prereqs]` to `[Prereqs / TestRequires]`. They had been runtime requires since 0.18, which made every user of `mean()` install two modules that are not core and that nothing under `lib/` or `LikeR.xs` loads.
    • They were flattened into the runtime phase by 0.18's "fix to dist.ini for dependencies", but the smoker failure that fix was chasing was not theirs. 0.17 had put `Devel::Confess` in `[Prereqs / DevelopRequires]`, a phase no CPAN client installs, while `t/col2col.t` and `t/transpose.t` still said `use Devel::Confess` -- so the test files would not compile on a smoker. 0.18 cured it by moving every phase back into runtime, which swept the three `Test::` modules along; 0.22 then re-declared `Devel::Confess` as a runtime dependency for the same reason, and 0.24 fixed it properly by dropping the module from the tests. Nothing under `t/` has loaded it since -- both files mention it only in comments -- and every other module the tests use is either core at 5.010 or declared.
    • Checked, not assumed. `dzil build` puts the three under `test.requires` in `META.json`; `ExtUtils::MakeMaker` carries them as `TEST_REQUIRES`; the generated `Makefile.PL`'s `%FallbackPrereqs` block folds them back into `PREREQ_PM` on any `ExtUtils::MakeMaker` older than 6.63_03, which was confirmed by running the configure step with `$ExtUtils::MakeMaker::VERSION` forced to 6.48 -- an old toolchain sees exactly what 0.316 gave it. `cpanm --showdeps` on the built tarball lists all three. The built dist configures, compiles and passes its 154 files clean.
    • The repo's own `Makefile.PL` is excluded from the dist and so is never regenerated by dzil; its prerequisite list still named `Devel::Confess` and `Digest::SHA`, neither of which has been used for many releases. It now mirrors `dist.ini`.
  • Tests
    • `t/read_table.xlsx.parser.t` is new: the tokenizer's own cases, which a well-formed workbook never reaches. Cells with no `r=` and with `r=` not first, self-closing `<c/>` and `<row/>`, blank rows, mid-row gaps, every entity form and the ones that must not be decoded, a shared-string index that is out of range or not a number, `t="str"`/`t="b"`/`t="e"`, a `>` inside an attribute value, `XFD` and the references past it, a repeated reference in one row, a cell with no `</c>`, an empty sheet, a header with no data rows, and the fast path and the callback path agreeing cell for cell. Leak checks on every path through the row emitter, with no `qr//` inside a measured block.
    • Those leak checks measure the parser with the worksheet part already decompressed. `read_table` itself is measured too, but only where decompressing is clean: `perl-5.10.1`'s bundled `Compress::Raw::Zlib` leaks 18 SVs per read in `crc32(undef)` (`IO/Uncompress/Unzip.pm` line 608), and a check around `read_table` counts those as this module's. The test asks this perl whether it leaks rather than naming versions, and skips only the four checks that need a file.
    • The whole suite passes on all seven perls in the local matrix -- 5.10.1, 5.12.5 (`long double`), 5.42.3-thr, 5.44.0, 5.44.0+x87, 5.44.0-quadmath (`__float128`) and 5.44.0-i686 (`ivsize=4`) -- and the generated `.c` compiles clean under strict `-std=c99`. Valgrind reports no errors over the truncation and mutation corpora.
  • aov() could abort the interpreter, or answer differently on every run
    • `aov` built the two halves of an `a:b` interaction in two 256-byte stack arrays, and filled the left one with `strncpy(left, term, colon - term)`. `strncpy` writes exactly the count it is given and knows nothing about the destination, so an interaction whose left component ran past 255 characters wrote off the end of the frame, and the `left[colon - term] = '\0'` after it stored past the end as well. glibc caught it as `*** buffer overflow detected ***` and aborted the interpreter, which no `eval` can catch. Both halves are now copied to the heap at the length they actually have. (`right` had been moved off `strcpy` already, but only as far as a truncating `snprintf`.)
    • Only the FIRST value of a hash-of-hashes was checked for being a reference. `SvRV()` on a later plain scalar reads a pointer out of a field that does not hold one, so `aov({r1 => {...}, bad => 42}, 'y ~ g')` segfaulted -- or did not, depending on hash order, since a non-reference that came first was rejected by the shape test. `lm` and `glm` have always checked this in `lm_read_rows`; `aov` now makes the same check with the same message.
    • The row count for a hash-of-arrays came from whichever column `hv_iternext()` returned first, so on a ragged frame which observations were fitted moved with perl's hash order from one run to the next. Refused now, naming the column, exactly as `lm_read_rows` refuses it.
    • `group.stats` had the same defect on the reporting side, and reached it by the documented no-formula form -- R's `stack()` -- whose columns are unequal by construction. `aov({short => [1..3], long => [1..20]})` reported `long`'s mean as 10.5 and its n as 20 on some runs and as 2 and 3 on others, on the same data in the same process image. Each column is now summarised over its own length, which is also what stops the evaluation running past the end of the short ones.
    • The groups of that no-formula form are stacked in sorted key order rather than hash order. The F statistic survived the shuffling, being invariant to the order of the rows, but only to within rounding -- `Pr(>F)` came back 0.0363396692989842 on one run and ...43 on the next -- and `fitted.values` named entirely different rows each time.
  • The `.` in a formula expanded in hash order
    • `get_all_columns()` returned the column names in hash-iteration order, which perl randomises per process and per hash. `.` is what fixes the term list, and a sequential (Type I) sum of squares is attributed in term order, so `y ~ .` produced a different ANOVA table on every run of the same script over the same data: over four runs of `aov({y, g, x, z}, 'y ~ .')` the `z` row came back with a sum of squares of 15 once and of 1.9e-30 another time, because `z` had been entered before or after the term it is orthogonal to. `lm` and `glm` take their `.` expansion through the same function, so their coefficient order moved the same way.
    • R has a column order to expand in and a Perl hash has none, so the names are sorted. That is the only order available that is the same twice.
  • aov()'s formula was truncated, and read its intercept markers with strstr()
    • The formula was copied into a 512-byte stack array and stopped at 511 characters, so a longer one was silently truncated and a different model was fitted than the one asked for. `.` expanded into a 2048-byte array and simply DROPPED every column that no longer fit, so `y ~ .` on a wide frame quietly fitted a smaller model.
    • `-1`, `+0`, `+1` and a leading `1+` were removed with `strstr()` over the whole right-hand side, so the `-1` inside `I(x-1)` was eaten and `y ~ I(x)`
      • a different model -- was fitted and reported under that name.
    • All of it now goes through `lm_formula_split()`, which is what `lm` and `glm` already parse with: it steps over `I(...)`, grows with the formula, and its buffer is on the save stack so a croak anywhere below releases it. The `.` expansion grows through `lm_append()`.
  • interpolate()'s splines solved a banded system densely
    • `ip_build_cubic()` and `ip_build_quad()` assembled an n x n dense matrix and ran a dense Gaussian elimination over it, where n is the number of numeric anchors in the column. Neither system is dense: the not-a-knot cubic spline is tridiagonal apart from its two boundary rows, and the degree-2 B-spline collocation matrix has three non-zero basis functions per row. Both are `kl = ku = 2`.
    • `ip_solve_band()` solves them in band storage, which is 7 NVs a row. A column of 12,800 anchors went from 0.835 s and 574 MB to 0.002 s and 13.6 MB; one of 400,000 -- which wanted 1.28 TB and could only ever fail -- now takes 0.071 s and 115 MB. `ip_eval_quad()` sums only the basis functions whose support reaches the point, rather than all n of them, so filling g gaps is O(g) and not O(n*g).
    • The band width is asserted rather than assumed: the degree-2 basis is a partition of unity, so every row of the collocation matrix sums to 1, and a window that had missed a non-zero basis function could not. Checked against SciPy's `CubicSpline(bc_type='not-a-knot')` and `make_interp_spline(k=2)` at 200, 2,000 and 20,000 anchors; worst relative disagreement 2.3e-16.
  • rbinom() cost O(size) per variate
    • `generate_binomial()` was the textbook Bernoulli loop: `size` draws from `Drand01()` per variate, counting successes. Exact, and unusable at any interesting size -- `rbinom(n => 10, size => 1e9)` asks it for ten billion uniforms, and `rbinom(n => 1e4, size => 1e5)` for a billion. 2000 variates at `size => 1e7` took 74 s.
    • Replaced by BTPE (Kachitvichyanukul and Schmeiser 1988, CACM 31, 216-222), transcribed from R 4.6.1 `src/nmath/rbinom.c`: the inverse-CDF walk below `n*p = 30` and the triangle / parallelogram / exponential-tail rejection above it. Neither draws a number of uniforms that grows with `size`. The same 2000 variates now take 0.0001 s, at any size.
    • Two departures from upstream. R caches the setup in file-static globals and its own comment there reads "FIXME: These should become THREAD_specific globals"; every variate of one `rbinom` call shares its n and p, so the setup is computed once per call into a caller-owned struct instead -- the same saving, no statics. And `unif_rand()` is `Drand01()`, which is what every other draw in this file uses and what makes `srand($seed)` govern the result.
    • THE SEEDED STREAM MOVED. BTPE consumes a different number of uniforms per variate than the Bernoulli loop did, so a given `srand` produces different numbers from those 0.315 gave. The distribution is unchanged and a run is still reproducible; only a script that hardcoded the values one seed used to give will see new ones.
  • ks_test(exact => 1) had no cost cap on the one-sample branch
    • `K2x()` builds an m x m matrix, `m = 2*floor(n*D) + 1`, and raises it to the n-th power: O(m^3 log n) time and O(m^2) memory, with nothing but the statistic bounding m. The default route cannot reach a large m, taking the exact branch only below n = 100, but `exact => 1` could. A sample of 800 whose D is 1 -- any badly-fitting reference distribution -- took 52 s and 69 MB, 3200 would have taken most of an hour, and past n ~ 23,000 the cell count overflowed the `int` it was computed in and went to `calloc()` wrapped. The two-sample branch has had `KS_EXACT_MAX_PRODUCT` for this all along.
    • m is now capped at `KS_EXACT_MAX_M` (500, about half a second and two megabytes) with the same warning and the same fall back to the asymptotic p-value the two-sample branch gives, and `K2x()`/`m_power()` compute in `size_t`.
  • dnorm() was computed at a double's width on whatever perl it ran on
    • `c_dnorm()` decides where the density has underflowed from the exponent range of the floating-point type, and asked `<float.h>` about a *double* -- `DBL_MAX`, `DBL_MIN_EXP`, `DBL_MANT_DIG` -- whatever perl's NV was. On a long-double or `__float128` build that cut the tail off at |x| ~ 38.57, where a double's subnormals run out, and returned a flat 0 beyond it: `dnorm(-100)` is 1.4e-2174, four thousand orders of magnitude inside a quadmath NV's range, and came back 0. `NV_MAX` / `NV_MIN_EXP` / `NV_MANT_DIG` now, which are the same constants on a double build.
    • `NV_MIN_EXP` needed a fallback: perl.h has defined it since 5.22, and both `perl-5.10.1` and `perl-5.12.5` are in the matrix. `ppport.h` does not backport a macro that is not an API function, so `LikeR.xs` derives it the way perl.h does -- from the same `<float.h>` constants and the same `USE_QUADMATH` / `USE_LONG_DOUBLE` tests the `nv_*` libm layer already switches on -- under an `#ifndef` a perl that has it never reaches. `NV_MANT_DIG` is guarded beside it so that a build with one and not the other cannot fail pointing at the wrong line. Caught by `./test.all.perls.pl`, which is what it is for: both builds failed at `make`, and nothing on the default perl would ever have shown it.
  • Convergence thresholds written against a double
    • The continued fraction for the incomplete beta, and the series and continued fraction for the incomplete gamma, stopped as soon as a term fell below a bare 1e-15 or 3e-15 -- so `pt`, `pf`, `pchisq`, `qchisq`, `qf` and everything built on them returned about sixteen digits on a perl carrying nineteen or thirty-four. `LIKER_EPS_SCALE` carries each of them to the build's own width; it is `NV_EPSILON / DBL_EPSILON`, which is exactly 1.0 when NV is a double, so every value produced on a double build is unchanged to the last bit. Thresholds that are part of an algorithm's definition -- `FT_TOL`, which is R's `uniroot()` default, and MASS's `double.eps^0.25` in `nb_theta_ml()` -- are numbers from the reference implementation and stay put; each says so where it is defined.
    • `FT_EPS` looks like one of those and is not: it is the other half of R's `uniroot()` stopping rule, `2*FT_EPS*|b| + FT_TOL/2`, and it is also the lower endpoint of the bracket `fisher_test`'s confidence interval is inverted over -- so it sets the largest odds ratio that interval can name, which R reports as `1/DBL_EPSILON`. Scaling it to the build's own epsilon moved the conditional odds ratio for SciPy's gh-3014 table 1.2e-9 off R on a `__float128` build and turned that table's upper limit from 4503599627370496 into 5.2e+33. It stays `DBL_EPSILON` on every build, and now says why.
    • `_qgamma()` in `LikeR.pm` inverted `1 - _igamc($shape, $x)`, which is the cancellation `igam()` was added to avoid: below a lower tail of about `NV_EPSILON` that difference can only be a multiple of `NV_EPSILON`, and below about 1e-16 it is exactly 0, so a bisection against it has nothing to bisect on. It now bisects against the new private `_pgamma_lower`, which is `igam()` itself. `age_standardize()`'s Fay-Feuer interval asks for the `alpha/2` quantile, so the lower limit is what this reaches: at `conf.level => 1 - 2e-12` it moved from 1.1e-6 off R to 3.3e-7 off. It cannot be pushed much further from the Perl side -- `1 - 2e-17` is already 1 in a double -- which is why `t/age_standardize.t` checks the primitive itself against R's `pgamma` as well as the interval.
  • The rank test in aov()'s QR was not scale invariant
    • `apply_householder_aov()` declared a column aliased on the absolute `max_val < 1e-10`, which is not a statement about collinearity but about units. A design whose columns are all smaller than 1e-10 -- a predictor in metres that wanted micrometres, a rate per person-year -- had every column declared aliased at step 0, and `aov` reported zero degrees of freedom and a zero sum of squares for every term on perfectly well-conditioned data: `aov` on `x` gave R's F of 1.9927680012954 and on `x * 1e-12` gave NaN. The test is now relative to each column's own scale, taken before the reduction starts, which is what `sweep_matrix_ols()` has always done for `lm`.
  • strtok() in the formula parsers
    • `lm_formula_terms()` and `aov` split the right-hand side with `strtok()`, whose position lives in a libc static. On a `-Dusethreads` perl every thread is its own interpreter inside one process and shares that static with all the others, so two threads fitting a model at the same moment could each be handed the other's term list. Replaced by `lm_tok()`, which takes the cursor from the caller and has no state of its own. It is the only routine with hidden state this file used.
  • Five symbols the shared object should not have exported
    • `approx_pnorm()`, `igamc()`, `get_p_value()` and the `cs_uninit_catcher` XSUB were compiled with external linkage, so the `.so` exported them alongside `boot_Stats__LikeR()`. `igamc` in particular is a name a numerical library might well define too, and the dynamic loader resolves the first definition it sees -- an interposed one would silently replace every chi-square tail this module computes. All four are static now. A fifth, `compare_doubles()`, had had no caller since the `qsort()` comparators were replaced by `LIKER_DEFINE_SORT()`, and is gone.
    • The generated `.c` is also clean under `-Wsign-compare`, which is not in `Makefile.PL`'s flags but is the check the mixed-sign comparisons this file's type rules can introduce would show up in. Seventeen of them were left; all were benign, and all are gone.
  • Smaller: memory and time that was being spent on nothing
    • `evaluate_term()` `savepv()`d the term string -- a malloc, a strcpy and a free -- for every cell of every design matrix `lm`, `glm`, `aov` and `anova` build, in order that two branches it was not going to take could write NULs into it. A bare column name, which is what nearly every cell asks for, now allocates nothing.
    • `melt()` built the whole long frame as one throwaway record per output row, with a nested arrayref of the id values, before materialising it: a melt of R rows over V value columns held R*V of them alive at once beside the result they were about to become. It emits into the requested shape as the loops go.
    • `table_one()` built its per-group row lists with one `grep` over the whole frame per group -- O(groups x rows), the same shape as the O(levels x groups x rows) counting beside it that 0.315 replaced with a single pass. One bucketing pass now.
    • `csort()`'s AoH -> HoA materialisation walked the sorted rows once per output column, fetching and type-checking each row `nk` times for the n*nk cells it produces. Row-outer now, one fetch per row, columns allocated at their final length and filled through `AvARRAY` as `filter()` and `mg_column()` do.
    • `colnames()`, `_present_keys()` and `_rename_inplace()` flattened a copy of every row reference in the frame (`my @rows = @$df`) to read it once; `assign()`'s HoA row view is filled with one hash slice rather than a keyed store per column.
  • Tests
    • `t/aov.regressions.t` pins every one of the `aov` items above: the long interaction component and the long formula (the answer must not depend on how long a column's name is), the malformed HoH, the ragged HoA, the determinism of `group.stats` and of `.`, the scale invariance of the rank test, and `I(x-1)`. Against the code these entries replace it fails in thirteen places, and the `ks_test`, `rbinom` and `interpolate` files below hang or exhaust memory on it rather than merely failing.
    • `t/interpolate.spline.banded.t` checks the banded solve against SciPy at 200, 2,000 and 20,000 anchors, with the generator committed beside it as `t/interpolate.spline.banded.py`, and interpolates a column of 200,000 -- which the dense build cannot allocate.
    • `t/rbinom.dist.t` runs a chi-square goodness-of-fit against the exact binomial CDF over thirteen `(size, prob)` pairs chosen to cross both branches of BTPE and the reflection, plus the moments, the support, the short circuits and reproducibility under `srand`. It deliberately pins no individual variate: pinning one would pin the algorithm rather than the distribution.
    • `t/ks_test.exact.guard.t` pins the one- and two-sample exact p-values against R and asserts that a forced exact run past the cap warns and falls back; `t/value_counts.utf8.t` checks `value_counts` against what a Perl hash makes of the same list, in both directions -- "\x{e9}" and "\xe9" are one value, "\x{263A}" and its three bytes are two; `t/dnorm.nv_width.t` pins `dnorm` against R and asserts the identity `dnorm(x) == exp(-x^2/2) / sqrt(2*pi)` at whatever width the perl running it carries -- which is one assertion that covers every build, since both sides underflow together on a double and neither does on a wider NV. Its R table stops at |x| = 29 on purpose: `perl-5.10.1` reads `2.1200065515246056e-298` as `1.999999999999999e-298` and anything below ~1e-308 as 0, so a table of R's far-tail values would be testing perl's own `atof`.
    • `t/age_standardize.t` gains the gamma quantile at a small tail probability, pinned against R's `qgamma` at four confidence levels, and checks `_pgamma_lower` directly against R's `pgamma` -- `conf.level` cannot reach far enough into the tail on its own to separate it from `1 - _igamc`, because `1 - 2e-17` is already 1 in a double. Where that subtraction starts losing the tail is a build property and the test asks the perl for it rather than assuming a double: its first draft pinned a literal 1e-20 and passed everywhere except `__float128`, where 1e-20 is still fourteen orders of magnitude above the point the subtraction fails at.
    • All of it passes on all seven perls in the local matrix -- 5.10.1, 5.12.5 (`long double`), 5.42.3-thr, 5.44.0, 5.44.0+x87, 5.44.0-quadmath (`__float128`) and 5.44.0-i686 (`ivsize=4`) -- 41,166 tests on the wider NV widths and 41,160 on the rest, the difference being the six `dnorm` assertions that only a build whose exponent range reaches past a double's has anything to check. The generated `.c` compiles clean under strict `-std=c99` and under `-Wall -Wextra -Wsign-compare`, on a quadmath CORE as well as a double one.

Documentation

CONTRIBUTING.md
README.md

Modules

Stats::LikeR
Get basic statistical functions, like in R, but with Perl using XS for performance

Other files