Branch Coverage

File:	LikeR.xs
Coverage:	66.9%

line	%	coverage		branch
30	50	T	F	if (df <= 0.0) return 0.0;
40	100	T	F	for (unsigned short i = 1; i < n_steps; i++) {
52	100	T	F	double weight = (i % 2 != 0) ? 4.0 : 2.0;
73	100	T	F	for (size_t i = 0; i < n; i++) {
79	100	T	F	for (size_t i = 0; i < n; ) {
81	100	T	F	while (j < n && items[j].val == items[i].val) j++;
	50	T	F	while (j < n && items[j].val == items[i].val) j++;
83	100	T	F	for (size_t k = i; k < j; k++) items[k].rank = avg_rank;
87	100	T	F	for (size_t i = 0; i < n; i++) ranks[i] = items[i].rank;
93	100	T	F	if (prob <= 0.0) return 0;
94	100	T	F	if (prob >= 1.0) return size;
97	100	T	F	for (size_t i = 0; i < size; i++) {
98	100	T	F	if (Drand01() <= prob) successes++;
111	100	T	F	for(size_t k = 0; k <= max_x - min_x; ++k) {
113	100	T	F	if (d_val > max_d) max_d = d_val;
117	100	T	F	for(size_t k = 0; k <= max_x - min_x; ++k) {
124	100	T	F	for(size_t k = 0; k <= max_x - min_x; ++k) {
126	100	T	F	if (min_x + k <= a) *lower_tail += p_prob;
127	100	T	F	if (min_x + k >= a) *upper_tail += p_prob;
135	100	T	F	size_t min_x = (r2 > c1) ? 0 : c1 - r2;
143	100	T	F	for(size_t x = min_x; x <= max_x; ++x) {
148	50	T	F	if (a == min_x && a == max_x) *mle_or = 1.0;
	0	T	F	if (a == min_x && a == max_x) *mle_or = 1.0;
149	50	T	F	else if (a == min_x) *mle_or = 0.0;
150	100	T	F	else if (a == max_x) *mle_or = INFINITY;
153	50	T	F	for (unsigned short int i = 0; i < 3000; i++) {
156	100	T	F	for(size_t k = 0; k <= max_x - min_x; ++k) {
158	100	T	F	if (d_val > max_d) max_d = d_val;
161	100	T	F	for(size_t k = 0; k <= max_x - min_x; ++k) {
168	100	T	F	if (exp_val > a) log_high = log_mid;
170	100	T	F	if (log_high - log_low < 1e-15) break;
179	100	T	F	if (!is_less) {
180	100	T	F	double target_alpha = is_greater ? alpha : alpha / 2.0;
181	50	T	F	if (a != min_x) {
183	50	T	F	for (unsigned short int i = 0; i < 1000; i++) {
188	100	T	F	if (err < best_err) { best_err = err; best = mid; }
189	100	T	F	if (ut > target_alpha) log_high = log_mid;
191	100	T	F	if (log_high - log_low < 1e-15) break;
198	100	T	F	if (!is_greater) {
199	100	T	F	double target_alpha = is_less ? alpha : alpha / 2.0;
200	100	T	F	if (a != max_x) {
202	50	T	F	for (unsigned short int i = 0; i < 1000; i++) {
207	100	T	F	if (err < best_err) { best_err = err; best = mid; }
208	100	T	F	if (lt > target_alpha) log_low = log_mid;
210	100	T	F	if (log_high - log_low < 1e-15) break;
221	100	T	F	size_t min_x = (r2 > c1) ? 0 : c1 - r2;
226	100	T	F	for(size_t x = min_x; x <= max_x; ++x) {
232	100	T	F	if (strcmp(alt, "less") == 0) {
233	100	T	F	for(size_t x = min_x; x <= a; ++x) p_val += exp(logdc[x - min_x]);
234	100	T	F	} else if (strcmp(alt, "greater") == 0) {
235	100	T	F	for(size_t x = a; x <= max_x; ++x) p_val += exp(logdc[x - min_x]);
239	100	T	F	for(size_t x = min_x; x <= max_x; ++x) {
241	100	T	F	if (p_cur <= p_obs * relErr) p_val += p_cur;
246	50	T	F	return (p_val > 1.0) ? 1.0 : p_val;
261	100	T	F	for (size_t k = 0; k < n; k++) {
266	100	T	F	for (size_t k = 0; k < n; k++) {
269	100	T	F	if (fabs(A[k * n + k]) <= 1e-10 * orig_diag[k] \|\| fabs(A[k * n + k]) < 1e-24) {
	50	T	F	if (fabs(A[k * n + k]) <= 1e-10 * orig_diag[k] \|\| fabs(A[k * n + k]) < 1e-24) {
272	100	T	F	for (size_t i = 0; i < n; i++) {
281	100	T	F	for (size_t j = 0; j < n; j++) A[k * n + j] *= pivot;
282	100	T	F	for (size_t i = 0; i < n; i++) {
283	100	T	F	if (i != k && A[i * n + k] != 0.0) {
	100	T	F	if (i != k && A[i * n + k] != 0.0) {
286	100	T	F	for (size_t j = 0; j < n; j++) {
299	100	T	F	if (row_hashes) {
301	50	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	50	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	50	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
305	50	T	F	} else if (data_hoa) {
307	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
312	50	T	F	if (val && SvOK(*val)) {
	100	T	F	if (val && SvOK(*val)) {
313	50	T	F	if (looks_like_number(val)) return SvNV(val);
322	50	T	F	if (data_hoa) {
325	100	T	F	while ((entry = hv_iternext(data_hoa))) {
328	0	T	F	} else if (row_hashes && n > 0 && row_hashes[0]) {
	0	T	F	} else if (row_hashes && n > 0 && row_hashes[0]) {
	0	T	F	} else if (row_hashes && n > 0 && row_hashes[0]) {
331	0	T	F	while ((entry = hv_iternext(row_hashes[0]))) {
340	50	T	F	if (!term \|\| term[0] == '\0') return NAN;
	50	T	F	if (!term \|\| term[0] == '\0') return NAN;
344	100	T	F	if (colon) {
350	50	T	F	if (isnan(left) \|\| isnan(right)) return NAN;
	50	T	F	if (isnan(left) \|\| isnan(right)) return NAN;
353	50	T	F	if (strncmp(term_cpy, "I(", 2) == 0) {
355	0	T	F	if (end) *end = '\0';
359	0	T	F	if (caret) {
366	0	T	F	if (isnan(v)) return NAN;
367	0	T	F	return power == 1 ? v : pow(v, power);
376	100	T	F	for (size_t i = 0; i < n; i++) {
378	100	T	F	if (row_hashes) {
380	100	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	50	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	50	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
384	50	T	F	} else if (data_hoa) {
386	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
391	100	T	F	if (val && SvOK(*val)) {
	50	T	F	if (val && SvOK(*val)) {
392	100	T	F	if (looks_like_number(*val)) return false; // First valid is number -> Numeric Column
402	50	T	F	if (row_hashes) {
404	0	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	0	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
	0	T	F	if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
408	50	T	F	} else if (data_hoa) {
410	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
	50	T	F	if (col && SvROK(col) && SvTYPE(SvRV(col)) == SVt_PVAV) {
415	50	T	F	if (val && SvOK(*val)) {
	50	T	F	if (val && SvOK(*val)) {
430	100	T	F	if (diff < 0) return -1;
431	50	T	F	if (diff > 0) return 1;
447	100	T	F	if (diff < 0) return -1;
448	100	T	F	if (diff > 0) return 1;
456	50	T	F	Newx(ri, n, RankItem);
457	100	T	F	for (size_t i = 0; i < n; i++) { ri[i].val = in[i]; ri[i].idx = i; }
461	100	T	F	while (i < n) {
464	100	T	F	while (j + 1 < n && ri[j + 1].val == ri[j].val) j++;
	100	T	F	while (j + 1 < n && ri[j + 1].val == ri[j].val) j++;
467	100	T	F	for (size_t k = i; k <= j; k++) out[ri[k].idx] = avg;
477	100	T	F	for (size_t i = 0; i < n; i++) {
483	50	T	F	if (den == 0.0) return NAN;
497	100	T	F	for (size_t i = 0; i < n - 1; i++) {
498	100	T	F	for (size_t j = i + 1; j < n; j++) {
501	100	T	F	if (sx == 0 && sy == 0) { /* joint tie â€” not counted */ }
	50	T	F	if (sx == 0 && sy == 0) { /* joint tie â€” not counted */ }
502	100	T	F	else if (sx == 0) tie_x++;
503	50	T	F	else if (sy == 0) tie_y++;
504	50	T	F	else if (sx == sy) C++;
509	50	T	F	if (denom == 0.0) return NAN;
517	100	T	F	if (strcmp(method, "spearman") == 0) {
519	50	T	F	Newx(rx, n, double); Newx(ry, n, double);
	50	T	F	Newx(rx, n, double); Newx(ry, n, double);
526	100	T	F	if (strcmp(method, "kendall") == 0)
542	50	T	F	if (fabs(d) < FPMIN) d = FPMIN;
544	50	T	F	for (m = 1; m <= MAX_ITER; m++) {
548	50	T	F	if (fabs(d) < FPMIN) d = FPMIN;
550	50	T	F	if (fabs(c) < FPMIN) c = FPMIN;
554	50	T	F	if (fabs(d) < FPMIN) d = FPMIN;
556	50	T	F	if (fabs(c) < FPMIN) c = FPMIN;
558	100	T	F	if (fabs(del - 1.0) < EPS) break;
564	50	T	F	if (x <= 0.0) return 0.0;
565	100	T	F	if (x >= 1.0) return 1.0;
567	100	T	F	if (x < (a + 1.0) / (a + b + 2.0)) return bt * _incbeta_cf(a, b, x) / a;
574	100	T	F	if (strcmp(alt, "less") == 0) return (t < 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
	100	T	F	if (strcmp(alt, "less") == 0) return (t < 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
575	100	T	F	if (strcmp(alt, "greater") == 0) return (t > 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
	50	T	F	if (strcmp(alt, "greater") == 0) return (t > 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
583	100	T	F	while (get_t_pvalue(high, df, "greater") > p_tail) {
586	50	T	F	if (high > 1000000.0) break; /* Fallback limit */
589	50	T	F	for (unsigned short int i = 0; i < 100; i++) {
592	100	T	F	if (p_mid > p_tail) {
597	100	T	F	if (high - low < 1e-8) break;
609	0	T	F	if (n == 0) return 1;
618	50	T	F	double step = (n_bins > 0) ? (breaks[1] - breaks[0]) : 0.0;
620	100	T	F	for (size_t i = 0; i < n_bins; i++) {
625	100	T	F	if (step > 0.0) {
626	100	T	F	for (size_t j = 0; j < n; j++) {
629	50	T	F	if (isnan(val) \|\| isinf(val) \|\| val < min_val) continue;
	50	T	F	if (isnan(val) \|\| isinf(val) \|\| val < min_val) continue;
	50	T	F	if (isnan(val) \|\| isinf(val) \|\| val < min_val) continue;
633	100	T	F	if (idx >= n_bins) {
639	100	T	F	while (idx > 0 && val <= breaks[idx]) {
	100	T	F	while (idx > 0 && val <= breaks[idx]) {
643	100	T	F	while (idx < n_bins - 1 && val > breaks[idx + 1]) {
	50	T	F	while (idx < n_bins - 1 && val > breaks[idx + 1]) {
648	50	T	F	} else if (n_bins > 0) {
653	100	T	F	for (size_t i = 0; i < n_bins; i++) {
655	100	T	F	if (bin_width > 0) {
658	50	T	F	density[i] = (n_bins == 1) ? 1.0 : 0.0;
681	100	T	F	if (fabs(y) < 0.42) {
687	100	T	F	if (y > 0) r = 1.0 - p;
691	100	T	F	if (y < 0) x = -x;
709	100	T	F	for (size_t i = 0; i < n; i++) { perm[i] = i + 1; c[i] = 0; }
724	100	T	F	TALLY_PERM(); /* initial permutation [1, 2, ..., n] */
	50	T	F	TALLY_PERM(); /* initial permutation [1, 2, ..., n] */
	50	T	F	TALLY_PERM(); /* initial permutation [1, 2, ..., n] */
727	100	T	F	while (k < n) {
728	100	T	F	if (c[k] < k) {
730	100	T	F	if (k % 2 == 0) {
735	100	T	F	TALLY_PERM();
	100	T	F	TALLY_PERM();
	100	T	F	TALLY_PERM();
751	50	T	F	if (strcmp(alt, "greater") == 0) return p_le;
752	50	T	F	if (strcmp(alt, "less") == 0) return p_ge;
754	50	T	F	double p = 2.0 * (p_le < p_ge ? p_le : p_ge);
755	50	T	F	return (p > 1.0) ? 1.0 : p;
764	100	T	F	for (long i = 0; i <= max_inv; i++) dp[i] = 0.0;
767	100	T	F	for (size_t i = 2; i <= n; i++) {
769	100	T	F	for (long k = 0; k <= max_inv; k++) next_dp[k] = 0.0;
771	100	T	F	for (int k = 0; k <= current_max_inv; k++) {
773	100	T	F	for (int j = 0; j <= i - 1 && k - j >= 0; j++) {
	100	T	F	for (int j = 0; j <= i - 1 && k - j >= 0; j++) {
784	50	T	F	if (i_obs < 0) i_obs = 0;
785	50	T	F	if (i_obs > max_inv) i_obs = max_inv;
787	100	T	F	for (long k = i_obs; k <= max_inv; k++) p_le += dp[k];
789	100	T	F	for (long k = 0; k <= i_obs; k++) p_ge += dp[k];
791	50	T	F	if (strcmp(alt, "greater") == 0) return p_ge;
792	100	T	F	if (strcmp(alt, "less") == 0) return p_le;
794	50	T	F	double p = 2.0 * (p_ge < p_le ? p_ge : p_le);
795	50	T	F	return p > 1.0 ? 1.0 : p;
799	50	T	F	if (f <= 0.0) return 0.0;
808	100	T	F	for (size_t k = 0; k < p; k++) {
810	50	T	F	if (r >= n) {
816	100	T	F	for (size_t i = r; i < n; i++) {
817	100	T	F	if (fabs(X[i][k]) > max_val) max_val = fabs(X[i][k]);
819	100	T	F	if (max_val < 1e-10) {
825	100	T	F	for (size_t i = r; i < n; i++) {
830	100	T	F	double s = (X[r][k] > 0) ? -norm : norm;
834	100	T	F	for (size_t j = k + 1; j < p; j++) {
836	100	T	F	for (size_t i = r + 1; i < n; i++) dot += X[i][j] * X[i][k];
839	100	T	F	for (size_t i = r + 1; i < n; i++) X[i][j] += tau * X[i][k];
844	100	T	F	for (size_t i = r + 1; i < n; i++) dot_y += y[i] * X[i][k];
847	100	T	F	for (size_t i = r + 1; i < n; i++) y[i] += tau_y * X[i][k];
863	50	T	F	if (!sv \|\| !SvOK(sv)) return 0;
	50	T	F	if (!sv \|\| !SvOK(sv)) return 0;
866	100	T	F	for (size_t i = 0; i < len; i++) {
867	50	T	F	if (!isdigit(s[i])) return 1;
874	50	T	F	if (!str) str = "";
876	100	T	F	if (strstr(str, sep) != NULL \|\| strchr(str, '"') != NULL \|\| strchr(str, '\r') != NULL \|\| strchr(str, '\n') != NULL) {
	100	T	F	if (strstr(str, sep) != NULL \|\| strchr(str, '"') != NULL \|\| strchr(str, '\r') != NULL \|\| strchr(str, '\n') != NULL) {
	50	T	F	if (strstr(str, sep) != NULL \|\| strchr(str, '"') != NULL \|\| strchr(str, '\r') != NULL \|\| strchr(str, '\n') != NULL) {
	50	T	F	if (strstr(str, sep) != NULL \|\| strchr(str, '"') != NULL \|\| strchr(str, '\r') != NULL \|\| strchr(str, '\n') != NULL) {
880	100	T	F	if (needs_quotes) {
882	100	T	F	for (const char restrict p = str; p; p++) {
883	100	T	F	if (*p == '"') {
899	100	T	F	for (size_t i = 0; i < len; i++) {
900	100	T	F	if (i > 0) PerlIO_write(fh, sep, sep_len);
901	50	T	F	if (row[i]) {
912	50	T	F	if (x < 0.0 \|\| a <= 0.0) return 1.0;
	50	T	F	if (x < 0.0 \|\| a <= 0.0) return 1.0;
913	50	T	F	if (x == 0.0) return 1.0;
916	100	T	F	if (x < a + 1.0) {
920	100	T	F	while (fabs(term) > 1e-15) {
933	50	T	F	while (i < 10000) { // Safety bound
937	50	T	F	if (fabs(d) < 1e-30) d = 1e-30;
939	50	T	F	if (fabs(c) < 1e-30) c = 1e-30;
943	100	T	F	if (fabs(del - 1.0) < 1e-15) break;
951	50	T	F	if (df <= 0) return 1.0;
952	50	T	F	if (stat <= 0.0) return 1.0;
963	50	T	F	if (k < 0 \|\| k > n) return 0.0L;
	50	T	F	if (k < 0 \|\| k > n) return 0.0L;
964	50	T	F	if (k > n / 2) k = n - k;
966	100	T	F	for (int i = 1; i <= k; i++) {
977	100	T	F	if (k < 0) return 0.0;
978	50	T	F	if (k >= max_u) return 1.0;
983	100	T	F	for (int j = 1; j <= n; j++) {
984	100	T	F	for (int i = j; i <= max_u; i++) w[i] += w[i - j];
985	100	T	F	for (int i = max_u; i >= j + m; i--) w[i] -= w[i - j - m];
989	100	T	F	for (int i = 0; i <= k; i++) cum_p += w[i];
1003	50	T	F	if (k < 0) return 0.0;
1004	100	T	F	if (k >= max_v) return 1.0;
1009	100	T	F	for (int i = 1; i <= n; i++) {
1010	100	T	F	for (int j = max_v; j >= i; j--) w[j] += w[j - i];
1014	100	T	F	for (int i = 0; i <= k; i++) cum_p += w[i];
1030	50	T	F	if (n == 0) return 0.0;
1035	100	T	F	while (i < n) {
1037	100	T	F	while (j < n && ri[j].val == ri[i].val) j++;
	100	T	F	while (j < n && ri[j].val == ri[i].val) j++;
1039	100	T	F	for (size_t k = i; k < j; k++) ri[k].rank = r;
1041	100	T	F	if (t > 1) { has_ties = 1; tie_adj += ((double)t t * t - t); }
1059	50	T	F	if (n == 0) return 1.0;
1060	50	T	F	if (n < 0) return 1.0 / r_pow_di(x, -n);
1062	100	T	F	for (int i = 0; i < n; i++) val *= x;
1070	0	T	F	if(x <= 0.) {
1071	0	T	F	if(lower) p = 0.;
1073	0	T	F	} else if(x < 1.) {
1078	0	T	F	for(k = 1; k < k_max; k += 2) {
1082	0	T	F	if(!lower) p = 1.0 - p;
1087	0	T	F	if(lower) {
1092	0	T	F	while(fabs(old_val - new_val) > tol) {
1105	100	T	F	for(unsigned int i = 0; i < m; i++) {
1106	100	T	F	for(unsigned int j = 0; j < m; j++) {
1108	100	T	F	for(unsigned int k = 0; k < m; k++) s += A[i * m + k] * B[k * m + j];
1115	100	T	F	if(n == 1) {
1116	100	T	F	for(int i = 0; i < m * m; i++) V[i] = A[i];
1124	100	T	F	if((n % 2) == 0) {
1125	100	T	F	for(int i = 0; i < m * m; i++) V[i] = B[i];
1131	50	T	F	if(V[(m / 2) * m + (m / 2)] > 1e140) {
1132	0	T	F	for(int i = 0; i < m * m; i++) V[i] = V[i] * 1e-140;
1146	100	T	F	for(int i = 0; i < m; i++) {
1147	100	T	F	for(int j = 0; j < m; j++) {
1148	100	T	F	if(i - j + 1 < 0) H[i * m + j] = 0;
1152	100	T	F	for(int i = 0; i < m; i++) {
1156	50	T	F	H[(m - 1) * m] += ((2 * h - 1 > 0) ? r_pow_di(2 * h - 1, m) : 0);
1158	100	T	F	for(int i = 0; i < m; i++) {
1159	100	T	F	for(int j = 0; j < m; j++) {
1160	100	T	F	if(i - j + 1 > 0) {
1161	100	T	F	for(int g = 1; g <= i - j + 1; g++) H[i * m + j] /= g;
1170	100	T	F	for(int i = 1; i <= n; i++) {
1172	50	T	F	if(s < 1e-140) {
1191	100	T	F	while(i < nx \|\| j < ny) {
	50	T	F	while(i < nx \|\| j < ny) {
1193	50	T	F	if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
	100	T	F	if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
	100	T	F	if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
1194	50	T	F	else if (i < nx) val = x[i];
1197	100	T	F	while(i < nx && x[i] <= val) i++;
	100	T	F	while(i < nx && x[i] <= val) i++;
1198	100	T	F	while(j < ny && y[j] <= val) j++;
	100	T	F	while(j < ny && y[j] <= val) j++;
1204	100	T	F	if (diff > max_d_plus) max_d_plus = diff;
1205	100	T	F	if (-diff > max_d_minus) max_d_minus = -diff;
1206	100	T	F	if (fabs(diff) > max_d) max_d = fabs(diff);
1215	100	T	F	if (two_sided) return (fabs(r - s) >= q);
1225	100	T	F	for(unsigned int j = 1; j <= n; j++) {
1226	100	T	F	if(psmirnov_exact_test(q, 0., j / nd, two_sided)) u[j] = 1.;
1229	100	T	F	for(unsigned int i = 1; i <= m; i++) {
1230	100	T	F	if(psmirnov_exact_test(q, i / md, 0., two_sided)) u[0] = 1.;
1231	100	T	F	for(int j = 1; j <= n; j++) {
1232	100	T	F	if(psmirnov_exact_test(q, i / md, j / nd, two_sided)) u[j] = 1.;
1247	50	T	F	if (nu < 1e-7) nu = 1e-7;
1255	50	T	F	if (strict && tside == 2) {
	0	T	F	if (strict && tside == 2) {
1276	50	T	F	}
	50	T	F	}
	50	T	F	}
1281	50	T	F	double p_mid = pf(mid, df1, df2);
1282	100	T	F
	50	T	F
1285	50	T	F	} else {
1292	100	T	F	// --- XS SECTION ---
1295	50	T	F	SV* ks_test(...)
1296	50	T	F	CODE:
1297	50	T	F	{
1298	0	T	F	SV restrict x_sv = NULL, restrict y_sv = NULL;
1301	50	T	F	int arg_idx = 0;
1305	50	T	F	x_sv = ST(arg_idx);
	50	T	F	x_sv = ST(arg_idx);
	50	T	F	x_sv = ST(arg_idx);
1313	100	T	F	} else if (SvPOK(ST(arg_idx))) {
	100	T	F	} else if (SvPOK(ST(arg_idx))) {
	50	T	F	} else if (SvPOK(ST(arg_idx))) {
1319	50	T	F	// Parse named arguments
1324	100	T	F	else if (strEQ(key, "y")) y_sv = val;
1326	50	T	F	if (!SvOK(val)) exact = -1;
	50	T	F	if (!SvOK(val)) exact = -1;
	50	T	F	if (!SvOK(val)) exact = -1;
1335	50	T	F	}
	100	T	F	}
	50	T	F	}
1341	100	T	F	if (!is_two_sided && !is_greater && !is_less) {
1343	50	T	F	}
	50	T	F	}
	50	T	F	}
1348	50	T	F
	50	T	F
1357	100	T	F	}
1358	100	T	F
1363	50	T	F	if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
1364	50	T	F	AV restrict y_av = (AV)SvRV(y_sv);
1370	100	T	F	SV**restrict el = av_fetch(y_av, i, 0);
1371	100	T	F	if (el && SvOK(el) && looks_like_number(el)) {
1375	100	T	F
1376	50	T	F	if (valid_nx < 1 \|\| valid_ny < 1) {
1379	50	T	F	}
	50	T	F	}
1383	50	T	F
1390	0	T	F	bool use_exact = false;
1399	50	T	F	for(size_t i=0; i<valid_ny; i++) comb[valid_nx+i] = y_data[i];
	50	T	F	for(size_t i=0; i<valid_ny; i++) comb[valid_nx+i] = y_data[i];
1401	50	T	F
1404	100	T	F	if(comb[i] == comb[i-1]) { has_ties = true; break; }
1412	50	T	F	method_desc = "Two-sample Kolmogorov-Smirnov exact test";
1413	100	T	F	double q = (0.5 + floor(statistic * valid_nx * valid_ny - 1e-7)) / ((double)valid_nx * valid_ny);
1414	100	T	F	p_value = psmirnov_exact_uniq_upper(q, valid_nx, valid_ny, is_two_sided);
1415	50	T	F	} else {
1417	100	T	F	double z = statistic * sqrt((double)(valid_nx * valid_ny) / (valid_nx + valid_ny));
1418	50	T	F	if (is_two_sided) {
1420	50	T	F	} else {
1421	50	T	F	p_value = exp(-2.0 * z * z); // One-sided limit distribution
1423	50	T	F	}
1424	50	T	F	Safefree(y_data);
1426	50	T	F	// --- ONE SAMPLE (e.g. against pnorm) ---
1436	0	T	F
1448	50	T	F	if (is_greater) statistic = max_d_plus;
1449	50	T	F	else if (is_less) statistic = max_d_minus;
1470	50	T	F	}
	100	T	F	}
	50	T	F	}
1475	50	T	F	Safefree(x_data);
	100	T	F	Safefree(x_data);
	50	T	F	Safefree(x_data);
1480	50	T	F	hv_stores(res, "p_value", newSVnv(p_value));
1484	100	T	F	}
1487	100	T	F
1488	100	T	F	SV* wilcox_test(...)
1489	100	T	F	CODE:
1490	50	T	F	{
1491	100	T	F	SV restrict x_sv = NULL, restrict y_sv = NULL;
1492	50	T	F	bool paired = false, correct = true;
1493	0	T	F	double mu = 0.0;
1496	50	T	F	int arg_idx = 0;
1500	100	T	F	arg_idx++;
	50	T	F	arg_idx++;
	50	T	F	arg_idx++;
1504	50	T	F	y_sv = ST(arg_idx);
1508	100	T	F	if ((items - arg_idx) % 2 != 0) {
	50	T	F	if ((items - arg_idx) % 2 != 0) {
	50	T	F	if ((items - arg_idx) % 2 != 0) {
1516	100	T	F	else if (strEQ(key, "y")) y_sv = val;
	100	T	F	else if (strEQ(key, "y")) y_sv = val;
1519	100	T	F	else if (strEQ(key, "mu")) mu = SvNV(val);
1521	50	T	F	if (!SvOK(val)) exact = -1;
	50	T	F	if (!SvOK(val)) exact = -1;
	50	T	F	if (!SvOK(val)) exact = -1;
1527	100	T	F	// --- Validate required / types ---
1529	50	T	F	croak("wilcox_test: 'x' is a required argument and must be an ARRAY reference");
	50	T	F	croak("wilcox_test: 'x' is a required argument and must be an ARRAY reference");
	50	T	F	croak("wilcox_test: 'x' is a required argument and must be an ARRAY reference");
1535	50	T	F	size_t ny = 0;
1536	50	T	F	if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
1541	100	T	F	const char *restrict method_desc = "";
	100	T	F	const char *restrict method_desc = "";
1544	50	T	F	if (ny > 0 && !paired) {
1545	50	T	F	RankInfo restrict ri = (RankInfo )safemalloc((nx + ny) * sizeof(RankInfo));
1546	50	T	F	size_t valid_nx = 0, valid_ny = 0;
	50	T	F	size_t valid_nx = 0, valid_ny = 0;
	100	T	F	size_t valid_nx = 0, valid_ny = 0;
1548	100	T	F	SV**restrict el = av_fetch(x_av, i, 0);
	50	T	F	SV**restrict el = av_fetch(x_av, i, 0);
1552	100	T	F	valid_nx++;
1557	100	T	F	if (el && SvOK(el) && looks_like_number(el)) {
1558	50	T	F	ri[valid_nx + valid_ny].val = SvNV(*el);
1560	0	T	F	valid_ny++;
1564	50	T	F	if (valid_ny == 0) { Safefree(ri); croak("not enough 'y' observations"); }
1570	50	T	F	statistic = w_rank_sum - (double)valid_nx * (valid_nx + 1.0) / 2.0;
1571	50	T	F
	100	T	F
1572	0	T	F	if (exact == 1) use_exact = true;
1573	0	T	F	else if (exact == 0) use_exact = false;
1577	50	T	F	warn("cannot compute exact p-value with ties; falling back to approximation");
1578	50	T	F	use_exact = false;
1583	100	T	F	double p_greater = 1.0 - exact_pwilcox(statistic - 1.0, valid_nx, valid_ny);
	50	T	F	double p_greater = 1.0 - exact_pwilcox(statistic - 1.0, valid_nx, valid_ny);
	100	T	F	double p_greater = 1.0 - exact_pwilcox(statistic - 1.0, valid_nx, valid_ny);
1587	100	T	F	else {
1589	50	T	F	p_value = 2.0 * p;
	50	T	F	p_value = 2.0 * p;
	50	T	F	p_value = 2.0 * p;
1592	100	T	F	method_desc = correct ? "Wilcoxon rank sum test with continuity correction" : "Wilcoxon rank sum test";
1594	50	T	F	double var = ((double)valid_nx * valid_ny / 12.0) * ((total_n + 1.0) - tie_adj / (total_n * (total_n - 1.0)));
	50	T	F	double var = ((double)valid_nx * valid_ny / 12.0) * ((total_n + 1.0) - tie_adj / (total_n * (total_n - 1.0)));
	50	T	F	double var = ((double)valid_nx * valid_ny / 12.0) * ((total_n + 1.0) - tie_adj / (total_n * (total_n - 1.0)));
1597	50	T	F	double CORRECTION = 0.0;
1601	50	T	F	else if (strcmp(alternative, "less") == 0) CORRECTION = -0.5;
1605	50	T	F	if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z);
1610	100	T	F	} else { // --- ONE SAMPLE / PAIRED ---
1617	100	T	F	if (!x_el \|\| !SvOK(x_el) \|\| !looks_like_number(x_el)) continue;
1618	100	T	F	double dx = SvNV(*x_el);
1620	50	T	F	if (paired) {
1621	50	T	F	SV**restrict y_el = av_fetch(y_av, i, 0);
1622	50	T	F	if (!y_el \|\| !SvOK(y_el) \|\| !looks_like_number(y_el)) continue;
	50	T	F	if (!y_el \|\| !SvOK(y_el) \|\| !looks_like_number(y_el)) continue;
1623	50	T	F	double dy = SvNV(*y_el);
	50	T	F	double dy = SvNV(*y_el);
1627	50	T	F	} else {
	50	T	F	} else {
1631	50	T	F	}
1636	50	T	F	}
1637	50	T	F	RankInfo ri = (RankInfo )safemalloc(n_nz * sizeof(RankInfo));
1639	50	T	F	ri[i].val = fabs(diffs[i]);
1643	0	T	F	double tie_adj = rank_and_count_ties(ri, n_nz, &has_ties);
1648	0	T	F	if (exact == 1) use_exact = true;
1649	0	T	F	else if (exact == 0) use_exact = false;
	0	T	F	else if (exact == 0) use_exact = false;
1650	0	T	F	else use_exact = (n_nz < 50 && !has_ties);
1651	0	T	F	if (use_exact && has_ties) {
1655	0	T	F	if (use_exact && has_zeroes) {
1656	0	T	F	warn("cannot compute exact p-value with zeroes; falling back to approximation");
1661	50	T	F	double p_less = exact_psignrank(statistic, n_nz);
1677	50	T	F	if (strcmp(alternative, "two.sided") == 0) CORRECTION = (z > 0 ? 0.5 : -0.5);
1680	50	T	F	}
	100	T	F	}
	50	T	F	}
1683	50	T	F	if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z);
1691	100	T	F	hv_stores(res, "statistic", newSVnv(statistic));
	50	T	F	hv_stores(res, "statistic", newSVnv(statistic));
	100	T	F	hv_stores(res, "statistic", newSVnv(statistic));
1694	100	T	F	hv_stores(res, "alternative", newSVpv(alternative, 0));
1697	100	T	F	OUTPUT:
1698	100	T	F	RETVAL
1699	100	T	F
1702	100	T	F	CODE:
1710	100	T	F	AVrestrict first_row = (AV)SvRV(*first_elem);
1714	100	T	F	r = 1;
1719	100	T	F	int yates = (is_2d && r == 2 && c == 2) ? 1 : 0;
1722	50	T	F	if (is_2d) {
1734	100	T	F	col_sum[j] += val;
1739	100	T	F	AV*restrict exp_row = newAV();
1753	100	T	F	} else {
1754	100	T	F	stat += ((O - E) * (O - E)) / E;
1777	50	T	F	hv_store(results, "statistic", 9, newSVnv(stat), 0);
	50	T	F	hv_store(results, "statistic", 9, newSVnv(stat), 0);
1779	100	T	F	hv_store(results, "p_value", 7, newSVnv(p_val), 0);
	50	T	F	hv_store(results, "p_value", 7, newSVnv(p_val), 0);
1784	50	T	F	} else {
1794	100	T	F
1795	50	T	F	PROTOTYPES: ENABLE
1798	50	T	F	PPCODE:
1799	100	T	F	{
1800	50	T	F	SV *restrict data_sv = NULL;
1801	100	T	F	SV *restrict file_sv = NULL;
1802	50	T	F	unsigned int arg_idx = 0;
1806	50	T	F	int type = SvTYPE(SvRV(ST(arg_idx)));
	50	T	F	int type = SvTYPE(SvRV(ST(arg_idx)));
1810	100	T	F	}
	50	T	F	}
1814	50	T	F	arg_idx++;
	50	T	F	arg_idx++;
1817	100	T	F	const char *restrict sep = ",";
	50	T	F	const char *restrict sep = ",";
1818	100	T	F	SV *restrict row_names_sv = sv_2mortal(newSViv(1));
	50	T	F	SV *restrict row_names_sv = sv_2mortal(newSViv(1));
1827	100	T	F	else if (strEQ(key, "col.names")) col_names_sv = val;
1829	50	T	F	else if (strEQ(key, "row.names")) row_names_sv = val;
1833	50	T	F
	50	T	F
1837	100	T	F	SV *restrict data_ref = SvRV(data_sv);
	50	T	F	SV *restrict data_ref = SvRV(data_sv);
1845	100	T	F	if (col_names_sv && SvOK(col_names_sv)) {
1847	50	T	F	croak("write_table: 'col.names' must be an ARRAY reference\n");
	50	T	F	croak("write_table: 'col.names' must be an ARRAY reference\n");
	50	T	F	croak("write_table: 'col.names' must be an ARRAY reference\n");
1848	0	T	F	}
1852	100	T	F	AV *restrict rows_av = NULL;
1855	100	T	F	if (SvTYPE(data_ref) == SVt_PVHV) {
1861	50	T	F	if (!first_val \|\| !SvROK(first_val)) {
1863	50	T	F	}
	50	T	F	}
	50	T	F	}
	50	T	F	}
1867	100	T	F	}
1869	50	T	F	is_hoa = (first_type == SVt_PVAV);
	50	T	F	is_hoa = (first_type == SVt_PVAV);
	50	T	F	is_hoa = (first_type == SVt_PVAV);
	50	T	F	is_hoa = (first_type == SVt_PVAV);
1877	50	T	F	if (is_hoh) {
1880	50	T	F	while ((entry = hv_iternext(hv))) {
	100	T	F	while ((entry = hv_iternext(hv))) {
1884	100	T	F	} else {
1885	100	T	F	AV restrict av = (AV)data_ref;
	50	T	F	AV restrict av = (AV)data_ref;
1887	100	T	F	SV **restrict first_ptr = av_fetch(av, 0, 0);
1889	50	T	F	croak("write_table: For ARRAY data, all elements must be HASH references (Array of Hashes)\n");
	50	T	F	croak("write_table: For ARRAY data, all elements must be HASH references (Array of Hashes)\n");
1895	100	T	F	croak("write_table: Mixed data types detected in Array of Hashes. All elements must be HASH references.\n");
1899	100	T	F	}
1905	100	T	F	bool inc_rownames = (row_names_sv && SvTRUE(row_names_sv)) ? 1 : 0;
1910	100	T	F	if (col_names_sv && SvOK(col_names_sv)) {
1918	50	T	F	hv_iterinit((HV*)data_ref);
1919	100	T	F	HE *restrict entry;
1921	50	T	F	HV restrict inner = (HV)SvRV(hv_iterval((HV*)data_ref, entry));
	50	T	F	HV restrict inner = (HV)SvRV(hv_iterval((HV*)data_ref, entry));
1928	100	T	F	unsigned num_cols = hv_iterinit(col_map);
1936	100	T	F	safefree(col_array);
1938	50	T	F	}
1941	50	T	F
1943	100	T	F	if (inc_rownames) header_row[h_idx++] = "";
1945	50	T	F	SV**restrict h_ptr = av_fetch(headers_av, i, 0);
	50	T	F	SV**restrict h_ptr = av_fetch(headers_av, i, 0);
1946	50	T	F	header_row[h_idx++] = (h_ptr && SvOK(h_ptr)) ? SvPV_nolen(h_ptr) : "";
1947	100	T	F	}
	100	T	F	}
1948	100	T	F	print_string_row(fh, header_row, h_idx, sep);
1952	50	T	F	const char *restrict row_array = safemalloc(num_rows sizeof(char*));
1953	50	T	F	for(size_t i=0; i<num_rows; i++) {
1965	100	T	F	SV **restrict inner_hv_ptr = hv_fetch(data_hv, row_array[i], strlen(row_array[i]), 0);
1970	100	T	F	const char restrict col_name = (h_ptr && SvOK(h_ptr)) ? SvPV_nolen(*h_ptr) : "";
1973	100	T	F	if (SvROK(*cell_ptr)) {
1976	100	T	F	safefree(row_data);
	50	T	F	safefree(row_data);
1978	100	T	F	if (rows_av) SvREFCNT_dec(rows_av);
1980	50	T	F	}
	50	T	F	}
1985	100	T	F	}
1990	100	T	F	} else if (is_hoa) { // ----- Hash of Arrays -----
1993	50	T	F	hv_iterinit(data_hv);
1994	100	T	F	HE *restrict entry;
	50	T	F	HE *restrict entry;
1998	0	T	F	if (len > max_rows) max_rows = len;
2000	0	T	F
	0	T	F
2002	0	T	F	AV restrict c_av = (AV)SvRV(col_names_sv);
2012	100	T	F	col_array[i] = SvPV_nolen(hv_iterkeysv(ce));
2013	100	T	F	}
2015	50	T	F	for(unsigned i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0));
	50	T	F	for(unsigned i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0));
2020	100	T	F	rownames_col = SvPV_nolen(row_names_sv);
2022	100	T	F
2023	50	T	F	for(size_t i=0; i<=av_len(headers_av); i++) {
2025	0	T	F	if (!h_ptr \|\| !*h_ptr) continue;
	0	T	F	if (!h_ptr \|\| !*h_ptr) continue;
2028	0	T	F	av_push(filtered_headers, newSVsv(h_sv));
	0	T	F	av_push(filtered_headers, newSVsv(h_sv));
2029	0	T	F	}
2032	0	T	F	headers_av = filtered_headers;
2048	100	T	F	if (rownames_col) {
2050	50	T	F	if (rn_arr_ptr && SvROK(*rn_arr_ptr)) {
	50	T	F	if (rn_arr_ptr && SvROK(*rn_arr_ptr)) {
2052	50	T	F	SV **restrict rn_val_ptr = av_fetch(rn_arr, i, 0);
	50	T	F	SV **restrict rn_val_ptr = av_fetch(rn_arr, i, 0);
2055	100	T	F	PerlIO_close(fh);
	100	T	F	PerlIO_close(fh);
2056	50	T	F	safefree(row_data);
2059	0	T	F	}
2071	100	T	F	}
	50	T	F	}
2074	50	T	F	SV**restrict h_ptr = av_fetch(headers_av, j, 0);
2077	100	T	F	if (arr_ptr && SvROK(*arr_ptr)) {
	50	T	F	if (arr_ptr && SvROK(*arr_ptr)) {
2079	100	T	F	SV **restrict cell_ptr = av_fetch(arr, i, 0);
2081	50	T	F	if (SvROK(*cell_ptr)) {
	50	T	F	if (SvROK(*cell_ptr)) {
2085	100	T	F	croak("write_table: Cannot write nested reference types to table\n");
2087	50	T	F	row_data[d_idx++] = SvPV_nolen(*cell_ptr);
	50	T	F	row_data[d_idx++] = SvPV_nolen(*cell_ptr);
2091	100	T	F	} else {
2098	100	T	F	safefree(row_data);
2103	100	T	F	AV restrict c_av = (AV)SvRV(col_names_sv);
2107	100	T	F	}
	50	T	F	}
2110	0	T	F	for(size_t i=0; i<num_rows; i++) {
2112	0	T	F	if (row_ptr && SvROK(*row_ptr)) {
	0	T	F	if (row_ptr && SvROK(*row_ptr)) {
2114	0	T	F	hv_iterinit(row_hv);
2124	100	T	F	HE *restrict ce = hv_iternext(col_map);
2125	100	T	F	col_array[i] = SvPV_nolen(hv_iterkeysv(ce));
2127	50	T	F	qsort(col_array, num_cols, sizeof(char*), cmp_string_wt);
	50	T	F	qsort(col_array, num_cols, sizeof(char*), cmp_string_wt);
2132	100	T	F	if (inc_rownames && contains_nondigit(row_names_sv)) {
2135	50	T	F	for(size_t i=0; i<=av_len(headers_av); i++) {
	50	T	F	for(size_t i=0; i<=av_len(headers_av); i++) {
2136	100	T	F	SV**restrict h_ptr = av_fetch(headers_av, i, 0);
2137	50	T	F	if (!h_ptr \|\| !*h_ptr) continue;
2138	0	T	F	SV restrict h_sv = h_ptr;
2139	0	T	F	if (strcmp(SvPV_nolen(h_sv), rownames_col) != 0) {
	0	T	F	if (strcmp(SvPV_nolen(h_sv), rownames_col) != 0) {
2140	0	T	F	av_push(filtered_headers, newSVsv(h_sv));
2143	0	T	F	SvREFCNT_dec(headers_av);
2157	100	T	F	for(size_t i=0; i<num_rows; i++) {
2159	50	T	F	SV **restrict row_ptr = av_fetch(data_av, i, 0);
	50	T	F	SV **restrict row_ptr = av_fetch(data_av, i, 0);
2160	50	T	F	HV restrict row_hv = (row_ptr && SvROK(row_ptr)) ? (HV)SvRV(row_ptr) : NULL;
2161	100	T	F	if (inc_rownames) {
	50	T	F	if (inc_rownames) {
2162	50	T	F	if (rownames_col) {
2165	0	T	F	if (SvROK(*rn_val_ptr)) {
2174	100	T	F	}
	50	T	F	}
2178	50	T	F	row_data[d_idx++] = savepv(buf);
2179	100	T	F	}
2196	50	T	F	}
	50	T	F	}
	50	T	F	}
2201	50	T	F	safefree(row_data);
2202	50	T	F	}
2205	50	T	F	PerlIO_close(fh);
2210	100	T	F	_parse_csv_file(char* file, const char* sep_str, const char* comment_str, SV* callback = &PL_sv_undef)
2214	50	T	F	AV *restrict current_row = newAV();
	100	T	F	AV *restrict current_row = newAV();
2216	50	T	F	bool in_quotes = 0, post_quote = 0;
	100	T	F	bool in_quotes = 0, post_quote = 0;
2220	100	T	F	CODE:
2223	50	T	F	} else {
2224	50	T	F	data = newAV();
	100	T	F	data = newAV();
2226	50	T	F	sep_len = sep_str ? strlen(sep_str) : 0;
2229	50	T	F	fp = PerlIO_open(file, "r");
	50	T	F	fp = PerlIO_open(file, "r");
	50	T	F	fp = PerlIO_open(file, "r");
2234	100	T	F	// Read line by line using PerlIO
2236	100	T	F	char *restrict line = SvPV_nolen(line_sv);
2237	100	T	F	size_t len = SvCUR(line_sv);
2238	100	T	F	// chomp \r\n (Handles Windows invisible \r natively)
	100	T	F	// chomp \r\n (Handles Windows invisible \r natively)
	100	T	F	// chomp \r\n (Handles Windows invisible \r natively)
2241	100	T	F	if (len > 0 && line[len-1] == '\r') {
2244	100	T	F	}
2247	100	T	F	bool is_empty = 1;
	50	T	F	bool is_empty = 1;
	50	T	F	bool is_empty = 1;
	100	T	F	bool is_empty = 1;
2256	100	T	F	}
2265	50	T	F	i++; // Skip the escaped second quote
2269	50	T	F	} else if (!post_quote) {
2270	50	T	F	in_quotes = 1; // Open quotes (only when not in post-quote state)
2273	50	T	F	av_push(current_row, newSVsv(field));
2285	100	T	F	post_quote = 0; // Reset post-quote state at row boundary
2287	50	T	F	av_push(current_row, newSVsv(field));
2291	50	T	F	dSP;
2292	50	T	F	ENTER;
2295	50	T	F	XPUSHs(sv_2mortal(newRV_inc((SV*)current_row)));
2306	50	T	F	}
2318	50	T	F	PUTBACK;
	50	T	F	PUTBACK;
2321	50	T	F	LEAVE;
	50	T	F	LEAVE;
2326	100	T	F	current_row = newAV();
2327	100	T	F	}
2328	50	T	F	SvREFCNT_dec(field);
2337	50	T	F
2348	100	T	F
2353	50	T	F	croak("cov: unknown method '%s' (use 'pearson', 'spearman', or 'kendall')", method);
	50	T	F	croak("cov: unknown method '%s' (use 'pearson', 'spearman', or 'kendall')", method);
	50	T	F	croak("cov: unknown method '%s' (use 'pearson', 'spearman', or 'kendall')", method);
2354	50	T	F	}
	50	T	F	}
	50	T	F	}
2357	50	T	F	AV restrict y_av = (AV)SvRV(y_sv);
	50	T	F	AV restrict y_av = (AV)SvRV(y_sv);
2365	50	T	F
2373	100	T	F	SV **restrict x_tv = av_fetch(x_av, i, 0);
2375	100	T	F
2376	100	T	F	// Extract numeric values, defaulting to NAN for missing/invalid data
2385	100	T	F	}
2394	100	T	F	// 5. Algorithm routing
2406	100	T	F	if (strcmp(method, "spearman") == 0) {
2462	50	T	F	size_t n = 0, valid_n = 0, i;
2464	100	T	F	int iter = 0, max_iter = 25, final_rank = 0, df_res = 0;
2467	100	T	F
2468	100	T	F	char **restrict row_names = NULL;
2469	50	T	F	char **restrict valid_row_names = NULL;
2472	50	T	F	SV *restrict ref = NULL;
2473	50	T	F
	50	T	F
2477	100	T	F	double restrict XtWX = NULL, restrict XtWZ = NULL;
	50	T	F	double restrict XtWX = NULL, restrict XtWZ = NULL;
2485	100	T	F	for (unsigned short i_arg = 0; i_arg < items; i_arg += 2) {
	100	T	F	for (unsigned short i_arg = 0; i_arg < items; i_arg += 2) {
	50	T	F	for (unsigned short i_arg = 0; i_arg < items; i_arg += 2) {
2489	50	T	F	else if (strEQ(key, "data")) data_sv = val;
2493	50	T	F	if (!formula) croak("glm: formula is required");
2494	50	T	F	if (!data_sv \|\| !SvROK(data_sv)) croak("glm: data is required and must be a reference");
2497	100	T	F	bool is_gaussian = (strcmp(family_str, "gaussian") == 0);
2498	50	T	F	if (!is_binomial && !is_gaussian) croak("glm: unsupported family '%s'", family_str);
2502	50	T	F	Newx(exp_terms, exp_cap, char*); Newx(is_dummy, exp_cap, bool);
	50	T	F	Newx(exp_terms, exp_cap, char*); Newx(is_dummy, exp_cap, bool);
2507	50	T	F	*dst = '\0';
2510	0	T	F	if (!tilde) croak("glm: invalid formula, missing '~'");
	0	T	F	if (!tilde) croak("glm: invalid formula, missing '~'");
2511	0	T	F	*tilde = '\0';
	0	T	F	*tilde = '\0';
2520	50	T	F	term_cap *= 2;
	0	T	F	term_cap *= 2;
2526	100	T	F	}
2528	100	T	F	if (star) {
2529	50	T	F	*star = '\0';
2531	50	T	F	char restrict c_l = strchr(left, '^'); if (c_l && strncmp(left, "I(", 2) != 0) c_l = '\0';
2537	50	T	F	terms[num_terms] = (char*)safemalloc(inter_len);
2539	50	T	F	} else {
2541	50	T	F	if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
2543	50	T	F	}
	100	T	F	}
2546	50	T	F
2547	100	T	F	for (i = 0; i < num_terms; i++) {
2551	50	T	F	}
	50	T	F	}
2553	50	T	F	}
	50	T	F	}
2555	100	T	F
2563	0	T	F	SV*restrict val = hv_iterval(hv, entry);
2566	0	T	F	n = av_len((AV*)SvRV(val)) + 1;
	0	T	F	n = av_len((AV*)SvRV(val)) + 1;
2567	0	T	F	Newx(row_names, n, char*);
2569	0	T	F	char buf[32]; snprintf(buf, sizeof(buf), "%lu", i+1);
	0	T	F	char buf[32]; snprintf(buf, sizeof(buf), "%lu", i+1);
	0	T	F	char buf[32]; snprintf(buf, sizeof(buf), "%lu", i+1);
2574	0	T	F	Newx(row_names, n, char); Newx(row_hashes, n, HV);
2582	100	T	F	} else croak("glm: Hash values must be ArrayRefs (HoA) or HashRefs (HoH)");
2583	50	T	F	}
2588	100	T	F	for (i = 0; i < n; i++) {
2591	100	T	F	row_hashes[i] = (HV)SvRV(val);
2593	50	T	F	row_names[i] = savepv(buf);
2594	100	T	F	} else {
2596	50	T	F	Safefree(row_names); Safefree(row_hashes);
2598	100	T	F	}
2599	100	T	F	}
2601	100	T	F
2602	50	T	F	// --- Categorical Expansion ---
	0	T	F	// --- Categorical Expansion ---
2608	50	T	F	}
2609	100	T	F	if (strcmp(uniq_terms[j], "Intercept") == 0) {
2610	100	T	F	exp_terms[p_exp] = savepv("Intercept"); is_dummy[p_exp] = false; p_exp++; continue;
2611	50	T	F	}
2616	100	T	F	char*restrict str_val = get_data_string_alloc(data_hoa, row_hashes, i, uniq_terms[j]);
2617	50	T	F	if (str_val) {
2628	100	T	F	}
2639	50	T	F	exp_cap *= 2;
	50	T	F	exp_cap *= 2;
2640	50	T	F	Renew(exp_terms, exp_cap, char*); Renew(is_dummy, exp_cap, bool);
2643	100	T	F	size_t t_len = strlen(uniq_terms[j]) + strlen(levels[l]) + 1;
2645	50	T	F	snprintf(exp_terms[p_exp], t_len, "%s%s", uniq_terms[j], levels[l]);
2649	100	T	F	for (size_t l = 0; l < num_levels; l++) Safefree(levels[l]);
2650	100	T	F	Safefree(levels);
2652	100	T	F	Safefree(levels); exp_terms[p_exp] = savepv(uniq_terms[j]); is_dummy[p_exp] = false; p_exp++;
2654	50	T	F	} else {
2655	100	T	F	exp_terms[p_exp] = savepv(uniq_terms[j]); is_dummy[p_exp] = false; p_exp++;
2660	50	T	F	Newx(X, n * p, double); Newx(Y, n, double);
2663	50	T	F	// --- Listwise Deletion ---
2665	100	T	F	double y_val = evaluate_term(data_hoa, row_hashes, i, lhs);
2671	50	T	F	if (strcmp(exp_terms[j], "Intercept") == 0) {
2672	0	T	F	row_x[j] = 1.0;
2681	100	T	F	if (isnan(row_x[j])) { row_ok = false; break; }
2684	100	T	F	if (!row_ok) { Safefree(row_names[i]); Safefree(row_x); continue; }
2686	100	T	F	for (size_t j = 0; j < p; j++) X[valid_n * p + j] = row_x[j];
2687	100	T	F	valid_row_names[valid_n] = row_names[i];
2688	50	T	F	valid_n++;
	50	T	F	valid_n++;
2692	100	T	F	if (valid_n <= p) {
2693	50	T	F	Safefree(X); Safefree(Y); Safefree(valid_row_names); if (row_hashes) Safefree(row_hashes);
2702	50	T	F	for (i = 0; i < p; i++) { beta[i] = 0.0; beta_old[i] = 0.0; }
2703	100	T	F	// Initialize (mustart / etastart equivalent)
2704	100	T	F	double sum_y = 0.0;
2707	50	T	F	for (i = 0; i < valid_n; i++) {
2716	100	T	F	deviance_old += dev;
	100	T	F	deviance_old += dev;
2717	100	T	F	} else {
2719	100	T	F	eta[i] = mu[i];
2722	100	T	F	// IRLS Loop
2726	100	T	F	double varmu = mu[i] * (1.0 - mu[i]);
2727	50	T	F	double mu_eta = varmu; // Link derivative for logit
2729	50	T	F	Z[i] = eta[i] + (Y[i] - mu[i]) / mu_eta;
	100	T	F	Z[i] = eta[i] + (Y[i] - mu[i]) / mu_eta;
2735	100	T	F	}
2737	100	T	F	for (i = 0; i < p; i++) { XtWZ[i] = 0.0; for (size_t j = 0; j < p; j++) XtWX[i * p + j] = 0.0; }
2739	50	T	F	double w = W[k], z = Z[k];
	100	T	F	double w = W[k], z = Z[k];
2741	100	T	F	XtWZ[i] += X[k * p + i] * w * z;
2744	50	T	F	}
2745	50	T	F	}
2748	100	T	F	if (aliased[i]) { beta[i] = NAN; } else {
2749	50	T	F	double sum = 0.0;
2759	100	T	F	double linear_pred = 0.0;
	100	T	F	double linear_pred = 0.0;
	50	T	F	double linear_pred = 0.0;
2764	100	T	F	// Boundary enforcement
2767	100	T	F
2771	100	T	F	else dev = 2.0 * (Y[i] * log(Y[i] / mu[i]) + (1.0 - Y[i]) * log((1.0 - Y[i]) / (1.0 - mu[i])));
2774	100	T	F	mu[i] = eta[i];
	100	T	F	mu[i] = eta[i];
2775	100	T	F	double res = Y[i] - mu[i];
2776	100	T	F	deviance_new += res * res;
2777	50	T	F	}
2778	100	T	F	}
2780	100	T	F	if (!is_binomial \|\| deviance_new <= deviance_old + 1e-7 \|\| !isfinite(deviance_new)) {
2786	100	T	F	}
2787	100	T	F	// Convergence Check
2788	100	T	F	if (fabs(deviance_new - deviance_old) / (0.1 + fabs(deviance_new)) < epsilon) {
2789	50	T	F	converged = true; break;
2797	100	T	F	double w = is_binomial ? (mu[k] * (1.0 - mu[k])) : 1.0;
2800	50	T	F	double xw = X[k * p + i] * w;
2806	100	T	F	double wtdmu = mean_y; // Since weights are 1.0 initially
	50	T	F	double wtdmu = mean_y; // Since weights are 1.0 initially
2807	100	T	F	for (i = 0; i < valid_n; i++) {
2809	100	T	F	if (Y[i] == 0.0) null_dev += -2.0 * log(1.0 - wtdmu);
2812	100	T	F	} else {
2813	50	T	F	double diff = Y[i] - wtdmu;
2815	100	T	F	}
2824	100	T	F	// --- Return Structures ---
2829	50	T	F	double res = Y[i] - mu[i];
2832	0	T	F	double d_res = 0.0;
2833	0	T	F	if (Y[i] == 0.0) d_res = sqrt(-2.0 * log(1.0 - mu[i]));
2837	100	T	F	}
2841	100	T	F	}
2842	100	T	F	Safefree(valid_row_names);
2864	100	T	F	}
2866	100	T	F	}
2868	100	T	F	hv_store(res_hv, "aic", 3, newSVnv(aic), 0);
2870	100	T	F	hv_store(res_hv, "converged", 9, newSVuv(converged ? 1 : 0), 0);
2877	100	T	F	hv_store(res_hv, "fitted.values", 13, newRV_noinc((SV*)fitted_hv), 0);
2887	50	T	F	for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]);
	50	T	F	for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]);
2899	100	T	F
2903	100	T	F	RETVAL
2904	100	T	F
2905	50	T	F	SV* cor_test(...)
2906	100	T	F	CODE:
	50	T	F	CODE:
2907	50	T	F	{
2920	50	T	F	for (unsigned short int i = 2; i < items; i += 2) {
	50	T	F	for (unsigned short int i = 2; i < items; i += 2) {
	50	T	F	for (unsigned short int i = 2; i < items; i += 2) {
2921	50	T	F	const char *restrict key = SvPV_nolen(ST(i));
	50	T	F	const char *restrict key = SvPV_nolen(ST(i));
	50	T	F	const char *restrict key = SvPV_nolen(ST(i));
2929	50	T	F	else croak("cor_test: unknown argument '%s'", key);
2935	100	T	F
2939	50	T	F	HV *restrict rhv;
	100	T	F	HV *restrict rhv;
	50	T	F	HV *restrict rhv;
2940	50	T	F
	100	T	F
	50	T	F
2943	100	T	F	croak("cor_test: x and y must be array references");
	100	T	F	croak("cor_test: x and y must be array references");
2950	50	T	F	if (n_raw != av_len(y_av) + 1) croak("incompatible dimensions");
2956	100	T	F	for (size_t i = 0; i < n_raw; i++) {
2959	100	T	F
2968	50	T	F	}
	50	T	F	}
2982	100	T	F	mean_x += dx / (i + 1);
2984	100	T	F	mean_y += dy / (i + 1);
2985	100	T	F	M2_x += dx * (x[i] - mean_x);
2989	50	T	F	estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
	0	T	F	estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
2990	50	T	F	df = n - 2;
2991	50	T	F	statistic = estimate * sqrt(df / (1.0 - estimate * estimate));
2992	100	T	F
2997	50	T	F	double q = inverse_normal_cdf(1.0 - alpha/2.0);
2999	50	T	F	ci_upper = tanh(z + q * se);
	50	T	F	ci_upper = tanh(z + q * se);
3003	50	T	F	} else if (is_kendall) {
	0	T	F	} else if (is_kendall) {
3004	50	T	F	int c = 0, d = 0, tie_x = 0, tie_y = 0;
	50	T	F	int c = 0, d = 0, tie_x = 0, tie_y = 0;
3009	50	T	F
	50	T	F
3011	50	T	F	else if (sign_x == 0) tie_x++;
3019	0	T	F
	0	T	F
3022	0	T	F
3024	0	T	F	if (!exact_sv \|\| !SvOK(exact_sv)) {
3030	50	T	F	if (do_exact && has_ties) do_exact = 0;
3038	100	T	F	double var_S = n * (n - 1) * (2.0 * n + 5.0) / 18.0;
3047	50	T	F	} else {
	50	T	F	} else {
3051	100	T	F	} else if (is_spearman) {
3058	100	T	F	double mean_x = 0.0, mean_y = 0.0, M2_x = 0.0, M2_y = 0.0, cov = 0.0;
3059	50	T	F	for (size_t i = 0; i < n; i++) {
	50	T	F	for (size_t i = 0; i < n; i++) {
3064	50	T	F	M2_x += dx * (rank_x[i] - mean_x);
	0	T	F	M2_x += dx * (rank_x[i] - mean_x);
3065	50	T	F	M2_y += dy * (rank_y[i] - mean_y);
	50	T	F	M2_y += dy * (rank_y[i] - mean_y);
3070	50	T	F	// S = sum of squared rank differences (R's reported statistic)
3075	0	T	F	}
3092	100	T	F	statistic = S_stat;
3114	50	T	F	hv_stores(rhv, "parameter", newSVnv(df));
	50	T	F	hv_stores(rhv, "parameter", newSVnv(df));
3121	50	T	F	RETVAL = newRV_noinc((SV*)rhv);
3124	100	T	F	RETVAL
3126	50	T	F	void shapiro_test(data)
	50	T	F	void shapiro_test(data)
3128	50	T	F	PREINIT:
3136	50	T	F	}
	50	T	F	}
3143	100	T	F	// Extract variables and calculate mean (skipping undefined/NaN values)
3146	50	T	F	if (elem && SvOK(*elem)) {
3153	50	T	F	}
3157	0	T	F	Safefree(x);
3163	50	T	F	for (size_t i = 0; i < n; i++) {
3164	50	T	F	ssq += (x[i] - mean) * (x[i] - mean);
3165	100	T	F	}
3173	50	T	F	if (n == 3) {
	100	T	F	if (n == 3) {
3175	100	T	F	double b_val = a_val * (x[2] - x[0]);
3183	100	T	F	Newx(m, n, double);
3187	100	T	F	sum_m2 += m[i] * m[i];
3197	100	T	F	}
3203	50	T	F	for (unsigned int i = 2; i < n-2; i++) {
3226	50	T	F	// Horner-form polynomials in n for mu and log(sigma)
3227	50	T	F	double mu = 0.544 + nn * (-0.39978 + nn * ( 0.025054 - nn * 0.0006714));
3237	50	T	F	double ln_n = log((double)n);
3247	100	T	F	if (p_val < 0.0) p_val = 0.0;
3249	100	T	F	Safefree(m); m = NULL; Safefree(a); a = NULL;
	50	T	F	Safefree(m); m = NULL; Safefree(a); a = NULL;
3252	100	T	F	ret_hash = newHV();
3254	50	T	F	hv_stores(ret_hash, "W", newSVnv(w));
	50	T	F	hv_stores(ret_hash, "W", newSVnv(w));
3256	100	T	F	hv_stores(ret_hash, "p.value", newSVnv(p_val));
	50	T	F	hv_stores(ret_hash, "p.value", newSVnv(p_val));
3263	50	T	F	double min_val = 0.0;
3265	100	T	F	bool first = TRUE;
	100	T	F	bool first = TRUE;
3272	100	T	F	for (size_t j = 0; j < len; j++) {
3273	50	T	F	SV** restrict tv = av_fetch(av, j, 0);
3284	100	T	F	}
3286	100	T	F	double val = SvNV(arg);
	50	T	F	double val = SvNV(arg);
3289	100	T	F	first = FALSE;
3291	50	T	F	count++;
	50	T	F	count++;
3293	100	T	F	croak("min: undefined value at argument index %d", (int)i);
	100	T	F	croak("min: undefined value at argument index %d", (int)i);
3300	50	T	F
3302	100	T	F	PROTOTYPE: @
	100	T	F	PROTOTYPE: @
3309	100	T	F	SV* restrict arg = ST(i);
3310	50	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3325	50	T	F	}
3329	100	T	F	max_val = val;
3331	100	T	F	}
	100	T	F	}
3333	100	T	F	} else {
3338	100	T	F	RETVAL = max_val;
3343	50	T	F	CODE:
3352	100	T	F
3355	100	T	F	}
3358	50	T	F	// 1. Check if the current argument is a string key for a named parameter
3366	50	T	F	} else if (strEQ(key, "min")) {
3370	50	T	F	continue;
3372	50	T	F	max = SvNV(ST(i+1));
3376	100	T	F	}
3378	50	T	F
3394	50	T	F	if (!n_set) {
3395	100	T	F	croak("runif() requires at least the 'n' parameter");
3403	100	T	F	const double range = max - min;
3407	100	T	F	r = NAN; // R behavior for inverted ranges
3408	100	T	F	} else {
3409	50	T	F	r = min + range * Drand01();
3414	100	T	F	}
	100	T	F	}
3415	100	T	F	OUTPUT:
	100	T	F	OUTPUT:
3418	50	T	F	SV* rbinom(...)
3420	100	T	F	{
3435	100	T	F	if (strEQ(key, "n")) n = (unsigned int)SvUV(val);
	100	T	F	if (strEQ(key, "n")) n = (unsigned int)SvUV(val);
3440	100	T	F
3444	50	T	F
3448	100	T	F	for (unsigned int i = 0; i < n; i++) {
3450	50	T	F	}
	50	T	F	}
3453	100	T	F	RETVAL = newRV_noinc((SV*)result_av);
3454	100	T	F	}
3457	50	T	F
3464	50	T	F	croak("hist: first argument must be an array reference");
3467	50	T	F	size_t n_raw = av_len(x_av) + 1;
3469	50	T	F
3471	50	T	F	double *restrict x;
	50	T	F	double *restrict x;
3477	50	T	F	SV**restrict tv = av_fetch(x_av, i, 0);
	0	T	F	SV**restrict tv = av_fetch(x_av, i, 0);
3481	50	T	F	if (val < min_val) min_val = val;
3485	50	T	F	if (n == 0) {
3486	50	T	F	Safefree(x);
3487	50	T	F	croak("hist: input contains no valid numeric data");
3488	50	T	F	}
3492	100	T	F	if (items == 2) {
3505	100	T	F	if (n_bins == 0 && looks_like_number(ST(1))) {
3507	100	T	F	}
3535	50	T	F	if (i < n_bins) {
	100	T	F	if (i < n_bins) {
	50	T	F	if (i < n_bins) {
3541	50	T	F	hv_stores(res_hv, "breaks", newRV_noinc((SV*)av_breaks));
3544	100	T	F	hv_stores(res_hv, "density", newRV_noinc((SV*)av_density));
3548	100	T	F	Safefree(density); Safefree(counts);
3549	50	T	F
3552	50	T	F	OUTPUT:
	50	T	F	OUTPUT:
	50	T	F	OUTPUT:
3556	50	T	F	CODE:
3560	50	T	F	int arg_idx = 0;
3562	100	T	F	/* --- 1. Consume first positional arg as 'x' if it's an array ref --- */
3564	50	T	F	x_sv = ST(arg_idx);
	50	T	F	x_sv = ST(arg_idx);
3568	50	T	F	/* --- 2. Remaining args must be key-value pairs --- */
3579	50	T	F	}
	50	T	F	}
	50	T	F	}
3583	100	T	F	size_t n_raw = av_len(x_av) + 1;
3585	50	T	F
	50	T	F
3586	50	T	F	/* --- Extract valid numeric data & drop NAs --- */
	50	T	F	/* --- Extract valid numeric data & drop NAs --- */
3593	0	T	F	x[n++] = SvNV(*tv);
3599	100	T	F	}
3602	100	T	F	// --- Parse Probabilities (Default matches R's c(0, .25, .5, .75, 1)) ---
3604	100	T	F	unsigned int n_probs = 5;
3606	100	T	F
3620	50	T	F	Newx(probs, n_probs, double);
3644	100	T	F	/* Format hash key to exactly match R's naming convention ("25%", "33.3%") */
3646	100	T	F	double pct = p * 100.0;
	50	T	F	double pct = p * 100.0;
3649	100	T	F	snprintf(key, sizeof(key), "%.0f%%", pct);
3651	50	T	F	snprintf(key, sizeof(key), "%.1f%%", pct);
	50	T	F	snprintf(key, sizeof(key), "%.1f%%", pct);
3653	50	T	F
3657	100	T	F	Safefree(x);
3658	100	T	F	Safefree(probs);
3669	100	T	F	double total = 0;
3671	100	T	F	CODE:
	50	T	F	CODE:
3674	100	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3676	50	T	F	size_t len = av_len(av) + 1;
	50	T	F	size_t len = av_len(av) + 1;
3684	50	T	F	}
3692	100	T	F	}
3693	50	T	F	if (count == 0) croak("mean needs >= 1 element");
3704	100	T	F	for (size_t i = 0; i < items; i++) {
3706	100	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
	50	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3709	100	T	F	for (size_t j = 0; j < len; j++) {
3711	50	T	F	if (tv && SvOK(*tv)) {
	50	T	F	if (tv && SvOK(*tv)) {
3719	50	T	F	total += SvNV(arg);
3727	100	T	F	OUTPUT:
3728	100	T	F	RETVAL
3744	50	T	F	if (tv && SvOK(*tv)) {
	100	T	F	if (tv && SvOK(*tv)) {
	50	T	F	if (tv && SvOK(*tv)) {
3750	50	T	F	} else {
	100	T	F	} else {
	50	T	F	} else {
3756	50	T	F	double val = SvNV(arg);
3761	100	T	F	croak("sd: undefined value at argument index %zu", i);
3765	100	T	F	RETVAL = sqrt(M2 / (count - 1));
3766	100	T	F	OUTPUT:
3767	100	T	F	RETVAL
3768	100	T	F
3769	100	T	F
3770	100	T	F	double var(...)
3771	50	T	F	PROTOTYPE: @
3776	100	T	F	/* Single Pass Variance via Welford's Algorithm */
	50	T	F	/* Single Pass Variance via Welford's Algorithm */
	50	T	F	/* Single Pass Variance via Welford's Algorithm */
3780	50	T	F	AV* restrict av = (AV*)SvRV(arg);
3782	100	T	F	for (size_t j = 0; j < len; j++) {
	50	T	F	for (size_t j = 0; j < len; j++) {
	50	T	F	for (size_t j = 0; j < len; j++) {
3785	50	T	F	count++;
	100	T	F	count++;
3790	100	T	F	} else {
3792	50	T	F	}
	50	T	F	}
3798	100	T	F	mean += delta / count;
	50	T	F	mean += delta / count;
3800	100	T	F	} else {
	100	T	F	} else {
3801	100	T	F	croak("var: undefined value at argument index %zu", i);
3803	100	T	F	}
	100	T	F	}
3805	100	T	F	RETVAL = M2 / (count - 1);
3807	50	T	F	RETVAL
	50	T	F	RETVAL
3813	100	T	F	SV*restrict y_sv = NULL;
3815	100	T	F	bool paired = FALSE, var_equal = FALSE;
3818	50	T	F	int arg_idx = 0;
	50	T	F	int arg_idx = 0;
3819	50	T	F
	50	T	F
3826	50	T	F	// 2. Shift second positional argument as 'y' if it's an array reference
3832	100	T	F	// Ensure the remaining arguments form complete key-value pairs
3833	50	T	F	if ((items - arg_idx) % 2 != 0) {
	0	T	F	if ((items - arg_idx) % 2 != 0) {
3842	50	T	F	if (strEQ(key, "x")) x_sv = val;
	0	T	F	if (strEQ(key, "x")) x_sv = val;
3862	100	T	F	if (conf_level <= 0.0 \|\| conf_level >= 1.0)
3866	100	T	F	HV*restrict results = newHV();
3889	100	T	F	var_y = M2_y / (ny - 1);
	50	T	F	var_y = M2_y / (ny - 1);
3895	100	T	F	double dx = (dx_ptr && SvOK(dx_ptr)) ? SvNV(dx_ptr) : 0.0;
3901	100	T	F	}
3903	100	T	F	if (var_d == 0.0) croak("t_test: data are essentially constant");
	100	T	F	if (var_d == 0.0) croak("t_test: data are essentially constant");
3904	100	T	F	cint_est = mean_d;
	50	T	F	cint_est = mean_d;
3905	50	T	F	std_err = sqrt(var_d / nx);
3909	50	T	F	} else if (var_equal) {
3910	50	T	F	if (var_x == 0.0 && var_y == 0.0) croak("t_test: data are essentially constant");
3912	100	T	F	cint_est = mean_x - mean_y;
3914	50	T	F	t_stat = (cint_est - mu) / std_err;
	50	T	F	t_stat = (cint_est - mu) / std_err;
3920	100	T	F	cint_est = mean_x - mean_y;
3921	100	T	F	double stderr_x2 = var_x / nx;
3923	100	T	F	std_err = sqrt(stderr_x2 + stderr_y2);
3925	100	T	F	df = pow(stderr_x2 + stderr_y2, 2) /
3927	100	T	F	hv_store(results, "estimate_x", 10, newSVnv(mean_x), 0);
3929	100	T	F	}
3930	100	T	F	} else {
3932	100	T	F	std_err = sqrt(var_x / nx);
3934	100	T	F	df = nx - 1;
3936	100	T	F	}
3937	50	T	F	p_val = get_t_pvalue(t_stat, df, alternative);
3939	100	T	F	if (strcmp(alternative, "less") == 0) {
3941	100	T	F	ci_lower = -INFINITY;
3943	100	T	F	} else if (strcmp(alternative, "greater") == 0) {
3944	50	T	F	t_crit = qt_tail(df, alpha);
3946	100	T	F	ci_upper = INFINITY;
3948	100	T	F	t_crit = qt_tail(df, alpha / 2.0);
3950	100	T	F	ci_upper = cint_est + t_crit * std_err;
3952	100	T	F	AV*restrict conf_int = newAV();
3953	100	T	F	av_push(conf_int, newSVnv(ci_lower));
3955	100	T	F	hv_store(results, "statistic", 9, newSVnv(t_stat), 0);
3957	50	T	F	hv_store(results, "p_value", 7, newSVnv(p_val), 0);
3958	50	T	F	hv_store(results, "conf_int", 8, newRV_noinc((SV*)conf_int), 0);
3961	100	T	F	OUTPUT:
3963	50	T	F
3968	100	T	F	}
3972	100	T	F	if (n == 0) {
3977	100	T	F	strncpy(meth, method, 63); meth[63] = '\0';
3979	100	T	F	// Resolve aliases
3984	100	T	F	PVal *restrict arr;
3986	100	T	F	Newx(arr, n, PVal);
3989	100	T	F	for (size_t i = 0; i < n; i++) {
3993	100	T	F	}
3994	100	T	F	// Sort ascending (Stable sort using original index)
4000	100	T	F	adj[arr[i].orig_idx] = (v < 1.0) ? v : 1.0;
4001	100	T	F	}
4002	50	T	F	} else if (strcmp(meth, "holm") == 0) {
4006	50	T	F	if (v > cummax) cummax = v;
4007	0	T	F	adj[arr[i].orig_idx] = (cummax < 1.0) ? cummax : 1.0;
4015	50	T	F	}
4016	100	T	F	} else if (strcmp(meth, "bh") == 0) {
4030	100	T	F	adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0;
4032	100	T	F	} else if (strcmp(meth, "hommel") == 0) {
	50	T	F	} else if (strcmp(meth, "hommel") == 0) {
4035	100	T	F	Newx(q_arr, n, double);
4037	50	T	F	double min_val = n * arr[0].p;
	50	T	F	double min_val = n * arr[0].p;
4039	50	T	F	double temp = (n * arr[i].p) / (i + 1.0);
4043	100	T	F	}
4045	50	T	F	for (size_t i = 0; i < n; i++) {
4047	100	T	F	q_arr[i] = min_val;
4049	100	T	F	for (size_t j = n - 1; j >= 2; j--) {
	50	T	F	for (size_t j = n - 1; j >= 2; j--) {
4052	100	T	F	// Calculate q1 = min(j * p[i2] / (2:j))
4054	50	T	F	for (size_t k = 1; k < i2_len; k++) {
	50	T	F	for (size_t k = 1; k < i2_len; k++) {
4058	50	T	F	}
4064	100	T	F	}
4070	100	T	F	for (size_t i = 0; i < n; i++) {
4077	100	T	F	for (size_t i = 0; i < n; i++) {
4078	100	T	F	double v = (pa[i] > arr[i].p) ? pa[i] : arr[i].p;
4079	100	T	F	if (v > 1.0) v = 1.0;
4084	50	T	F	for (size_t i = 0; i < n; i++) {
	50	T	F	for (size_t i = 0; i < n; i++) {
4089	50	T	F	croak("Unknown p-value adjustment method: %s", method);
4095	50	T	F	}
	100	T	F	}
	50	T	F	}
4100	50	T	F	PROTOTYPE: @
	50	T	F	PROTOTYPE: @
4101	50	T	F	INIT:
4103	50	T	F	double* restrict nums;
4104	50	T	F	double median_val = 0.0;
4107	50	T	F	for (size_t i = 0; i < items; i++) {
	50	T	F	for (size_t i = 0; i < items; i++) {
4109	50	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
	100	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
	50	T	F	if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
4115	100	T	F	total_count++;
	50	T	F	total_count++;
4116	50	T	F	} else {
4120	100	T	F	} else if (SvOK(arg)) {
4124	50	T	F	}
4128	50	T	F	/* Allocate C array now that we know the exact size */
4129	50	T	F	Newx(nums, total_count, double);
4131	100	T	F	/* Pass 2: Populate the C array â€” Safefree before any croak */
4133	50	T	F	SV* restrict arg = ST(i);
	50	T	F	SV* restrict arg = ST(i);
	50	T	F	SV* restrict arg = ST(i);
4135	100	T	F	AV* restrict av = (AV*)SvRV(arg);
4137	50	T	F	for (size_t j = 0; j < len; j++) {
	50	T	F	for (size_t j = 0; j < len; j++) {
	50	T	F	for (size_t j = 0; j < len; j++) {
4146	50	T	F	} else if (SvOK(arg)) {
4151	50	T	F	}
	50	T	F	}
	50	T	F	}
4155	50	T	F	qsort(nums, total_count, sizeof(double), compare_doubles);
4160	100	T	F	}
4162	50	T	F	nums = NULL;
	50	T	F	nums = NULL;
	50	T	F	nums = NULL;
4166	50	T	F
	50	T	F
4167	50	T	F	SV* cor(SV* x_sv, SV* y_sv = &PL_sv_undef, const char* method = "pearson")
4168	100	T	F	INIT:
4170	50	T	F	if (strcmp(method, "pearson") != 0 &&
	50	T	F	if (strcmp(method, "pearson") != 0 &&
	50	T	F	if (strcmp(method, "pearson") != 0 &&
4177	50	T	F	if (!SvROK(x_sv) \|\| SvTYPE(SvRV(x_sv)) != SVt_PVAV)
4179	100	T	F
4180	50	T	F	AVrestrict x_av = (AV)SvRV(x_sv);
4181	100	T	F	size_t nx = av_len(x_av) + 1;
4185	50	T	F	bool x_is_matrix = 0;
	50	T	F	bool x_is_matrix = 0;
	50	T	F	bool x_is_matrix = 0;
4195	50	T	F
	50	T	F
4199	50	T	F	bool y_is_matrix = 0;
4201	50	T	F	SV**restrict fp = av_fetch(y_av, 0, 0);
4202	100	T	F	if (fp && SvROK(fp) && SvTYPE(SvRV(fp)) == SVt_PVAV)
4203	50	T	F	y_is_matrix = 1;
4204	100	T	F	}
4208	50	T	F	if (!x_is_matrix && !y_is_matrix) {
	50	T	F	if (!x_is_matrix && !y_is_matrix) {
	50	T	F	if (!x_is_matrix && !y_is_matrix) {
4216	50	T	F
4223	50	T	F
4224	100	T	F	for (size_t i = 0; i < nx; i++) {
4228	50	T	F	for (size_t i = 0; i < ny; i++) {
4231	0	T	F	}
4232	0	T	F
4233	0	T	F	double r = compute_cor(xd, yd, nx, method);
4235	0	T	F	RETVAL = newSVnv(r);
4237	0	T	F	} else {//Branch 2: x is a matrix (or y is a matrix) â†’ AoA result
4244	0	T	F	if (!xr0 \|\| !SvROK(xr0) \|\| SvTYPE(SvRV(xr0)) != SVt_PVAV)
4245	0	T	F	croak("cor: each row of x must be an ARRAY reference");
4248	0	T	F	if (ncols_x == 0) croak("cor: x matrix has zero columns");
4252	100	T	F	// PRE-VALIDATION PASS: Ensure all rows are arrays to prevent memory leaks on croak
4253	100	T	F	for (size_t i = 0; i < nrows; i++) {
4257	100	T	F	}
4261	100	T	F	for (size_t i = 0; i < nrows; i++) {
4263	50	T	F	if (!rv \|\| !SvROK(rv) \|\| SvTYPE(SvRV(rv)) != SVt_PVAV)
4264	100	T	F	croak("cor: y row %lu is not an array ref", i);
4280	50	T	F	}
4282	100	T	F	// -- resolve y: separate matrix or re-use x (symmetric)
	100	T	F	// -- resolve y: separate matrix or re-use x (symmetric)
4287	50	T	F	// 1 = cor(X) â€” result is symmetric
4289	50	T	F	// cross-correlation: X (nrows Ã— p) vs Y (nrows Ã— q)
4294	50	T	F	Newx(col_y, ncols_y, double*);
	50	T	F	Newx(col_y, ncols_y, double*);
	100	T	F	Newx(col_y, ncols_y, double*);
4296	50	T	F	Newx(col_y[j], nrows, double);
	50	T	F	Newx(col_y[j], nrows, double);
	50	T	F	Newx(col_y[j], nrows, double);
	50	T	F	Newx(col_y[j], nrows, double);
4298	0	T	F	SV**restrict rv = av_fetch(y_av, i, 0);
4300	0	T	F	SV**restrict cv = av_fetch(row, j, 0);
4309	100	T	F	if (nrows < 2)
4311	50	T	F	// -- build cache for symmetric case: compute upper triangle, store results, mirror to lower triangle
4315	50	T	F	AV **restrict rows_out;
	50	T	F	AV **restrict rows_out;
	50	T	F	AV **restrict rows_out;
4317	50	T	F	for (size_t i = 0; i < ncols_x; i++) {
	50	T	F	for (size_t i = 0; i < ncols_x; i++) {
	50	T	F	for (size_t i = 0; i < ncols_x; i++) {
	50	T	F	for (size_t i = 0; i < ncols_x; i++) {
4319	0	T	F	av_extend(rows_out[i], ncols_y - 1);
4321	0	T	F	if (symmetric) {
4322	0	T	F	/* Upper triangle + diagonal, then mirror. r_cache[i][j] (j >= i) holds the computed value. */
4333	100	T	F	r_cache[j][i] = r; // symmetry
4335	100	T	F	}
	50	T	F	}
4337	50	T	F	for (size_t i = 0; i < ncols_x; i++)
4339	50	T	F	av_store(rows_out[i], j, newSVnv(r_cache[i][j]));
	50	T	F	av_store(rows_out[i], j, newSVnv(r_cache[i][j]));
	50	T	F	av_store(rows_out[i], j, newSVnv(r_cache[i][j]));
4345	100	T	F	for (size_t i = 0; i < ncols_x; i++)
4355	50	T	F	Safefree(col_x); col_x = NULL;
	50	T	F	Safefree(col_x); col_x = NULL;
4361	100	T	F	}
4367	100	T	F	PPCODE:
4370	50	T	F	double center_val = 0.0, scale_val = 1.0;
4372	100	T	F	// 1. Parse Options Hash (if it exists as the last argument)
4374	50	T	F	SV*restrict last_arg = ST(items - 1);
	50	T	F	SV*restrict last_arg = ST(items - 1);
4377	50	T	F	HVrestrict opt_hv = (HV)SvRV(last_arg);
	50	T	F	HVrestrict opt_hv = (HV)SvRV(last_arg);
4384	50	T	F	} else {
4387	50	T	F	if (strcasecmp(str, "mean") == 0 \|\| strcasecmp(str, "true") == 0 \|\| strcmp(str, "1") == 0) {
4388	50	T	F	do_center_mean = true;
4394	100	T	F	do_center_mean = true;
4401	100	T	F	SV**restrict scale_sv = hv_fetch(opt_hv, "scale", 5, 0);
4403	50	T	F	SVrestrict val_sv = scale_sv;
4410	50	T	F	} else if (strcasecmp(str, "none") == 0 \|\| strcasecmp(str, "false") == 0 \|\| strcmp(str, "0") == 0 \|\| strcmp(str, "") == 0) {
4419	100	T	F	}
4421	50	T	F	}
	0	T	F	}
4424	0	T	F	// 2. Detect if the input is a Matrix (Array of Arrays)
4426	0	T	F	if (data_items == 1) {
	0	T	F	if (data_items == 1) {
4428	50	T	F	if (SvROK(first_arg) && SvTYPE(SvRV(first_arg)) == SVt_PVAV) {
4432	50	T	F	if (first_elem && SvROK(first_elem) && SvTYPE(SvRV(first_elem)) == SVt_PVAV) {
4433	50	T	F	is_matrix = true;
4434	100	T	F	}
4436	50	T	F	}
	0	T	F	}
4439	0	T	F	//=========================================================
4441	0	T	F	//=========================================================
	0	T	F	//=========================================================
4446	50	T	F	ncol = av_len((AV)SvRV(first_row)) + 1;
4452	100	T	F	av_extend(result_av, nrow - 1);
4454	100	T	F	for (size_t r = 0; r < nrow; r++) {
4455	100	T	F	row_ptrs[r] = newAV();
4460	100	T	F	for (size_t c = 0; c < ncol; c++) {
4466	50	T	F	SV**restrict row_sv = av_fetch(mat_av, r, 0);
4467	100	T	F	if (row_sv && SvROK(*row_sv)) {
4469	50	T	F	SV**restrict cell_sv = av_fetch(row_av, c, 0);
4479	50	T	F	// Calculate Standard Deviation for this specific column if needed
4486	100	T	F	double sum_sq = 0.0;
4489	100	T	F	sum_sq += diff * diff;
4491	100	T	F	col_scale = sqrt(sum_sq / (nrow - 1));
4494	50	T	F	for (size_t r = 0; r < nrow; r++) {
4497	0	T	F	av_store(row_ptrs[r], c, newSVnv(final_val));
4504	50	T	F	PUSHs(sv_2mortal(newRV_noinc((SV*)result_av)));
	100	T	F	PUSHs(sv_2mortal(newRV_noinc((SV*)result_av)));
	50	T	F	PUSHs(sv_2mortal(newRV_noinc((SV*)result_av)));
4508	50	T	F	// ======================================
4509	100	T	F	size_t total_count = 0, k = 0;
4513	50	T	F	SV*restrict arg = ST(i);
	0	T	F	SV*restrict arg = ST(i);
4516	50	T	F	size_t len = av_len(av) + 1;
	100	T	F	size_t len = av_len(av) + 1;
4518	50	T	F	SV**restrict tv = av_fetch(av, j, 0);
	0	T	F	SV**restrict tv = av_fetch(av, j, 0);
4522	100	T	F	total_count++;
	50	T	F	total_count++;
4530	100	T	F	AVrestrict av = (AV)SvRV(arg);
4537	100	T	F	}
4540	50	T	F	double val = SvNV(arg);
4541	50	T	F	nums[k++] = val; sum += val;
4584	50	T	F	nrow_set = TRUE;
4586	100	T	F	ncol = (size_t)SvUV(val);
4589	100	T	F	byrow = SvTRUE(val);
4590	50	T	F	} else {
4593	100	T	F	}
4594	100	T	F	// Validate data input
	100	T	F	// Validate data input
4600	50	T	F	if (data_len == 0) {
4602	50	T	F	}
4604	50	T	F	if (!nrow_set && !ncol_set) {
4606	50	T	F	ncol = 1;
	100	T	F	ncol = 1;
4609	50	T	F	} else if (!nrow_set && ncol_set) {
4610	100	T	F	nrow = (data_len + ncol - 1) / ncol;
4615	50	T	F	}
	50	T	F	}
4617	50	T	F	AV*restrict result_av = newAV();
	50	T	F	AV*restrict result_av = newAV();
4619	100	T	F	size_t r, c;// Use unsigned types for counters to prevent negative indexing
4627	0	T	F	size_t total_cells = nrow * ncol;
4629	0	T	F	// Vector recycling logic
4630	0	T	F	SV**restrict fetched = av_fetch(data_av, i % data_len, 0);
4631	0	T	F	SVrestrict val = fetched ? newSVsv(fetched) : newSV(0);
4633	0	T	F	r = i / ncol;
	0	T	F	r = i / ncol;
	0	T	F	r = i / ncol;
4638	0	T	F	}
4649	100	T	F	const char *restrict formula = NULL;
	100	T	F	const char *restrict formula = NULL;
	50	T	F	const char *restrict formula = NULL;
4653	100	T	F
4654	100	T	F	char restrict terms = NULL, restrict uniq_terms = NULL, **restrict exp_terms = NULL;
4655	50	T	F	bool *restrict is_dummy = NULL;
4667	100	T	F	bool *restrict aliased = NULL;
4669	50	T	F	int final_rank = 0, df_res = 0;
	0	T	F	int final_rank = 0, df_res = 0;
4671	0	T	F	AV *restrict terms_av;
	0	T	F	AV *restrict terms_av;
	0	T	F	AV *restrict terms_av;
	0	T	F	AV *restrict terms_av;
4675	100	T	F	if (items % 2 != 0) croak("Usage: lm(formula => 'mpg ~ wt * hp', data => \\%%mtcars)");
	50	T	F	if (items % 2 != 0) croak("Usage: lm(formula => 'mpg ~ wt * hp', data => \\%%mtcars)");
4676	50	T	F
	0	T	F
	0	T	F
4682	100	T	F	else croak("lm: unknown argument '%s'", key);
	50	T	F	else croak("lm: unknown argument '%s'", key);
4683	0	T	F	}
	0	T	F	}
	0	T	F	}
4689	100	T	F	// ========================================================================
	50	T	F	// ========================================================================
	0	T	F	// ========================================================================
4695	100	T	F	if (entry) {
	50	T	F	if (entry) {
	0	T	F	if (entry) {
4699	100	T	F	n = av_len((AV*)SvRV(val)) + 1;
	50	T	F	n = av_len((AV*)SvRV(val)) + 1;
4700	0	T	F	Newx(row_names, n, char*);
	0	T	F	Newx(row_names, n, char*);
	0	T	F	Newx(row_names, n, char*);
4705	100	T	F	}
4706	50	T	F	} else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
	0	T	F	} else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
4707	50	T	F	n = hv_iterinit(hv);
	0	T	F	n = hv_iterinit(hv);
4716	50	T	F	} else croak("lm: Hash values must be ArrayRefs (HoA) or HashRefs (HoH)");
4718	50	T	F	} else if (SvTYPE(ref) == SVt_PVAV) {
4720	50	T	F	Newx(row_names, n, char*);
	50	T	F	Newx(row_names, n, char*);
4727	100	T	F	row_names[i] = savepv(buf);
4728	100	T	F	} else {
4730	100	T	F	Safefree(row_names); Safefree(row_hashes);
4732	50	T	F	}
	50	T	F	}
4734	100	T	F	} else croak("lm: Data must be an Array or Hash reference");
4736	50	T	F	// ========================================================================
4737	100	T	F	// PHASE 2: Formula Parsing & `.` Expansion
4747	50	T	F	croak("lm: invalid formula, missing '~'");
4748	100	T	F	}
4760	100	T	F	if (p_idx[0] == 'I' && p_idx[1] == '(') {
4762	50	T	F	while (p_idx) { if (p_idx == '(') depth++; else if (*p_idx == ')') { depth--; if (depth == 0) { p_idx++; break; } } p_idx++; }
4764	100	T	F	}
4765	50	T	F	// Match bare -1
4770	100	T	F	continue; // re-examine same position
4775	50	T	F	has_intercept = false;
	0	T	F	has_intercept = false;
4777	50	T	F	continue;
	50	T	F	continue;
4785	50	T	F	// Match bare 0 (entire rhs)
	0	T	F	// Match bare 0 (entire rhs)
4792	100	T	F	memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1);
4794	50	T	F	}
	100	T	F	}
4795	50	T	F	// Strip leading bare 1 or 1+
4802	100	T	F	}
4803	50	T	F
4808	100	T	F	memmove(p_idx, p_idx + 1, strlen(p_idx + 1) + 1);
4811	100	T	F	if (len_rhs > 0 && rhs[len_rhs - 1] == '+') rhs[len_rhs - 1] = '\0';
4815	100	T	F	char rhs_expanded[2048] = "";
4817	50	T	F	chunk = strtok(rhs, "+");
4819	100	T	F	if (strcmp(chunk, ".") == 0) {
	100	T	F	if (strcmp(chunk, ".") == 0) {
4820	100	T	F	AV *cols = get_all_columns(data_hoa, row_hashes, n);
4821	50	T	F	for (size_t c = 0; c <= (size_t)av_len(cols); c++) {
4827	50	T	F	if (rhs_len + slen + 2 < sizeof(rhs_expanded)) {
4828	100	T	F	if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
4829	100	T	F	strcat(rhs_expanded, col_name);
4830	100	T	F	rhs_len += slen;
4831	100	T	F	}
4832	50	T	F	}
4845	100	T	F	}
4856	50	T	F	if (num_terms >= term_cap - 3) {
	50	T	F	if (num_terms >= term_cap - 3) {
4857	50	T	F	term_cap *= 2;
4862	100	T	F	*star = '\0';
4864	100	T	F	char *restrict right = star + 1;
4868	100	T	F	if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
4869	100	T	F	terms[num_terms++] = savepv(left);
4871	100	T	F	size_t inter_len = strlen(left) + strlen(right) + 2;
4873	50	T	F	snprintf(terms[num_terms++], inter_len, "%s:%s", left, right);
4874	100	T	F	} else {
4879	50	T	F	chunk = strtok(NULL, "+");
4882	50	T	F
4885	100	T	F	for (j = 0; j < num_uniq; j++) { if (strcmp(terms[i], uniq_terms[j]) == 0) { found = true; break; } }
4892	100	T	F	// ========================================================================
4893	100	T	F	for (j = 0; j < p; j++) {
4894	100	T	F	if (p_exp + 32 >= exp_cap) {
4895	100	T	F	exp_cap *= 2;
4897	50	T	F	Renew(dummy_base, exp_cap, char); Renew(dummy_level, exp_cap, char);
4901	50	T	F	}
4909	100	T	F	bool found = false;
4910	100	T	F	for (l = 0; l < num_levels; l++) { if (strcmp(levels[l], str_val) == 0) { found = true; break; } }
4912	100	T	F	if (num_levels >= levels_cap) { levels_cap = 2; Renew(levels, levels_cap, char); }
4916	100	T	F	}
4918	100	T	F	if (num_levels > 0) {
4924	100	T	F	exp_cap *= 2;
4925	100	T	F	Renew(exp_terms, exp_cap, char*); Renew(is_dummy, exp_cap, bool);
4928	100	T	F	size_t t_len = strlen(uniq_terms[j]) + strlen(levels[l]) + 1;
	100	T	F	size_t t_len = strlen(uniq_terms[j]) + strlen(levels[l]) + 1;
4943	100	T	F	exp_terms[p_exp] = savepv(uniq_terms[j]); is_dummy[p_exp] = false; p_exp++;
4946	100	T	F	p = p_exp;
4948	100	T	F	Newx(valid_row_names, n, char*);
	100	T	F	Newx(valid_row_names, n, char*);
4951	100	T	F	// PHASE 4: Matrix Construction & Listwise Deletion
4960	50	T	F	if (strcmp(exp_terms[j], "Intercept") == 0) {
4966	50	T	F	Safefree(str_val);
	50	T	F	Safefree(str_val);
4969	50	T	F	row_x[j] = evaluate_term(data_hoa, row_hashes, i, exp_terms[j]);
	50	T	F	row_x[j] = evaluate_term(data_hoa, row_hashes, i, exp_terms[j]);
4972	0	T	F	}
4976	0	T	F	for (j = 0; j < p; j++) X[valid_n * p + j] = row_x[j];
4980	100	T	F	}
4984	100	T	F	for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
4991	50	T	F	Safefree(X); Safefree(Y); Safefree(valid_row_names);
	0	T	F	Safefree(X); Safefree(Y); Safefree(valid_row_names);
5011	50	T	F	}
5021	100	T	F	}
5022	100	T	F	}
5023	100	T	F
5025	100	T	F	// PHASE 6: Metrics & Cleanup
5030	100	T	F	df_res = (int)valid_n - final_rank;
5044	50	T	F	hv_store(fitted_hv, valid_row_names[i], strlen(valid_row_names[i]), newSVnv(y_hat), 0);
5045	0	T	F	hv_store(resid_hv, valid_row_names[i], strlen(valid_row_names[i]), newSVnv(res), 0);
5046	0	T	F	Safefree(valid_row_names[i]);
5054	100	T	F	double r_squared = 0.0, adj_r_squared = 0.0, f_stat = NAN, f_pvalue = NAN;
	50	T	F	double r_squared = 0.0, adj_r_squared = 0.0, f_stat = NAN, f_pvalue = NAN;
	100	T	F	double r_squared = 0.0, adj_r_squared = 0.0, f_stat = NAN, f_pvalue = NAN;
	50	T	F	double r_squared = 0.0, adj_r_squared = 0.0, f_stat = NAN, f_pvalue = NAN;
5062	50	T	F	f_pvalue = 1.0 - pf(f_stat, (double)numdf, (double)df_res);
5065	50	T	F	f_pvalue = 0.0;
5066	100	T	F	}
5076	100	T	F	hv_store(row_hv, "Estimate", 8, newSVpv("NaN", 0), 0);
5083	50	T	F	double p_val = get_t_pvalue(t_val, df_res, "two.sided");
	50	T	F	double p_val = get_t_pvalue(t_val, df_res, "two.sided");
	0	T	F	double p_val = get_t_pvalue(t_val, df_res, "two.sided");
	0	T	F	double p_val = get_t_pvalue(t_val, df_res, "two.sided");
5089	50	T	F	hv_store(summary_hv, exp_terms[j], strlen(exp_terms[j]), newRV_noinc((SV*)row_hv), 0);
5093	100	T	F	hv_store(res_hv, "fitted.values", 13, newRV_noinc((SV*)fitted_hv), 0);
5097	100	T	F	hv_store(res_hv, "rss", 3, newSVnv(rss), 0);
5098	100	T	F	hv_store(res_hv, "summary", 7, newRV_noinc((SV*)summary_hv),0);
5099	50	T	F	hv_store(res_hv, "terms", 5, newRV_noinc((SV*)terms_av), 0);
5103	100	T	F	AV *fstat_av = newAV();
5106	50	T	F	av_push(fstat_av, newSViv(df_res));
5109	100	T	F	}
5116	100	T	F	if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
	50	T	F	if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
5121	100	T	F	if (row_hashes) Safefree(row_hashes);
5156	50	T	F	EXTEND(SP, n_elements);
5162	50	T	F
5164	50	T	F	CODE:
5166	50	T	F	// Auto-seed the PRNG if the Perl script hasn't done so yet
5168	50	T	F
	50	T	F
5171	50	T	F	int arg_start = 0;
5172	100	T	F
5176	0	T	F	arg_start = 1; // Start parsing named arguments from the second element
	0	T	F	arg_start = 1; // Start parsing named arguments from the second element
5178	0	T	F
	0	T	F
5180	0	T	F	if ((items - arg_start) % 2 != 0) {
5188	0	T	F	if (strEQ(key, "n")) n = (unsigned int)SvUV(val);
5191	0	T	F	else croak("rnorm: unknown argument '%s'", key);
5192	0	T	F	}
5193	0	T	F
5195	0	T	F
	0	T	F
	0	T	F
5201	0	T	F	double u, v, s;
5212	100	T	F	if (i < n) {
	100	T	F	if (i < n) {
	50	T	F	if (i < n) {
5216	100	T	F	}
5217	100	T	F	RETVAL = newRV_noinc((SV*)result_av);
5218	50	T	F	}
5226	50	T	F	{
5227	50	T	F	const char *restrict formula = SvPV_nolen(formula_sv);
5228	50	T	F	char f_cpy[512];
5229	50	T	F	char restrict src, restrict dst, restrict tilde, restrict lhs, restrict rhs, restrict chunk;
	0	T	F	char restrict src, restrict dst, restrict tilde, restrict lhs, restrict rhs, restrict chunk;
5230	50	T	F
5231	50	T	F	char restrict terms = NULL, restrict uniq_terms = NULL, restrict exp_terms = NULL, restrict parent_term = NULL;
	0	T	F	char restrict terms = NULL, restrict uniq_terms = NULL, restrict exp_terms = NULL, restrict parent_term = NULL;
5232	50	T	F	bool restrict is_dummy = NULL, is_interact = NULL;
	0	T	F	bool restrict is_dummy = NULL, is_interact = NULL;
5234	50	T	F	int restrict term_map = NULL, restrict left_idx = NULL, *restrict right_idx = NULL;
5235	50	T	F	unsigned int term_cap = 64, exp_cap = 64, num_terms = 0, num_uniq = 0, p = 0, p_exp = 0;
5237	50	T	F	bool has_intercept = true;
	50	T	F	bool has_intercept = true;
5242	100	T	F	SV *restrict ref = NULL;
5243	100	T	F	HE *restrict entry;
5245	100	T	F	double *restrict Y = NULL;
5247	50	T	F	char *restrict term_base_level = NULL; / reference level for each uniq_term (NULL if not categorical) */
	50	T	F	char *restrict term_base_level = NULL; / reference level for each uniq_term (NULL if not categorical) */
5249	100	T	F	// ========================================================================
5251	50	T	F	// ========================================================================
5252	100	T	F	ref = SvRV(data_sv);
5262	50	T	F	Newx(row_names, n, char*);
5263	100	T	F	for(i = 0; i < n; i++) {
5279	50	T	F	} else if (SvTYPE(ref) == SVt_PVAV) {
5281	50	T	F	n = av_len(av) + 1;
5283	100	T	F	Newx(row_hashes, n, HV*);
5284	50	T	F	for (i = 0; i < n; i++) {
5289	100	T	F	snprintf(buf, sizeof(buf), "%lu", (unsigned long)(i + 1));
5294	50	T	F	croak("aov: Array values must be HashRefs (AoH)");
	0	T	F	croak("aov: Array values must be HashRefs (AoH)");
5295	50	T	F	}
	0	T	F	}
5303	50	T	F	while (src && (dst - f_cpy < 511)) { if (!isspace(src)) { dst++ = src; } src++; }
	0	T	F	while (src && (dst - f_cpy < 511)) { if (!isspace(src)) { dst++ = src; } src++; }
5310	100	T	F	croak("aov: invalid formula, missing '~'");
5312	100	T	F	*tilde = '\0';
5313	50	T	F	lhs = f_cpy;
5315	50	T	F
5322	100	T	F	if (rhs[0] == '1' && rhs[1] == '\0') { rhs[0] = '\0'; }
5323	100	T	F	else if (rhs[0] == '1' && rhs[1] == '+') { memmove(rhs, rhs + 2, strlen(rhs + 2) + 1); }
5331	100	T	F	size_t rhs_len = 0;
5341	100	T	F	size_t slen = strlen(col_name);
5349	100	T	F	}
5350	100	T	F	SvREFCNT_dec(cols);
5351	100	T	F	} else {
5354	100	T	F	if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
	50	T	F	if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
5358	100	T	F	}
5359	100	T	F	chunk = strtok(NULL, "+");
5360	50	T	F	}
5381	100	T	F	*star = '\0';
5385	100	T	F	if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
5387	50	T	F	terms[num_terms++] = savepv(left);
5389	100	T	F	size_t inter_len = strlen(left) + strlen(right) + 2;
5390	100	T	F	terms[num_terms] = (char*)safemalloc(inter_len);
5392	100	T	F	} else {
5393	50	T	F	char *restrict c_chunk = strchr(chunk, '^');
5400	50	T	F
5401	100	T	F	for (i = 0; i < num_terms; i++) {
5402	100	T	F	bool found = false;
5403	100	T	F	for (size_t k = 0; k < num_uniq; k++) {
5408	100	T	F	p = num_uniq;
5409	50	T	F
5426	100	T	F
5446	100	T	F	if (strcmp(parent_term[e], left) == 0) l_indices[l_count++] = e;
5447	50	T	F	if (strcmp(parent_term[e], right) == 0) r_indices[r_count++] = e;
5451	100	T	F	Safefree(l_indices); Safefree(r_indices);
5453	50	T	F	} else {
5456	100	T	F	if (p_exp >= exp_cap) {
5457	100	T	F	exp_cap *= 2;
5459	100	T	F	Renew(is_dummy, exp_cap, bool); Renew(is_interact, exp_cap, bool);
5461	100	T	F	Renew(term_map, exp_cap, int); Renew(left_idx, exp_cap, int); Renew(right_idx, exp_cap, int);
5463	50	T	F	size_t t_len = strlen(exp_terms[l_indices[li]]) + strlen(exp_terms[r_indices[ri]]) + 2;
5464	100	T	F	exp_terms[p_exp] = (char*)safemalloc(t_len);
5469	50	T	F	right_idx[p_exp] = r_indices[ri];
5472	50	T	F	}
5475	100	T	F	Safefree(l_indices); Safefree(r_indices);
5481	100	T	F	for (i = 0; i < n; i++) {
5483	100	T	F	if (str_val) {
5484	100	T	F	bool found = false;
5485	100	T	F	for (size_t l = 0; l < num_levels; l++) {
5487	50	T	F	}
5493	100	T	F	}
5495	50	T	F
5509	100	T	F	for (size_t l = 1; l < num_levels; l++) {
5510	100	T	F	if (p_exp >= exp_cap) {
5511	100	T	F	exp_cap *= 2;
5518	100	T	F	exp_terms[p_exp] = (char*)safemalloc(t_len);
5519	100	T	F	snprintf(exp_terms[p_exp], t_len, "%s%s", uniq_terms[j], levels[l]);
5522	100	T	F	dummy_base[p_exp] = savepv(uniq_terms[j]);
5526	50	T	F	}
5529	100	T	F	} else {
5530	100	T	F	Safefree(levels);
5534	100	T	F	term_map[p_exp] = j;
5539	50	T	F	parent_term[p_exp] = savepv(uniq_terms[j]);
	100	T	F	parent_term[p_exp] = savepv(uniq_terms[j]);
5557	100	T	F	for (j = 0; j < p_exp; j++) {
5558	100	T	F	if (strcmp(exp_terms[j], "Intercept") == 0) {
5559	100	T	F	row_x[j] = 1.0;
5561	100	T	F	row_x[j] = row_x[left_idx[j]] * row_x[right_idx[j]];
5568	100	T	F	} else {
5571	50	T	F	}
5582	100	T	F	if (valid_n <= p_exp) {
5589	100	T	F	}
5590	50	T	F	Safefree(exp_terms); Safefree(parent_term);
5593	50	T	F	Safefree(term_map); Safefree(left_idx); Safefree(right_idx);
	50	T	F	Safefree(term_map); Safefree(left_idx); Safefree(right_idx);
5595	50	T	F	Safefree(X_mat); Safefree(Y);
5600	50	T	F	/* ------------ */
5604	100	T	F	// PHASE 5: Math & Output Formatting
5606	50	T	F	bool restrict aliased_qr = (bool)safemalloc(p_exp * sizeof(bool));
5609	50	T	F	double *restrict term_ss;
	50	T	F	double *restrict term_ss;
	50	T	F	double *restrict term_ss;
	50	T	F	double *restrict term_ss;
5616	50	T	F	if (aliased_qr[i]) continue;
	50	T	F	if (aliased_qr[i]) continue;
5617	50	T	F	int t_idx = term_map[i];
	50	T	F	int t_idx = term_map[i];
5618	50	T	F	size_t r_k = rank_map[i];
	50	T	F	size_t r_k = rank_map[i];
5619	50	T	F	term_ss[t_idx] += Y[r_k] * Y[r_k];
	50	T	F	term_ss[t_idx] += Y[r_k] * Y[r_k];
5623	50	T	F	for (i = 0; i < p_exp; i++) {
5626	50	T	F	double rss_prev = 0.0;
5629	50	T	F	}
	50	T	F	}
5632	100	T	F
5633	100	T	F	HV*restrict ret_hash = newHV();
5636	50	T	F	HV*restrict term_stats = newHV();
	50	T	F	HV*restrict term_stats = newHV();
5637	50	T	F	double ss = term_ss[j];
	50	T	F	double ss = term_ss[j];
5642	50	T	F	hv_stores(term_stats, "Sum Sq", newSVnv(ss));
	50	T	F	hv_stores(term_stats, "Sum Sq", newSVnv(ss));
5647	100	T	F	hv_stores(term_stats, "Pr(>F)", newSVnv(1.0 - pf(f_val, (double)df, (double)res_df)));
5648	100	T	F	} else {
5653	100	T	F	}
5654	100	T	F
5655	50	T	F	HV*restrict res_stats = newHV();
	50	T	F	HV*restrict res_stats = newHV();
5656	50	T	F	hv_stores(res_stats, "Df", newSViv(res_df));
	50	T	F	hv_stores(res_stats, "Df", newSViv(res_df));
5657	50	T	F	hv_stores(res_stats, "Sum Sq", newSVnv(rss_prev));
	50	T	F	hv_stores(res_stats, "Sum Sq", newSVnv(rss_prev));
5658	50	T	F	hv_stores(res_stats, "Mean Sq", newSVnv(ms_res));
	50	T	F	hv_stores(res_stats, "Mean Sq", newSVnv(ms_res));
5700	50	T	F	for (i = 0; i < n; i++) Safefree(X_mat[i]);
5701	100	T	F	Safefree(X_mat); Safefree(Y);
5705	100	T	F	}
5706	100	T	F	OUTPUT:
5707	100	T	F	RETVAL
5708	50	T	F
	100	T	F
5709	100	T	F	PROTOTYPES: DISABLE
5710	100	T	F
5711	50	T	F	SV* fisher_test(...)
5712	0	T	F	CODE:
5713	0	T	F	{
5717	100	T	F	double conf_level = 0.95;
	50	T	F	double conf_level = 0.95;
5718	50	T	F	const char*restrict alternative = "two.sided";
	50	T	F	const char*restrict alternative = "two.sided";
5719	100	T	F
	50	T	F
5720	50	T	F	// Parse named arguments
	50	T	F	// Parse named arguments
5721	100	T	F	for (unsigned short int i = 1; i < items; i += 2) {
	50	T	F	for (unsigned short int i = 1; i < items; i += 2) {
5724	100	T	F	SV*restrict val = ST(i + 1);
5725	50	T	F	if (strEQ(key, "conf_level") \|\| strEQ(key, "conf.level")) {
5726	100	T	F	conf_level = SvNV(val);
5727	50	T	F	} else if (strEQ(key, "alternative")) {
5728	50	T	F	alternative = SvPV_nolen(val);
5730	50	T	F	}
5734	100	T	F	size_t a = 0, b = 0, c = 0, d = 0;
5735	50	T	F	// Extract Data
5736	50	T	F	if (SvTYPE(deref) == SVt_PVAV) {
	50	T	F	if (SvTYPE(deref) == SVt_PVAV) {
5737	100	T	F	AVrestrict outer = (AV)deref;
	50	T	F	AVrestrict outer = (AV)deref;
5738	100	T	F	if (av_len(outer) != 1) croak("Outer array must have exactly 2 rows");
5739	100	T	F	SV**restrict row1_ptr = av_fetch(outer, 0, 0);
	100	T	F	SV**restrict row1_ptr = av_fetch(outer, 0, 0);
5740	100	T	F	SV**restrict row2_ptr = av_fetch(outer, 1, 0);
	50	T	F	SV**restrict row2_ptr = av_fetch(outer, 1, 0);
	50	T	F	SV**restrict row2_ptr = av_fetch(outer, 1, 0);
5741	100	T	F	if (row1_ptr && row2_ptr && SvROK(row1_ptr) && SvROK(row2_ptr)) {
	50	T	F	if (row1_ptr && row2_ptr && SvROK(row1_ptr) && SvROK(row2_ptr)) {
5742	100	T	F	AVrestrict row1 = (AV)SvRV(*row1_ptr);
5744	50	T	F	SV**restrict a_ptr = av_fetch(row1, 0, 0);
5746	50	T	F	SV**restrict c_ptr = av_fetch(row2, 0, 0);
	0	T	F	SV**restrict c_ptr = av_fetch(row2, 0, 0);
5747	100	T	F	SV**restrict d_ptr = av_fetch(row2, 1, 0);
5749	100	T	F	b = (b_ptr && SvOK(b_ptr)) ? SvIV(b_ptr) : 0;
5753	0	T	F	croak("Invalid 2D Array structure");
5755	0	T	F	} else if (SvTYPE(deref) == SVt_PVHV) {
5757	0	T	F	HVrestrict outer = (HV)deref;
5761	0	T	F	if (!he1 \|\| !he2) croak("Invalid outer hash");
5763	0	T	F	const char*restrict k2 = SvPV_nolen(hv_iterkeysv(he2));
	0	T	F	const char*restrict k2 = SvPV_nolen(hv_iterkeysv(he2));
5764	0	T	F	HE*restrict row1_he = (strcmp(k1, k2) < 0) ? he1 : he2;
5766	0	T	F	SV*restrict row1_sv = hv_iterval(outer, row1_he);
5770	0	T	F	croak("Inner elements must be hashes");
5772	0	T	F	HVrestrict in1 = (HV)SvRV(row1_sv);
5774	0	T	F	if (hv_iterinit(in1) != 2 \|\| hv_iterinit(in2) != 2) croak("Inner hashes must have exactly 2 keys");
5786	100	T	F	HE*restrict in2_c2 = (strcmp(in2_k1, in2_k2) < 0) ? in2_he2 : in2_he1;
	100	T	F	HE*restrict in2_c2 = (strcmp(in2_k1, in2_k2) < 0) ? in2_he2 : in2_he1;
5788	100	T	F	b = (hv_iterval(in1, in1_c2) && SvOK(hv_iterval(in1, in1_c2))) ? SvIV(hv_iterval(in1, in1_c2)) : 0;
	100	T	F	b = (hv_iterval(in1, in1_c2) && SvOK(hv_iterval(in1, in1_c2))) ? SvIV(hv_iterval(in1, in1_c2)) : 0;
5789	100	T	F	c = (hv_iterval(in2, in2_c1) && SvOK(hv_iterval(in2, in2_c1))) ? SvIV(hv_iterval(in2, in2_c1)) : 0;
5803	50	T	F	hv_stores(ret_hash, "alternative", newSVpv(alternative, 0));
	100	T	F	hv_stores(ret_hash, "alternative", newSVpv(alternative, 0));
5805	100	T	F	av_push(ci_array, newSVnv(ci_low));
5807	50	T	F	hv_stores(ret_hash, "conf_int", newRV_noinc((SV*)ci_array));
5811	100	T	F	hv_stores(ret_hash, "p_value", newSVnv(p_val));
	50	T	F	hv_stores(ret_hash, "p_value", newSVnv(p_val));
5812	100	T	F	// Return the HashRef
5813	50	T	F	RETVAL = newRV_noinc((SV*)ret_hash);
5817	100	T	F
5820	100	T	F	{
5821	50	T	F	SV* sv_n = NULL;
5822	0	T	F	SV* sv_delta = NULL;
5826	100	T	F
	50	T	F
	50	T	F
5838	100	T	F	else if (strEQ(key, "delta")) sv_delta = val;
5839	50	T	F	else if (strEQ(key, "sd")) sv_sd = val;
	50	T	F	else if (strEQ(key, "sd")) sv_sd = val;
5846	100	T	F	else croak("power_t_test: unknown argument '%s'", key);
5848	50	T	F
	50	T	F
5852	50	T	F	bool is_null_sd = (sv_sd && !SvOK(sv_sd));
5859	100	T	F	if (is_null_sd) missing_count++;
5862	100	T	F	if (missing_count != 1) {
5864	50	T	F	}
	50	T	F	}
	50	T	F	}
5878	50	T	F	while (p_body(high, delta, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
	50	T	F	while (p_body(high, delta, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
	50	T	F	while (p_body(high, delta, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
5880	50	T	F	double mid = low + (high - low) / 2.0;
	50	T	F	double mid = low + (high - low) / 2.0;
	50	T	F	double mid = low + (high - low) / 2.0;
5887	50	T	F	while (high - low > tol) {
5888	50	T	F	double mid = low + (high - low) / 2.0;
5896	100	T	F	while (high - low > tol) {
5899	50	T	F	else high = mid;
	50	T	F	else high = mid;
	50	T	F	else high = mid;
5900	50	T	F	}
	50	T	F	}
5905	100	T	F	double mid = low + (high - low) / 2.0;
5923	50	T	F	}
5924	50	T	F	OUTPUT:
5933	100	T	F	// 1. Shift positional arguments
5941	100	T	F	}
5942	50	T	F	}
5950	50	T	F	const char *restrict key = SvPV_nolen(ST(arg_idx));