Branch Coverage

File:LikeR.xs
Coverage:66.9%

line%coveragebranch
3050TF        if (df <= 0.0) return 0.0;
40100TF        for (unsigned short i = 1; i < n_steps; i++) {
52100TF                double weight = (i % 2 != 0) ? 4.0 : 2.0;
73100TF        for (size_t i = 0; i < n; i++) {
79100TF        for (size_t i = 0; i < n; ) {
81100TF                while (j < n && items[j].val == items[i].val) j++;
50TF                while (j < n && items[j].val == items[i].val) j++;
83100TF                for (size_t k = i; k < j; k++) items[k].rank = avg_rank;
87100TF        for (size_t i = 0; i < n; i++) ranks[i] = items[i].rank;
93100TF        if (prob <= 0.0) return 0;
94100TF        if (prob >= 1.0) return size;
97100TF        for (size_t i = 0; i < size; i++) {
98100TF                if (Drand01() <= prob) successes++;
111100TF        for(size_t k = 0; k <= max_x - min_x; ++k) {
113100TF          if (d_val > max_d) max_d = d_val;
117100TF        for(size_t k = 0; k <= max_x - min_x; ++k) {
124100TF        for(size_t k = 0; k <= max_x - min_x; ++k) {
126100TF          if (min_x + k <= a) *lower_tail += p_prob;
127100TF          if (min_x + k >= a) *upper_tail += p_prob;
135100TF    size_t min_x = (r2 > c1) ? 0 : c1 - r2;
143100TF    for(size_t x = min_x; x <= max_x; ++x) {
14850TF    if (a == min_x && a == max_x) *mle_or = 1.0;
0TF    if (a == min_x && a == max_x) *mle_or = 1.0;
14950TF    else if (a == min_x) *mle_or = 0.0;
150100TF    else if (a == max_x) *mle_or = INFINITY;
15350TF        for (unsigned short int i = 0; i < 3000; i++) {
156100TF            for(size_t k = 0; k <= max_x - min_x; ++k) {
158100TF                if (d_val > max_d) max_d = d_val;
161100TF            for(size_t k = 0; k <= max_x - min_x; ++k) {
168100TF            if (exp_val > a) log_high = log_mid;
170100TF            if (log_high - log_low < 1e-15) break;
179100TF    if (!is_less) {
180100TF        double target_alpha = is_greater ? alpha : alpha / 2.0;
18150TF        if (a != min_x) {
18350TF            for (unsigned short int i = 0; i < 1000; i++) {
188100TF                if (err < best_err) { best_err = err; best = mid; }
189100TF                if (ut > target_alpha) log_high = log_mid;
191100TF                if (log_high - log_low < 1e-15) break;
198100TF    if (!is_greater) {
199100TF        double target_alpha = is_less ? alpha : alpha / 2.0;
200100TF        if (a != max_x) {
20250TF            for (unsigned short int i = 0; i < 1000; i++) {
207100TF                if (err < best_err) { best_err = err; best = mid; }
208100TF                if (lt > target_alpha) log_low = log_mid;
210100TF                if (log_high - log_low < 1e-15) break;
221100TF    size_t min_x = (r2 > c1) ? 0 : c1 - r2;
226100TF    for(size_t x = min_x; x <= max_x; ++x) {
232100TF    if (strcmp(alt, "less") == 0) {
233100TF        for(size_t x = min_x; x <= a; ++x) p_val += exp(logdc[x - min_x]);
234100TF    } else if (strcmp(alt, "greater") == 0) {
235100TF        for(size_t x = a; x <= max_x; ++x) p_val += exp(logdc[x - min_x]);
239100TF        for(size_t x = min_x; x <= max_x; ++x) {
241100TF            if (p_cur <= p_obs * relErr) p_val += p_cur;
24650TF    return (p_val > 1.0) ? 1.0 : p_val;
261100TF        for (size_t k = 0; k < n; k++) {
266100TF        for (size_t k = 0; k < n; k++) {
269100TF                if (fabs(A[k * n + k]) <= 1e-10 * orig_diag[k] || fabs(A[k * n + k]) < 1e-24) {
50TF                if (fabs(A[k * n + k]) <= 1e-10 * orig_diag[k] || fabs(A[k * n + k]) < 1e-24) {
272100TF                        for (size_t i = 0; i < n; i++) {
281100TF                for (size_t j = 0; j < n; j++) A[k * n + j] *= pivot;
282100TF                for (size_t i = 0; i < n; i++) {
283100TF                        if (i != k && A[i * n + k] != 0.0) {
100TF                        if (i != k && A[i * n + k] != 0.0) {
286100TF                                  for (size_t j = 0; j < n; j++) {
299100TF    if (row_hashes) {
30150TF        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
50TF        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
50TF        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
30550TF    } else if (data_hoa) {
30750TF        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
31250TF    if (val && SvOK(*val)) {
100TF    if (val && SvOK(*val)) {
31350TF        if (looks_like_number(*val)) return SvNV(*val);
32250TF    if (data_hoa) {
325100TF        while ((entry = hv_iternext(data_hoa))) {
3280TF    } else if (row_hashes && n > 0 && row_hashes[0]) {
0TF    } else if (row_hashes && n > 0 && row_hashes[0]) {
0TF    } else if (row_hashes && n > 0 && row_hashes[0]) {
3310TF        while ((entry = hv_iternext(row_hashes[0]))) {
34050TF    if (!term || term[0] == '\0') return NAN;
50TF    if (!term || term[0] == '\0') return NAN;
344100TF    if (colon) {
35050TF        if (isnan(left) || isnan(right)) return NAN;
50TF        if (isnan(left) || isnan(right)) return NAN;
35350TF    if (strncmp(term_cpy, "I(", 2) == 0) {
3550TF        if (end) *end = '\0';
3590TF        if (caret) {
3660TF        if (isnan(v)) return NAN;
3670TF        return power == 1 ? v : pow(v, power);
376100TF        for (size_t i = 0; i < n; i++) {
378100TF                if (row_hashes) {
380100TF                        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
50TF                        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
50TF                        if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
38450TF                } else if (data_hoa) {
38650TF                        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF                        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF                        if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
391100TF                if (val && SvOK(*val)) {
50TF                if (val && SvOK(*val)) {
392100TF                        if (looks_like_number(*val)) return false; // First valid is number -> Numeric Column
40250TF        if (row_hashes) {
4040TF                if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
0TF                if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
0TF                if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) {
40850TF        } else if (data_hoa) {
41050TF                if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF                if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
50TF                if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV) {
41550TF        if (val && SvOK(*val)) {
50TF        if (val && SvOK(*val)) {
430100TF        if (diff < 0) return -1;
43150TF        if (diff > 0) return 1;
447100TF        if (diff < 0) return -1;
448100TF        if (diff > 0) return  1;
45650TF        Newx(ri, n, RankItem);
457100TF        for (size_t i = 0; i < n; i++) { ri[i].val = in[i]; ri[i].idx = i; }
461100TF        while (i < n) {
464100TF                while (j + 1 < n && ri[j + 1].val == ri[j].val) j++;
100TF                while (j + 1 < n && ri[j + 1].val == ri[j].val) j++;
467100TF                for (size_t k = i; k <= j; k++) out[ri[k].idx] = avg;
477100TF        for (size_t i = 0; i < n; i++) {
48350TF        if (den == 0.0) return NAN;
497100TF        for (size_t i = 0; i < n - 1; i++) {
498100TF          for (size_t j = i + 1; j < n; j++) {
501100TF                   if      (sx == 0 && sy == 0) { /* joint tie — not counted */ }
50TF                   if      (sx == 0 && sy == 0) { /* joint tie — not counted */ }
502100TF                   else if (sx == 0)            tie_x++;
50350TF                   else if (sy == 0)            tie_y++;
50450TF                   else if (sx == sy)           C++;
50950TF        if (denom == 0.0) return NAN;
517100TF        if (strcmp(method, "spearman") == 0) {
51950TF          Newx(rx, n, double); Newx(ry, n, double);
50TF          Newx(rx, n, double); Newx(ry, n, double);
526100TF        if (strcmp(method, "kendall") == 0)
54250TF        if (fabs(d) < FPMIN) d = FPMIN;
54450TF        for (m = 1; m <= MAX_ITER; m++) {
54850TF          if (fabs(d) < FPMIN) d = FPMIN;
55050TF          if (fabs(c) < FPMIN) c = FPMIN;
55450TF          if (fabs(d) < FPMIN) d = FPMIN;
55650TF          if (fabs(c) < FPMIN) c = FPMIN;
558100TF          if (fabs(del - 1.0) < EPS) break;
56450TF        if (x <= 0.0) return 0.0;
565100TF        if (x >= 1.0) return 1.0;
567100TF        if (x < (a + 1.0) / (a + b + 2.0)) return bt * _incbeta_cf(a, b, x) / a;
574100TF        if (strcmp(alt, "less") == 0) return (t < 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
100TF        if (strcmp(alt, "less") == 0) return (t < 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
575100TF        if (strcmp(alt, "greater") == 0) return (t > 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
50TF        if (strcmp(alt, "greater") == 0) return (t > 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail;
583100TF        while (get_t_pvalue(high, df, "greater") > p_tail) {
58650TF          if (high > 1000000.0) break; /* Fallback limit */
58950TF        for (unsigned short int i = 0; i < 100; i++) {
592100TF          if (p_mid > p_tail) {
597100TF          if (high - low < 1e-8) break;
6090TF        if (n == 0) return 1;
61850TF        double step = (n_bins > 0) ? (breaks[1] - breaks[0]) : 0.0;
620100TF        for (size_t i = 0; i < n_bins; i++) {
625100TF        if (step > 0.0) {
626100TF          for (size_t j = 0; j < n; j++) {
62950TF                   if (isnan(val) || isinf(val) || val < min_val) continue;
50TF                   if (isnan(val) || isinf(val) || val < min_val) continue;
50TF                   if (isnan(val) || isinf(val) || val < min_val) continue;
633100TF                   if (idx >= n_bins) {
639100TF                   while (idx > 0 && val <= breaks[idx]) {
100TF                   while (idx > 0 && val <= breaks[idx]) {
643100TF                   while (idx < n_bins - 1 && val > breaks[idx + 1]) {
50TF                   while (idx < n_bins - 1 && val > breaks[idx + 1]) {
64850TF        } else if (n_bins > 0) {
653100TF        for (size_t i = 0; i < n_bins; i++) {
655100TF          if (bin_width > 0) {
65850TF                   density[i] = (n_bins == 1) ? 1.0 : 0.0;
681100TF        if (fabs(y) < 0.42) {
687100TF          if (y > 0) r = 1.0 - p;
691100TF          if (y < 0) x = -x;
709100TF        for (size_t i = 0; i < n; i++) { perm[i] = i + 1; c[i] = 0; }
724100TF        TALLY_PERM();   /* initial permutation [1, 2, ..., n] */
50TF        TALLY_PERM();   /* initial permutation [1, 2, ..., n] */
50TF        TALLY_PERM();   /* initial permutation [1, 2, ..., n] */
727100TF        while (k < n) {
728100TF          if (c[k] < k) {
730100TF                   if (k % 2 == 0) {
735100TF                   TALLY_PERM();
100TF                   TALLY_PERM();
100TF                   TALLY_PERM();
75150TF        if (strcmp(alt, "greater") == 0) return p_le;
75250TF        if (strcmp(alt, "less")    == 0) return p_ge;
75450TF        double p = 2.0 * (p_le < p_ge ? p_le : p_ge);
75550TF        return (p > 1.0) ? 1.0 : p;
764100TF        for (long i = 0; i <= max_inv; i++) dp[i] = 0.0;
767100TF        for (size_t i = 2; i <= n; i++) {
769100TF          for (long k = 0; k <= max_inv; k++) next_dp[k] = 0.0;
771100TF          for (int k = 0; k <= current_max_inv; k++) {
773100TF                   for (int j = 0; j <= i - 1 && k - j >= 0; j++) {
100TF                   for (int j = 0; j <= i - 1 && k - j >= 0; j++) {
78450TF        if (i_obs < 0) i_obs = 0;
78550TF        if (i_obs > max_inv) i_obs = max_inv;
787100TF        for (long k = i_obs; k <= max_inv; k++) p_le += dp[k];
789100TF        for (long k = 0; k <= i_obs; k++) p_ge += dp[k];
79150TF        if (strcmp(alt, "greater") == 0) return p_ge;
792100TF        if (strcmp(alt, "less") == 0) return p_le;
79450TF        double p = 2.0 * (p_ge < p_le ? p_ge : p_le);
79550TF        return p > 1.0 ? 1.0 : p;
79950TF        if (f <= 0.0) return 0.0;
808100TF        for (size_t k = 0; k < p; k++) {
81050TF                if (r >= n) {
816100TF                for (size_t i = r; i < n; i++) {
817100TF                        if (fabs(X[i][k]) > max_val) max_val = fabs(X[i][k]);
819100TF                if (max_val < 1e-10) {
825100TF                for (size_t i = r; i < n; i++) {
830100TF                double s = (X[r][k] > 0) ? -norm : norm;
834100TF                for (size_t j = k + 1; j < p; j++) {
836100TF                        for (size_t i = r + 1; i < n; i++) dot += X[i][j] * X[i][k];
839100TF                        for (size_t i = r + 1; i < n; i++) X[i][j] += tau * X[i][k];
844100TF                for (size_t i = r + 1; i < n; i++) dot_y += y[i] * X[i][k];
847100TF                for (size_t i = r + 1; i < n; i++) y[i] += tau_y * X[i][k];
86350TF        if (!sv || !SvOK(sv)) return 0;
50TF        if (!sv || !SvOK(sv)) return 0;
866100TF        for (size_t i = 0; i < len; i++) {
86750TF          if (!isdigit(s[i])) return 1;
87450TF        if (!str) str = "";
876100TF        if (strstr(str, sep) != NULL || strchr(str, '"') != NULL || strchr(str, '\r') != NULL || strchr(str, '\n') != NULL) {
100TF        if (strstr(str, sep) != NULL || strchr(str, '"') != NULL || strchr(str, '\r') != NULL || strchr(str, '\n') != NULL) {
50TF        if (strstr(str, sep) != NULL || strchr(str, '"') != NULL || strchr(str, '\r') != NULL || strchr(str, '\n') != NULL) {
50TF        if (strstr(str, sep) != NULL || strchr(str, '"') != NULL || strchr(str, '\r') != NULL || strchr(str, '\n') != NULL) {
880100TF        if (needs_quotes) {
882100TF          for (const char *restrict p = str; *p; p++) {
883100TF                   if (*p == '"') {
899100TF        for (size_t i = 0; i < len; i++) {
900100TF          if (i > 0) PerlIO_write(fh, sep, sep_len);
90150TF          if (row[i]) {
91250TF        if (x < 0.0 || a <= 0.0) return 1.0;
50TF        if (x < 0.0 || a <= 0.0) return 1.0;
91350TF        if (x == 0.0) return 1.0;
916100TF        if (x < a + 1.0) {
920100TF                while (fabs(term) > 1e-15) {
93350TF        while (i < 10000) { // Safety bound
93750TF                if (fabs(d) < 1e-30) d = 1e-30;
93950TF                if (fabs(c) < 1e-30) c = 1e-30;
943100TF                if (fabs(del - 1.0) < 1e-15) break;
95150TF        if (df <= 0) return 1.0;
95250TF        if (stat <= 0.0) return 1.0;
96350TF        if (k < 0 || k > n) return 0.0L;
50TF        if (k < 0 || k > n) return 0.0L;
96450TF        if (k > n / 2) k = n - k;
966100TF        for (int i = 1; i <= k; i++) {
977100TF    if (k < 0) return 0.0;
97850TF    if (k >= max_u) return 1.0;
983100TF    for (int j = 1; j <= n; j++) {
984100TF        for (int i = j; i <= max_u; i++) w[i] += w[i - j];
985100TF        for (int i = max_u; i >= j + m; i--) w[i] -= w[i - j - m];
989100TF    for (int i = 0; i <= k; i++) cum_p += w[i];
100350TF        if (k < 0) return 0.0;
1004100TF        if (k >= max_v) return 1.0;
1009100TF        for (int i = 1; i <= n; i++) {
1010100TF          for (int j = max_v; j >= i; j--) w[j] += w[j - i];
1014100TF        for (int i = 0; i <= k; i++) cum_p += w[i];
103050TF        if (n == 0) return 0.0;
1035100TF        while (i < n) {
1037100TF                while (j < n && ri[j].val == ri[i].val) j++;
100TF                while (j < n && ri[j].val == ri[i].val) j++;
1039100TF                for (size_t k = i; k < j; k++) ri[k].rank = r;
1041100TF                if (t > 1) { *has_ties = 1; tie_adj += ((double)t * t * t - t); }
105950TF        if (n == 0) return 1.0;
106050TF        if (n < 0) return 1.0 / r_pow_di(x, -n);
1062100TF        for (int i = 0; i < n; i++) val *= x;
10700TF        if(x <= 0.) {
10710TF          if(lower) p = 0.;
10730TF        } else if(x < 1.) {
10780TF          for(k = 1; k < k_max; k += 2) {
10820TF          if(!lower) p = 1.0 - p;
10870TF          if(lower) {
10920TF          while(fabs(old_val - new_val) > tol) {
1105100TF        for(unsigned int i = 0; i < m; i++) {
1106100TF          for(unsigned int j = 0; j < m; j++) {
1108100TF                   for(unsigned int k = 0; k < m; k++) s += A[i * m + k] * B[k * m + j];
1115100TF    if(n == 1) {
1116100TF        for(int i = 0; i < m * m; i++) V[i] = A[i];
1124100TF    if((n % 2) == 0) {
1125100TF        for(int i = 0; i < m * m; i++) V[i] = B[i];
113150TF    if(V[(m / 2) * m + (m / 2)] > 1e140) {
11320TF        for(int i = 0; i < m * m; i++) V[i] = V[i] * 1e-140;
1146100TF        for(int i = 0; i < m; i++) {
1147100TF          for(int j = 0; j < m; j++) {
1148100TF                   if(i - j + 1 < 0) H[i * m + j] = 0;
1152100TF        for(int i = 0; i < m; i++) {
115650TF        H[(m - 1) * m] += ((2 * h - 1 > 0) ? r_pow_di(2 * h - 1, m) : 0);
1158100TF        for(int i = 0; i < m; i++) {
1159100TF          for(int j = 0; j < m; j++) {
1160100TF                   if(i - j + 1 > 0) {
1161100TF                       for(int g = 1; g <= i - j + 1; g++) H[i * m + j] /= g;
1170100TF        for(int i = 1; i <= n; i++) {
117250TF          if(s < 1e-140) {
1191100TF        while(i < nx || j < ny) {
50TF        while(i < nx || j < ny) {
119350TF          if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
100TF          if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
100TF          if (i < nx && j < ny) val = (x[i] < y[j]) ? x[i] : y[j];
119450TF          else if (i < nx) val = x[i];
1197100TF          while(i < nx && x[i] <= val) i++;
100TF          while(i < nx && x[i] <= val) i++;
1198100TF          while(j < ny && y[j] <= val) j++;
100TF          while(j < ny && y[j] <= val) j++;
1204100TF          if (diff > max_d_plus) max_d_plus = diff;
1205100TF          if (-diff > max_d_minus) max_d_minus = -diff;
1206100TF          if (fabs(diff) > max_d) max_d = fabs(diff);
1215100TF        if (two_sided) return (fabs(r - s) >= q);
1225100TF        for(unsigned int j = 1; j <= n; j++) {
1226100TF          if(psmirnov_exact_test(q, 0., j / nd, two_sided)) u[j] = 1.;
1229100TF        for(unsigned int i = 1; i <= m; i++) {
1230100TF          if(psmirnov_exact_test(q, i / md, 0., two_sided)) u[0] = 1.;
1231100TF          for(int j = 1; j <= n; j++) {
1232100TF                   if(psmirnov_exact_test(q, i / md, j / nd, two_sided)) u[j] = 1.;
124750TF        if (nu < 1e-7) nu = 1e-7;
125550TF        if (strict && tside == 2) {
0TF        if (strict && tside == 2) {
127650TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
128150TF        if (arg_idx < items) {
1282100TF                if (SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF                if (SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
128550TF                } else if (SvPOK(ST(arg_idx))) {
1292100TF        for (; arg_idx < items; arg_idx += 2) {
129550TF          if      (strEQ(key, "x"))           x_sv = val;
129650TF          else if (strEQ(key, "y"))           y_sv = val;
129750TF          else if (strEQ(key, "exact"))       {
12980TF                   if (!SvOK(val)) exact = -1;
130150TF          else if (strEQ(key, "alternative")) alternative = SvPV_nolen(val);
130550TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) {
50TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) {
50TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) {
1313100TF        if (!is_two_sided && !is_greater && !is_less) {
100TF        if (!is_two_sided && !is_greater && !is_less) {
50TF        if (!is_two_sided && !is_greater && !is_less) {
131950TF        if (nx == 0) croak("Not enough 'x' observations");
1324100TF        for (size_t i = 0; i < nx; i++) {
132650TF          if (el && SvOK(*el) && looks_like_number(*el)) {
50TF          if (el && SvOK(*el) && looks_like_number(*el)) {
50TF          if (el && SvOK(*el) && looks_like_number(*el)) {
133550TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
100TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
50TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
1341100TF          for (size_t i = 0; i < ny; i++) {
134350TF                   if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                   if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                   if (el && SvOK(*el) && looks_like_number(*el)) {
134850TF          if (valid_nx < 1 || valid_ny < 1) {
50TF          if (valid_nx < 1 || valid_ny < 1) {
1357100TF          if (is_greater) statistic = d_plus;
1358100TF          else if (is_less) statistic = d_minus;
136350TF          if (exact == 1) use_exact = true;
136450TF          else if (exact == 0) use_exact = false;
1370100TF          for(size_t i=0; i<valid_nx; i++) comb[i] = x_data[i];
1371100TF          for(size_t i=0; i<valid_ny; i++) comb[valid_nx+i] = y_data[i];
1375100TF          for(size_t i = 1; i < total_n; i++) {
137650TF                   if(comb[i] == comb[i-1]) { has_ties = true; break; }
137950TF          if (use_exact && has_ties) {
50TF          if (use_exact && has_ties) {
138350TF          if (use_exact) {
13900TF                   if (is_two_sided) {
139950TF        else if (y_sv && SvPOK(y_sv)) {
50TF        else if (y_sv && SvPOK(y_sv)) {
140150TF          if (strEQ(dist, "pnorm")) {
1404100TF                   for(size_t i = 0; i < valid_nx; i++) {
141250TF                       if (diff1 > max_d_plus) max_d_plus = diff1;
1413100TF                       if (diff2 > max_d_plus) max_d_plus = diff2;
1414100TF                       if (-diff1 > max_d_minus) max_d_minus = -diff1;
141550TF                       if (-diff2 > max_d_minus) max_d_minus = -diff2;
1417100TF                       if (fabs(diff1) > max_d) max_d = fabs(diff1);
141850TF                       if (fabs(diff2) > max_d) max_d = fabs(diff2);
142050TF                   if (is_greater) statistic = max_d_plus;
142150TF                   else if (is_less) statistic = max_d_minus;
142350TF                   bool use_exact = (exact == -1) ? (valid_nx < 100) : (exact == 1);
142450TF                   if (use_exact) {
142650TF                       if (is_two_sided) {
14360TF                       if (is_two_sided) p_value = K2l(z, 0, 1e-6);
144850TF        if (p_value > 1.0) p_value = 1.0;
144950TF        if (p_value < 0.0) p_value = 0.0;
147050TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
100TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
147550TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
100TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF        if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
148050TF        if ((items - arg_idx) % 2 != 0) {
1484100TF        for (; arg_idx < items; arg_idx += 2) {
1487100TF                if      (strEQ(key, "x"))           x_sv = val;
1488100TF                else if (strEQ(key, "y"))           y_sv = val;
1489100TF                else if (strEQ(key, "paired"))      paired = SvTRUE(val);
149050TF                else if (strEQ(key, "correct"))     correct = SvTRUE(val);
1491100TF                else if (strEQ(key, "mu"))          mu = SvNV(val);
149250TF                else if (strEQ(key, "exact"))       {
14930TF                        if (!SvOK(val)) exact = -1;
149650TF                else if (strEQ(key, "alternative")) alternative = SvPV_nolen(val);
1500100TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF        if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
150450TF        if (nx == 0) croak("Not enough 'x' observations");
1508100TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
50TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
50TF        if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) {
1516100TF        if (ny > 0 && !paired) {
100TF        if (ny > 0 && !paired) {
1519100TF                for (size_t i = 0; i < nx; i++) {
152150TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
1527100TF                for (size_t i = 0; i < ny; i++) {
152950TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                        if (el && SvOK(*el) && looks_like_number(*el)) {
153550TF                if (valid_nx == 0) { Safefree(ri); croak("not enough (non-missing) 'x' observations"); }
153650TF                if (valid_ny == 0) { Safefree(ri); croak("not enough 'y' observations"); }
1541100TF                for (size_t i = 0; i < total_n; i++) if (ri[i].idx == 1) w_rank_sum += ri[i].rank;
100TF                for (size_t i = 0; i < total_n; i++) if (ri[i].idx == 1) w_rank_sum += ri[i].rank;
154450TF                if (exact == 1) use_exact = true;
154550TF                else if (exact == 0) use_exact = false;
154650TF                else use_exact = (valid_nx < 50 && valid_ny < 50 && !has_ties);
50TF                else use_exact = (valid_nx < 50 && valid_ny < 50 && !has_ties);
100TF                else use_exact = (valid_nx < 50 && valid_ny < 50 && !has_ties);
1548100TF                if (use_exact && has_ties) {
50TF                if (use_exact && has_ties) {
1552100TF                if (use_exact) {
1557100TF                        if (strcmp(alternative, "less") == 0) p_value = p_less;
155850TF                        else if (strcmp(alternative, "greater") == 0) p_value = p_greater;
15600TF                                double p = (p_less < p_greater) ? p_less : p_greater;
156450TF                        method_desc = correct ? "Wilcoxon rank sum test with continuity correction" : "Wilcoxon rank sum test";
157050TF                        if (correct) {
157150TF                                if (strcmp(alternative, "two.sided") == 0) CORRECTION = (z > 0 ? 0.5 : -0.5);
100TF                                if (strcmp(alternative, "two.sided") == 0) CORRECTION = (z > 0 ? 0.5 : -0.5);
15720TF                                else if (strcmp(alternative, "greater") == 0) CORRECTION = 0.5;
15730TF                                else if (strcmp(alternative, "less") == 0) CORRECTION = -0.5;
157750TF                        if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z);
157850TF                        else if (strcmp(alternative, "greater") == 0) p_value = 1.0 - approx_pnorm(z);
1583100TF                if (paired && (!y_av || nx != ny)) croak("'x' and 'y' must have the same length for paired test");
50TF                if (paired && (!y_av || nx != ny)) croak("'x' and 'y' must have the same length for paired test");
100TF                if (paired && (!y_av || nx != ny)) croak("'x' and 'y' must have the same length for paired test");
1587100TF                for (size_t i = 0; i < nx; i++) {
158950TF                        if (!x_el || !SvOK(*x_el) || !looks_like_number(*x_el)) continue;
50TF                        if (!x_el || !SvOK(*x_el) || !looks_like_number(*x_el)) continue;
50TF                        if (!x_el || !SvOK(*x_el) || !looks_like_number(*x_el)) continue;
1592100TF                        if (paired) {
159450TF                                if (!y_el || !SvOK(*y_el) || !looks_like_number(*y_el)) continue;
50TF                                if (!y_el || !SvOK(*y_el) || !looks_like_number(*y_el)) continue;
50TF                                if (!y_el || !SvOK(*y_el) || !looks_like_number(*y_el)) continue;
159750TF                                if (d == 0.0) has_zeroes = true; // Drop exact zeroes
160150TF                                if (d == 0.0) has_zeroes = true;
160550TF                if (n_nz == 0) {
1610100TF                for (size_t i = 0; i < n_nz; i++) {
1617100TF                for (size_t i = 0; i < n_nz; i++) {
1618100TF                        if (ri[i].idx) statistic += ri[i].rank;
162050TF                if (exact == 1) use_exact = true;
162150TF                else if (exact == 0) use_exact = false;
162250TF                else use_exact = (n_nz < 50 && !has_ties);
50TF                else use_exact = (n_nz < 50 && !has_ties);
162350TF                if (use_exact && has_ties) {
50TF                if (use_exact && has_ties) {
162750TF                if (use_exact && has_zeroes) {
50TF                if (use_exact && has_zeroes) {
163150TF                if (use_exact) {
163650TF                        if (strcmp(alternative, "less") == 0) p_value = p_less;
163750TF                        else if (strcmp(alternative, "greater") == 0) p_value = p_greater;
163950TF                                double p = (p_less < p_greater) ? p_less : p_greater;
16430TF                        method_desc = correct ? "Wilcoxon signed rank test with continuity correction" : "Wilcoxon signed rank test";
16480TF                        if (correct) {
16490TF                                if (strcmp(alternative, "two.sided") == 0) CORRECTION = (z > 0 ? 0.5 : -0.5);
0TF                                if (strcmp(alternative, "two.sided") == 0) CORRECTION = (z > 0 ? 0.5 : -0.5);
16500TF                                else if (strcmp(alternative, "greater") == 0) CORRECTION = 0.5;
16510TF                                else if (strcmp(alternative, "less") == 0) CORRECTION = -0.5;
16550TF                        if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z);
16560TF                        else if (strcmp(alternative, "greater") == 0) p_value = 1.0 - approx_pnorm(z);
166150TF        if (p_value > 1.0) p_value = 1.0;
167750TF        int r = av_top_index(obs_av) + 1, c = 0;
168050TF        if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
100TF        if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
50TF        if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
168350TF                c = av_top_index(first_row) + 1;
1691100TF        int yates = (is_2d && r == 2 && c == 2) ? 1 : 0;
50TF        int yates = (is_2d && r == 2 && c == 2) ? 1 : 0;
100TF        int yates = (is_2d && r == 2 && c == 2) ? 1 : 0;
1694100TF        if (is_2d) {
1697100TF                for(unsigned int i=0; i<r; i++) row_sum[i] = 0.0;
1698100TF                for(unsigned int j=0; j<c; j++) col_sum[j] = 0.0;
1699100TF                for (unsigned int i = 0; i < r; i++) {
1702100TF                        for (unsigned int j = 0; j < c; j++) {
1710100TF                for (unsigned int i = 0; i < r; i++) {
1714100TF                        for (unsigned int j = 0; j < c; j++) {
1719100TF                                if (yates) {
172250TF                                  double y_corr = (abs_diff > 0.5) ? 0.5 : abs_diff;
1734100TF          for (unsigned int j = 0; j < c; j++) {
1739100TF          for (unsigned int j = 0; j < c; j++) {
1753100TF        if (is_2d) {
1754100TF                if (yates) {
177750TF        if (arg_idx < items && SvROK(ST(arg_idx))) {
50TF        if (arg_idx < items && SvROK(ST(arg_idx))) {
1779100TF          if (type == SVt_PVHV || type == SVt_PVAV) {
50TF          if (type == SVt_PVHV || type == SVt_PVAV) {
178450TF        if (arg_idx < items) {
1794100TF        for (; arg_idx < items; arg_idx += 2) {
179550TF          if (arg_idx + 1 >= items) croak("write_table: Odd number of arguments passed");
179850TF          if (strEQ(key, "data")) data_sv = val;
1799100TF          else if (strEQ(key, "col.names")) col_names_sv = val;
180050TF          else if (strEQ(key, "file")) file_sv = val;
1801100TF          else if (strEQ(key, "row.names")) row_names_sv = val;
180250TF          else if (strEQ(key, "sep")) sep = SvPV_nolen(val);
180650TF        if (!data_sv || !SvROK(data_sv)) {
50TF        if (!data_sv || !SvROK(data_sv)) {
1810100TF        if (SvTYPE(data_ref) != SVt_PVHV && SvTYPE(data_ref) != SVt_PVAV) {
50TF        if (SvTYPE(data_ref) != SVt_PVHV && SvTYPE(data_ref) != SVt_PVAV) {
181450TF        if (!file_sv || !SvOK(file_sv)) croak("write_table: file name missing\n");
50TF        if (!file_sv || !SvOK(file_sv)) croak("write_table: file name missing\n");
1817100TF        if (col_names_sv && SvOK(col_names_sv)) {
50TF        if (col_names_sv && SvOK(col_names_sv)) {
1818100TF          if (!SvROK(col_names_sv) || SvTYPE(SvRV(col_names_sv)) != SVt_PVAV) {
50TF          if (!SvROK(col_names_sv) || SvTYPE(SvRV(col_names_sv)) != SVt_PVAV) {
1827100TF        if (SvTYPE(data_ref) == SVt_PVHV) {
182950TF                if (hv_iterinit(hv) == 0) XSRETURN_EMPTY;
183350TF                if (!first_val || !SvROK(first_val)) {
50TF                if (!first_val || !SvROK(first_val)) {
1837100TF                if (first_type != SVt_PVHV && first_type != SVt_PVAV) {
50TF                if (first_type != SVt_PVHV && first_type != SVt_PVAV) {
1845100TF                while ((entry = hv_iternext(hv))) {
184750TF                        if (!val || !SvROK(val) || SvTYPE(SvRV(val)) != first_type) {
50TF                        if (!val || !SvROK(val) || SvTYPE(SvRV(val)) != first_type) {
50TF                        if (!val || !SvROK(val) || SvTYPE(SvRV(val)) != first_type) {
18480TF                                 croak("write_table: Mixed data types detected. Ensure all values are %s references.\n", is_hoh ? "HASH" : "ARRAY");
1852100TF                if (is_hoh) {
1855100TF                        while ((entry = hv_iternext(hv))) {
186150TF                if (av_len(av) < 0) XSRETURN_EMPTY;
186350TF                if (!first_ptr || !*first_ptr || !SvROK(*first_ptr) || SvTYPE(SvRV(*first_ptr)) != SVt_PVHV) {
50TF                if (!first_ptr || !*first_ptr || !SvROK(*first_ptr) || SvTYPE(SvRV(*first_ptr)) != SVt_PVHV) {
50TF                if (!first_ptr || !*first_ptr || !SvROK(*first_ptr) || SvTYPE(SvRV(*first_ptr)) != SVt_PVHV) {
50TF                if (!first_ptr || !*first_ptr || !SvROK(*first_ptr) || SvTYPE(SvRV(*first_ptr)) != SVt_PVHV) {
1867100TF                for (size_t i = 0; i <= av_len(av); i++) {
186950TF                        if (!ptr || !*ptr || !SvROK(*ptr) || SvTYPE(SvRV(*ptr)) != SVt_PVHV) {
50TF                        if (!ptr || !*ptr || !SvROK(*ptr) || SvTYPE(SvRV(*ptr)) != SVt_PVHV) {
50TF                        if (!ptr || !*ptr || !SvROK(*ptr) || SvTYPE(SvRV(*ptr)) != SVt_PVHV) {
50TF                        if (!ptr || !*ptr || !SvROK(*ptr) || SvTYPE(SvRV(*ptr)) != SVt_PVHV) {
187750TF        if (!fh) croak("write_table: Could not open '%s' for writing", file);
188050TF        int inc_rownames = (row_names_sv && SvTRUE(row_names_sv)) ? 1 : 0;
100TF        int inc_rownames = (row_names_sv && SvTRUE(row_names_sv)) ? 1 : 0;
1884100TF        if (is_hoh) {
1885100TF                if (col_names_sv && SvOK(col_names_sv)) {
50TF                if (col_names_sv && SvOK(col_names_sv)) {
1887100TF                        for(size_t i=0; i<=av_len(c_av); i++) {
188950TF                                if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
50TF                                if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
1895100TF                        while((entry = hv_iternext((HV*)data_ref))) {
1899100TF                                 while((inner_entry = hv_iternext(inner))) {
1905100TF                        for(unsigned i=0; i<num_cols; i++) {
1910100TF                        for(unsigned i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0));
191850TF          if (inc_rownames) header_row[h_idx++] = "";
1919100TF          for(unsigned short int i=0; i<num_headers; i++) {
192150TF                   header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                   header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
1928100TF          for(size_t i=0; i<num_rows; i++) {
1936100TF          for(size_t i=0; i<num_rows; i++) {
193850TF                   if (inc_rownames) row_data[d_idx++] = row_array[i];
194150TF                   HV *restrict inner_hv = inner_hv_ptr ? (HV*)SvRV(*inner_hv_ptr) : NULL;
1943100TF                   for(size_t j=0; j<num_headers; j++) {
194550TF                       const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                       const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
194650TF                       SV **restrict cell_ptr = inner_hv ? hv_fetch(inner_hv, col_name, strlen(col_name), 0) : NULL;
1947100TF                       if (cell_ptr && SvOK(*cell_ptr)) {
100TF                       if (cell_ptr && SvOK(*cell_ptr)) {
1948100TF                           if (SvROK(*cell_ptr)) {
195250TF                       if (headers_av) SvREFCNT_dec(headers_av);
195350TF                       if (rows_av) SvREFCNT_dec(rows_av);
1965100TF        } else if (is_hoa) { // ----- Hash of Arrays -----
1970100TF                while((entry = hv_iternext(data_hv))) {
1973100TF                        if (len > max_rows) max_rows = len;
1976100TF                if (col_names_sv && SvOK(col_names_sv)) {
50TF                if (col_names_sv && SvOK(col_names_sv)) {
1978100TF                        for(size_t i=0; i<=av_len(c_av); i++) {
198050TF                                 if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
50TF                                 if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
1985100TF                        for(unsigned int i=0; i<num_cols; i++) {
1990100TF                        for(unsigned i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0));
199350TF                if (av_len(headers_av) < 0) croak("Could not get headers in write_table");
1994100TF                if (inc_rownames && contains_nondigit(row_names_sv)) {
50TF                if (inc_rownames && contains_nondigit(row_names_sv)) {
19980TF                        for(size_t i=0; i<=av_len(headers_av); i++) {
20000TF                                 if (!h_ptr || !*h_ptr) continue;
0TF                                 if (!h_ptr || !*h_ptr) continue;
20020TF                                 if (strcmp(SvPV_nolen(h_sv), rownames_col) != 0) {
2012100TF                if (inc_rownames) header_row[h_idx++] = "";
2013100TF                for(size_t i=0; i<num_headers; i++) {
201550TF                        header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                        header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
2020100TF                for(size_t i=0; i<max_rows; i++) {
2022100TF                        if (inc_rownames) {
202350TF                                 if (rownames_col) {
20250TF                                     if (rn_arr_ptr && SvROK(*rn_arr_ptr)) {
0TF                                     if (rn_arr_ptr && SvROK(*rn_arr_ptr)) {
20280TF                                         if (rn_val_ptr && SvOK(*rn_val_ptr)) {
0TF                                         if (rn_val_ptr && SvOK(*rn_val_ptr)) {
20290TF                                             if (SvROK(*rn_val_ptr)) {
20320TF                                                 if (headers_av) SvREFCNT_dec(headers_av);
2048100TF                        for(size_t j=0; j<num_headers; j++) {
205050TF                                 const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                                 const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
205250TF                                 if (arr_ptr && SvROK(*arr_ptr)) {
50TF                                 if (arr_ptr && SvROK(*arr_ptr)) {
2055100TF                                     if (cell_ptr && SvOK(*cell_ptr)) {
100TF                                     if (cell_ptr && SvOK(*cell_ptr)) {
205650TF                                         if (SvROK(*cell_ptr)) {
20590TF                                             if (headers_av) SvREFCNT_dec(headers_av);
2071100TF                        if (inc_rownames && !rownames_col) safefree((char*)row_data[0]);
50TF                        if (inc_rownames && !rownames_col) safefree((char*)row_data[0]);
207450TF        } else if (is_aoh) {// ----- Array of Hashes -----
2077100TF                if (col_names_sv && SvOK(col_names_sv)) {
50TF                if (col_names_sv && SvOK(col_names_sv)) {
2079100TF                        for(size_t i=0; i<=av_len(c_av); i++) {
208150TF                                 if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
50TF                                 if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c));
2085100TF                        for(size_t i=0; i<num_rows; i++) {
208750TF                                 if (row_ptr && SvROK(*row_ptr)) {
50TF                                 if (row_ptr && SvROK(*row_ptr)) {
2091100TF                                     while((entry = hv_iternext(row_hv))) {
2098100TF                        for(unsigned int i=0; i<num_cols; i++) {
2103100TF                        for(unsigned int i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0));
2107100TF                if (inc_rownames && contains_nondigit(row_names_sv)) {
50TF                if (inc_rownames && contains_nondigit(row_names_sv)) {
21100TF                        for(size_t i=0; i<=av_len(headers_av); i++) {
21120TF                                 if (!h_ptr || !*h_ptr) continue;
0TF                                 if (!h_ptr || !*h_ptr) continue;
21140TF                                 if (strcmp(SvPV_nolen(h_sv), rownames_col) != 0) {
2124100TF                if (inc_rownames) header_row[h_idx++] = "";
2125100TF                for(size_t i=0; i<num_headers; i++) {
212750TF                        header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                        header_row[h_idx++] = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
2132100TF                for(size_t i=0; i<num_rows; i++) {
213550TF                        HV *restrict row_hv = (row_ptr && SvROK(*row_ptr)) ? (HV*)SvRV(*row_ptr) : NULL;
50TF                        HV *restrict row_hv = (row_ptr && SvROK(*row_ptr)) ? (HV*)SvRV(*row_ptr) : NULL;
2136100TF                        if (inc_rownames) {
213750TF                                if (rownames_col) {
21380TF                                  SV **restrict rn_val_ptr = row_hv ? hv_fetch(row_hv, rownames_col, strlen(rownames_col), 0) : NULL;
21390TF                                  if (rn_val_ptr && SvOK(*rn_val_ptr)) {
0TF                                  if (rn_val_ptr && SvOK(*rn_val_ptr)) {
21400TF                                                if (SvROK(*rn_val_ptr)) {
21430TF                                                                  if (headers_av) SvREFCNT_dec(headers_av);
2157100TF                        for(size_t j=0; j<num_headers; j++) {
215950TF                                 const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
50TF                                 const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : "";
216050TF                                 SV **restrict cell_ptr = row_hv ? hv_fetch(row_hv, col_name, strlen(col_name), 0) : NULL;
2161100TF                                 if (cell_ptr && SvOK(*cell_ptr)) {
50TF                                 if (cell_ptr && SvOK(*cell_ptr)) {
216250TF                                     if (SvROK(*cell_ptr)) {
21650TF                         if (headers_av) SvREFCNT_dec(headers_av);
2174100TF                        if (inc_rownames && !rownames_col) safefree((char*)row_data[0]);
50TF                        if (inc_rownames && !rownames_col) safefree((char*)row_data[0]);
217850TF        if (headers_av) SvREFCNT_dec(headers_av);
2179100TF        if (rows_av) SvREFCNT_dec(rows_av);
219650TF        if (SvOK(callback) && SvROK(callback) && SvTYPE(SvRV(callback)) == SVt_PVCV) {
50TF        if (SvOK(callback) && SvROK(callback) && SvTYPE(SvRV(callback)) == SVt_PVCV) {
50TF        if (SvOK(callback) && SvROK(callback) && SvTYPE(SvRV(callback)) == SVt_PVCV) {
220150TF        sep_len = sep_str ? strlen(sep_str) : 0;
220250TF        comment_len = comment_str ? strlen(comment_str) : 0;
220550TF        if (!fp) {
2210100TF        while (sv_gets(line_sv, fp, 0) != NULL) {
221450TF                if (len > 0 && line[len-1] == '\n') {
100TF                if (len > 0 && line[len-1] == '\n') {
221650TF                        if (len > 0 && line[len-1] == '\r') {
100TF                        if (len > 0 && line[len-1] == '\r') {
2220100TF                if (!in_quotes) {
222350TF                        for (size_t i = 0; i < len; i++) {
222450TF                                if (line[i] != ' ' && line[i] != '\t') { is_empty = 0; break; }
100TF                                if (line[i] != ' ' && line[i] != '\t') { is_empty = 0; break; }
222650TF                        if (is_empty) continue;
222950TF                        if (comment_len > 0 && len >= comment_len && strncmp(line, comment_str, comment_len) == 0) {
50TF                        if (comment_len > 0 && len >= comment_len && strncmp(line, comment_str, comment_len) == 0) {
50TF                        if (comment_len > 0 && len >= comment_len && strncmp(line, comment_str, comment_len) == 0) {
2234100TF                for (size_t i = 0; i < len; i++) {
2236100TF                        if (ch == '\r') continue;
2237100TF                        if (ch == '"') {
2238100TF                                if (in_quotes && (i + 1 < len) && line[i+1] == '"') {
100TF                                if (in_quotes && (i + 1 < len) && line[i+1] == '"') {
100TF                                if (in_quotes && (i + 1 < len) && line[i+1] == '"') {
2241100TF                                } else if (in_quotes) {
2244100TF                                } else if (!post_quote) {
2247100TF                        } else if (!in_quotes && sep_len > 0 && (len - i) >= sep_len && strncmp(line + i, sep_str, sep_len) == 0) {
50TF                        } else if (!in_quotes && sep_len > 0 && (len - i) >= sep_len && strncmp(line + i, sep_str, sep_len) == 0) {
50TF                        } else if (!in_quotes && sep_len > 0 && (len - i) >= sep_len && strncmp(line + i, sep_str, sep_len) == 0) {
100TF                        } else if (!in_quotes && sep_len > 0 && (len - i) >= sep_len && strncmp(line + i, sep_str, sep_len) == 0) {
2256100TF                if (in_quotes) {
226550TF                        if (use_cb) {
226950TF                                PUSHMARK(SP);
227050TF                                XPUSHs(sv_2mortal(newRV_inc((SV*)current_row)));
227350TF                                FREETMPS;
2285100TF        if (in_quotes) {
228750TF                if (use_cb) {
229150TF                        PUSHMARK(SP);
229250TF                        XPUSHs(sv_2mortal(newRV_inc((SV*)current_row)));
229550TF                        FREETMPS;
230650TF        if (use_cb) {
231850TF                if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) {
50TF                if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) {
232150TF                if (!SvROK(y_sv) || SvTYPE(SvRV(y_sv)) != SVt_PVAV) {
50TF                if (!SvROK(y_sv) || SvTYPE(SvRV(y_sv)) != SVt_PVAV) {
2326100TF                if (strcmp(method, "pearson") != 0 &&
2327100TF                        strcmp(method, "spearman") != 0 &&
232850TF                        strcmp(method, "kendall") != 0) {
233750TF                if (nx != ny) {
2348100TF                for (size_t i = 0; i < nx; i++) {
235350TF                        double xv = (x_tv && SvOK(*x_tv) && looks_like_number(*x_tv)) ? SvNV(*x_tv) : NAN;
50TF                        double xv = (x_tv && SvOK(*x_tv) && looks_like_number(*x_tv)) ? SvNV(*x_tv) : NAN;
50TF                        double xv = (x_tv && SvOK(*x_tv) && looks_like_number(*x_tv)) ? SvNV(*x_tv) : NAN;
235450TF                        double yv = (y_tv && SvOK(*y_tv) && looks_like_number(*y_tv)) ? SvNV(*y_tv) : NAN;
50TF                        double yv = (y_tv && SvOK(*y_tv) && looks_like_number(*y_tv)) ? SvNV(*y_tv) : NAN;
50TF                        double yv = (y_tv && SvOK(*y_tv) && looks_like_number(*y_tv)) ? SvNV(*y_tv) : NAN;
235750TF                        if (!isnan(xv) && !isnan(yv)) {
50TF                        if (!isnan(xv) && !isnan(yv)) {
236550TF                if (n < 2) {
2373100TF                        if (strcmp(method, "kendall") == 0) {
2375100TF                                 for (size_t i = 0; i < n; i++) {
2376100TF                                     for (size_t j = 0; j < n; j++) {
2385100TF                                 if (strcmp(method, "spearman") == 0) {
2394100TF                                     for (size_t i = 0; i < n; i++) {
2406100TF                                     for (size_t i = 0; i < n; i++) {
246250TF        if (items % 2 != 0) croak("Usage: glm(formula => 'am ~ wt + hp', data => \\%mtcars)");
2464100TF        for (unsigned short i_arg = 0; i_arg < items; i_arg += 2) {
2467100TF          if      (strEQ(key, "formula")) formula = SvPV_nolen(val);
2468100TF          else if (strEQ(key, "data"))    data_sv = val;
246950TF          else if (strEQ(key, "family"))  family_str = SvPV_nolen(val);
247250TF        if (!formula) croak("glm: formula is required");
247350TF        if (!data_sv || !SvROK(data_sv)) croak("glm: data is required and must be a reference");
50TF        if (!data_sv || !SvROK(data_sv)) croak("glm: data is required and must be a reference");
2477100TF        if (!is_binomial && !is_gaussian) croak("glm: unsupported family '%s'", family_str);
50TF        if (!is_binomial && !is_gaussian) croak("glm: unsupported family '%s'", family_str);
2485100TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
100TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
50TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
248950TF        if (!tilde) croak("glm: invalid formula, missing '~'");
249350TF        if (strstr(rhs, "-1")) has_intercept = false;
249450TF        if (has_intercept) terms[num_terms++] = savepv("Intercept");
2497100TF        while (chunk != NULL) {
249850TF          if (num_terms >= term_cap - 3) {
250250TF          if (strcmp(chunk, "1") == 0 || strcmp(chunk, "-1") == 0) {
50TF          if (strcmp(chunk, "1") == 0 || strcmp(chunk, "-1") == 0) {
250750TF          if (star) {
25100TF                   char *restrict c_l = strchr(left, '^'); if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
0TF                   char *restrict c_l = strchr(left, '^'); if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
25110TF                   char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
0TF                   char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
252050TF                   if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
0TF                   if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
2526100TF        for (i = 0; i < num_terms; i++) {
2528100TF          for (size_t j = 0; j < num_uniq; j++) {
252950TF                   if (strcmp(terms[i], uniq_terms[j]) == 0) { found = true; break; }
253150TF          if (!found) uniq_terms[num_uniq++] = savepv(terms[i]);
253750TF        if (SvTYPE(ref) == SVt_PVHV) {
253950TF                if (hv_iterinit(hv) == 0) croak("glm: Data hash is empty");
254150TF                if (entry) {
254350TF                        if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
100TF                        if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
254650TF                                 Newx(row_names, n, char*);
2547100TF                                 for(i = 0; i < n; i++) {
255150TF                        } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
50TF                        } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
255350TF                                 Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
50TF                                 Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
2555100TF                                 while ((entry = hv_iternext(hv))) {
25630TF        } else if (SvTYPE(ref) == SVt_PVAV) {
25660TF          Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
0TF          Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
25670TF          for (i = 0; i < n; i++) {
25690TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
25740TF                       for (size_t k = 0; k < i; k++) Safefree(row_names[k]);
2582100TF        for (size_t j = 0; j < p; j++) {
258350TF          if (p_exp + 32 >= exp_cap) {
2588100TF          if (strcmp(uniq_terms[j], "Intercept") == 0) {
2591100TF          if (is_column_categorical(data_hoa, row_hashes, n, uniq_terms[j])) {
259350TF                   Newx(levels, levels_cap, char*);
2594100TF                   for (i = 0; i < n; i++) {
259650TF                       if (str_val) {
2598100TF                           for (size_t l = 0; l < num_levels; l++) {
2599100TF                               if (strcmp(levels[l], str_val) == 0) { found = true; break; }
2601100TF                           if (!found) {
260250TF                               if (num_levels >= levels_cap) { levels_cap *= 2; Renew(levels, levels_cap, char*); }
0TF                               if (num_levels >= levels_cap) { levels_cap *= 2; Renew(levels, levels_cap, char*); }
260850TF                   if (num_levels > 0) {
2609100TF                       for (size_t l1 = 0; l1 < num_levels - 1; l1++) {
2610100TF                           for (size_t l2 = l1 + 1; l2 < num_levels; l2++) {
261150TF                               if (strcmp(levels[l1], levels[l2]) > 0) {
2616100TF                       for (size_t l = 1; l < num_levels; l++) {
261750TF                           if (p_exp >= exp_cap) {
2628100TF                       for (size_t l = 0; l < num_levels; l++) Safefree(levels[l]);
263950TF        Newx(X, n * p, double); Newx(Y, n, double);
50TF        Newx(X, n * p, double); Newx(Y, n, double);
264050TF        Newx(valid_row_names, n, char*);
2643100TF        for (size_t i = 0; i < n; i++) {
264550TF                if (isnan(y_val)) { Safefree(row_names[i]); continue; }
2649100TF                for (size_t j = 0; j < p; j++) {
2650100TF                        if (strcmp(exp_terms[j], "Intercept") == 0) {
2652100TF                        } else if (is_dummy[j]) {
265450TF                                 if (str_val) {
2655100TF                                     row_x[j] = (strcmp(str_val, dummy_level[j]) == 0) ? 1.0 : 0.0;
266050TF                                 if (isnan(row_x[j])) { row_ok = false; break; }
266350TF                if (!row_ok) { Safefree(row_names[i]); Safefree(row_x); continue; }
2665100TF                for (size_t j = 0; j < p; j++) X[valid_n * p + j] = row_x[j];
267150TF        if (valid_n <= p) {
26720TF          Safefree(X); Safefree(Y); Safefree(valid_row_names); if (row_hashes) Safefree(row_hashes);
2681100TF        for (i = 0; i < p; i++) { beta[i] = 0.0; beta_old[i] = 0.0; }
2684100TF        for (i = 0; i < valid_n; i++) sum_y += Y[i];
2686100TF        for (i = 0; i < valid_n; i++) {
2687100TF          if (is_binomial) {
268850TF                   if (Y[i] < 0.0 || Y[i] > 1.0) croak("glm: binomial family requires response between 0 and 1");
50TF                   if (Y[i] < 0.0 || Y[i] > 1.0) croak("glm: binomial family requires response between 0 and 1");
2692100TF                   if (Y[i] == 0.0)      dev = -2.0 * log(1.0 - mu[i]);
269350TF                   else if (Y[i] == 1.0) dev = -2.0 * log(mu[i]);
270250TF        for (iter = 1; iter <= max_iter; iter++) {
2703100TF                for (i = 0; i < valid_n; i++) {
2704100TF                        if (is_binomial) {
270750TF                                 if (varmu < 1e-10) varmu = 1e-10;
2716100TF                for (i = 0; i < p; i++) { XtWZ[i] = 0.0; for (size_t j = 0; j < p; j++) XtWX[i * p + j] = 0.0; }
100TF                for (i = 0; i < p; i++) { XtWZ[i] = 0.0; for (size_t j = 0; j < p; j++) XtWX[i * p + j] = 0.0; }
2717100TF                for (size_t k = 0; k < valid_n; k++) {
2719100TF                        for (i = 0; i < p; i++) {
2722100TF                                 for (size_t j = 0; j < p; j++) XtWX[i * p + j] += xw * X[k * p + j];
2726100TF                for (i = 0; i < p; i++) {
272750TF                        if (aliased[i]) { beta[i] = NAN; } else {
272950TF                                 for (size_t j = 0; j < p; j++) if (!aliased[j]) sum += XtWX[i * p + j] * XtWZ[j];
100TF                                 for (size_t j = 0; j < p; j++) if (!aliased[j]) sum += XtWX[i * p + j] * XtWZ[j];
2735100TF                for (unsigned short int half = 0; half < 10; half++) {
2737100TF                        for (i = 0; i < valid_n; i++) {
273950TF                                 for (size_t j = 0; j < p; j++) if (!aliased[j]) linear_pred += X[i * p + j] * beta[j];
100TF                                 for (size_t j = 0; j < p; j++) if (!aliased[j]) linear_pred += X[i * p + j] * beta[j];
2741100TF                                 if (is_binomial) {
274450TF                                     if (mu[i] < 10 * DBL_EPSILON) mu[i] = 10 * DBL_EPSILON;
274550TF                                     if (mu[i] > 1.0 - 10 * DBL_EPSILON) mu[i] = 1.0 - 10 * DBL_EPSILON;
2748100TF                                     if (Y[i] == 0.0)      dev = -2.0 * log(1.0 - mu[i]);
274950TF                                     else if (Y[i] == 1.0) dev = -2.0 * log(mu[i]);
2759100TF                        if (!is_binomial || deviance_new <= deviance_old + 1e-7 || !isfinite(deviance_new)) {
100TF                        if (!is_binomial || deviance_new <= deviance_old + 1e-7 || !isfinite(deviance_new)) {
50TF                        if (!is_binomial || deviance_new <= deviance_old + 1e-7 || !isfinite(deviance_new)) {
2764100TF                        for (size_t j = 0; j < p; j++) beta[j] = (beta[j] + beta_old[j]) / 2.0;
2767100TF                if (fabs(deviance_new - deviance_old) / (0.1 + fabs(deviance_new)) < epsilon) {
2771100TF                for (size_t j = 0; j < p; j++) beta_old[j] = beta[j];
2774100TF        for (i = 0; i < p; i++) { for (size_t j = 0; j < p; j++) XtWX[i * p + j] = 0.0; }
100TF        for (i = 0; i < p; i++) { for (size_t j = 0; j < p; j++) XtWX[i * p + j] = 0.0; }
2775100TF        for (size_t k = 0; k < valid_n; k++) {
2776100TF          double w = is_binomial ? (mu[k] * (1.0 - mu[k])) : 1.0;
277750TF          if (w < 1e-10) w = 1e-10;
2778100TF          for (i = 0; i < p; i++) {
2780100TF                   for (size_t j = 0; j < p; j++) XtWX[i * p + j] += xw * X[k * p + j];
2786100TF        for (i = 0; i < valid_n; i++) {
2787100TF          if (is_binomial) {
2788100TF                   if (Y[i] == 0.0)      null_dev += -2.0 * log(1.0 - wtdmu);
278950TF                   else if (Y[i] == 1.0) null_dev += -2.0 * log(wtdmu);
2797100TF        if (is_gaussian) {
280050TF        } else if (is_binomial) {
2806100TF        dispersion = is_binomial ? 1.0 : ((df_res > 0) ? (deviance_new / df_res) : NAN);
50TF        dispersion = is_binomial ? 1.0 : ((df_res > 0) ? (deviance_new / df_res) : NAN);
2807100TF        for (size_t i = 0; i < valid_n; i++) {
2809100TF                if (is_binomial) {
2812100TF                        if (Y[i] == 0.0)      d_res = sqrt(-2.0 * log(1.0 - mu[i]));
281350TF                        else if (Y[i] == 1.0) d_res = sqrt(-2.0 * log(mu[i]));
2815100TF                        res = (Y[i] > mu[i]) ? d_res : -d_res;
2824100TF        for (size_t j = 0; j < p; j++) {
282950TF                if (aliased[j]) {
28320TF                        hv_store(row_hv, is_binomial ? "z value" : "t value", 7, newSVpv("NaN", 0), 0);
28330TF                        hv_store(row_hv, is_binomial ? "Pr(>|z|)" : "Pr(>|t|)", 8, newSVpv("NaN", 0), 0);
2837100TF                        double p_val = is_binomial ? 2.0 * (1.0 - approx_pnorm(fabs(val_stat))) : get_t_pvalue(val_stat, df_res, "two.sided");
2841100TF                        hv_store(row_hv, is_binomial ? "z value" : "t value", 7, newSVnv(val_stat), 0);
2842100TF                        hv_store(row_hv, is_binomial ? "Pr(>|z|)" : "Pr(>|t|)", 8, newSVnv(p_val), 0);
2864100TF        for (i = 0; i < num_terms; i++) Safefree(terms[i]);
2866100TF        for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]);
2868100TF        for (size_t j = 0; j < p_exp; j++) {
2870100TF                if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
2877100TF        if (row_hashes) Safefree(row_hashes);
288750TF        if (items < 2 || items % 2 != 0)
50TF        if (items < 2 || items % 2 != 0)
2899100TF        for (unsigned short int i = 2; i < items; i += 2) {
2903100TF          if      (strEQ(key, "alternative")) alternative = SvPV_nolen(val);
2904100TF          else if (strEQ(key, "method"))      method = SvPV_nolen(val);
290550TF          else if (strEQ(key, "exact"))       exact_sv = val;
2906100TF          else if (strEQ(key, "conf.level") || strEQ(key, "conf_level")) conf_level = SvNV(val);
50TF          else if (strEQ(key, "conf.level") || strEQ(key, "conf_level")) conf_level = SvNV(val);
290750TF          else if (strEQ(key, "continuity"))  continuity = SvTRUE(val);
292050TF        if (!SvOK(x_ref) || !SvROK(x_ref) || SvTYPE(SvRV(x_ref)) != SVt_PVAV ||
50TF        if (!SvOK(x_ref) || !SvROK(x_ref) || SvTYPE(SvRV(x_ref)) != SVt_PVAV ||
50TF        if (!SvOK(x_ref) || !SvROK(x_ref) || SvTYPE(SvRV(x_ref)) != SVt_PVAV ||
292150TF            !SvOK(y_ref) || !SvROK(y_ref) || SvTYPE(SvRV(y_ref)) != SVt_PVAV) {
50TF            !SvOK(y_ref) || !SvROK(y_ref) || SvTYPE(SvRV(y_ref)) != SVt_PVAV) {
50TF            !SvOK(y_ref) || !SvROK(y_ref) || SvTYPE(SvRV(y_ref)) != SVt_PVAV) {
292950TF        if (n_raw != av_len(y_av) + 1) croak("incompatible dimensions");
2935100TF        for (size_t i = 0; i < n_raw; i++) {
293950TF          double xv = (x_val && SvOK(*x_val) && looks_like_number(*x_val)) ? SvNV(*x_val) : NAN;
100TF          double xv = (x_val && SvOK(*x_val) && looks_like_number(*x_val)) ? SvNV(*x_val) : NAN;
50TF          double xv = (x_val && SvOK(*x_val) && looks_like_number(*x_val)) ? SvNV(*x_val) : NAN;
294050TF          double yv = (y_val && SvOK(*y_val) && looks_like_number(*y_val)) ? SvNV(*y_val) : NAN;
100TF          double yv = (y_val && SvOK(*y_val) && looks_like_number(*y_val)) ? SvNV(*y_val) : NAN;
50TF          double yv = (y_val && SvOK(*y_val) && looks_like_number(*y_val)) ? SvNV(*y_val) : NAN;
2943100TF          if (!isnan(xv) && !isnan(yv)) {
100TF          if (!isnan(xv) && !isnan(yv)) {
295050TF        if (n < 3) {
2956100TF        if (is_pearson) {
2959100TF          for (size_t i = 0; i < n; i++) {
296850TF          estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
50TF          estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
2982100TF        } else if (is_kendall) {
2984100TF          for (size_t i = 0; i < n - 1; i++) {
2985100TF                   for (size_t j = i + 1; j < n; j++) {
298950TF                       if (sign_x == 0 && sign_y == 0) { /* Joint tie, ignore */ }
0TF                       if (sign_x == 0 && sign_y == 0) { /* Joint tie, ignore */ }
299050TF                       else if (sign_x == 0) tie_x++;
299150TF                       else if (sign_y == 0) tie_y++;
2992100TF                       else if (sign_x * sign_y > 0) c++;
299750TF          estimate = (denom == 0.0) ? (0.0/0.0) : (double)(c - d) / denom;
299950TF          bool has_ties = (tie_x > 0 || tie_y > 0);
50TF          bool has_ties = (tie_x > 0 || tie_y > 0);
300350TF          if (!exact_sv || !SvOK(exact_sv)) {
0TF          if (!exact_sv || !SvOK(exact_sv)) {
300450TF                   do_exact = (n < 50) && !has_ties;
50TF                   do_exact = (n < 50) && !has_ties;
300950TF          if (do_exact && has_ties) do_exact = 0;
50TF          if (do_exact && has_ties) do_exact = 0;
301150TF          if (do_exact) {
30190TF                   if (continuity) S -= (S > 0 ? 1 : -1);
0TF                   if (continuity) S -= (S > 0 ? 1 : -1);
30220TF                   if (strcmp(alternative, "two.sided") == 0) {
30240TF                   } else if (strcmp(alternative, "less") == 0) {
303050TF        } else if (is_spearman) {
3038100TF          for (size_t i = 0; i < n; i++) {
304750TF          estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
50TF          estimate = (M2_x > 0.0 && M2_y > 0.0) ? cov / sqrt(M2_x * M2_y) : 0.0;
3051100TF          for (size_t i = 0; i < n; i++) {
3058100TF          for (size_t i = 0; i < n; i++) {
305950TF                   if (rank_x[i] != floor(rank_x[i]) || rank_y[i] != floor(rank_y[i])) {
50TF                   if (rank_x[i] != floor(rank_x[i]) || rank_y[i] != floor(rank_y[i])) {
306450TF          if (!exact_sv || !SvOK(exact_sv)) {
0TF          if (!exact_sv || !SvOK(exact_sv)) {
306550TF                   do_exact = (n < 10) && !has_ties;
50TF                   do_exact = (n < 10) && !has_ties;
307050TF          if (do_exact) {
30750TF                   if (continuity)
3092100TF        if (is_pearson) {
311450TF        if (!SvROK(data) || SvTYPE(SvRV(data)) != SVt_PVAV) {
50TF        if (!SvROK(data) || SvTYPE(SvRV(data)) != SVt_PVAV) {
312150TF        Newx(x, n_raw, double);
3124100TF        for (size_t i = 0; i < n_raw; i++) {
312650TF          if (elem && SvOK(*elem)) {
50TF          if (elem && SvOK(*elem)) {
312850TF                   if (!isnan(val)) {
313650TF        if (n < 3 || n > 5000) {
50TF        if (n < 3 || n > 5000) {
3143100TF        for (size_t i = 0; i < n; i++) {
314650TF        if (ssq == 0.0) {
315350TF        if (n == 3) {
31570TF          if (w < 0.75) w = 0.75;
316350TF          Newx(m, n, double);
316450TF          Newx(a, n, double);
3165100TF          for (size_t i = 0; i < n; i++) {
317350TF          if (n == 4 || n == 5) {
100TF          if (n == 4 || n == 5) {
3175100TF                   for (unsigned int i = 1; i < n-1; i++) {
3183100TF                   for (unsigned int i = 2; i < n-2; i++) {
3187100TF          for (size_t i = 0; i < n; i++) {
3197100TF          if (n <= 11) {
320350TF                   if (y >= gamma) {
322650TF          if (p_val > 1.0) p_val = 1.0;
322750TF          if (p_val < 0.0) p_val = 0.0;
323750TF        EXTEND(SP, 1);
3247100TF                for (unsigned short int i = 0; i < items; i++) {
3249100TF                        if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                        if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3252100TF                                for (size_t j = 0; j < len; j++) {
325450TF                                  if (tv && SvOK(*tv)) {
50TF                                  if (tv && SvOK(*tv)) {
3256100TF                                                if (first || val < min_val) {
50TF                                                if (first || val < min_val) {
326350TF                        } else if (SvOK(arg)) {
3265100TF                                if (first || val < min_val) {
100TF                                if (first || val < min_val) {
3272100TF                if (count == 0) croak("min needs >= 1 numeric element");
327350TF                RETVAL = min_val;
3284100TF          for (size_t i = 0; i < items; i++) {
3286100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3289100TF                       for (size_t j = 0; j < len; j++) {
329150TF                           if (tv && SvOK(*tv)) {
50TF                           if (tv && SvOK(*tv)) {
3293100TF                               if (first || val > max_val) {
100TF                               if (first || val > max_val) {
330050TF                   } else if (SvOK(arg)) {
3302100TF                       if (first || val > max_val) {
100TF                       if (first || val > max_val) {
3309100TF          if (count == 0) croak("max needs >= 1 numeric element");
331050TF          RETVAL = max_val;
332550TF        if (items == 0) {
3329100TF        while (i < items) {
3331100TF                if (i + 1 < items && SvPOK(ST(i))) {
100TF                if (i + 1 < items && SvPOK(ST(i))) {
3333100TF                        if (strEQ(key, "n")) {
3338100TF                        } else if (strEQ(key, "min")) {
334350TF                        } else if (strEQ(key, "max")) {
3352100TF                if (!n_set) {
3355100TF                } else if (!min_set) {
335850TF                } else if (!max_set) {
336650TF        if (!n_set) {
337050TF        AUTO_SEED_PRNG();
337250TF        if (n > 0) {
3376100TF        for (size_t j = 0; j < n; j++) {
337850TF                if (max < min) {
339450TF          AUTO_SEED_PRNG();
3395100TF          if (items % 2 != 0)
3403100TF          for (unsigned short i = 0; i < items; i += 2) {
3407100TF                   if      (strEQ(key, "n"))      n    = (unsigned int)SvUV(val);
3408100TF                   else if (strEQ(key, "size")) { size = (unsigned int)SvUV(val); size_set = true; }
340950TF                   else if (strEQ(key, "prob")) { prob = SvNV(val); prob_set = true; }
3414100TF          if (!size_set || !prob_set) croak("rbinom: 'size' and 'prob' are required arguments");
100TF          if (!size_set || !prob_set) croak("rbinom: 'size' and 'prob' are required arguments");
3415100TF          if (prob < 0.0 || prob > 1.0) croak("rbinom: prob must be between 0 and 1");
100TF          if (prob < 0.0 || prob > 1.0) croak("rbinom: prob must be between 0 and 1");
341850TF          if (n > 0) {
3420100TF                   for (unsigned int i = 0; i < n; i++) {
3435100TF                if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
100TF                if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
3440100TF                if (n_raw == 0) croak("hist: input array is empty");
344450TF                Newx(x, n_raw, double);
3448100TF                for (size_t i = 0; i < n_raw; i++) {
345050TF                        if (tv && SvOK(*tv)) {
50TF                        if (tv && SvOK(*tv)) {
3453100TF                                 if (val < min_val) min_val = val;
3454100TF                                 if (val > max_val) max_val = val;
345750TF                if (n == 0) {
346450TF                if (items == 2) {
346750TF                } else if (items > 2) {
346950TF                        for (unsigned short i = 1; i < items - 1; i++) {
347150TF                                 if (SvPOK(ST(i)) && strEQ(SvPV_nolen(ST(i)), "breaks")) {
50TF                                 if (SvPOK(ST(i)) && strEQ(SvPV_nolen(ST(i)), "breaks")) {
347750TF                        if (n_bins == 0 && looks_like_number(ST(1))) {
0TF                        if (n_bins == 0 && looks_like_number(ST(1))) {
348150TF                if (n_bins == 0) n_bins = calculate_sturges_bins(n);
348550TF                Newx(breaks,  n_bins + 1, double);
348650TF                Newx(mids,    n_bins,     double);
348750TF                Newx(density, n_bins,     double);
348850TF                Newx(counts,  n_bins,     size_t);
3492100TF                for (size_t i = 0; i <= n_bins; i++) {
3505100TF                for (size_t i = 0; i <= n_bins; i++) {
3507100TF                        if (i < n_bins) {
353550TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
100TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
354150TF                if ((items - arg_idx) % 2 != 0)
3544100TF                for (; arg_idx < items; arg_idx += 2) {
3548100TF                         if      (strEQ(key, "x"))     x_sv     = val;
354950TF                         else if (strEQ(key, "probs")) probs_sv = val;
355250TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
355650TF                if (n_raw == 0) croak("quantile: 'x' is empty");
356050TF                Newx(x, n_raw, double);
3562100TF                for (size_t i = 0; i < n_raw; i++) {
356450TF                        if (tv && SvOK(*tv)) {
50TF                        if (tv && SvOK(*tv)) {
356850TF                if (n == 0) {
357950TF                if (probs_sv && SvROK(probs_sv) && SvTYPE(SvRV(probs_sv)) == SVt_PVAV) {
50TF                if (probs_sv && SvROK(probs_sv) && SvTYPE(SvRV(probs_sv)) == SVt_PVAV) {
50TF                if (probs_sv && SvROK(probs_sv) && SvTYPE(SvRV(probs_sv)) == SVt_PVAV) {
3583100TF                        for (unsigned int i = 0; i < n_probs; i++) {
358550TF                                 probs[i] = (tv && SvOK(*tv)) ? SvNV(*tv) : 0.0;
50TF                                 probs[i] = (tv && SvOK(*tv)) ? SvNV(*tv) : 0.0;
358650TF                                 if (probs[i] < 0.0 || probs[i] > 1.0) {
50TF                                 if (probs[i] < 0.0 || probs[i] > 1.0) {
35930TF                        for (unsigned int i = 0; i < n_probs; i++) probs[i] = default_probs[i];
3599100TF                for (size_t i = 0; i < n_probs; i++) {
3602100TF                        if (n == 1) {
3604100TF                        } else if (p == 1.0) {
3606100TF                        } else if (p == 0.0) {
362050TF                        if (pct == (unsigned int)pct) {
3644100TF          for (size_t i = 0; i < items; i++) {
3646100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3649100TF                       for (size_t j = 0; j < len; j++) {
365150TF                           if (tv && SvOK(*tv)) { total += SvNV(*tv); count++; }
50TF                           if (tv && SvOK(*tv)) { total += SvNV(*tv); count++; }
365350TF                   } else if (SvOK(arg)) {
3657100TF          if (count == 0) croak("mean needs >= 1 element");
3658100TF          RETVAL = total / count;
3669100TF          for (size_t i = 0; i < items; i++) {
3671100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3674100TF                       for (size_t j = 0; j < len; j++) {
367650TF                           if (tv && SvOK(*tv)) {
50TF                           if (tv && SvOK(*tv)) {
368450TF                   } else if (SvOK(arg)) {
3692100TF          if (count < 2) croak("stdev needs >= 2 elements");
369350TF          RETVAL = sqrt(M2 / (count - 1));
3704100TF          for (size_t i = 0; i < items; i++) {
3706100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
3709100TF                       for (size_t j = 0; j < len; j++) {
371150TF                           if (tv && SvOK(*tv)) {
50TF                           if (tv && SvOK(*tv)) {
371950TF                   } else if (SvOK(arg)) {
3727100TF          if (count < 2) croak("var needs >= 2 elements");
3728100TF          RETVAL = M2 / (count - 1);
374450TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
100TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
375050TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
100TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
50TF                if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
375650TF                if ((items - arg_idx) % 2 != 0) {
3761100TF                for (; arg_idx < items; arg_idx += 2) {
3765100TF                        if      (strEQ(key, "x"))           x_sv        = val;
3766100TF                        else if (strEQ(key, "y"))           y_sv        = val;
3767100TF                        else if (strEQ(key, "mu"))          mu          = SvNV(val);
3768100TF                        else if (strEQ(key, "paired"))      paired      = SvTRUE(val);
3769100TF                        else if (strEQ(key, "var_equal"))   var_equal   = SvTRUE(val);
3770100TF                        else if (strEQ(key, "conf_level"))  conf_level  = SvNV(val);
377150TF                        else if (strEQ(key, "alternative")) alternative = SvPV_nolen(val);
3776100TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
378050TF                if (nx < 2) croak("t_test: 'x' needs at least 2 elements");
3782100TF                if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV)
50TF                if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV)
50TF                if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV)
378550TF                if (conf_level <= 0.0 || conf_level >= 1.0)
100TF                if (conf_level <= 0.0 || conf_level >= 1.0)
3790100TF                for (size_t i = 0; i < nx; i++) {
379250TF                        double val = (tv && SvOK(*tv)) ? SvNV(*tv) : 0;
50TF                        double val = (tv && SvOK(*tv)) ? SvNV(*tv) : 0;
3798100TF                if (var_x == 0.0 && !y_av) croak("t_test: data are essentially constant");
50TF                if (var_x == 0.0 && !y_av) croak("t_test: data are essentially constant");
3800100TF                if (paired || y_av) {
100TF                if (paired || y_av) {
3801100TF                        if (!y_av) croak("t_test: 'y' must be provided for paired or two-sample tests");
3803100TF                        if (paired && ny != nx) croak("t_test: Paired arrays must be same length");
100TF                        if (paired && ny != nx) croak("t_test: Paired arrays must be same length");
3805100TF                        for (size_t i = 0; i < ny; i++) {
380750TF                                 double val = (tv && SvOK(*tv)) ? SvNV(*tv) : 0;
50TF                                 double val = (tv && SvOK(*tv)) ? SvNV(*tv) : 0;
3813100TF                        if (paired) {
3815100TF                                 for (size_t i = 0; i < nx; i++) {
381850TF                                     double dx = (dx_ptr && SvOK(*dx_ptr)) ? SvNV(*dx_ptr) : 0.0;
50TF                                     double dx = (dx_ptr && SvOK(*dx_ptr)) ? SvNV(*dx_ptr) : 0.0;
381950TF                                     double dy = (dy_ptr && SvOK(*dy_ptr)) ? SvNV(*dy_ptr) : 0.0;
50TF                                     double dy = (dy_ptr && SvOK(*dy_ptr)) ? SvNV(*dy_ptr) : 0.0;
382650TF                                 if (var_d == 0.0) croak("t_test: data are essentially constant");
3832100TF                        } else if (var_equal) {
383350TF                                 if (var_x == 0.0 && var_y == 0.0) croak("t_test: data are essentially constant");
0TF                                 if (var_x == 0.0 && var_y == 0.0) croak("t_test: data are essentially constant");
384250TF                                 if (var_x == 0.0 && var_y == 0.0) croak("t_test: data are essentially constant");
0TF                                 if (var_x == 0.0 && var_y == 0.0) croak("t_test: data are essentially constant");
3862100TF                if (strcmp(alternative, "less") == 0) {
3866100TF                } else if (strcmp(alternative, "greater") == 0) {
3889100TF                if (!SvROK(p_sv) || SvTYPE(SvRV(p_sv)) != SVt_PVAV) {
50TF                if (!SvROK(p_sv) || SvTYPE(SvRV(p_sv)) != SVt_PVAV) {
3895100TF                if (n == 0) {
3901100TF                for(unsigned short int i = 0; meth[i]; i++) meth[i] = tolower(meth[i]);
3903100TF                if (strstr(meth, "benjamini") && strstr(meth, "hochberg")) strcpy(meth, "bh");
100TF                if (strstr(meth, "benjamini") && strstr(meth, "hochberg")) strcpy(meth, "bh");
3904100TF                if (strstr(meth, "benjamini") && strstr(meth, "yekutieli")) strcpy(meth, "by");
50TF                if (strstr(meth, "benjamini") && strstr(meth, "yekutieli")) strcpy(meth, "by");
390550TF                if (strcmp(meth, "fdr") == 0) strcpy(meth, "bh");
390950TF                Newx(arr, n, PVal);
391050TF                Newx(adj, n, double);
3912100TF                for (size_t i = 0; i < n; i++) {
391450TF                        arr[i].p = (tv && SvOK(*tv)) ? SvNV(*tv) : 1.0;
50TF                        arr[i].p = (tv && SvOK(*tv)) ? SvNV(*tv) : 1.0;
3920100TF                if (strcmp(meth, "bonferroni") == 0) {
3921100TF                        for (size_t i = 0; i < n; i++) {
3923100TF                                adj[arr[i].orig_idx] = (v < 1.0) ? v : 1.0;
3925100TF                } else if (strcmp(meth, "holm") == 0) {
3927100TF                        for (size_t i = 0; i < n; i++) {
3929100TF                                 if (v > cummax) cummax = v;
3930100TF                                 adj[arr[i].orig_idx] = (cummax < 1.0) ? cummax : 1.0;
3932100TF                } else if (strcmp(meth, "hochberg") == 0) {
3934100TF                        for (ssize_t i = n - 1; i >= 0; i--) {
3936100TF                                 if (v < cummin) cummin = v;
393750TF                                 adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0;
3939100TF                } else if (strcmp(meth, "bh") == 0) {
3941100TF                        for (ssize_t i = n - 1; i >= 0; i--) {
3943100TF                                if (v < cummin) cummin = v;
394450TF                                adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0;
3946100TF                } else if (strcmp(meth, "by") == 0) {
3948100TF                        for (size_t i = 1; i <= n; i++) q += 1.0 / i;
3950100TF                        for (ssize_t i = n - 1; i >= 0; i--) {
3952100TF                                if (v < cummin) cummin = v;
3953100TF                                adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0;
3955100TF                } else if (strcmp(meth, "hommel") == 0) {
395750TF                        Newx(pa, n, double);
395850TF                        Newx(q_arr, n, double);
3961100TF                        for (size_t i = 1; i < n; i++) {
396350TF                                if (temp < min_val) {
3968100TF                        for (size_t i = 0; i < n; i++) {
3972100TF                        for (size_t j = n - 1; j >= 2; j--) {
3977100TF                                 for (size_t k = 1; k < i2_len; k++) {
3979100TF                                     if (temp_q1 < q1) {
3984100TF                                 for (size_t i = 0; i <= n_mj; i++) {
3986100TF                                     q_arr[i] = (v < q1) ? v : q1;
3989100TF                                 for (size_t i = 0; i < i2_len; i++) {
3993100TF                                for (size_t i = 0; i < n; i++) {
3994100TF                                    if (pa[i] < q_arr[i]) {
4000100TF                        for (size_t i = 0; i < n; i++) {
4001100TF                                double v = (pa[i] > arr[i].p) ? pa[i] : arr[i].p;
400250TF                                if (v > 1.0) v = 1.0;
400650TF                } else if (strcmp(meth, "none") == 0) {
40070TF                        for (size_t i = 0; i < n; i++) {
401550TF                EXTEND(SP, n);
4016100TF                for (size_t i = 0; i < n; i++) {
4030100TF          for (size_t i = 0; i < items; i++) {
4032100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
4035100TF                       for (size_t j = 0; j < len; j++) {
403750TF                           if (tv && SvOK(*tv)) { total_count++; }
50TF                           if (tv && SvOK(*tv)) { total_count++; }
403950TF                   } else if (SvOK(arg)) {
4043100TF          if (total_count == 0) croak("median needs >= 1 element");
404550TF          Newx(nums, total_count, double);
4047100TF          for (size_t i = 0; i < items; i++) {
4049100TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
50TF                   if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
4052100TF                       for (size_t j = 0; j < len; j++) {
405450TF                           if (tv && SvOK(*tv)) {
50TF                           if (tv && SvOK(*tv)) {
405850TF                   } else if (SvOK(arg)) {
4064100TF          if (total_count % 2 == 0) {
4070100TF          RETVAL = median_val;
4077100TF          if (strcmp(method, "pearson")  != 0 &&
4078100TF                   strcmp(method, "spearman") != 0 &&
4079100TF                   strcmp(method, "kendall")  != 0)
408450TF          if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF          if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
408950TF          if (nx == 0) croak("cor: x is empty");
409550TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
100TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
50TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
410050TF          bool has_y = (SvOK(y_sv) && SvROK(y_sv) &&
50TF          bool has_y = (SvOK(y_sv) && SvROK(y_sv) &&
410150TF                            SvTYPE(SvRV(y_sv)) == SVt_PVAV);
410350TF          AV*restrict y_av = has_y ? (AV*)SvRV(y_sv) : NULL;
410450TF          size_t ny = has_y ? av_len(y_av) + 1 : 0;
410750TF          if (has_y && ny > 0) {
50TF          if (has_y && ny > 0) {
410950TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
100TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
50TF                   if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
4115100TF          if (!x_is_matrix && !y_is_matrix) {
50TF          if (!x_is_matrix && !y_is_matrix) {
411650TF                   if (!has_y) {
4120100TF                       if (nx != ny)
412450TF                       if (nx < 2)
412850TF                       Newx(xd, nx, double);
412950TF                       Newx(yd, ny, double);
4131100TF                       for (size_t i = 0; i < nx; i++) {
413350TF                           xd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
50TF                           xd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
50TF                           xd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
4135100TF                       for (size_t i = 0; i < ny; i++) {
413750TF                           yd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
50TF                           yd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
50TF                           yd[i] = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN;
414650TF                   if (!x_is_matrix)
415150TF                   if (!xr0 || !SvROK(*xr0) || SvTYPE(SvRV(*xr0)) != SVt_PVAV)
50TF                   if (!xr0 || !SvROK(*xr0) || SvTYPE(SvRV(*xr0)) != SVt_PVAV)
50TF                   if (!xr0 || !SvROK(*xr0) || SvTYPE(SvRV(*xr0)) != SVt_PVAV)
415550TF                   if (ncols_x == 0) croak("cor: x matrix has zero columns");
4160100TF                   for (size_t i = 0; i < nrows; i++) {
416250TF                       if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
50TF                       if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
50TF                       if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
416650TF           if (has_y && y_is_matrix) {
50TF           if (has_y && y_is_matrix) {
416750TF                       if (ny != nrows) croak("cor: x and y must have the same number of rows (%lu vs %lu)", nrows, ny);
4168100TF               for (size_t i = 0; i < nrows; i++) {
417050TF                           if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
50TF                           if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
50TF                           if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV)
417750TF                   Newx(col_x, ncols_x, double*);
4179100TF                   for (size_t j = 0; j < ncols_x; j++) {
418050TF                       Newx(col_x[j], nrows, double);
4181100TF                       for (size_t i = 0; i < nrows; i++) {
418550TF                           col_x[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
50TF                           col_x[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
50TF                           col_x[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
419550TF                   if (has_y && y_is_matrix) {
50TF                   if (has_y && y_is_matrix) {
419950TF                       if (ncols_y == 0) croak("cor: y matrix has zero columns");
420150TF                       Newx(col_y, ncols_y, double*);
4202100TF                       for (size_t j = 0; j < ncols_y; j++) {
420350TF                           Newx(col_y[j], nrows, double);
4204100TF                           for (size_t i = 0; i < nrows; i++) {
420850TF                               col_y[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
50TF                               col_y[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
50TF                               col_y[j][i] = (cv && SvOK(*cv) && looks_like_number(*cv)) ? SvNV(*cv) : NAN;
421650TF                   if (nrows < 2)
422350TF                   Newx(rows_out, ncols_x, AV*);
4224100TF                   for (size_t i = 0; i < ncols_x; i++) {
422850TF                   if (symmetric) {
42310TF                       Newx(r_cache, ncols_x, double*);
42320TF                       for (size_t i = 0; i < ncols_x; i++)
42330TF                           Newx(r_cache[i], ncols_x, double);
42350TF                       for (size_t i = 0; i < ncols_x; i++) {
42370TF                           for (size_t j = i + 1; j < ncols_x; j++) {
42440TF                       for (size_t i = 0; i < ncols_x; i++)
42450TF                           for (size_t j = 0; j < ncols_x; j++)
42480TF                       for (size_t i = 0; i < ncols_x; i++) Safefree(r_cache[i]);
4252100TF                       for (size_t i = 0; i < ncols_x; i++)
4253100TF                           for (size_t j = 0; j < ncols_y; j++)
4257100TF                   for (size_t i = 0; i < ncols_x; i++)
4261100TF                   for (size_t j = 0; j < ncols_x; j++) Safefree(col_x[j]);
426350TF                   if (!symmetric) {
4264100TF                       for (size_t j = 0; j < ncols_y; j++) Safefree(col_y[j]);
428050TF                if (items > 0) {
4282100TF                        if (SvROK(last_arg) && SvTYPE(SvRV(last_arg)) == SVt_PVHV) {
100TF                        if (SvROK(last_arg) && SvTYPE(SvRV(last_arg)) == SVt_PVHV) {
428750TF                                 if (center_sv) {
428950TF                                     if (!SvOK(val_sv)) {
429450TF                                         if (strcasecmp(str, "mean") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
50TF                                         if (strcasecmp(str, "mean") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
100TF                                         if (strcasecmp(str, "mean") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
429650TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
42980TF                                         } else if (looks_like_number(val_sv)) {
43000TF                                         } else if (SvTRUE(val_sv)) {
4309100TF                                 if (scale_sv) {
431150TF                                     if (!SvOK(val_sv)) {
431550TF                                         if (strcasecmp(str, "sd") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
50TF                                         if (strcasecmp(str, "sd") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
50TF                                         if (strcasecmp(str, "sd") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
431750TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
50TF                                         } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
43190TF                                         } else if (looks_like_number(val_sv)) {
43210TF                                             if (scale_val == 0.0) scale_val = 1.0; /* Prevent Division By Zero */
43220TF                                         } else if (SvTRUE(val_sv)) {
4333100TF                if (data_items == 1) {
4335100TF                        if (SvROK(first_arg) && SvTYPE(SvRV(first_arg)) == SVt_PVAV) {
50TF                        if (SvROK(first_arg) && SvTYPE(SvRV(first_arg)) == SVt_PVAV) {
433750TF                                 if (av_len(av) >= 0) {
433950TF                                     if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
50TF                                     if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
50TF                                     if (first_elem && SvROK(*first_elem) && SvTYPE(SvRV(*first_elem)) == SVt_PVAV) {
4345100TF                if (is_matrix) {
435550TF                        if (nrow == 0 || ncol == 0) croak("scale requires non-empty matrix");
50TF                        if (nrow == 0 || ncol == 0) croak("scale requires non-empty matrix");
4361100TF                        for (size_t r = 0; r < nrow; r++) {
4367100TF                        for (size_t c = 0; c < ncol; c++) {
437050TF                                 Newx(col_data, nrow, double);
4372100TF                                 for (size_t r = 0; r < nrow; r++) {
437450TF                                     if (row_sv && SvROK(*row_sv)) {
50TF                                     if (row_sv && SvROK(*row_sv)) {
437750TF                                         col_data[r] = (cell_sv && SvOK(*cell_sv)) ? SvNV(*cell_sv) : 0.0;
50TF                                         col_data[r] = (cell_sv && SvOK(*cell_sv)) ? SvNV(*cell_sv) : 0.0;
438450TF                                 double col_center = do_center_mean ? (col_sum / nrow) : center_val;
438750TF                                 if (do_scale_sd) {
438850TF                                     if (nrow <= 1) {
4394100TF                                     for (size_t r = 0; r < nrow; r++) {
4401100TF                                 for (size_t r = 0; r < nrow; r++) {
440350TF                                     double final_val = (col_scale == 0.0) ? (0.0 / 0.0) : (centered / col_scale);
441050TF                        EXTEND(SP, 1);
4419100TF                        for (size_t i = 0; i < data_items; i++) {
442150TF                                if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
0TF                                if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
44240TF                                  for (unsigned int j = 0; j < len; j++) {
44260TF                                                if (tv && SvOK(*tv)) { total_count++; }
0TF                                                if (tv && SvOK(*tv)) { total_count++; }
442850TF                                } else if (SvOK(arg)) {
443250TF                        if (total_count == 0) croak("scale requires at least 1 numeric element");
443350TF                        Newx(nums, total_count, double);
4434100TF                        for (size_t i = 0; i < data_items; i++) {
443650TF                                 if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
0TF                                 if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
44390TF                                     for (size_t j = 0; j < len; j++) {
44410TF                                         if (tv && SvOK(*tv)) {
0TF                                         if (tv && SvOK(*tv)) {
444650TF                                 } else if (SvOK(arg)) {
4452100TF                        if (do_center_mean) center_val = sum / total_count;
4454100TF                        if (do_scale_sd) {
4455100TF                                 if (total_count <= 1) {
4460100TF                                 for (size_t i = 0; i < total_count; i++) {
446650TF                        EXTEND(SP, total_count);
4467100TF                        for (size_t i = 0; i < total_count; i++) {
446950TF                                double final_val = (scale_val == 0.0) ? (0.0 / 0.0) : (centered / scale_val);
447950TF        if (items % 2 != 0) {
4486100TF        for (size_t i = 0; i < items; i += 2) {
4489100TF          if (strEQ(key, "data")) {
4491100TF          } else if (strEQ(key, "nrow")) {
449450TF          } else if (strEQ(key, "ncol")) {
44970TF          } else if (strEQ(key, "byrow")) {
450450TF        if (!data_sv || !SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
100TF        if (!data_sv || !SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
50TF        if (!data_sv || !SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
450850TF        size_t data_len = (UV)(av_top_index(data_av) + 1);
4509100TF        if (data_len == 0) {
451350TF        if (!nrow_set && !ncol_set) {
0TF        if (!nrow_set && !ncol_set) {
451650TF        } else if (nrow_set && !ncol_set) {
100TF        } else if (nrow_set && !ncol_set) {
451850TF        } else if (!nrow_set && ncol_set) {
0TF        } else if (!nrow_set && ncol_set) {
4522100TF        if (nrow == 0 || ncol == 0) {
50TF        if (nrow == 0 || ncol == 0) {
4530100TF        for (r = 0; r < nrow; r++) {
4537100TF        for (size_t i = 0; i < total_cells; i++) {
454050TF          SV*restrict val = fetched ? newSVsv(*fetched) : newSV(0);
454150TF          if (byrow) {
458450TF        if (items % 2 != 0) croak("Usage: lm(formula => 'mpg ~ wt * hp', data => \\%%mtcars)");
4586100TF        for (unsigned short i_arg = 0; i_arg < items; i_arg += 2) {
4589100TF          if      (strEQ(key, "formula")) formula = SvPV_nolen(val);
459050TF          else if (strEQ(key, "data"))    data_sv = val;
4593100TF        if (!formula) croak("lm: formula is required");
4594100TF        if (!data_sv || !SvROK(data_sv)) croak("lm: data is required and must be a reference");
100TF        if (!data_sv || !SvROK(data_sv)) croak("lm: data is required and must be a reference");
460050TF        if (SvTYPE(ref) == SVt_PVHV) {
460250TF          if (hv_iterinit(hv) == 0) croak("lm: Data hash is empty");
460450TF          if (entry) {
460650TF                   if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
100TF                   if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
460950TF                       Newx(row_names, n, char*);
4610100TF                       for (i = 0; i < n; i++) {
461550TF                   } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
50TF                   } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
461750TF                       Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
50TF                       Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
4619100TF                       while ((entry = hv_iternext(hv))) {
46270TF        } else if (SvTYPE(ref) == SVt_PVAV) {
46290TF          Newx(row_names, n, char*);
46300TF          Newx(row_hashes, n, HV*);
46310TF          for (i = 0; i < n; i++) {
46330TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                   if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
46380TF                       for (k = 0; k < i; k++) Safefree(row_names[k]);
4649100TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
100TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
50TF        while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
4653100TF        if (!tilde) {
4654100TF          for (i = 0; i < n; i++) Safefree(row_names[i]);
465550TF          Safefree(row_names); if (row_hashes) Safefree(row_hashes);
4667100TF          while (*p_idx) {
466950TF                   if (p_idx[0] == 'I' && p_idx[1] == '(') {
0TF                   if (p_idx[0] == 'I' && p_idx[1] == '(') {
46710TF                       while (*p_idx) { if (*p_idx == '(') depth++; else if (*p_idx == ')') { depth--; if (depth == 0) { p_idx++; break; } } p_idx++; }
0TF                       while (*p_idx) { if (*p_idx == '(') depth++; else if (*p_idx == ')') { depth--; if (depth == 0) { p_idx++; break; } } p_idx++; }
0TF                       while (*p_idx) { if (*p_idx == '(') depth++; else if (*p_idx == ')') { depth--; if (depth == 0) { p_idx++; break; } } p_idx++; }
0TF                       while (*p_idx) { if (*p_idx == '(') depth++; else if (*p_idx == ')') { depth--; if (depth == 0) { p_idx++; break; } } p_idx++; }
4675100TF                   if (p_idx[0] == '-' && p_idx[1] == '1' &&
50TF                   if (p_idx[0] == '-' && p_idx[1] == '1' &&
467650TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
4682100TF                   if (p_idx[0] == '+' && p_idx[1] == '0' &&
50TF                   if (p_idx[0] == '+' && p_idx[1] == '0' &&
46830TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
4689100TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '+') {
50TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '+') {
0TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '+') {
4695100TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '\0') {
50TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '\0') {
0TF                   if (p_idx == rhs && p_idx[0] == '0' && p_idx[1] == '\0') {
4699100TF                   if (p_idx[0] == '+' && p_idx[1] == '1' &&
50TF                   if (p_idx[0] == '+' && p_idx[1] == '1' &&
47000TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
0TF                       (p_idx[2] == '\0' || p_idx[2] == '+' || p_idx[2] == '-')) {
4705100TF                   if (p_idx == rhs) {
470650TF                       if (p_idx[0] == '1' && p_idx[1] == '\0') { p_idx[0] = '\0'; break; }
0TF                       if (p_idx[0] == '1' && p_idx[1] == '\0') { p_idx[0] = '\0'; break; }
470750TF                       if (p_idx[0] == '1' && p_idx[1] == '+') { memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); continue; }
0TF                       if (p_idx[0] == '1' && p_idx[1] == '+') { memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); continue; }
471650TF          while ((p_idx = strstr(rhs, "++")) != NULL)
471850TF          if (rhs[0] == '+') memmove(rhs, rhs + 1, strlen(rhs + 1) + 1);
472050TF          if (len_rhs > 0 && rhs[len_rhs - 1] == '+') rhs[len_rhs - 1] = '\0';
50TF          if (len_rhs > 0 && rhs[len_rhs - 1] == '+') rhs[len_rhs - 1] = '\0';
4727100TF        while (chunk != NULL) {
4728100TF          if (strcmp(chunk, ".") == 0) {
4730100TF                   for (size_t c = 0; c <= (size_t)av_len(cols); c++) {
473250TF                       if (col_sv && SvOK(*col_sv)) {
50TF                       if (col_sv && SvOK(*col_sv)) {
4734100TF                           if (strcmp(col_name, lhs) != 0) {
473650TF                               if (rhs_len + slen + 2 < sizeof(rhs_expanded)) {
4737100TF                                   if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
474750TF                   if (rhs_len + slen + 2 < sizeof(rhs_expanded)) {
4748100TF                       if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
4760100TF        if (has_intercept) { terms[num_terms++] = savepv("Intercept"); }
476250TF        if (strlen(rhs_expanded) > 0) {
4764100TF          while (chunk != NULL) {
476550TF                   if (num_terms >= term_cap - 3) {
4770100TF                   if (star) {
477550TF                       if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
0TF                       if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
477750TF                       if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
50TF                       if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
478550TF                       if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
0TF                       if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
4792100TF        for (i = 0; i < num_terms; i++) {
479450TF          for (j = 0; j < num_uniq; j++) { if (strcmp(terms[i], uniq_terms[j]) == 0) { found = true; break; } }
100TF          for (j = 0; j < num_uniq; j++) { if (strcmp(terms[i], uniq_terms[j]) == 0) { found = true; break; } }
479550TF          if (!found) uniq_terms[num_uniq++] = savepv(terms[i]);
4802100TF        for (j = 0; j < p; j++) {
480350TF          if (p_exp + 32 >= exp_cap) {
4808100TF          if (strcmp(uniq_terms[j], "Intercept") == 0) {
4811100TF          if (is_column_categorical(data_hoa, row_hashes, n, uniq_terms[j])) {
4815100TF                   for (i = 0; i < n; i++) {
481750TF                       if (str_val) {
4819100TF                           for (l = 0; l < num_levels; l++) { if (strcmp(levels[l], str_val) == 0) { found = true; break; } }
100TF                           for (l = 0; l < num_levels; l++) { if (strcmp(levels[l], str_val) == 0) { found = true; break; } }
4820100TF                           if (!found) {
482150TF                               if (num_levels >= levels_cap) { levels_cap *= 2; Renew(levels, levels_cap, char*); }
482750TF                   if (num_levels > 0) {
4828100TF                       for (l1 = 0; l1 < num_levels - 1; l1++)
4829100TF                           for (l2 = l1 + 1; l2 < num_levels; l2++)
4830100TF                               if (strcmp(levels[l1], levels[l2]) > 0) { char *tmp = levels[l1]; levels[l1] = levels[l2]; levels[l2] = tmp; }
4831100TF                       for (l = 1; l < num_levels; l++) {
483250TF                           if (p_exp >= exp_cap) {
4845100TF                       for (l = 0; l < num_levels; l++) Safefree(levels[l]);
485650TF        Newx(X, n * p, double); Newx(Y, n, double);
50TF        Newx(X, n * p, double); Newx(Y, n, double);
485750TF        Newx(valid_row_names, n, char*);
4862100TF        for (i = 0; i < n; i++) {
4864100TF          if (isnan(y_val)) { Safefree(row_names[i]); continue; }
4868100TF          for (j = 0; j < p; j++) {
4869100TF                   if (strcmp(exp_terms[j], "Intercept") == 0) {
4871100TF                   } else if (is_dummy[j]) {
487350TF                       if (str_val) {
4874100TF                           row_x[j] = (strcmp(str_val, dummy_level[j]) == 0) ? 1.0 : 0.0;
487950TF                       if (isnan(row_x[j])) { row_ok = false; break; }
488250TF          if (!row_ok) { Safefree(row_names[i]); Safefree(row_x); continue; }
4885100TF          for (j = 0; j < p; j++) X[valid_n * p + j] = row_x[j];
4892100TF        if (valid_n <= p) {
4893100TF          for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
4894100TF          for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms);
4895100TF          for (j = 0; j < p_exp; j++) {
489750TF                   if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
490150TF          if (row_hashes) Safefree(row_hashes);
4909100TF        for (i = 0; i < p; i++)
4910100TF          for (j = 0; j < p; j++) {
4912100TF                   for (k = 0; k < valid_n; k++) sum += X[k * p + i] * X[k * p + j];
4916100TF        for (i = 0; i < p; i++) {
4918100TF          for (k = 0; k < valid_n; k++) sum += X[k * p + i] * Y[k];
4924100TF        for (i = 0; i < p; i++) {
4925100TF          if (aliased[i]) { beta[i] = NAN; }
4928100TF                   for (j = 0; j < p; j++) if (!aliased[j]) sum += XtX[i * p + j] * XtY[j];
100TF                   for (j = 0; j < p; j++) if (!aliased[j]) sum += XtX[i * p + j] * XtY[j];
4943100TF        for (i = 0; i < valid_n; i++) sum_y += Y[i];
4946100TF        for (i = 0; i < valid_n; i++) {
4948100TF          for (j = 0; j < p; j++) if (!aliased[j]) y_hat += X[i * p + j] * beta[j];
100TF          for (j = 0; j < p; j++) if (!aliased[j]) y_hat += X[i * p + j] * beta[j];
4951100TF          double diff_m = has_intercept ? (y_hat - mean_y) : y_hat;
496050TF        rse_sq = (df_res > 0) ? (rss / (double)df_res) : NAN;
496650TF        if (final_rank != df_int && (mss + rss) > 0.0) {
50TF        if (final_rank != df_int && (mss + rss) > 0.0) {
496950TF          if (rse_sq > 0.0 && numdf > 0) {
50TF          if (rse_sq > 0.0 && numdf > 0) {
49720TF          } else if (rse_sq == 0.0) {
49760TF        } else if (final_rank == df_int) {
4980100TF        for (j = 0; j < p; j++) {
4984100TF          if (aliased[j]) {
499150TF                   double t_val = (se > 0.0) ? (beta[j] / se) : (INFINITY * (beta[j] >= 0.0 ? 1.0 : -1.0));
0TF                   double t_val = (se > 0.0) ? (beta[j] / se) : (INFINITY * (beta[j] >= 0.0 ? 1.0 : -1.0));
501150TF        if (!isnan(f_stat)) {
5021100TF        for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
5022100TF        for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms);
5023100TF        for (j = 0; j < p_exp; j++) {
5025100TF          if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
5030100TF        if (row_hashes) Safefree(row_hashes);
504450TF                if (by == 0.0) {
50450TF                        if (from == to) {
50460TF                                 EXTEND(SP, 1);
5054100TF                if ((from < to && by < 0.0) || (from > to && by > 0.0)) {
50TF                if ((from < to && by < 0.0) || (from > to && by > 0.0)) {
100TF                if ((from < to && by < 0.0) || (from > to && by > 0.0)) {
50TF                if ((from < to && by < 0.0) || (from > to && by > 0.0)) {
506250TF                if (n_elements_d < 0.0) n_elements_d = 0.0;
506550TF                EXTEND(SP, n_elements);
5066100TF                for (size_t i = 0; i < n_elements; i++) {
5076100TF        AUTO_SEED_PRNG();
508350TF        if (items > 0 && SvIOK(ST(0)) && (items == 1 || items % 2 != 0)) {
50TF        if (items > 0 && SvIOK(ST(0)) && (items == 1 || items % 2 != 0)) {
0TF        if (items > 0 && SvIOK(ST(0)) && (items == 1 || items % 2 != 0)) {
0TF        if (items > 0 && SvIOK(ST(0)) && (items == 1 || items % 2 != 0)) {
508950TF        if ((items - arg_start) % 2 != 0) {
5093100TF        for (int i = arg_start; i < items; i += 2) {
5097100TF            if      (strEQ(key, "n"))    n    = (unsigned int)SvUV(val);
5098100TF            else if (strEQ(key, "mean")) mean = SvNV(val);
509950TF            else if (strEQ(key, "sd"))   sd   = SvNV(val);
5103100TF        if (sd < 0.0) croak("rnorm: standard deviation must be non-negative");
510650TF        if (n > 0) {
5109100TF            for (size_t i = 0; i < n; ) {
5116100TF                 } while (s >= 1.0 || s == 0.0);
50TF                 } while (s >= 1.0 || s == 0.0);
5121100TF                 if (i < n) {
515650TF                if (!SvROK(data_sv)) croak("aov: data is required and must be a reference");
516250TF                if (SvTYPE(ref) == SVt_PVHV) {
516450TF                        if (hv_iterinit(hv) == 0) croak("aov: Data hash is empty");
516650TF                        if (entry) {
516850TF                                 if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
50TF                                 if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVAV) {
517150TF                                     Newx(row_names, n, char*);
5172100TF                                     for(i = 0; i < n; i++) {
51760TF                                 } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
0TF                                 } else if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV) {
51780TF                                     Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
0TF                                     Newx(row_names, n, char*); Newx(row_hashes, n, HV*);
51800TF                                     while ((entry = hv_iternext(hv))) {
51880TF                } else if (SvTYPE(ref) == SVt_PVAV) {
51910TF                        Newx(row_names, n, char*);
51920TF                        Newx(row_hashes, n, HV*);
51930TF                        for (i = 0; i < n; i++) {
51950TF                                 if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                                 if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
0TF                                 if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) {
52010TF                                     for (size_t k = 0; k < i; k++) Safefree(row_names[k]);
5212100TF                while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
100TF                while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
50TF                while (*src && (dst - f_cpy < 511)) { if (!isspace(*src)) { *dst++ = *src; } src++; }
5216100TF                if (!tilde) {
5217100TF                        for (i = 0; i < n; i++) Safefree(row_names[i]);
521850TF                        Safefree(row_names); if (row_hashes) Safefree(row_hashes);
522650TF                while ((p_idx = strstr(rhs, "-1")) != NULL) { has_intercept = false; memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); }
522750TF                while ((p_idx = strstr(rhs, "+0")) != NULL) { has_intercept = false; memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); }
522850TF                while ((p_idx = strstr(rhs, "0+")) != NULL) { has_intercept = false; memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); }
522950TF                if (rhs[0] == '0' && rhs[1] == '\0')        { has_intercept = false; rhs[0] = '\0'; }
0TF                if (rhs[0] == '0' && rhs[1] == '\0')        { has_intercept = false; rhs[0] = '\0'; }
523050TF                while ((p_idx = strstr(rhs, "+1")) != NULL) { memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); }
523150TF                if (rhs[0] == '1' && rhs[1] == '\0')        { rhs[0] = '\0'; }
0TF                if (rhs[0] == '1' && rhs[1] == '\0')        { rhs[0] = '\0'; }
523250TF                else if (rhs[0] == '1' && rhs[1] == '+')    { memmove(rhs, rhs + 2, strlen(rhs + 2) + 1); }
0TF                else if (rhs[0] == '1' && rhs[1] == '+')    { memmove(rhs, rhs + 2, strlen(rhs + 2) + 1); }
523450TF                while ((p_idx = strstr(rhs, "++")) != NULL) memmove(p_idx, p_idx + 1, strlen(p_idx + 1) + 1);
523550TF                if (rhs[0] == '+') memmove(rhs, rhs + 1, strlen(rhs + 1) + 1);
523750TF                if (len_rhs > 0 && rhs[len_rhs - 1] == '+') rhs[len_rhs - 1] = '\0';
50TF                if (len_rhs > 0 && rhs[len_rhs - 1] == '+') rhs[len_rhs - 1] = '\0';
5242100TF                while (chunk != NULL) {
5243100TF                        if (strcmp(chunk, ".") == 0) {
5245100TF                                 for (size_t c = 0; c <= av_len(cols); c++) {
524750TF                                     if (col_sv && SvOK(*col_sv)) {
50TF                                     if (col_sv && SvOK(*col_sv)) {
5249100TF                                         if (strcmp(col_name, lhs) != 0) {
525150TF                                             if (rhs_len + slen + 2 < sizeof(rhs_expanded)) {
5252100TF                                                 if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
526250TF                                 if (rhs_len + slen + 2 < sizeof(rhs_expanded)) {
5263100TF                                     if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; }
527950TF                if (has_intercept) { terms[num_terms++] = savepv("Intercept"); }
528150TF                if (strlen(rhs_expanded) > 0) {
5283100TF                        while (chunk != NULL) {
528450TF                                 if (num_terms >= term_cap - 3) {
5289100TF                                 if (star) {
529450TF                                     if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
0TF                                     if (c_l && strncmp(left, "I(", 2) != 0) *c_l = '\0';
529550TF                                     char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
0TF                                     char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0';
530350TF                                     if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
0TF                                     if (c_chunk && strncmp(chunk, "I(", 2) != 0) *c_chunk = '\0';
5310100TF                for (i = 0; i < num_terms; i++) {
5312100TF                        for (size_t k = 0; k < num_uniq; k++) {
531350TF                                  if (strcmp(terms[i], uniq_terms[k]) == 0) { found = true; break; }
531550TF                        if (!found) uniq_terms[num_uniq++] = savepv(terms[i]);
5322100TF                for (j = 0; j < p; j++) {
5323100TF                        if (p_exp + 64 >= exp_cap) {
5331100TF                        if (strcmp(uniq_terms[j], "Intercept") == 0) {
5341100TF                        if (colon) {
5349100TF                                  for (size_t e = 0; e < p_exp; e++) {
5350100TF                                      if (strcmp(parent_term[e], left) == 0) l_indices[l_count++] = e;
5351100TF                                      if (strcmp(parent_term[e], right) == 0) r_indices[r_count++] = e;
5354100TF                                  if (l_count == 0 || r_count == 0) {
50TF                                  if (l_count == 0 || r_count == 0) {
5358100TF                                      for (int li = 0; li < l_count; li++) {
5359100TF                                          for (int ri = 0; ri < r_count; ri++) {
536050TF                                              if (p_exp >= exp_cap) {
5381100TF                                  if (is_column_categorical(data_hoa, row_hashes, n, uniq_terms[j])) {
5385100TF                                      for (i = 0; i < n; i++) {
538750TF                                          if (str_val) {
5389100TF                                              for (size_t l = 0; l < num_levels; l++) {
5390100TF                                                  if (strcmp(levels[l], str_val) == 0) { found = true; break; }
5392100TF                                              if (!found) {
539350TF                                                  if (num_levels >= levels_cap) { levels_cap *= 2; Renew(levels, levels_cap, char*); }
540050TF                                      if (num_levels > 0) {
5401100TF                                          for (size_t l1 = 0; l1 < num_levels - 1; l1++) {
5402100TF                                              for (size_t l2 = l1 + 1; l2 < num_levels; l2++) {
5403100TF                                                  if (strcmp(levels[l1], levels[l2]) > 0) {
5408100TF                                          for (size_t l = 1; l < num_levels; l++) {
540950TF                                              if (p_exp >= exp_cap) {
5426100TF                                          for (size_t l = 0; l < num_levels; l++) Safefree(levels[l]);
5446100TF                for(i = 0; i < n; i++) X_mat[i] = (double*)safemalloc(p_exp * sizeof(double));
544750TF                Newx(Y, n, double);
5451100TF                for (i = 0; i < n; i++) {
545350TF                        if (isnan(y_val)) { Safefree(row_names[i]); continue; }
5456100TF                        for (j = 0; j < p_exp; j++) {
5457100TF                                  if (strcmp(exp_terms[j], "Intercept") == 0) {
5459100TF                                  } else if (is_interact[j]) {
5461100TF                                  } else if (is_dummy[j]) {
546350TF                                      if (str_val) {
5464100TF                                          row_x[j] = (strcmp(str_val, dummy_level[j]) == 0) ? 1.0 : 0.0;
546950TF                                      if (isnan(row_x[j])) { row_ok = false; break; }
547250TF                        if (!row_ok) { Safefree(row_names[i]); Safefree(row_x); continue; }
5475100TF                        for (j = 0; j < p_exp; j++) X_mat[valid_n][j] = row_x[j];
5481100TF                if (valid_n <= p_exp) {
5483100TF                        for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
5484100TF                        for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms);
5485100TF                        for (j = 0; j < p_exp; j++) {
548750TF                                 if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
5493100TF                        for(i = 0; i < n; i++) Safefree(X_mat[i]);
549550TF                        if (row_hashes) Safefree(row_hashes);
5509100TF                for (i = 0; i < p_exp; i++) {
5510100TF                        if (strcmp(exp_terms[i], "Intercept") == 0) continue;
5511100TF                        if (aliased_qr[i]) continue;
5518100TF                for (i = 0; i < p_exp; i++) {
5519100TF                        if (!aliased_qr[i]) rank++;
5522100TF                for (i = rank; i < valid_n; i++) {
552650TF                double ms_res = (res_df > 0) ? rss_prev / res_df : 0.0;
5529100TF                for (j = 0; j < num_uniq; j++) {
5530100TF                        if (strcmp(uniq_terms[j], "Intercept") == 0) continue;
5534100TF                        double ms = (df > 0) ? ss / df : 0.0;
553950TF                        if (ms_res > 0.0 && df > 0) {
100TF                        if (ms_res > 0.0 && df > 0) {
5557100TF                for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
5558100TF                for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms);
5559100TF                for (j = 0; j < p_exp; j++) {
5561100TF                        if (is_dummy[j]) { Safefree(dummy_base[j]); Safefree(dummy_level[j]); }
5568100TF                for (i = 0; i < n; i++) Safefree(X_mat[i]);
557150TF                if (row_hashes) Safefree(row_hashes);
5582100TF        if (items < 1) croak("fisher_test requires at least a data reference");
5589100TF        for (unsigned short int i = 1; i < items; i += 2) {
559050TF                if (i + 1 >= items) croak("fisher_test: odd number of arguments");
559350TF                if (strEQ(key, "conf_level") || strEQ(key, "conf.level")) {
50TF                if (strEQ(key, "conf_level") || strEQ(key, "conf.level")) {
559550TF                } else if (strEQ(key, "alternative")) {
560050TF        if (!SvROK(data_ref)) croak("fisher_test requires a reference to an Array or Hash");
5604100TF        if (SvTYPE(deref) == SVt_PVAV) {
560650TF          if (av_len(outer) != 1) croak("Outer array must have exactly 2 rows");
560950TF          if (row1_ptr && row2_ptr && SvROK(*row1_ptr) && SvROK(*row2_ptr)) {
50TF          if (row1_ptr && row2_ptr && SvROK(*row1_ptr) && SvROK(*row2_ptr)) {
50TF          if (row1_ptr && row2_ptr && SvROK(*row1_ptr) && SvROK(*row2_ptr)) {
50TF          if (row1_ptr && row2_ptr && SvROK(*row1_ptr) && SvROK(*row2_ptr)) {
561650TF                   a = (a_ptr && SvOK(*a_ptr)) ? SvIV(*a_ptr) : 0;
50TF                   a = (a_ptr && SvOK(*a_ptr)) ? SvIV(*a_ptr) : 0;
561750TF                   b = (b_ptr && SvOK(*b_ptr)) ? SvIV(*b_ptr) : 0;
50TF                   b = (b_ptr && SvOK(*b_ptr)) ? SvIV(*b_ptr) : 0;
561850TF                   c = (c_ptr && SvOK(*c_ptr)) ? SvIV(*c_ptr) : 0;
50TF                   c = (c_ptr && SvOK(*c_ptr)) ? SvIV(*c_ptr) : 0;
561950TF                   d = (d_ptr && SvOK(*d_ptr)) ? SvIV(*d_ptr) : 0;
50TF                   d = (d_ptr && SvOK(*d_ptr)) ? SvIV(*d_ptr) : 0;
562350TF        } else if (SvTYPE(deref) == SVt_PVHV) {
562650TF                if (hv_iterinit(outer) != 2) croak("Outer hash must have exactly 2 keys");
562950TF                if (!he1 || !he2) croak("Invalid outer hash");
50TF                if (!he1 || !he2) croak("Invalid outer hash");
5632100TF                HE*restrict row1_he = (strcmp(k1, k2) < 0) ? he1 : he2;
5633100TF                HE*restrict row2_he = (strcmp(k1, k2) < 0) ? he2 : he1;
563650TF                if (!SvROK(row1_sv) || SvTYPE(SvRV(row1_sv)) != SVt_PVHV ||
50TF                if (!SvROK(row1_sv) || SvTYPE(SvRV(row1_sv)) != SVt_PVHV ||
563750TF                        !SvROK(row2_sv) || SvTYPE(SvRV(row2_sv)) != SVt_PVHV) {
50TF                        !SvROK(row2_sv) || SvTYPE(SvRV(row2_sv)) != SVt_PVHV) {
564250TF                if (hv_iterinit(in1) != 2 || hv_iterinit(in2) != 2) croak("Inner hashes must have exactly 2 keys");
50TF                if (hv_iterinit(in1) != 2 || hv_iterinit(in2) != 2) croak("Inner hashes must have exactly 2 keys");
5647100TF                HE*restrict in1_c1 = (strcmp(in1_k1, in1_k2) < 0) ? in1_he1 : in1_he2;
5648100TF                HE*restrict in1_c2 = (strcmp(in1_k1, in1_k2) < 0) ? in1_he2 : in1_he1;
5653100TF                HE*restrict in2_c1 = (strcmp(in2_k1, in2_k2) < 0) ? in2_he1 : in2_he2;
5654100TF                HE*restrict in2_c2 = (strcmp(in2_k1, in2_k2) < 0) ? in2_he2 : in2_he1;
565550TF                a = (hv_iterval(in1, in1_c1) && SvOK(hv_iterval(in1, in1_c1))) ? SvIV(hv_iterval(in1, in1_c1)) : 0;
50TF                a = (hv_iterval(in1, in1_c1) && SvOK(hv_iterval(in1, in1_c1))) ? SvIV(hv_iterval(in1, in1_c1)) : 0;
565650TF                b = (hv_iterval(in1, in1_c2) && SvOK(hv_iterval(in1, in1_c2))) ? SvIV(hv_iterval(in1, in1_c2)) : 0;
50TF                b = (hv_iterval(in1, in1_c2) && SvOK(hv_iterval(in1, in1_c2))) ? SvIV(hv_iterval(in1, in1_c2)) : 0;
565750TF                c = (hv_iterval(in2, in2_c1) && SvOK(hv_iterval(in2, in2_c1))) ? SvIV(hv_iterval(in2, in2_c1)) : 0;
50TF                c = (hv_iterval(in2, in2_c1) && SvOK(hv_iterval(in2, in2_c1))) ? SvIV(hv_iterval(in2, in2_c1)) : 0;
565850TF                d = (hv_iterval(in2, in2_c2) && SvOK(hv_iterval(in2, in2_c2))) ? SvIV(hv_iterval(in2, in2_c2)) : 0;
50TF                d = (hv_iterval(in2, in2_c2) && SvOK(hv_iterval(in2, in2_c2))) ? SvIV(hv_iterval(in2, in2_c2)) : 0;
570050TF    if (items % 2 != 0) croak("Usage: power_t_test(n => 30, delta => 0.5, sd => 1.0, ...)");
5701100TF    for (unsigned short int i = 0; i < items; i += 2) {
5705100TF        if      (strEQ(key, "n"))           sv_n = val;
5706100TF        else if (strEQ(key, "delta"))       sv_delta = val;
5707100TF        else if (strEQ(key, "sd"))          sv_sd = val;
570850TF        else if (strEQ(key, "sig.level") || strEQ(key, "sig_level")) sv_sig_level = val;
100TF        else if (strEQ(key, "sig.level") || strEQ(key, "sig_level")) sv_sig_level = val;
5709100TF        else if (strEQ(key, "power"))       sv_power = val;
5710100TF        else if (strEQ(key, "type"))        type = SvPV_nolen(val);
571150TF        else if (strEQ(key, "alternative")) alternative = SvPV_nolen(val);
57120TF        else if (strEQ(key, "strict"))      strict = SvTRUE(val);
57130TF        else if (strEQ(key, "tol"))         tol = SvNV(val);
5717100TF    bool is_null_n = (!sv_n || !SvOK(sv_n));
50TF    bool is_null_n = (!sv_n || !SvOK(sv_n));
571850TF    bool is_null_delta = (!sv_delta || !SvOK(sv_delta));
50TF    bool is_null_delta = (!sv_delta || !SvOK(sv_delta));
5719100TF    bool is_null_power = (!sv_power || !SvOK(sv_power));
50TF    bool is_null_power = (!sv_power || !SvOK(sv_power));
572050TF    bool is_null_sd = (sv_sd && !SvOK(sv_sd));
50TF    bool is_null_sd = (sv_sd && !SvOK(sv_sd));
5721100TF    bool is_null_sig_level = (sv_sig_level && !SvOK(sv_sig_level));
50TF    bool is_null_sig_level = (sv_sig_level && !SvOK(sv_sig_level));
5724100TF    if (is_null_n) missing_count++;
572550TF    if (is_null_delta) missing_count++;
5726100TF    if (is_null_power) missing_count++;
572750TF    if (is_null_sd) missing_count++;
572850TF    if (is_null_sig_level) missing_count++;
573050TF    if (missing_count != 1) {
5734100TF    double n = is_null_n ? 0.0 : SvNV(sv_n);
573550TF    double delta = is_null_delta ? 0.0 : SvNV(sv_delta);
573650TF    double sd = (!sv_sd || is_null_sd) ? 1.0 : SvNV(sv_sd);
50TF    double sd = (!sv_sd || is_null_sd) ? 1.0 : SvNV(sv_sd);
5737100TF    double sig_level = (!sv_sig_level || is_null_sig_level) ? 0.05 : SvNV(sv_sig_level);
50TF    double sig_level = (!sv_sig_level || is_null_sig_level) ? 0.05 : SvNV(sv_sig_level);
5738100TF    double power = is_null_power ? 0.0 : SvNV(sv_power);
5739100TF    short int tsample = (strEQ(type, "one.sample") || strEQ(type, "paired")) ? 1 : 2;
100TF    short int tsample = (strEQ(type, "one.sample") || strEQ(type, "paired")) ? 1 : 2;
5740100TF    short int tside = (strEQ(alternative, "one.sided") || strEQ(alternative, "greater") || strEQ(alternative, "less")) ? 1 : 2;
50TF    short int tside = (strEQ(alternative, "one.sided") || strEQ(alternative, "greater") || strEQ(alternative, "less")) ? 1 : 2;
50TF    short int tside = (strEQ(alternative, "one.sided") || strEQ(alternative, "greater") || strEQ(alternative, "less")) ? 1 : 2;
5741100TF    if (tside == 2 && !is_null_delta) delta = fabs(delta);
50TF    if (tside == 2 && !is_null_delta) delta = fabs(delta);
5742100TF    if (is_null_power) {
574450TF    } else if (is_null_n) {
574650TF        while (p_body(high, delta, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
0TF        while (p_body(high, delta, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
5747100TF        while (high - low > tol) {
5749100TF            if (p_body(mid, delta, sd, sig_level, tsample, tside, strict) < power) low = mid;
57530TF    } else if (is_null_sd) {
57550TF        while (high - low > tol) {
57570TF            if (p_body(n, delta, mid, sig_level, tsample, tside, strict) > power) low = mid;
57610TF    } else if (is_null_delta) {
57630TF        while (p_body(n, high, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
0TF        while (p_body(n, high, sd, sig_level, tsample, tside, strict) < power && high < 1e12) high *= 2.0;
57640TF        while (high - low > tol) {
57660TF            if (p_body(n, mid, sd, sig_level, tsample, tside, strict) < power) low = mid;
57700TF    } else if (is_null_sig_level) {
57720TF        while (high - low > tol) {
57740TF            if (p_body(n, delta, sd, mid, tsample, tside, strict) < power) low = mid;
5786100TF    const char*restrict m_str = (tsample == 1) ? (strEQ(type, "paired") ? "Paired t test power calculation" : "One-sample t test power calculation") : "Two-sample t test power calculation";
100TF    const char*restrict m_str = (tsample == 1) ? (strEQ(type, "paired") ? "Paired t test power calculation" : "One-sample t test power calculation") : "Two-sample t test power calculation";
5788100TF    const char*restrict n_str = (tsample == 2) ? "n is number in *each* group" : (strEQ(type, "paired") ? "n is number of *pairs*, sd is std.dev. of *differences* within pairs" : "");
100TF    const char*restrict n_str = (tsample == 2) ? "n is number in *each* group" : (strEQ(type, "paired") ? "n is number of *pairs*, sd is std.dev. of *differences* within pairs" : "");
5789100TF    if (n_str[0] != '\0') hv_stores(ret, "note", newSVpv(n_str, 0));
580350TF        if (arg_idx < items && SvROK(ST(arg_idx))) {
100TF        if (arg_idx < items && SvROK(ST(arg_idx))) {
5805100TF                if (t == SVt_PVAV) {
580750TF                } else if (t == SVt_PVHV) {
5811100TF        if (!h_sv && arg_idx < items
50TF        if (!h_sv && arg_idx < items
5812100TF                  && SvROK(ST(arg_idx))
581350TF                  && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) {
5817100TF        for (; arg_idx < items; arg_idx += 2) {
5820100TF                if      (strEQ(key, "x")) x_sv = val;
582150TF                else if (strEQ(key, "g")) g_sv = val;
58220TF                else if (strEQ(key, "h")) h_sv = val;
5826100TF        if (h_sv && (x_sv || g_sv))
50TF        if (h_sv && (x_sv || g_sv))
50TF        if (h_sv && (x_sv || g_sv))
5838100TF        if (h_sv) {
583950TF                if (!SvROK(h_sv) || SvTYPE(SvRV(h_sv)) != SVt_PVHV)
50TF                if (!SvROK(h_sv) || SvTYPE(SvRV(h_sv)) != SVt_PVHV)
5846100TF                while ((he = hv_iternext(h_hv))) {
584850TF                        if (!SvROK(val) || SvTYPE(SvRV(val)) != SVt_PVAV)
50TF                        if (!SvROK(val) || SvTYPE(SvRV(val)) != SVt_PVAV)
585250TF                if (total < 2) croak("not enough observations");
5859100TF                while ((he = hv_iternext(h_hv))) {
5862100TF                        for (size_t i = 0; i < n_g; i++) {
586450TF                                if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                                if (el && SvOK(*el) && looks_like_number(*el)) {
50TF                                if (el && SvOK(*el) && looks_like_number(*el)) {
587850TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
50TF                if (!x_sv || !SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
588050TF                if (!g_sv || !SvROK(g_sv) || SvTYPE(SvRV(g_sv)) != SVt_PVAV)
50TF                if (!g_sv || !SvROK(g_sv) || SvTYPE(SvRV(g_sv)) != SVt_PVAV)
50TF                if (!g_sv || !SvROK(g_sv) || SvTYPE(SvRV(g_sv)) != SVt_PVAV)
588750TF                if (nx != ng) croak("kruskal_test: 'x' and 'g' must have the same length");
588850TF                if (nx < 2)   croak("not enough observations");
5896100TF                for (size_t i = 0; i < nx; i++) {
589950TF                        if (x_el && SvOK(*x_el) && looks_like_number(*x_el)
50TF                        if (x_el && SvOK(*x_el) && looks_like_number(*x_el)
50TF                        if (x_el && SvOK(*x_el) && looks_like_number(*x_el)
590050TF                                 && g_el && SvOK(*g_el)) {
50TF                                 && g_el && SvOK(*g_el)) {
5905100TF                                if (id_sv) {
592350TF        if (valid_n < 2) { Safefree(ri); croak("not enough observations");            }
592450TF        if (k       < 2) { Safefree(ri); croak("all observations are in the same group"); }
5933100TF        for (size_t i = 0; i < valid_n; i++) {
5941100TF        for (size_t i = 0; i < k; i++) {
594250TF                if (group_counts[i] > 0)
595050TF        if (tie_adj > 0.0) {