formualizer_eval/builtins/stats/
mod.rs

1//! Statistical basic functions (Sprint 6)
2//!
3//! Implementations target Excel semantic parity for:
4//! LARGE, SMALL, RANK.EQ, RANK.AVG, MEDIAN, STDEV.S, STDEV.P, VAR.S, VAR.P,
5//! PERCENTILE.INC, PERCENTILE.EXC, QUARTILE.INC, QUARTILE.EXC.
6//!
7//! Notes:
8//! - We currently materialize numeric values into a Vec<f64>. For large ranges this could be
9//!   optimized with streaming selection algorithms (nth_element / partial sort). TODO(perf).
10//! - Text/boolean coercion nuance: For Excel statistical functions, values coming from range
11//!   references should ignore text and logical values (they are skipped), while direct scalar
12//!   arguments still coerce (e.g. =STDEV(1,TRUE) treats TRUE as 1). This file now implements that
13//!   distinction. TODO(excel-nuance): refine numeric text literal vs non‑numeric text handling.
14//! - Errors encountered in any argument propagate immediately.
15//! - Empty numeric sets produce Excel-specific errors (#NUM! for LARGE/SMALL, #N/A for rank target
16//!   out of range, #DIV/0! for STDEV/VAR sample with n < 2, etc.).
17
18use super::super::builtins::utils::{ARG_RANGE_NUM_LENIENT_ONE, coerce_num};
19use crate::args::ArgSchema;
20use crate::function::Function;
21use crate::traits::{ArgumentHandle, FunctionContext};
22use formualizer_common::{ExcelError, LiteralValue};
23// use std::collections::BTreeMap; // removed unused import
24use formualizer_macros::func_caps;
25
26fn scalar_like_value(arg: &ArgumentHandle<'_, '_>) -> Result<LiteralValue, ExcelError> {
27    Ok(match arg.value()? {
28        crate::traits::CalcValue::Scalar(v) => v,
29        crate::traits::CalcValue::Range(rv) => rv.get_cell(0, 0),
30    })
31}
32
33/// Collect numeric inputs applying Excel statistical semantics:
34/// - Range references: include only numeric cells; skip text, logical, blank. Errors propagate.
35/// - Direct scalar arguments: attempt numeric coercion (so TRUE/FALSE, numeric text are included if
36///   coerce_num succeeds). Non-numeric text is ignored (Excel would treat a direct non-numeric text
37///   argument as #VALUE! in some contexts; covered by TODO for finer parity).
38fn collect_numeric_stats(args: &[ArgumentHandle]) -> Result<Vec<f64>, ExcelError> {
39    let mut out = Vec::new();
40    for a in args {
41        // Special-case: inline array literal argument should be treated like a list of direct scalar
42        // arguments (not a by-ref range). This allows boolean/text coercion per element akin to
43        // passing multiple scalars to the function.
44        if let Some(arr) = a.inline_array_literal()? {
45            for row in arr.into_iter() {
46                for cell in row.into_iter() {
47                    match cell {
48                        LiteralValue::Error(e) => return Err(e),
49                        other => {
50                            if let Ok(n) = coerce_num(&other) {
51                                out.push(n);
52                            }
53                        }
54                    }
55                }
56            }
57            continue;
58        }
59
60        if let Ok(view) = a.range_view() {
61            view.for_each_cell(&mut |v| {
62                match v {
63                    LiteralValue::Error(e) => return Err(e.clone()),
64                    LiteralValue::Number(n) => out.push(*n),
65                    LiteralValue::Int(i) => out.push(*i as f64),
66                    _ => {}
67                }
68                Ok(())
69            })?;
70        } else {
71            let v = scalar_like_value(a)?;
72            match v {
73                LiteralValue::Error(e) => return Err(e),
74                other => {
75                    if let Ok(n) = coerce_num(&other) {
76                        out.push(n);
77                    }
78                }
79            }
80        }
81    }
82    Ok(out)
83}
84
85fn percentile_inc(sorted: &[f64], p: f64) -> Result<f64, ExcelError> {
86    if sorted.is_empty() {
87        return Err(ExcelError::new_num());
88    }
89    if !(0.0..=1.0).contains(&p) {
90        return Err(ExcelError::new_num());
91    }
92    if sorted.len() == 1 {
93        return Ok(sorted[0]);
94    }
95    let n = sorted.len() as f64;
96    let rank = p * (n - 1.0); // 0-based rank
97    let lo = rank.floor() as usize;
98    let hi = rank.ceil() as usize;
99    if lo == hi {
100        return Ok(sorted[lo]);
101    }
102    let frac = rank - (lo as f64);
103    Ok(sorted[lo] + (sorted[hi] - sorted[lo]) * frac)
104}
105
106fn percentile_exc(sorted: &[f64], p: f64) -> Result<f64, ExcelError> {
107    // Excel PERCENTILE.EXC requires 0 < p < 1 and uses (n+1) basis; invalid if rank<1 or >n
108    if sorted.is_empty() {
109        return Err(ExcelError::new_num());
110    }
111    if !(0.0..=1.0).contains(&p) || p <= 0.0 || p >= 1.0 {
112        return Err(ExcelError::new_num());
113    }
114    let n = sorted.len() as f64;
115    let rank = p * (n + 1.0); // 1..n domain
116    if rank < 1.0 || rank > n {
117        return Err(ExcelError::new_num());
118    }
119    let lo = rank.floor();
120    let hi = rank.ceil();
121    if (lo - hi).abs() < f64::EPSILON {
122        return Ok(sorted[(lo as usize) - 1]);
123    }
124    let frac = rank - lo;
125    let lo_idx = (lo as usize) - 1;
126    let hi_idx = (hi as usize) - 1;
127    Ok(sorted[lo_idx] + (sorted[hi_idx] - sorted[lo_idx]) * frac)
128}
129
130/// RANK.EQ(number, ref, [order]) Excel semantics:
131/// - order omitted or 0 => descending (largest value rank 1)
132/// - order non-zero => ascending (smallest value rank 1)
133/// - ties return same rank (position of first in ordering)
134#[derive(Debug)]
135pub struct RankEqFn;
136impl Function for RankEqFn {
137    func_caps!(PURE, NUMERIC_ONLY);
138    fn name(&self) -> &'static str {
139        "RANK.EQ"
140    }
141    fn aliases(&self) -> &'static [&'static str] {
142        &["RANK"]
143    }
144    fn min_args(&self) -> usize {
145        2
146    }
147    fn variadic(&self) -> bool {
148        true
149    } // allow optional order
150    fn arg_schema(&self) -> &'static [ArgSchema] {
151        &ARG_RANGE_NUM_LENIENT_ONE[..]
152    }
153    fn eval<'a, 'b, 'c>(
154        &self,
155        args: &'c [ArgumentHandle<'a, 'b>],
156        _ctx: &dyn FunctionContext<'b>,
157    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
158        if args.len() < 2 {
159            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
160                ExcelError::new_na(),
161            )));
162        }
163        let target = match coerce_num(&args[0].value()?.into_literal()) {
164            Ok(n) => n,
165            Err(_) => {
166                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
167                    ExcelError::new_na(),
168                )));
169            }
170        };
171        // optional order arg at end if 3 args
172        let order = if args.len() >= 3 {
173            coerce_num(&args[2].value()?.into_literal()).unwrap_or(0.0)
174        } else {
175            0.0
176        };
177        let nums = collect_numeric_stats(&args[1..2])?; // only one ref range per Excel spec
178        if nums.is_empty() {
179            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
180                ExcelError::new_na(),
181            )));
182        }
183        let mut sorted = nums; // copy
184        if order.abs() < 1e-12 {
185            // descending
186            sorted.sort_by(|a, b| b.partial_cmp(a).unwrap());
187        } else {
188            // ascending
189            sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
190        }
191        for (i, &v) in sorted.iter().enumerate() {
192            if (v - target).abs() < 1e-12 {
193                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
194                    (i + 1) as f64,
195                )));
196            }
197        }
198        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
199            ExcelError::new_na(),
200        )))
201    }
202}
203
204/// RANK.AVG(number, ref, [order]) ties return average of ranks
205#[derive(Debug)]
206pub struct RankAvgFn;
207impl Function for RankAvgFn {
208    func_caps!(PURE, NUMERIC_ONLY);
209    fn name(&self) -> &'static str {
210        "RANK.AVG"
211    }
212    fn min_args(&self) -> usize {
213        2
214    }
215    fn variadic(&self) -> bool {
216        true
217    }
218    fn arg_schema(&self) -> &'static [ArgSchema] {
219        &ARG_RANGE_NUM_LENIENT_ONE[..]
220    }
221    fn eval<'a, 'b, 'c>(
222        &self,
223        args: &'c [ArgumentHandle<'a, 'b>],
224        _ctx: &dyn FunctionContext<'b>,
225    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
226        if args.len() < 2 {
227            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
228                ExcelError::new_na(),
229            )));
230        }
231        let t0 = scalar_like_value(&args[0])?;
232        let target = match coerce_num(&t0) {
233            Ok(n) => n,
234            Err(_) => {
235                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
236                    ExcelError::new_na(),
237                )));
238            }
239        };
240        let order = if args.len() >= 3 {
241            let ord = scalar_like_value(&args[2])?;
242            coerce_num(&ord).unwrap_or(0.0)
243        } else {
244            0.0
245        };
246        let nums = collect_numeric_stats(&args[1..2])?;
247        if nums.is_empty() {
248            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
249                ExcelError::new_na(),
250            )));
251        }
252        let mut sorted = nums;
253        if order.abs() < 1e-12 {
254            sorted.sort_by(|a, b| b.partial_cmp(a).unwrap());
255        } else {
256            sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
257        }
258        let mut positions = Vec::new();
259        for (i, &v) in sorted.iter().enumerate() {
260            if (v - target).abs() < 1e-12 {
261                positions.push(i + 1);
262            }
263        }
264        if positions.is_empty() {
265            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
266                ExcelError::new_na(),
267            )));
268        }
269        let avg = positions.iter().copied().sum::<usize>() as f64 / positions.len() as f64;
270        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(avg)))
271    }
272}
273
274#[derive(Debug)]
275pub struct LARGE;
276impl Function for LARGE {
277    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
278    fn name(&self) -> &'static str {
279        "LARGE"
280    }
281    fn min_args(&self) -> usize {
282        2
283    }
284    fn variadic(&self) -> bool {
285        true
286    }
287    fn arg_schema(&self) -> &'static [ArgSchema] {
288        &ARG_RANGE_NUM_LENIENT_ONE[..]
289    }
290    fn eval<'a, 'b, 'c>(
291        &self,
292        args: &'c [ArgumentHandle<'a, 'b>],
293        _ctx: &dyn FunctionContext<'b>,
294    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
295        if args.len() < 2 {
296            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
297                ExcelError::new_num(),
298            )));
299        }
300        let k = match coerce_num(&args.last().unwrap().value()?.into_literal()) {
301            Ok(n) => n,
302            Err(_) => {
303                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
304                    ExcelError::new_num(),
305                )));
306            }
307        };
308        let k = k as i64;
309        if k < 1 {
310            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
311                ExcelError::new_num(),
312            )));
313        }
314        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
315        if nums.is_empty() || k as usize > nums.len() {
316            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
317                ExcelError::new_num(),
318            )));
319        }
320        nums.sort_by(|a, b| b.partial_cmp(a).unwrap());
321        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
322            nums[(k as usize) - 1],
323        )))
324    }
325}
326
327#[derive(Debug)]
328pub struct SMALL;
329impl Function for SMALL {
330    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
331    fn name(&self) -> &'static str {
332        "SMALL"
333    }
334    fn min_args(&self) -> usize {
335        2
336    }
337    fn variadic(&self) -> bool {
338        true
339    }
340    fn arg_schema(&self) -> &'static [ArgSchema] {
341        &ARG_RANGE_NUM_LENIENT_ONE[..]
342    }
343    fn eval<'a, 'b, 'c>(
344        &self,
345        args: &'c [ArgumentHandle<'a, 'b>],
346        _ctx: &dyn FunctionContext<'b>,
347    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
348        if args.len() < 2 {
349            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
350                ExcelError::new_num(),
351            )));
352        }
353        let k = match coerce_num(&args.last().unwrap().value()?.into_literal()) {
354            Ok(n) => n,
355            Err(_) => {
356                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
357                    ExcelError::new_num(),
358                )));
359            }
360        };
361        let k = k as i64;
362        if k < 1 {
363            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
364                ExcelError::new_num(),
365            )));
366        }
367        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
368        if nums.is_empty() || k as usize > nums.len() {
369            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
370                ExcelError::new_num(),
371            )));
372        }
373        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
374        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
375            nums[(k as usize) - 1],
376        )))
377    }
378}
379
380#[derive(Debug)]
381pub struct MEDIAN;
382impl Function for MEDIAN {
383    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
384    fn name(&self) -> &'static str {
385        "MEDIAN"
386    }
387    fn min_args(&self) -> usize {
388        1
389    }
390    fn variadic(&self) -> bool {
391        true
392    }
393    fn arg_schema(&self) -> &'static [ArgSchema] {
394        &ARG_RANGE_NUM_LENIENT_ONE[..]
395    }
396    fn eval<'a, 'b, 'c>(
397        &self,
398        args: &'c [ArgumentHandle<'a, 'b>],
399        _ctx: &dyn FunctionContext<'b>,
400    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
401        let mut nums = collect_numeric_stats(args)?;
402        if nums.is_empty() {
403            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
404                ExcelError::new_num(),
405            )));
406        }
407        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
408        let n = nums.len();
409        let mid = n / 2;
410        let med = if n % 2 == 1 {
411            nums[mid]
412        } else {
413            (nums[mid - 1] + nums[mid]) / 2.0
414        };
415        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(med)))
416    }
417}
418
419#[derive(Debug)]
420pub struct StdevSample; // sample
421impl Function for StdevSample {
422    func_caps!(PURE, NUMERIC_ONLY, REDUCTION, STREAM_OK);
423    fn name(&self) -> &'static str {
424        "STDEV.S"
425    }
426    fn aliases(&self) -> &'static [&'static str] {
427        &["STDEV"]
428    }
429    fn min_args(&self) -> usize {
430        1
431    }
432    fn variadic(&self) -> bool {
433        true
434    }
435    fn arg_schema(&self) -> &'static [ArgSchema] {
436        &ARG_RANGE_NUM_LENIENT_ONE[..]
437    }
438    fn eval<'a, 'b, 'c>(
439        &self,
440        args: &'c [ArgumentHandle<'a, 'b>],
441        _ctx: &dyn FunctionContext<'b>,
442    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
443        let nums = collect_numeric_stats(args)?;
444        let n = nums.len();
445        if n < 2 {
446            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
447                ExcelError::from_error_string("#DIV/0!"),
448            )));
449        }
450        let mean = nums.iter().sum::<f64>() / (n as f64);
451        let mut ss = 0.0;
452        for &v in &nums {
453            let d = v - mean;
454            ss += d * d;
455        }
456        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
457            (ss / ((n - 1) as f64)).sqrt(),
458        )))
459    }
460}
461
462#[derive(Debug)]
463pub struct StdevPop; // population
464impl Function for StdevPop {
465    func_caps!(PURE, NUMERIC_ONLY, REDUCTION, STREAM_OK);
466    fn name(&self) -> &'static str {
467        "STDEV.P"
468    }
469    fn aliases(&self) -> &'static [&'static str] {
470        &["STDEVP"]
471    }
472    fn min_args(&self) -> usize {
473        1
474    }
475    fn variadic(&self) -> bool {
476        true
477    }
478    fn arg_schema(&self) -> &'static [ArgSchema] {
479        &ARG_RANGE_NUM_LENIENT_ONE[..]
480    }
481    fn eval<'a, 'b, 'c>(
482        &self,
483        args: &'c [ArgumentHandle<'a, 'b>],
484        _ctx: &dyn FunctionContext<'b>,
485    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
486        let nums = collect_numeric_stats(args)?;
487        let n = nums.len();
488        if n == 0 {
489            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
490                ExcelError::from_error_string("#DIV/0!"),
491            )));
492        }
493        let mean = nums.iter().sum::<f64>() / (n as f64);
494        let mut ss = 0.0;
495        for &v in &nums {
496            let d = v - mean;
497            ss += d * d;
498        }
499        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
500            (ss / (n as f64)).sqrt(),
501        )))
502    }
503}
504
505#[derive(Debug)]
506pub struct VarSample; // sample variance
507impl Function for VarSample {
508    func_caps!(PURE, NUMERIC_ONLY, REDUCTION, STREAM_OK);
509    fn name(&self) -> &'static str {
510        "VAR.S"
511    }
512    fn aliases(&self) -> &'static [&'static str] {
513        &["VAR"]
514    }
515    fn min_args(&self) -> usize {
516        1
517    }
518    fn variadic(&self) -> bool {
519        true
520    }
521    fn arg_schema(&self) -> &'static [ArgSchema] {
522        &ARG_RANGE_NUM_LENIENT_ONE[..]
523    }
524    fn eval<'a, 'b, 'c>(
525        &self,
526        args: &'c [ArgumentHandle<'a, 'b>],
527        _ctx: &dyn FunctionContext<'b>,
528    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
529        let nums = collect_numeric_stats(args)?;
530        let n = nums.len();
531        if n < 2 {
532            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
533                ExcelError::from_error_string("#DIV/0!"),
534            )));
535        }
536        let mean = nums.iter().sum::<f64>() / (n as f64);
537        let mut ss = 0.0;
538        for &v in &nums {
539            let d = v - mean;
540            ss += d * d;
541        }
542        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
543            ss / ((n - 1) as f64),
544        )))
545    }
546}
547
548#[derive(Debug)]
549pub struct VarPop; // population variance
550impl Function for VarPop {
551    func_caps!(PURE, NUMERIC_ONLY, REDUCTION, STREAM_OK);
552    fn name(&self) -> &'static str {
553        "VAR.P"
554    }
555    fn aliases(&self) -> &'static [&'static str] {
556        &["VARP"]
557    }
558    fn min_args(&self) -> usize {
559        1
560    }
561    fn variadic(&self) -> bool {
562        true
563    }
564    fn arg_schema(&self) -> &'static [ArgSchema] {
565        &ARG_RANGE_NUM_LENIENT_ONE[..]
566    }
567    fn eval<'a, 'b, 'c>(
568        &self,
569        args: &'c [ArgumentHandle<'a, 'b>],
570        _ctx: &dyn FunctionContext<'b>,
571    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
572        let nums = collect_numeric_stats(args)?;
573        let n = nums.len();
574        if n == 0 {
575            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
576                ExcelError::from_error_string("#DIV/0!"),
577            )));
578        }
579        let mean = nums.iter().sum::<f64>() / (n as f64);
580        let mut ss = 0.0;
581        for &v in &nums {
582            let d = v - mean;
583            ss += d * d;
584        }
585        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
586            ss / (n as f64),
587        )))
588    }
589}
590
591// MODE.SNGL (alias MODE) and MODE.MULT
592#[derive(Debug)]
593pub struct ModeSingleFn;
594impl Function for ModeSingleFn {
595    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
596    fn name(&self) -> &'static str {
597        "MODE.SNGL"
598    }
599    fn aliases(&self) -> &'static [&'static str] {
600        &["MODE"]
601    }
602    fn min_args(&self) -> usize {
603        1
604    }
605    fn variadic(&self) -> bool {
606        true
607    }
608    fn arg_schema(&self) -> &'static [ArgSchema] {
609        &ARG_RANGE_NUM_LENIENT_ONE[..]
610    }
611    fn eval<'a, 'b, 'c>(
612        &self,
613        args: &'c [ArgumentHandle<'a, 'b>],
614        _ctx: &dyn FunctionContext<'b>,
615    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
616        let mut nums = collect_numeric_stats(args)?;
617        if nums.is_empty() {
618            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
619                ExcelError::new_na(),
620            )));
621        }
622        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
623        let mut best_val = nums[0];
624        let mut best_cnt = 1usize;
625        let mut cur_val = nums[0];
626        let mut cur_cnt = 1usize;
627        for &v in &nums[1..] {
628            if (v - cur_val).abs() < 1e-12 {
629                cur_cnt += 1;
630            } else {
631                if cur_cnt > best_cnt {
632                    best_cnt = cur_cnt;
633                    best_val = cur_val;
634                }
635                cur_val = v;
636                cur_cnt = 1;
637            }
638        }
639        if cur_cnt > best_cnt {
640            best_cnt = cur_cnt;
641            best_val = cur_val;
642        }
643        if best_cnt <= 1 {
644            Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
645                ExcelError::new_na(),
646            )))
647        } else {
648            Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
649                best_val,
650            )))
651        }
652    }
653}
654
655#[derive(Debug)]
656pub struct ModeMultiFn;
657impl Function for ModeMultiFn {
658    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
659    fn name(&self) -> &'static str {
660        "MODE.MULT"
661    }
662    fn min_args(&self) -> usize {
663        1
664    }
665    fn variadic(&self) -> bool {
666        true
667    }
668    fn arg_schema(&self) -> &'static [ArgSchema] {
669        &ARG_RANGE_NUM_LENIENT_ONE[..]
670    }
671    fn eval<'a, 'b, 'c>(
672        &self,
673        args: &'c [ArgumentHandle<'a, 'b>],
674        _ctx: &dyn FunctionContext<'b>,
675    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
676        let mut nums = collect_numeric_stats(args)?;
677        if nums.is_empty() {
678            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
679                ExcelError::new_na(),
680            )));
681        }
682        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
683        let mut runs: Vec<(f64, usize)> = Vec::new();
684        let mut cur_val = nums[0];
685        let mut cur_cnt = 1usize;
686        for &v in &nums[1..] {
687            if (v - cur_val).abs() < 1e-12 {
688                cur_cnt += 1;
689            } else {
690                runs.push((cur_val, cur_cnt));
691                cur_val = v;
692                cur_cnt = 1;
693            }
694        }
695        runs.push((cur_val, cur_cnt));
696        let max_freq = runs.iter().map(|r| r.1).max().unwrap_or(0);
697        if max_freq <= 1 {
698            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
699                ExcelError::new_na(),
700            )));
701        }
702        let rows: Vec<Vec<LiteralValue>> = runs
703            .into_iter()
704            .filter(|&(_, c)| c == max_freq)
705            .map(|(v, _)| vec![LiteralValue::Number(v)])
706            .collect();
707        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(rows)))
708    }
709}
710
711#[derive(Debug)]
712pub struct PercentileInc; // inclusive
713impl Function for PercentileInc {
714    func_caps!(PURE, NUMERIC_ONLY);
715    fn name(&self) -> &'static str {
716        "PERCENTILE.INC"
717    }
718    fn aliases(&self) -> &'static [&'static str] {
719        &["PERCENTILE"]
720    }
721    fn min_args(&self) -> usize {
722        2
723    }
724    fn variadic(&self) -> bool {
725        true
726    }
727    fn arg_schema(&self) -> &'static [ArgSchema] {
728        &ARG_RANGE_NUM_LENIENT_ONE[..]
729    }
730    fn eval<'a, 'b, 'c>(
731        &self,
732        args: &'c [ArgumentHandle<'a, 'b>],
733        _ctx: &dyn FunctionContext<'b>,
734    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
735        if args.len() < 2 {
736            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
737                ExcelError::new_num(),
738            )));
739        }
740        let pv = scalar_like_value(args.last().unwrap())?;
741        let p = match coerce_num(&pv) {
742            Ok(n) => n,
743            Err(_) => {
744                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
745                    ExcelError::new_num(),
746                )));
747            }
748        };
749        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
750        if nums.is_empty() {
751            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
752                ExcelError::new_num(),
753            )));
754        }
755        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
756        match percentile_inc(&nums, p) {
757            Ok(v) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(v))),
758            Err(e) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
759        }
760    }
761}
762
763#[derive(Debug)]
764pub struct PercentileExc; // exclusive
765impl Function for PercentileExc {
766    func_caps!(PURE, NUMERIC_ONLY);
767    fn name(&self) -> &'static str {
768        "PERCENTILE.EXC"
769    }
770    fn min_args(&self) -> usize {
771        2
772    }
773    fn variadic(&self) -> bool {
774        true
775    }
776    fn arg_schema(&self) -> &'static [ArgSchema] {
777        &ARG_RANGE_NUM_LENIENT_ONE[..]
778    }
779    fn eval<'a, 'b, 'c>(
780        &self,
781        args: &'c [ArgumentHandle<'a, 'b>],
782        _ctx: &dyn FunctionContext<'b>,
783    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
784        if args.len() < 2 {
785            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
786                ExcelError::new_num(),
787            )));
788        }
789        let pv = scalar_like_value(args.last().unwrap())?;
790        let p = match coerce_num(&pv) {
791            Ok(n) => n,
792            Err(_) => {
793                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
794                    ExcelError::new_num(),
795                )));
796            }
797        };
798        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
799        if nums.is_empty() {
800            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
801                ExcelError::new_num(),
802            )));
803        }
804        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
805        match percentile_exc(&nums, p) {
806            Ok(v) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(v))),
807            Err(e) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
808        }
809    }
810}
811
812#[derive(Debug)]
813pub struct QuartileInc; // quartile inclusive
814impl Function for QuartileInc {
815    func_caps!(PURE, NUMERIC_ONLY);
816    fn name(&self) -> &'static str {
817        "QUARTILE.INC"
818    }
819    fn aliases(&self) -> &'static [&'static str] {
820        &["QUARTILE"]
821    }
822    fn min_args(&self) -> usize {
823        2
824    }
825    fn variadic(&self) -> bool {
826        true
827    }
828    fn arg_schema(&self) -> &'static [ArgSchema] {
829        &ARG_RANGE_NUM_LENIENT_ONE[..]
830    }
831    fn eval<'a, 'b, 'c>(
832        &self,
833        args: &'c [ArgumentHandle<'a, 'b>],
834        _ctx: &dyn FunctionContext<'b>,
835    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
836        if args.len() < 2 {
837            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
838                ExcelError::new_num(),
839            )));
840        }
841        let qv = scalar_like_value(args.last().unwrap())?;
842        let q = match coerce_num(&qv) {
843            Ok(n) => n,
844            Err(_) => {
845                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
846                    ExcelError::new_num(),
847                )));
848            }
849        };
850        let q_i = q as i64;
851        if !(0..=4).contains(&q_i) {
852            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
853                ExcelError::new_num(),
854            )));
855        }
856        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
857        if nums.is_empty() {
858            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
859                ExcelError::new_num(),
860            )));
861        }
862        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
863        let p = match q_i {
864            0 => {
865                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
866                    nums[0],
867                )));
868            }
869            4 => {
870                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
871                    nums[nums.len() - 1],
872                )));
873            }
874            1 => 0.25,
875            2 => 0.5,
876            3 => 0.75,
877            _ => {
878                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
879                    ExcelError::new_num(),
880                )));
881            }
882        };
883        match percentile_inc(&nums, p) {
884            Ok(v) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(v))),
885            Err(e) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
886        }
887    }
888}
889
890#[derive(Debug)]
891pub struct QuartileExc; // quartile exclusive
892impl Function for QuartileExc {
893    func_caps!(PURE, NUMERIC_ONLY);
894    fn name(&self) -> &'static str {
895        "QUARTILE.EXC"
896    }
897    fn min_args(&self) -> usize {
898        2
899    }
900    fn variadic(&self) -> bool {
901        true
902    }
903    fn arg_schema(&self) -> &'static [ArgSchema] {
904        &ARG_RANGE_NUM_LENIENT_ONE[..]
905    }
906    fn eval<'a, 'b, 'c>(
907        &self,
908        args: &'c [ArgumentHandle<'a, 'b>],
909        _ctx: &dyn FunctionContext<'b>,
910    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
911        if args.len() < 2 {
912            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
913                ExcelError::new_num(),
914            )));
915        }
916        let qv = scalar_like_value(args.last().unwrap())?;
917        let q = match coerce_num(&qv) {
918            Ok(n) => n,
919            Err(_) => {
920                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
921                    ExcelError::new_num(),
922                )));
923            }
924        };
925        let q_i = q as i64;
926        if !(1..=3).contains(&q_i) {
927            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
928                ExcelError::new_num(),
929            )));
930        }
931        let mut nums = collect_numeric_stats(&args[..args.len() - 1])?;
932        if nums.len() < 2 {
933            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
934                ExcelError::new_num(),
935            )));
936        }
937        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
938        let p = match q_i {
939            1 => 0.25,
940            2 => 0.5,
941            3 => 0.75,
942            _ => {
943                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
944                    ExcelError::new_num(),
945                )));
946            }
947        };
948        match percentile_exc(&nums, p) {
949            Ok(v) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(v))),
950            Err(e) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
951        }
952    }
953}
954
955/// PRODUCT(number1, [number2], ...) - Multiplies all arguments
956#[derive(Debug)]
957pub struct ProductFn;
958impl Function for ProductFn {
959    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
960    fn name(&self) -> &'static str {
961        "PRODUCT"
962    }
963    fn min_args(&self) -> usize {
964        1
965    }
966    fn variadic(&self) -> bool {
967        true
968    }
969    fn arg_schema(&self) -> &'static [ArgSchema] {
970        &ARG_RANGE_NUM_LENIENT_ONE[..]
971    }
972    fn eval<'a, 'b, 'c>(
973        &self,
974        args: &'c [ArgumentHandle<'a, 'b>],
975        _ctx: &dyn FunctionContext<'b>,
976    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
977        let nums = collect_numeric_stats(args)?;
978        if nums.is_empty() {
979            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0)));
980        }
981        let result = nums.iter().product::<f64>();
982        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
983            result,
984        )))
985    }
986}
987
988/// GEOMEAN(number1, [number2], ...) - Returns the geometric mean
989#[derive(Debug)]
990pub struct GeomeanFn;
991impl Function for GeomeanFn {
992    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
993    fn name(&self) -> &'static str {
994        "GEOMEAN"
995    }
996    fn min_args(&self) -> usize {
997        1
998    }
999    fn variadic(&self) -> bool {
1000        true
1001    }
1002    fn arg_schema(&self) -> &'static [ArgSchema] {
1003        &ARG_RANGE_NUM_LENIENT_ONE[..]
1004    }
1005    fn eval<'a, 'b, 'c>(
1006        &self,
1007        args: &'c [ArgumentHandle<'a, 'b>],
1008        _ctx: &dyn FunctionContext<'b>,
1009    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1010        let nums = collect_numeric_stats(args)?;
1011        if nums.is_empty() {
1012            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1013                ExcelError::new_num(),
1014            )));
1015        }
1016        // All values must be positive
1017        if nums.iter().any(|&n| n <= 0.0) {
1018            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1019                ExcelError::new_num(),
1020            )));
1021        }
1022        // Geometric mean = (x1 * x2 * ... * xn)^(1/n)
1023        // Use log to avoid overflow: exp(mean(ln(x)))
1024        let log_sum: f64 = nums.iter().map(|x| x.ln()).sum();
1025        let result = (log_sum / nums.len() as f64).exp();
1026        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1027            result,
1028        )))
1029    }
1030}
1031
1032/// HARMEAN(number1, [number2], ...) - Returns the harmonic mean
1033#[derive(Debug)]
1034pub struct HarmeanFn;
1035impl Function for HarmeanFn {
1036    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1037    fn name(&self) -> &'static str {
1038        "HARMEAN"
1039    }
1040    fn min_args(&self) -> usize {
1041        1
1042    }
1043    fn variadic(&self) -> bool {
1044        true
1045    }
1046    fn arg_schema(&self) -> &'static [ArgSchema] {
1047        &ARG_RANGE_NUM_LENIENT_ONE[..]
1048    }
1049    fn eval<'a, 'b, 'c>(
1050        &self,
1051        args: &'c [ArgumentHandle<'a, 'b>],
1052        _ctx: &dyn FunctionContext<'b>,
1053    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1054        let nums = collect_numeric_stats(args)?;
1055        if nums.is_empty() {
1056            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1057                ExcelError::new_num(),
1058            )));
1059        }
1060        // All values must be positive
1061        if nums.iter().any(|&n| n <= 0.0) {
1062            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1063                ExcelError::new_num(),
1064            )));
1065        }
1066        // Harmonic mean = n / sum(1/x)
1067        let sum_reciprocals: f64 = nums.iter().map(|x| 1.0 / x).sum();
1068        let result = nums.len() as f64 / sum_reciprocals;
1069        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1070            result,
1071        )))
1072    }
1073}
1074
1075/// AVEDEV(number1, [number2], ...) - Returns the average of absolute deviations from mean
1076#[derive(Debug)]
1077pub struct AvedevFn;
1078impl Function for AvedevFn {
1079    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1080    fn name(&self) -> &'static str {
1081        "AVEDEV"
1082    }
1083    fn min_args(&self) -> usize {
1084        1
1085    }
1086    fn variadic(&self) -> bool {
1087        true
1088    }
1089    fn arg_schema(&self) -> &'static [ArgSchema] {
1090        &ARG_RANGE_NUM_LENIENT_ONE[..]
1091    }
1092    fn eval<'a, 'b, 'c>(
1093        &self,
1094        args: &'c [ArgumentHandle<'a, 'b>],
1095        _ctx: &dyn FunctionContext<'b>,
1096    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1097        let nums = collect_numeric_stats(args)?;
1098        if nums.is_empty() {
1099            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1100                ExcelError::new_num(),
1101            )));
1102        }
1103        let mean = nums.iter().sum::<f64>() / nums.len() as f64;
1104        let avedev = nums.iter().map(|x| (x - mean).abs()).sum::<f64>() / nums.len() as f64;
1105        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1106            avedev,
1107        )))
1108    }
1109}
1110
1111/// DEVSQ(number1, [number2], ...) - Returns the sum of squared deviations from mean
1112#[derive(Debug)]
1113pub struct DevsqFn;
1114
1115/* ─────────────────────────── MAXIFS / MINIFS ──────────────────────────── */
1116
1117use super::utils::{ARG_ANY_ONE, criteria_match};
1118
1119/// MAXIFS(max_range, criteria_range1, criteria1, [criteria_range2, criteria2], ...)
1120/// Returns the maximum value among cells specified by given conditions.
1121#[derive(Debug)]
1122pub struct MaxIfsFn;
1123impl Function for MaxIfsFn {
1124    func_caps!(PURE, REDUCTION);
1125    fn name(&self) -> &'static str {
1126        "MAXIFS"
1127    }
1128    fn min_args(&self) -> usize {
1129        3
1130    }
1131    fn variadic(&self) -> bool {
1132        true
1133    }
1134    fn arg_schema(&self) -> &'static [ArgSchema] {
1135        &ARG_ANY_ONE[..]
1136    }
1137    fn eval<'a, 'b, 'c>(
1138        &self,
1139        args: &'c [ArgumentHandle<'a, 'b>],
1140        _ctx: &dyn FunctionContext<'b>,
1141    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1142        eval_maxminifs(args, true)
1143    }
1144}
1145
1146/// MINIFS(min_range, criteria_range1, criteria1, [criteria_range2, criteria2], ...)
1147/// Returns the minimum value among cells specified by given conditions.
1148#[derive(Debug)]
1149pub struct MinIfsFn;
1150impl Function for MinIfsFn {
1151    func_caps!(PURE, REDUCTION);
1152    fn name(&self) -> &'static str {
1153        "MINIFS"
1154    }
1155    fn min_args(&self) -> usize {
1156        3
1157    }
1158    fn variadic(&self) -> bool {
1159        true
1160    }
1161    fn arg_schema(&self) -> &'static [ArgSchema] {
1162        &ARG_ANY_ONE[..]
1163    }
1164    fn eval<'a, 'b, 'c>(
1165        &self,
1166        args: &'c [ArgumentHandle<'a, 'b>],
1167        _ctx: &dyn FunctionContext<'b>,
1168    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1169        eval_maxminifs(args, false)
1170    }
1171}
1172
1173/// Shared implementation for MAXIFS and MINIFS
1174fn eval_maxminifs<'a, 'b>(
1175    args: &[ArgumentHandle<'a, 'b>],
1176    is_max: bool,
1177) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1178    // Validate argument count: must be target_range + N pairs
1179    if args.len() < 3 || !(args.len() - 1).is_multiple_of(2) {
1180        return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1181            ExcelError::new_value().with_message(format!(
1182                "Function expects 1 target_range followed by N pairs (criteria_range, criteria); got {} args",
1183                args.len()
1184            )),
1185        )));
1186    }
1187
1188    // Get target range
1189    let target_view = match args[0].range_view() {
1190        Ok(v) => v,
1191        Err(_) => {
1192            // Single value case - if criteria match, return that value
1193            let target_val = args[0].value()?.into_literal();
1194            if let LiteralValue::Error(e) = target_val {
1195                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e)));
1196            }
1197            // Check all criteria against empty/scalar
1198            let mut all_match = true;
1199            for i in (1..args.len()).step_by(2) {
1200                let crit_val = args[i].value()?.into_literal();
1201                let pred = crate::args::parse_criteria(&args[i + 1].value()?.into_literal())?;
1202                if !criteria_match(&pred, &crit_val) {
1203                    all_match = false;
1204                    break;
1205                }
1206            }
1207            if all_match {
1208                return match coerce_num(&target_val) {
1209                    Ok(n) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(n))),
1210                    Err(_) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0))),
1211                };
1212            }
1213            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0)));
1214        }
1215    };
1216
1217    let (rows, cols) = target_view.dims();
1218
1219    // Parse all criteria
1220    let mut criteria_ranges = Vec::new();
1221    let mut predicates = Vec::new();
1222    for i in (1..args.len()).step_by(2) {
1223        let crit_view = args[i].range_view().ok();
1224        let pred = crate::args::parse_criteria(&args[i + 1].value()?.into_literal())?;
1225        criteria_ranges.push(crit_view);
1226        predicates.push(pred);
1227    }
1228
1229    // Iterate through all cells and find max/min where all criteria match
1230    let mut result: Option<f64> = None;
1231
1232    for r in 0..rows {
1233        for c in 0..cols {
1234            // Check all criteria
1235            let mut all_match = true;
1236            for (crit_idx, pred) in predicates.iter().enumerate() {
1237                let crit_val = match &criteria_ranges[crit_idx] {
1238                    Some(view) => {
1239                        let (cr, cc) = view.dims();
1240                        if r < cr && c < cc {
1241                            view.get_cell(r, c)
1242                        } else {
1243                            LiteralValue::Empty
1244                        }
1245                    }
1246                    None => LiteralValue::Empty,
1247                };
1248                if !criteria_match(pred, &crit_val) {
1249                    all_match = false;
1250                    break;
1251                }
1252            }
1253
1254            if all_match {
1255                let target_val = target_view.get_cell(r, c);
1256                match target_val {
1257                    LiteralValue::Error(e) => {
1258                        return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e)));
1259                    }
1260                    LiteralValue::Number(n) => {
1261                        result = Some(match result {
1262                            None => n,
1263                            Some(curr) => {
1264                                if is_max {
1265                                    curr.max(n)
1266                                } else {
1267                                    curr.min(n)
1268                                }
1269                            }
1270                        });
1271                    }
1272                    LiteralValue::Int(i) => {
1273                        let n = i as f64;
1274                        result = Some(match result {
1275                            None => n,
1276                            Some(curr) => {
1277                                if is_max {
1278                                    curr.max(n)
1279                                } else {
1280                                    curr.min(n)
1281                                }
1282                            }
1283                        });
1284                    }
1285                    _ => {} // Skip non-numeric
1286                }
1287            }
1288        }
1289    }
1290
1291    // Excel returns 0 if no matches found
1292    Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1293        result.unwrap_or(0.0),
1294    )))
1295}
1296
1297/* ─────────────────────────── TRIMMEAN ──────────────────────────── */
1298
1299/// TRIMMEAN(array, percent) - Returns the mean of the interior of a data set
1300/// Excludes a percentage of data points from both ends
1301#[derive(Debug)]
1302pub struct TrimmeanFn;
1303impl Function for TrimmeanFn {
1304    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1305    fn name(&self) -> &'static str {
1306        "TRIMMEAN"
1307    }
1308    fn min_args(&self) -> usize {
1309        2
1310    }
1311    fn arg_schema(&self) -> &'static [ArgSchema] {
1312        &ARG_RANGE_NUM_LENIENT_ONE[..]
1313    }
1314    fn eval<'a, 'b, 'c>(
1315        &self,
1316        args: &'c [ArgumentHandle<'a, 'b>],
1317        _ctx: &dyn FunctionContext<'b>,
1318    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1319        let mut nums = collect_numeric_stats(&args[0..1])?;
1320        if nums.is_empty() {
1321            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1322                ExcelError::new_num(),
1323            )));
1324        }
1325
1326        let percent = match args[1].value()?.into_literal() {
1327            LiteralValue::Error(e) => {
1328                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e)));
1329            }
1330            other => coerce_num(&other)?,
1331        };
1332
1333        // Percent must be between 0 and 1 (exclusive of 1)
1334        if !(0.0..1.0).contains(&percent) {
1335            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1336                ExcelError::new_num(),
1337            )));
1338        }
1339
1340        nums.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1341
1342        let n = nums.len();
1343        // Number of values to exclude from each end
1344        let exclude = ((n as f64 * percent) / 2.0).floor() as usize;
1345
1346        if 2 * exclude >= n {
1347            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1348                ExcelError::new_num(),
1349            )));
1350        }
1351
1352        let trimmed = &nums[exclude..n - exclude];
1353        let sum: f64 = trimmed.iter().sum();
1354        let mean = sum / trimmed.len() as f64;
1355
1356        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(mean)))
1357    }
1358}
1359
1360/* ─────────────────────────── CORREL ──────────────────────────── */
1361
1362/// Helper to collect two paired arrays for regression/correlation functions
1363fn collect_paired_arrays(args: &[ArgumentHandle]) -> Result<(Vec<f64>, Vec<f64>), ExcelError> {
1364    let y_nums = collect_numeric_stats(&args[0..1])?;
1365    let x_nums = collect_numeric_stats(&args[1..2])?;
1366
1367    // Arrays must have same length
1368    if y_nums.len() != x_nums.len() {
1369        return Err(ExcelError::new_na());
1370    }
1371
1372    if y_nums.is_empty() {
1373        return Err(ExcelError::new_div());
1374    }
1375
1376    Ok((y_nums, x_nums))
1377}
1378
1379/// CORREL(array1, array2) - Returns the correlation coefficient between two data sets
1380#[derive(Debug)]
1381pub struct CorrelFn;
1382impl Function for CorrelFn {
1383    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1384    fn name(&self) -> &'static str {
1385        "CORREL"
1386    }
1387    fn min_args(&self) -> usize {
1388        2
1389    }
1390    fn arg_schema(&self) -> &'static [ArgSchema] {
1391        &ARG_RANGE_NUM_LENIENT_ONE[..]
1392    }
1393    fn eval<'a, 'b, 'c>(
1394        &self,
1395        args: &'c [ArgumentHandle<'a, 'b>],
1396        _ctx: &dyn FunctionContext<'b>,
1397    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1398        let (y, x) = match collect_paired_arrays(args) {
1399            Ok(v) => v,
1400            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
1401        };
1402
1403        let n = x.len() as f64;
1404        let mean_x = x.iter().sum::<f64>() / n;
1405        let mean_y = y.iter().sum::<f64>() / n;
1406
1407        let mut sum_xy = 0.0;
1408        let mut sum_x2 = 0.0;
1409        let mut sum_y2 = 0.0;
1410
1411        for i in 0..x.len() {
1412            let dx = x[i] - mean_x;
1413            let dy = y[i] - mean_y;
1414            sum_xy += dx * dy;
1415            sum_x2 += dx * dx;
1416            sum_y2 += dy * dy;
1417        }
1418
1419        let denom = (sum_x2 * sum_y2).sqrt();
1420        if denom == 0.0 {
1421            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1422                ExcelError::new_div(),
1423            )));
1424        }
1425
1426        let correl = sum_xy / denom;
1427        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1428            correl,
1429        )))
1430    }
1431}
1432
1433/* ─────────────────────────── SLOPE ──────────────────────────── */
1434
1435/// SLOPE(known_y's, known_x's) - Returns the slope of the linear regression line
1436#[derive(Debug)]
1437pub struct SlopeFn;
1438impl Function for SlopeFn {
1439    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1440    fn name(&self) -> &'static str {
1441        "SLOPE"
1442    }
1443    fn min_args(&self) -> usize {
1444        2
1445    }
1446    fn arg_schema(&self) -> &'static [ArgSchema] {
1447        &ARG_RANGE_NUM_LENIENT_ONE[..]
1448    }
1449    fn eval<'a, 'b, 'c>(
1450        &self,
1451        args: &'c [ArgumentHandle<'a, 'b>],
1452        _ctx: &dyn FunctionContext<'b>,
1453    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1454        let (y, x) = match collect_paired_arrays(args) {
1455            Ok(v) => v,
1456            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
1457        };
1458
1459        let n = x.len() as f64;
1460        let mean_x = x.iter().sum::<f64>() / n;
1461        let mean_y = y.iter().sum::<f64>() / n;
1462
1463        let mut sum_xy = 0.0;
1464        let mut sum_x2 = 0.0;
1465
1466        for i in 0..x.len() {
1467            let dx = x[i] - mean_x;
1468            let dy = y[i] - mean_y;
1469            sum_xy += dx * dy;
1470            sum_x2 += dx * dx;
1471        }
1472
1473        if sum_x2 == 0.0 {
1474            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1475                ExcelError::new_div(),
1476            )));
1477        }
1478
1479        let slope = sum_xy / sum_x2;
1480        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1481            slope,
1482        )))
1483    }
1484}
1485
1486/* ─────────────────────────── INTERCEPT ──────────────────────────── */
1487
1488/// INTERCEPT(known_y's, known_x's) - Returns the y-intercept of the linear regression line
1489#[derive(Debug)]
1490pub struct InterceptFn;
1491impl Function for InterceptFn {
1492    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1493    fn name(&self) -> &'static str {
1494        "INTERCEPT"
1495    }
1496    fn min_args(&self) -> usize {
1497        2
1498    }
1499    fn arg_schema(&self) -> &'static [ArgSchema] {
1500        &ARG_RANGE_NUM_LENIENT_ONE[..]
1501    }
1502    fn eval<'a, 'b, 'c>(
1503        &self,
1504        args: &'c [ArgumentHandle<'a, 'b>],
1505        _ctx: &dyn FunctionContext<'b>,
1506    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1507        let (y, x) = match collect_paired_arrays(args) {
1508            Ok(v) => v,
1509            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
1510        };
1511
1512        let n = x.len() as f64;
1513        let mean_x = x.iter().sum::<f64>() / n;
1514        let mean_y = y.iter().sum::<f64>() / n;
1515
1516        let mut sum_xy = 0.0;
1517        let mut sum_x2 = 0.0;
1518
1519        for i in 0..x.len() {
1520            let dx = x[i] - mean_x;
1521            let dy = y[i] - mean_y;
1522            sum_xy += dx * dy;
1523            sum_x2 += dx * dx;
1524        }
1525
1526        if sum_x2 == 0.0 {
1527            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1528                ExcelError::new_div(),
1529            )));
1530        }
1531
1532        let slope = sum_xy / sum_x2;
1533        let intercept = mean_y - slope * mean_x;
1534        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1535            intercept,
1536        )))
1537    }
1538}
1539
1540impl Function for DevsqFn {
1541    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
1542    fn name(&self) -> &'static str {
1543        "DEVSQ"
1544    }
1545    fn min_args(&self) -> usize {
1546        1
1547    }
1548    fn variadic(&self) -> bool {
1549        true
1550    }
1551    fn arg_schema(&self) -> &'static [ArgSchema] {
1552        &ARG_RANGE_NUM_LENIENT_ONE[..]
1553    }
1554    fn eval<'a, 'b, 'c>(
1555        &self,
1556        args: &'c [ArgumentHandle<'a, 'b>],
1557        _ctx: &dyn FunctionContext<'b>,
1558    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1559        let nums = collect_numeric_stats(args)?;
1560        if nums.is_empty() {
1561            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1562                ExcelError::new_num(),
1563            )));
1564        }
1565        let mean = nums.iter().sum::<f64>() / nums.len() as f64;
1566        let devsq = nums.iter().map(|x| (x - mean).powi(2)).sum::<f64>();
1567        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1568            devsq,
1569        )))
1570    }
1571}
1572
1573/* ═══════════════════════════════════════════════════════════════════════════
1574STATISTICAL DISTRIBUTION FUNCTIONS
1575═══════════════════════════════════════════════════════════════════════════ */
1576
1577/// Helper: Standard normal CDF using error function approximation
1578fn std_norm_cdf(z: f64) -> f64 {
1579    // Use the complementary error function: Φ(z) = 0.5 * erfc(-z / sqrt(2))
1580    // Approximation using Abramowitz and Stegun formula 7.1.26
1581    let a1 = 0.254829592;
1582    let a2 = -0.284496736;
1583    let a3 = 1.421413741;
1584    let a4 = -1.453152027;
1585    let a5 = 1.061405429;
1586    let p = 0.3275911;
1587
1588    let sign = if z < 0.0 { -1.0 } else { 1.0 };
1589    let z_abs = z.abs() / std::f64::consts::SQRT_2;
1590
1591    let t = 1.0 / (1.0 + p * z_abs);
1592    let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-z_abs * z_abs).exp();
1593
1594    0.5 * (1.0 + sign * y)
1595}
1596
1597/// Helper: Standard normal PDF
1598fn std_norm_pdf(z: f64) -> f64 {
1599    let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
1600    inv_sqrt_2pi * (-0.5 * z * z).exp()
1601}
1602
1603/// Helper: Inverse standard normal CDF (probit function)
1604/// Uses Rational approximation from Abramowitz and Stegun
1605#[allow(clippy::excessive_precision)]
1606fn std_norm_inv(p: f64) -> Option<f64> {
1607    if p <= 0.0 || p >= 1.0 {
1608        return None;
1609    }
1610
1611    // Coefficients for rational approximation
1612    const A: [f64; 6] = [
1613        -3.969683028665376e+01,
1614        2.209460984245205e+02,
1615        -2.759285104469687e+02,
1616        1.383577518672690e+02,
1617        -3.066479806614716e+01,
1618        2.506628277459239e+00,
1619    ];
1620    const B: [f64; 5] = [
1621        -5.447609879822406e+01,
1622        1.615858368580409e+02,
1623        -1.556989798598866e+02,
1624        6.680131188771972e+01,
1625        -1.328068155288572e+01,
1626    ];
1627    const C: [f64; 6] = [
1628        -7.784894002430293e-03,
1629        -3.223964580411365e-01,
1630        -2.400758277161838e+00,
1631        -2.549732539343734e+00,
1632        4.374664141464968e+00,
1633        2.938163982698783e+00,
1634    ];
1635    const D: [f64; 4] = [
1636        7.784695709041462e-03,
1637        3.224671290700398e-01,
1638        2.445134137142996e+00,
1639        3.754408661907416e+00,
1640    ];
1641
1642    const P_LOW: f64 = 0.02425;
1643    const P_HIGH: f64 = 1.0 - P_LOW;
1644
1645    let q = p - 0.5;
1646
1647    if p < P_LOW {
1648        // Lower tail
1649        let r = (-2.0 * p.ln()).sqrt();
1650        let num = ((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5];
1651        let den = (((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0;
1652        Some(num / den)
1653    } else if p <= P_HIGH {
1654        // Central region
1655        let r = q * q;
1656        let num = ((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5];
1657        let den = ((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0;
1658        Some(q * num / den)
1659    } else {
1660        // Upper tail
1661        let r = (-2.0 * (1.0 - p).ln()).sqrt();
1662        let num = ((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5];
1663        let den = (((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0;
1664        Some(-num / den)
1665    }
1666}
1667
1668/// NORM.S.DIST(z, cumulative) - Standard normal distribution
1669#[derive(Debug)]
1670pub struct NormSDistFn;
1671impl Function for NormSDistFn {
1672    func_caps!(PURE);
1673    fn name(&self) -> &'static str {
1674        "NORM.S.DIST"
1675    }
1676    fn min_args(&self) -> usize {
1677        2
1678    }
1679    fn arg_schema(&self) -> &'static [ArgSchema] {
1680        use std::sync::LazyLock;
1681        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
1682            vec![
1683                ArgSchema::number_lenient_scalar(),
1684                ArgSchema::number_lenient_scalar(),
1685            ]
1686        });
1687        &SCHEMA[..]
1688    }
1689    fn eval<'a, 'b, 'c>(
1690        &self,
1691        args: &'c [ArgumentHandle<'a, 'b>],
1692        _ctx: &dyn FunctionContext<'b>,
1693    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1694        let z = coerce_num(&scalar_like_value(&args[0])?)?;
1695        let cumulative = coerce_num(&scalar_like_value(&args[1])?)? != 0.0;
1696
1697        let result = if cumulative {
1698            std_norm_cdf(z)
1699        } else {
1700            std_norm_pdf(z)
1701        };
1702        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1703            result,
1704        )))
1705    }
1706}
1707
1708/// NORM.S.INV(probability) - Inverse standard normal distribution
1709#[derive(Debug)]
1710pub struct NormSInvFn;
1711impl Function for NormSInvFn {
1712    func_caps!(PURE);
1713    fn name(&self) -> &'static str {
1714        "NORM.S.INV"
1715    }
1716    fn min_args(&self) -> usize {
1717        1
1718    }
1719    fn arg_schema(&self) -> &'static [ArgSchema] {
1720        use std::sync::LazyLock;
1721        static SCHEMA: LazyLock<Vec<ArgSchema>> =
1722            LazyLock::new(|| vec![ArgSchema::number_lenient_scalar()]);
1723        &SCHEMA[..]
1724    }
1725    fn eval<'a, 'b, 'c>(
1726        &self,
1727        args: &'c [ArgumentHandle<'a, 'b>],
1728        _ctx: &dyn FunctionContext<'b>,
1729    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1730        let p = coerce_num(&scalar_like_value(&args[0])?)?;
1731
1732        match std_norm_inv(p) {
1733            Some(z) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(z))),
1734            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1735                ExcelError::new_num(),
1736            ))),
1737        }
1738    }
1739}
1740
1741/// NORM.DIST(x, mean, standard_dev, cumulative) - Normal distribution
1742#[derive(Debug)]
1743pub struct NormDistFn;
1744impl Function for NormDistFn {
1745    func_caps!(PURE);
1746    fn name(&self) -> &'static str {
1747        "NORM.DIST"
1748    }
1749    fn min_args(&self) -> usize {
1750        4
1751    }
1752    fn arg_schema(&self) -> &'static [ArgSchema] {
1753        use std::sync::LazyLock;
1754        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
1755            vec![
1756                ArgSchema::number_lenient_scalar(),
1757                ArgSchema::number_lenient_scalar(),
1758                ArgSchema::number_lenient_scalar(),
1759                ArgSchema::number_lenient_scalar(),
1760            ]
1761        });
1762        &SCHEMA[..]
1763    }
1764    fn eval<'a, 'b, 'c>(
1765        &self,
1766        args: &'c [ArgumentHandle<'a, 'b>],
1767        _ctx: &dyn FunctionContext<'b>,
1768    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1769        let x = coerce_num(&scalar_like_value(&args[0])?)?;
1770        let mean = coerce_num(&scalar_like_value(&args[1])?)?;
1771        let std_dev = coerce_num(&scalar_like_value(&args[2])?)?;
1772        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
1773
1774        if std_dev <= 0.0 {
1775            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1776                ExcelError::new_num(),
1777            )));
1778        }
1779
1780        let z = (x - mean) / std_dev;
1781
1782        let result = if cumulative {
1783            std_norm_cdf(z)
1784        } else {
1785            std_norm_pdf(z) / std_dev
1786        };
1787        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1788            result,
1789        )))
1790    }
1791}
1792
1793/// NORM.INV(probability, mean, standard_dev) - Inverse normal distribution
1794#[derive(Debug)]
1795pub struct NormInvFn;
1796impl Function for NormInvFn {
1797    func_caps!(PURE);
1798    fn name(&self) -> &'static str {
1799        "NORM.INV"
1800    }
1801    fn min_args(&self) -> usize {
1802        3
1803    }
1804    fn arg_schema(&self) -> &'static [ArgSchema] {
1805        use std::sync::LazyLock;
1806        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
1807            vec![
1808                ArgSchema::number_lenient_scalar(),
1809                ArgSchema::number_lenient_scalar(),
1810                ArgSchema::number_lenient_scalar(),
1811            ]
1812        });
1813        &SCHEMA[..]
1814    }
1815    fn eval<'a, 'b, 'c>(
1816        &self,
1817        args: &'c [ArgumentHandle<'a, 'b>],
1818        _ctx: &dyn FunctionContext<'b>,
1819    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1820        let p = coerce_num(&scalar_like_value(&args[0])?)?;
1821        let mean = coerce_num(&scalar_like_value(&args[1])?)?;
1822        let std_dev = coerce_num(&scalar_like_value(&args[2])?)?;
1823
1824        if std_dev <= 0.0 {
1825            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1826                ExcelError::new_num(),
1827            )));
1828        }
1829
1830        match std_norm_inv(p) {
1831            Some(z) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1832                mean + z * std_dev,
1833            ))),
1834            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1835                ExcelError::new_num(),
1836            ))),
1837        }
1838    }
1839}
1840
1841/// LOGNORM.DIST(x, mean, standard_dev, cumulative) - Log-normal distribution
1842#[derive(Debug)]
1843pub struct LognormDistFn;
1844impl Function for LognormDistFn {
1845    func_caps!(PURE);
1846    fn name(&self) -> &'static str {
1847        "LOGNORM.DIST"
1848    }
1849    fn min_args(&self) -> usize {
1850        4
1851    }
1852    fn arg_schema(&self) -> &'static [ArgSchema] {
1853        use std::sync::LazyLock;
1854        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
1855            vec![
1856                ArgSchema::number_lenient_scalar(),
1857                ArgSchema::number_lenient_scalar(),
1858                ArgSchema::number_lenient_scalar(),
1859                ArgSchema::number_lenient_scalar(),
1860            ]
1861        });
1862        &SCHEMA[..]
1863    }
1864    fn eval<'a, 'b, 'c>(
1865        &self,
1866        args: &'c [ArgumentHandle<'a, 'b>],
1867        _ctx: &dyn FunctionContext<'b>,
1868    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1869        let x = coerce_num(&scalar_like_value(&args[0])?)?;
1870        let mean = coerce_num(&scalar_like_value(&args[1])?)?;
1871        let std_dev = coerce_num(&scalar_like_value(&args[2])?)?;
1872        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
1873
1874        if x <= 0.0 || std_dev <= 0.0 {
1875            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1876                ExcelError::new_num(),
1877            )));
1878        }
1879
1880        let z = (x.ln() - mean) / std_dev;
1881
1882        let result = if cumulative {
1883            std_norm_cdf(z)
1884        } else {
1885            std_norm_pdf(z) / (x * std_dev)
1886        };
1887        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1888            result,
1889        )))
1890    }
1891}
1892
1893/// LOGNORM.INV(probability, mean, standard_dev) - Inverse log-normal distribution
1894#[derive(Debug)]
1895pub struct LognormInvFn;
1896impl Function for LognormInvFn {
1897    func_caps!(PURE);
1898    fn name(&self) -> &'static str {
1899        "LOGNORM.INV"
1900    }
1901    fn min_args(&self) -> usize {
1902        3
1903    }
1904    fn arg_schema(&self) -> &'static [ArgSchema] {
1905        use std::sync::LazyLock;
1906        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
1907            vec![
1908                ArgSchema::number_lenient_scalar(),
1909                ArgSchema::number_lenient_scalar(),
1910                ArgSchema::number_lenient_scalar(),
1911            ]
1912        });
1913        &SCHEMA[..]
1914    }
1915    fn eval<'a, 'b, 'c>(
1916        &self,
1917        args: &'c [ArgumentHandle<'a, 'b>],
1918        _ctx: &dyn FunctionContext<'b>,
1919    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1920        let p = coerce_num(&scalar_like_value(&args[0])?)?;
1921        let mean = coerce_num(&scalar_like_value(&args[1])?)?;
1922        let std_dev = coerce_num(&scalar_like_value(&args[2])?)?;
1923
1924        if std_dev <= 0.0 {
1925            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1926                ExcelError::new_num(),
1927            )));
1928        }
1929
1930        match std_norm_inv(p) {
1931            Some(z) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1932                (mean + z * std_dev).exp(),
1933            ))),
1934            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
1935                ExcelError::new_num(),
1936            ))),
1937        }
1938    }
1939}
1940
1941/// PHI(x) - Standard normal distribution density function (alias for NORM.S.DIST PDF)
1942#[derive(Debug)]
1943pub struct PhiFn;
1944impl Function for PhiFn {
1945    func_caps!(PURE);
1946    fn name(&self) -> &'static str {
1947        "PHI"
1948    }
1949    fn min_args(&self) -> usize {
1950        1
1951    }
1952    fn arg_schema(&self) -> &'static [ArgSchema] {
1953        use std::sync::LazyLock;
1954        static SCHEMA: LazyLock<Vec<ArgSchema>> =
1955            LazyLock::new(|| vec![ArgSchema::number_lenient_scalar()]);
1956        &SCHEMA[..]
1957    }
1958    fn eval<'a, 'b, 'c>(
1959        &self,
1960        args: &'c [ArgumentHandle<'a, 'b>],
1961        _ctx: &dyn FunctionContext<'b>,
1962    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1963        let z = coerce_num(&scalar_like_value(&args[0])?)?;
1964        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1965            std_norm_pdf(z),
1966        )))
1967    }
1968}
1969
1970/// GAUSS(z) - Returns the probability that a member of a standard normal population will fall between the mean and z standard deviations from the mean
1971#[derive(Debug)]
1972pub struct GaussFn;
1973impl Function for GaussFn {
1974    func_caps!(PURE);
1975    fn name(&self) -> &'static str {
1976        "GAUSS"
1977    }
1978    fn min_args(&self) -> usize {
1979        1
1980    }
1981    fn arg_schema(&self) -> &'static [ArgSchema] {
1982        use std::sync::LazyLock;
1983        static SCHEMA: LazyLock<Vec<ArgSchema>> =
1984            LazyLock::new(|| vec![ArgSchema::number_lenient_scalar()]);
1985        &SCHEMA[..]
1986    }
1987    fn eval<'a, 'b, 'c>(
1988        &self,
1989        args: &'c [ArgumentHandle<'a, 'b>],
1990        _ctx: &dyn FunctionContext<'b>,
1991    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
1992        let z = coerce_num(&scalar_like_value(&args[0])?)?;
1993        // GAUSS(z) = Φ(z) - 0.5
1994        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
1995            std_norm_cdf(z) - 0.5,
1996        )))
1997    }
1998}
1999
2000/// Helper: Log-gamma function
2001#[allow(clippy::excessive_precision)]
2002fn ln_gamma(x: f64) -> f64 {
2003    // Lanczos approximation
2004    const G: f64 = 7.0;
2005    const C: [f64; 9] = [
2006        0.99999999999980993,
2007        676.5203681218851,
2008        -1259.1392167224028,
2009        771.32342877765313,
2010        -176.61502916214059,
2011        12.507343278686905,
2012        -0.13857109526572012,
2013        9.9843695780195716e-6,
2014        1.5056327351493116e-7,
2015    ];
2016
2017    if x < 0.5 {
2018        // Reflection formula
2019        let pi = std::f64::consts::PI;
2020        pi.ln() - (pi * x).sin().ln() - ln_gamma(1.0 - x)
2021    } else {
2022        let x = x - 1.0;
2023        let mut ag = C[0];
2024        for (i, c) in C.iter().enumerate().skip(1) {
2025            ag += c / (x + i as f64);
2026        }
2027        let tmp = x + G + 0.5;
2028        0.5 * (2.0 * std::f64::consts::PI).ln() + (tmp).ln() * (x + 0.5) - tmp + ag.ln()
2029    }
2030}
2031
2032/// Helper: Regularized lower incomplete gamma function P(a, x)
2033fn gamma_p(a: f64, x: f64) -> f64 {
2034    if x < 0.0 || a <= 0.0 {
2035        return 0.0;
2036    }
2037    if x == 0.0 {
2038        return 0.0;
2039    }
2040
2041    // Use series expansion for x < a+1
2042    if x < a + 1.0 {
2043        gamma_series(a, x)
2044    } else {
2045        // Use continued fraction for x >= a+1
2046        1.0 - gamma_cf(a, x)
2047    }
2048}
2049
2050/// Helper: Series expansion for incomplete gamma
2051fn gamma_series(a: f64, x: f64) -> f64 {
2052    let ln_ga = ln_gamma(a);
2053    let mut sum = 1.0 / a;
2054    let mut term = sum;
2055    for n in 1..200 {
2056        term *= x / (a + n as f64);
2057        sum += term;
2058        if term.abs() < sum.abs() * 1e-15 {
2059            break;
2060        }
2061    }
2062    sum * (-x + a * x.ln() - ln_ga).exp()
2063}
2064
2065/// Helper: Continued fraction for upper incomplete gamma Q(a,x)
2066/// Using modified Lentz's algorithm (Numerical Recipes formulation)
2067fn gamma_cf(a: f64, x: f64) -> f64 {
2068    let ln_ga = ln_gamma(a);
2069    const TINY: f64 = 1e-30;
2070    const EPS: f64 = 1e-14;
2071
2072    // Set up for evaluating continued fraction by modified Lentz's method
2073    let mut b = x + 1.0 - a;
2074    let mut c = 1.0 / TINY;
2075    let mut d = 1.0 / b;
2076    let mut h = d;
2077
2078    for i in 1..=200 {
2079        let an = -(i as f64) * (i as f64 - a);
2080        b += 2.0;
2081        d = an * d + b;
2082        if d.abs() < TINY {
2083            d = TINY;
2084        }
2085        c = b + an / c;
2086        if c.abs() < TINY {
2087            c = TINY;
2088        }
2089        d = 1.0 / d;
2090        let delta = d * c;
2091        h *= delta;
2092        if (delta - 1.0).abs() <= EPS {
2093            break;
2094        }
2095    }
2096
2097    h * (-x + a * x.ln() - ln_ga).exp()
2098}
2099
2100/// Helper: Regularized incomplete beta function I_x(a,b)
2101/// Uses the continued fraction representation (NIST DLMF 8.17.22)
2102fn beta_i(x: f64, a: f64, b: f64) -> f64 {
2103    if x <= 0.0 {
2104        return 0.0;
2105    }
2106    if x >= 1.0 {
2107        return 1.0;
2108    }
2109    if a <= 0.0 || b <= 0.0 {
2110        return f64::NAN;
2111    }
2112
2113    // Use symmetry for better convergence: I_x(a,b) = 1 - I_{1-x}(b,a)
2114    if x > (a + 1.0) / (a + b + 2.0) {
2115        return 1.0 - beta_i(1.0 - x, b, a);
2116    }
2117
2118    // Compute the prefactor: x^a * (1-x)^b / (a * B(a,b))
2119    let ln_beta = ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b);
2120    let ln_prefactor = a * x.ln() + b * (1.0 - x).ln() - ln_beta - a.ln();
2121    let prefactor = ln_prefactor.exp();
2122
2123    // Evaluate the continued fraction using modified Lentz algorithm
2124    // The CF is: 1 / (1 + d1/(1 + d2/(1 + ...)))
2125    // where d_{2m+1} = -(a+m)(a+b+m)x / ((a+2m)(a+2m+1))
2126    //       d_{2m}   = m(b-m)x / ((a+2m-1)(a+2m))
2127    const EPS: f64 = 1e-14;
2128    const TINY: f64 = 1e-30;
2129
2130    let qab = a + b;
2131    let qap = a + 1.0;
2132    let qam = a - 1.0;
2133    let mut c = 1.0;
2134    let mut d = 1.0 - qab * x / qap;
2135    if d.abs() < TINY {
2136        d = TINY;
2137    }
2138    d = 1.0 / d;
2139    let mut h = d;
2140
2141    for m in 1..=200 {
2142        let m_f64 = m as f64;
2143        let m2 = 2.0 * m_f64;
2144
2145        // Even step: d_{2m} = m(b-m)x / ((a+2m-1)(a+2m))
2146        let aa = m_f64 * (b - m_f64) * x / ((qam + m2) * (a + m2));
2147        d = 1.0 + aa * d;
2148        if d.abs() < TINY {
2149            d = TINY;
2150        }
2151        c = 1.0 + aa / c;
2152        if c.abs() < TINY {
2153            c = TINY;
2154        }
2155        d = 1.0 / d;
2156        h *= d * c;
2157
2158        // Odd step: d_{2m+1} = -(a+m)(a+b+m)x / ((a+2m)(a+2m+1))
2159        let aa = -((a + m_f64) * (qab + m_f64) * x) / ((a + m2) * (qap + m2));
2160        d = 1.0 + aa * d;
2161        if d.abs() < TINY {
2162            d = TINY;
2163        }
2164        c = 1.0 + aa / c;
2165        if c.abs() < TINY {
2166            c = TINY;
2167        }
2168        d = 1.0 / d;
2169        let delta = d * c;
2170        h *= delta;
2171
2172        if (delta - 1.0).abs() <= EPS {
2173            break;
2174        }
2175    }
2176
2177    prefactor * h
2178}
2179
2180/// Helper: T distribution CDF
2181fn t_cdf(t: f64, df: f64) -> f64 {
2182    let x = df / (df + t * t);
2183    0.5 * (1.0 + t.signum() * (1.0 - beta_i(x, df / 2.0, 0.5)))
2184}
2185
2186/// Helper: T distribution inverse CDF using Newton-Raphson
2187fn t_inv(p: f64, df: f64) -> Option<f64> {
2188    if p <= 0.0 || p >= 1.0 {
2189        return None;
2190    }
2191
2192    // Initial guess using normal approximation
2193    let mut t = std_norm_inv(p)?;
2194
2195    // Newton-Raphson iteration
2196    for _ in 0..50 {
2197        let cdf = t_cdf(t, df);
2198        let pdf = t_pdf(t, df);
2199        if pdf.abs() < 1e-30 {
2200            break;
2201        }
2202        let delta = (cdf - p) / pdf;
2203        t -= delta;
2204        if delta.abs() < 1e-12 {
2205            break;
2206        }
2207    }
2208
2209    Some(t)
2210}
2211
2212/// Helper: T distribution PDF
2213fn t_pdf(t: f64, df: f64) -> f64 {
2214    let coef =
2215        (ln_gamma((df + 1.0) / 2.0) - ln_gamma(df / 2.0) - 0.5 * (df * std::f64::consts::PI).ln())
2216            .exp();
2217    coef * (1.0 + t * t / df).powf(-(df + 1.0) / 2.0)
2218}
2219
2220/// Helper: Chi-square CDF
2221fn chisq_cdf(x: f64, df: f64) -> f64 {
2222    if x <= 0.0 {
2223        return 0.0;
2224    }
2225    gamma_p(df / 2.0, x / 2.0)
2226}
2227
2228/// Helper: Chi-square inverse CDF using Newton-Raphson
2229fn chisq_inv(p: f64, df: f64) -> Option<f64> {
2230    if p <= 0.0 || p >= 1.0 {
2231        return None;
2232    }
2233
2234    // Initial guess
2235    let mut x = df.max(1.0);
2236    if p < 0.5 {
2237        x = x.min(1.0);
2238    }
2239
2240    // Newton-Raphson iteration
2241    for _ in 0..100 {
2242        let cdf = chisq_cdf(x, df);
2243        let pdf = chisq_pdf(x, df);
2244        if pdf.abs() < 1e-30 {
2245            break;
2246        }
2247        let delta = (cdf - p) / pdf;
2248        let new_x = (x - delta).max(1e-15);
2249        if (new_x - x).abs() < 1e-12 * x {
2250            x = new_x;
2251            break;
2252        }
2253        x = new_x;
2254    }
2255
2256    Some(x)
2257}
2258
2259/// Helper: Chi-square PDF
2260fn chisq_pdf(x: f64, df: f64) -> f64 {
2261    if x <= 0.0 {
2262        return 0.0;
2263    }
2264    let k = df / 2.0;
2265    ((k - 1.0) * x.ln() - x / 2.0 - k * 2.0_f64.ln() - ln_gamma(k)).exp()
2266}
2267
2268/// Helper: F distribution CDF
2269fn f_cdf(f: f64, d1: f64, d2: f64) -> f64 {
2270    if f <= 0.0 {
2271        return 0.0;
2272    }
2273    let x = d1 * f / (d1 * f + d2);
2274    beta_i(x, d1 / 2.0, d2 / 2.0)
2275}
2276
2277/// Helper: F distribution inverse CDF using Newton-Raphson
2278fn f_inv(p: f64, d1: f64, d2: f64) -> Option<f64> {
2279    if p <= 0.0 || p >= 1.0 {
2280        return None;
2281    }
2282
2283    // Initial guess
2284    let mut f = 1.0;
2285
2286    // Newton-Raphson iteration
2287    for _ in 0..100 {
2288        let cdf = f_cdf(f, d1, d2);
2289        let pdf = f_pdf(f, d1, d2);
2290        if pdf.abs() < 1e-30 {
2291            break;
2292        }
2293        let delta = (cdf - p) / pdf;
2294        let new_f = (f - delta).max(1e-15);
2295        if (new_f - f).abs() < 1e-12 * f {
2296            f = new_f;
2297            break;
2298        }
2299        f = new_f;
2300    }
2301
2302    Some(f)
2303}
2304
2305/// Helper: F distribution PDF
2306fn f_pdf(f: f64, d1: f64, d2: f64) -> f64 {
2307    if f <= 0.0 {
2308        return 0.0;
2309    }
2310    let ln_beta = ln_gamma(d1 / 2.0) + ln_gamma(d2 / 2.0) - ln_gamma((d1 + d2) / 2.0);
2311    let coef = (d1 / 2.0) * (d1 / d2).ln() + (d1 / 2.0 - 1.0) * f.ln()
2312        - ((d1 + d2) / 2.0) * (1.0 + d1 * f / d2).ln()
2313        - ln_beta;
2314    coef.exp()
2315}
2316
2317/// T.DIST(x, deg_freedom, cumulative) - Student's t-distribution
2318#[derive(Debug)]
2319pub struct TDistFn;
2320impl Function for TDistFn {
2321    func_caps!(PURE);
2322    fn name(&self) -> &'static str {
2323        "T.DIST"
2324    }
2325    fn min_args(&self) -> usize {
2326        3
2327    }
2328    fn arg_schema(&self) -> &'static [ArgSchema] {
2329        use std::sync::LazyLock;
2330        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2331            vec![
2332                ArgSchema::number_lenient_scalar(),
2333                ArgSchema::number_lenient_scalar(),
2334                ArgSchema::number_lenient_scalar(),
2335            ]
2336        });
2337        &SCHEMA[..]
2338    }
2339    fn eval<'a, 'b, 'c>(
2340        &self,
2341        args: &'c [ArgumentHandle<'a, 'b>],
2342        _ctx: &dyn FunctionContext<'b>,
2343    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2344        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2345        let df = coerce_num(&scalar_like_value(&args[1])?)?;
2346        let cumulative = coerce_num(&scalar_like_value(&args[2])?)? != 0.0;
2347
2348        if df < 1.0 {
2349            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2350                ExcelError::new_num(),
2351            )));
2352        }
2353
2354        let result = if cumulative {
2355            t_cdf(x, df)
2356        } else {
2357            t_pdf(x, df)
2358        };
2359        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2360            result,
2361        )))
2362    }
2363}
2364
2365/// T.INV(probability, deg_freedom) - Inverse of Student's t-distribution
2366#[derive(Debug)]
2367pub struct TInvFn;
2368impl Function for TInvFn {
2369    func_caps!(PURE);
2370    fn name(&self) -> &'static str {
2371        "T.INV"
2372    }
2373    fn min_args(&self) -> usize {
2374        2
2375    }
2376    fn arg_schema(&self) -> &'static [ArgSchema] {
2377        use std::sync::LazyLock;
2378        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2379            vec![
2380                ArgSchema::number_lenient_scalar(),
2381                ArgSchema::number_lenient_scalar(),
2382            ]
2383        });
2384        &SCHEMA[..]
2385    }
2386    fn eval<'a, 'b, 'c>(
2387        &self,
2388        args: &'c [ArgumentHandle<'a, 'b>],
2389        _ctx: &dyn FunctionContext<'b>,
2390    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2391        let p = coerce_num(&scalar_like_value(&args[0])?)?;
2392        let df = coerce_num(&scalar_like_value(&args[1])?)?;
2393
2394        if df < 1.0 {
2395            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2396                ExcelError::new_num(),
2397            )));
2398        }
2399
2400        match t_inv(p, df) {
2401            Some(result) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2402                result,
2403            ))),
2404            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2405                ExcelError::new_num(),
2406            ))),
2407        }
2408    }
2409}
2410
2411/// CHISQ.DIST(x, deg_freedom, cumulative) - Chi-squared distribution
2412#[derive(Debug)]
2413pub struct ChisqDistFn;
2414impl Function for ChisqDistFn {
2415    func_caps!(PURE);
2416    fn name(&self) -> &'static str {
2417        "CHISQ.DIST"
2418    }
2419    fn min_args(&self) -> usize {
2420        3
2421    }
2422    fn arg_schema(&self) -> &'static [ArgSchema] {
2423        use std::sync::LazyLock;
2424        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2425            vec![
2426                ArgSchema::number_lenient_scalar(),
2427                ArgSchema::number_lenient_scalar(),
2428                ArgSchema::number_lenient_scalar(),
2429            ]
2430        });
2431        &SCHEMA[..]
2432    }
2433    fn eval<'a, 'b, 'c>(
2434        &self,
2435        args: &'c [ArgumentHandle<'a, 'b>],
2436        _ctx: &dyn FunctionContext<'b>,
2437    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2438        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2439        let df = coerce_num(&scalar_like_value(&args[1])?)?;
2440        let cumulative = coerce_num(&scalar_like_value(&args[2])?)? != 0.0;
2441
2442        if df < 1.0 || x < 0.0 {
2443            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2444                ExcelError::new_num(),
2445            )));
2446        }
2447
2448        let result = if cumulative {
2449            chisq_cdf(x, df)
2450        } else {
2451            chisq_pdf(x, df)
2452        };
2453        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2454            result,
2455        )))
2456    }
2457}
2458
2459/// CHISQ.INV(probability, deg_freedom) - Inverse of chi-squared distribution
2460#[derive(Debug)]
2461pub struct ChisqInvFn;
2462impl Function for ChisqInvFn {
2463    func_caps!(PURE);
2464    fn name(&self) -> &'static str {
2465        "CHISQ.INV"
2466    }
2467    fn min_args(&self) -> usize {
2468        2
2469    }
2470    fn arg_schema(&self) -> &'static [ArgSchema] {
2471        use std::sync::LazyLock;
2472        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2473            vec![
2474                ArgSchema::number_lenient_scalar(),
2475                ArgSchema::number_lenient_scalar(),
2476            ]
2477        });
2478        &SCHEMA[..]
2479    }
2480    fn eval<'a, 'b, 'c>(
2481        &self,
2482        args: &'c [ArgumentHandle<'a, 'b>],
2483        _ctx: &dyn FunctionContext<'b>,
2484    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2485        let p = coerce_num(&scalar_like_value(&args[0])?)?;
2486        let df = coerce_num(&scalar_like_value(&args[1])?)?;
2487
2488        if df < 1.0 {
2489            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2490                ExcelError::new_num(),
2491            )));
2492        }
2493
2494        match chisq_inv(p, df) {
2495            Some(result) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2496                result,
2497            ))),
2498            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2499                ExcelError::new_num(),
2500            ))),
2501        }
2502    }
2503}
2504
2505/// F.DIST(x, deg_freedom1, deg_freedom2, cumulative) - F distribution
2506#[derive(Debug)]
2507pub struct FDistFn;
2508impl Function for FDistFn {
2509    func_caps!(PURE);
2510    fn name(&self) -> &'static str {
2511        "F.DIST"
2512    }
2513    fn min_args(&self) -> usize {
2514        4
2515    }
2516    fn arg_schema(&self) -> &'static [ArgSchema] {
2517        use std::sync::LazyLock;
2518        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2519            vec![
2520                ArgSchema::number_lenient_scalar(),
2521                ArgSchema::number_lenient_scalar(),
2522                ArgSchema::number_lenient_scalar(),
2523                ArgSchema::number_lenient_scalar(),
2524            ]
2525        });
2526        &SCHEMA[..]
2527    }
2528    fn eval<'a, 'b, 'c>(
2529        &self,
2530        args: &'c [ArgumentHandle<'a, 'b>],
2531        _ctx: &dyn FunctionContext<'b>,
2532    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2533        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2534        let d1 = coerce_num(&scalar_like_value(&args[1])?)?;
2535        let d2 = coerce_num(&scalar_like_value(&args[2])?)?;
2536        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
2537
2538        if d1 < 1.0 || d2 < 1.0 || x < 0.0 {
2539            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2540                ExcelError::new_num(),
2541            )));
2542        }
2543
2544        let result = if cumulative {
2545            f_cdf(x, d1, d2)
2546        } else {
2547            f_pdf(x, d1, d2)
2548        };
2549        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2550            result,
2551        )))
2552    }
2553}
2554
2555/// F.INV(probability, deg_freedom1, deg_freedom2) - Inverse of F distribution
2556#[derive(Debug)]
2557pub struct FInvFn;
2558impl Function for FInvFn {
2559    func_caps!(PURE);
2560    fn name(&self) -> &'static str {
2561        "F.INV"
2562    }
2563    fn min_args(&self) -> usize {
2564        3
2565    }
2566    fn arg_schema(&self) -> &'static [ArgSchema] {
2567        use std::sync::LazyLock;
2568        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2569            vec![
2570                ArgSchema::number_lenient_scalar(),
2571                ArgSchema::number_lenient_scalar(),
2572                ArgSchema::number_lenient_scalar(),
2573            ]
2574        });
2575        &SCHEMA[..]
2576    }
2577    fn eval<'a, 'b, 'c>(
2578        &self,
2579        args: &'c [ArgumentHandle<'a, 'b>],
2580        _ctx: &dyn FunctionContext<'b>,
2581    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2582        let p = coerce_num(&scalar_like_value(&args[0])?)?;
2583        let d1 = coerce_num(&scalar_like_value(&args[1])?)?;
2584        let d2 = coerce_num(&scalar_like_value(&args[2])?)?;
2585
2586        if d1 < 1.0 || d2 < 1.0 {
2587            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2588                ExcelError::new_num(),
2589            )));
2590        }
2591
2592        match f_inv(p, d1, d2) {
2593            Some(result) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2594                result,
2595            ))),
2596            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2597                ExcelError::new_num(),
2598            ))),
2599        }
2600    }
2601}
2602
2603/// STANDARDIZE(x, mean, standard_dev) - Returns the normalized value
2604#[derive(Debug)]
2605pub struct StandardizeFn;
2606impl Function for StandardizeFn {
2607    func_caps!(PURE);
2608    fn name(&self) -> &'static str {
2609        "STANDARDIZE"
2610    }
2611    fn min_args(&self) -> usize {
2612        3
2613    }
2614    fn arg_schema(&self) -> &'static [ArgSchema] {
2615        use std::sync::LazyLock;
2616        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2617            vec![
2618                ArgSchema::number_lenient_scalar(),
2619                ArgSchema::number_lenient_scalar(),
2620                ArgSchema::number_lenient_scalar(),
2621            ]
2622        });
2623        &SCHEMA[..]
2624    }
2625    fn eval<'a, 'b, 'c>(
2626        &self,
2627        args: &'c [ArgumentHandle<'a, 'b>],
2628        _ctx: &dyn FunctionContext<'b>,
2629    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2630        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2631        let mean = coerce_num(&scalar_like_value(&args[1])?)?;
2632        let std_dev = coerce_num(&scalar_like_value(&args[2])?)?;
2633
2634        if std_dev <= 0.0 {
2635            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2636                ExcelError::new_num(),
2637            )));
2638        }
2639
2640        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2641            (x - mean) / std_dev,
2642        )))
2643    }
2644}
2645
2646/// Helper: Factorial function
2647fn factorial(n: i64) -> f64 {
2648    if n < 0 {
2649        return f64::NAN;
2650    }
2651    if n <= 1 {
2652        return 1.0;
2653    }
2654    // For large n, use gamma function: n! = Gamma(n+1)
2655    if n > 20 {
2656        return ln_gamma((n + 1) as f64).exp();
2657    }
2658    let mut result = 1.0;
2659    for i in 2..=n {
2660        result *= i as f64;
2661    }
2662    result
2663}
2664
2665/// Helper: Log of binomial coefficient (n choose k)
2666fn ln_binom(n: i64, k: i64) -> f64 {
2667    if k < 0 || k > n {
2668        return f64::NEG_INFINITY;
2669    }
2670    if k == 0 || k == n {
2671        return 0.0;
2672    }
2673    ln_gamma((n + 1) as f64) - ln_gamma((k + 1) as f64) - ln_gamma((n - k + 1) as f64)
2674}
2675
2676/// BINOM.DIST(number_s, trials, probability_s, cumulative) - Binomial distribution
2677#[derive(Debug)]
2678pub struct BinomDistFn;
2679impl Function for BinomDistFn {
2680    func_caps!(PURE);
2681    fn name(&self) -> &'static str {
2682        "BINOM.DIST"
2683    }
2684    fn min_args(&self) -> usize {
2685        4
2686    }
2687    fn arg_schema(&self) -> &'static [ArgSchema] {
2688        use std::sync::LazyLock;
2689        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2690            vec![
2691                ArgSchema::number_lenient_scalar(),
2692                ArgSchema::number_lenient_scalar(),
2693                ArgSchema::number_lenient_scalar(),
2694                ArgSchema::number_lenient_scalar(),
2695            ]
2696        });
2697        &SCHEMA[..]
2698    }
2699    fn eval<'a, 'b, 'c>(
2700        &self,
2701        args: &'c [ArgumentHandle<'a, 'b>],
2702        _ctx: &dyn FunctionContext<'b>,
2703    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2704        let k = coerce_num(&scalar_like_value(&args[0])?)?.trunc() as i64;
2705        let n = coerce_num(&scalar_like_value(&args[1])?)?.trunc() as i64;
2706        let p = coerce_num(&scalar_like_value(&args[2])?)?;
2707        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
2708
2709        if n < 0 || k < 0 || k > n || !(0.0..=1.0).contains(&p) {
2710            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2711                ExcelError::new_num(),
2712            )));
2713        }
2714
2715        let result = if cumulative {
2716            // CDF: sum from i=0 to k of P(X=i)
2717            let mut sum = 0.0;
2718            for i in 0..=k {
2719                let ln_prob =
2720                    ln_binom(n, i) + (i as f64) * p.ln() + ((n - i) as f64) * (1.0 - p).ln();
2721                sum += ln_prob.exp();
2722            }
2723            sum
2724        } else {
2725            // PMF: P(X=k)
2726            let ln_prob = ln_binom(n, k) + (k as f64) * p.ln() + ((n - k) as f64) * (1.0 - p).ln();
2727            ln_prob.exp()
2728        };
2729
2730        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2731            result,
2732        )))
2733    }
2734}
2735
2736/// POISSON.DIST(x, mean, cumulative) - Poisson distribution
2737#[derive(Debug)]
2738pub struct PoissonDistFn;
2739impl Function for PoissonDistFn {
2740    func_caps!(PURE);
2741    fn name(&self) -> &'static str {
2742        "POISSON.DIST"
2743    }
2744    fn min_args(&self) -> usize {
2745        3
2746    }
2747    fn arg_schema(&self) -> &'static [ArgSchema] {
2748        use std::sync::LazyLock;
2749        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2750            vec![
2751                ArgSchema::number_lenient_scalar(),
2752                ArgSchema::number_lenient_scalar(),
2753                ArgSchema::number_lenient_scalar(),
2754            ]
2755        });
2756        &SCHEMA[..]
2757    }
2758    fn eval<'a, 'b, 'c>(
2759        &self,
2760        args: &'c [ArgumentHandle<'a, 'b>],
2761        _ctx: &dyn FunctionContext<'b>,
2762    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2763        let k = coerce_num(&scalar_like_value(&args[0])?)?.trunc() as i64;
2764        let lambda = coerce_num(&scalar_like_value(&args[1])?)?;
2765        let cumulative = coerce_num(&scalar_like_value(&args[2])?)? != 0.0;
2766
2767        if k < 0 || lambda < 0.0 {
2768            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2769                ExcelError::new_num(),
2770            )));
2771        }
2772
2773        let result = if cumulative {
2774            // CDF: sum from i=0 to k of P(X=i) = 1 - Q(k+1, lambda)
2775            // Using the regularized incomplete gamma function
2776            1.0 - gamma_p((k + 1) as f64, lambda)
2777        } else {
2778            // PMF: P(X=k) = lambda^k * e^(-lambda) / k!
2779            // Use log to avoid overflow
2780            let ln_prob = (k as f64) * lambda.ln() - lambda - ln_gamma((k + 1) as f64);
2781            ln_prob.exp()
2782        };
2783
2784        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2785            result,
2786        )))
2787    }
2788}
2789
2790/// EXPON.DIST(x, lambda, cumulative) - Exponential distribution
2791#[derive(Debug)]
2792pub struct ExponDistFn;
2793impl Function for ExponDistFn {
2794    func_caps!(PURE);
2795    fn name(&self) -> &'static str {
2796        "EXPON.DIST"
2797    }
2798    fn min_args(&self) -> usize {
2799        3
2800    }
2801    fn arg_schema(&self) -> &'static [ArgSchema] {
2802        use std::sync::LazyLock;
2803        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2804            vec![
2805                ArgSchema::number_lenient_scalar(),
2806                ArgSchema::number_lenient_scalar(),
2807                ArgSchema::number_lenient_scalar(),
2808            ]
2809        });
2810        &SCHEMA[..]
2811    }
2812    fn eval<'a, 'b, 'c>(
2813        &self,
2814        args: &'c [ArgumentHandle<'a, 'b>],
2815        _ctx: &dyn FunctionContext<'b>,
2816    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2817        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2818        let lambda = coerce_num(&scalar_like_value(&args[1])?)?;
2819        let cumulative = coerce_num(&scalar_like_value(&args[2])?)? != 0.0;
2820
2821        if x < 0.0 || lambda <= 0.0 {
2822            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2823                ExcelError::new_num(),
2824            )));
2825        }
2826
2827        let result = if cumulative {
2828            // CDF: 1 - e^(-lambda*x)
2829            1.0 - (-lambda * x).exp()
2830        } else {
2831            // PDF: lambda * e^(-lambda*x)
2832            lambda * (-lambda * x).exp()
2833        };
2834
2835        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2836            result,
2837        )))
2838    }
2839}
2840
2841/// GAMMA.DIST(x, alpha, beta, cumulative) - Gamma distribution
2842#[derive(Debug)]
2843pub struct GammaDistFn;
2844impl Function for GammaDistFn {
2845    func_caps!(PURE);
2846    fn name(&self) -> &'static str {
2847        "GAMMA.DIST"
2848    }
2849    fn min_args(&self) -> usize {
2850        4
2851    }
2852    fn arg_schema(&self) -> &'static [ArgSchema] {
2853        use std::sync::LazyLock;
2854        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2855            vec![
2856                ArgSchema::number_lenient_scalar(),
2857                ArgSchema::number_lenient_scalar(),
2858                ArgSchema::number_lenient_scalar(),
2859                ArgSchema::number_lenient_scalar(),
2860            ]
2861        });
2862        &SCHEMA[..]
2863    }
2864    fn eval<'a, 'b, 'c>(
2865        &self,
2866        args: &'c [ArgumentHandle<'a, 'b>],
2867        _ctx: &dyn FunctionContext<'b>,
2868    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2869        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2870        let alpha = coerce_num(&scalar_like_value(&args[1])?)?; // shape
2871        let beta = coerce_num(&scalar_like_value(&args[2])?)?; // scale
2872        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
2873
2874        if x < 0.0 || alpha <= 0.0 || beta <= 0.0 {
2875            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2876                ExcelError::new_num(),
2877            )));
2878        }
2879
2880        let result = if cumulative {
2881            // CDF: P(alpha, x/beta) where P is the regularized lower incomplete gamma
2882            gamma_p(alpha, x / beta)
2883        } else {
2884            // PDF: x^(alpha-1) * e^(-x/beta) / (beta^alpha * Gamma(alpha))
2885            let ln_pdf = (alpha - 1.0) * x.ln() - x / beta - alpha * beta.ln() - ln_gamma(alpha);
2886            ln_pdf.exp()
2887        };
2888
2889        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2890            result,
2891        )))
2892    }
2893}
2894
2895/// WEIBULL.DIST(x, alpha, beta, cumulative) - Weibull distribution
2896/// alpha = shape parameter, beta = scale parameter
2897#[derive(Debug)]
2898pub struct WeibullDistFn;
2899impl Function for WeibullDistFn {
2900    func_caps!(PURE);
2901    fn name(&self) -> &'static str {
2902        "WEIBULL.DIST"
2903    }
2904    fn min_args(&self) -> usize {
2905        4
2906    }
2907    fn arg_schema(&self) -> &'static [ArgSchema] {
2908        use std::sync::LazyLock;
2909        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2910            vec![
2911                ArgSchema::number_lenient_scalar(),
2912                ArgSchema::number_lenient_scalar(),
2913                ArgSchema::number_lenient_scalar(),
2914                ArgSchema::number_lenient_scalar(),
2915            ]
2916        });
2917        &SCHEMA[..]
2918    }
2919    fn eval<'a, 'b, 'c>(
2920        &self,
2921        args: &'c [ArgumentHandle<'a, 'b>],
2922        _ctx: &dyn FunctionContext<'b>,
2923    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2924        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2925        let alpha = coerce_num(&scalar_like_value(&args[1])?)?; // shape
2926        let beta = coerce_num(&scalar_like_value(&args[2])?)?; // scale
2927        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
2928
2929        if x < 0.0 || alpha <= 0.0 || beta <= 0.0 {
2930            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
2931                ExcelError::new_num(),
2932            )));
2933        }
2934
2935        let result = if cumulative {
2936            // CDF: 1 - e^(-(x/beta)^alpha)
2937            1.0 - (-(x / beta).powf(alpha)).exp()
2938        } else {
2939            // PDF: (alpha/beta) * (x/beta)^(alpha-1) * e^(-(x/beta)^alpha)
2940            if x == 0.0 {
2941                if alpha < 1.0 {
2942                    f64::INFINITY
2943                } else if alpha == 1.0 {
2944                    alpha / beta
2945                } else {
2946                    0.0
2947                }
2948            } else {
2949                (alpha / beta) * (x / beta).powf(alpha - 1.0) * (-(x / beta).powf(alpha)).exp()
2950            }
2951        };
2952
2953        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
2954            result,
2955        )))
2956    }
2957}
2958
2959/// BETA.DIST(x, alpha, beta, cumulative, [A], [B]) - Beta distribution
2960/// A and B are optional bounds, defaults to 0 and 1
2961#[derive(Debug)]
2962pub struct BetaDistFn;
2963impl Function for BetaDistFn {
2964    func_caps!(PURE);
2965    fn name(&self) -> &'static str {
2966        "BETA.DIST"
2967    }
2968    fn min_args(&self) -> usize {
2969        4
2970    }
2971    fn variadic(&self) -> bool {
2972        true
2973    }
2974    fn arg_schema(&self) -> &'static [ArgSchema] {
2975        use std::sync::LazyLock;
2976        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
2977            vec![
2978                ArgSchema::number_lenient_scalar(),
2979                ArgSchema::number_lenient_scalar(),
2980                ArgSchema::number_lenient_scalar(),
2981                ArgSchema::number_lenient_scalar(),
2982                ArgSchema::number_lenient_scalar(),
2983                ArgSchema::number_lenient_scalar(),
2984            ]
2985        });
2986        &SCHEMA[..]
2987    }
2988    fn eval<'a, 'b, 'c>(
2989        &self,
2990        args: &'c [ArgumentHandle<'a, 'b>],
2991        _ctx: &dyn FunctionContext<'b>,
2992    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
2993        let x = coerce_num(&scalar_like_value(&args[0])?)?;
2994        let alpha = coerce_num(&scalar_like_value(&args[1])?)?;
2995        let beta_param = coerce_num(&scalar_like_value(&args[2])?)?;
2996        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
2997
2998        // Optional bounds A and B (default 0 and 1)
2999        let a = if args.len() > 4 {
3000            coerce_num(&scalar_like_value(&args[4])?)?
3001        } else {
3002            0.0
3003        };
3004        let b = if args.len() > 5 {
3005            coerce_num(&scalar_like_value(&args[5])?)?
3006        } else {
3007            1.0
3008        };
3009
3010        if alpha <= 0.0 || beta_param <= 0.0 || a >= b {
3011            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3012                ExcelError::new_num(),
3013            )));
3014        }
3015
3016        // x must be in [a, b]
3017        if x < a || x > b {
3018            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3019                ExcelError::new_num(),
3020            )));
3021        }
3022
3023        // Transform x to standard [0,1] interval
3024        let x_std = (x - a) / (b - a);
3025
3026        let result = if cumulative {
3027            // CDF: I_x(alpha, beta) - regularized incomplete beta function
3028            beta_i(x_std, alpha, beta_param)
3029        } else {
3030            // PDF: (x-A)^(alpha-1) * (B-x)^(beta-1) / ((B-A)^(alpha+beta-1) * B(alpha, beta))
3031            let ln_beta = ln_gamma(alpha) + ln_gamma(beta_param) - ln_gamma(alpha + beta_param);
3032            let scale = b - a;
3033            if (x_std == 0.0 && alpha < 1.0) || (x_std == 1.0 && beta_param < 1.0) {
3034                f64::INFINITY
3035            } else if x_std == 0.0 {
3036                if alpha == 1.0 {
3037                    (1.0 - x_std).powf(beta_param - 1.0) / (scale * ln_beta.exp())
3038                } else {
3039                    0.0
3040                }
3041            } else if x_std == 1.0 {
3042                if beta_param == 1.0 {
3043                    x_std.powf(alpha - 1.0) / (scale * ln_beta.exp())
3044                } else {
3045                    0.0
3046                }
3047            } else {
3048                let ln_pdf =
3049                    (alpha - 1.0) * x_std.ln() + (beta_param - 1.0) * (1.0 - x_std).ln() - ln_beta;
3050                ln_pdf.exp() / scale
3051            }
3052        };
3053
3054        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3055            result,
3056        )))
3057    }
3058}
3059
3060/// NEGBINOM.DIST(number_f, number_s, probability_s, cumulative) - Negative binomial distribution
3061/// number_f = number of failures, number_s = threshold number of successes, probability_s = probability of success
3062#[derive(Debug)]
3063pub struct NegbinomDistFn;
3064impl Function for NegbinomDistFn {
3065    func_caps!(PURE);
3066    fn name(&self) -> &'static str {
3067        "NEGBINOM.DIST"
3068    }
3069    fn min_args(&self) -> usize {
3070        4
3071    }
3072    fn arg_schema(&self) -> &'static [ArgSchema] {
3073        use std::sync::LazyLock;
3074        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
3075            vec![
3076                ArgSchema::number_lenient_scalar(),
3077                ArgSchema::number_lenient_scalar(),
3078                ArgSchema::number_lenient_scalar(),
3079                ArgSchema::number_lenient_scalar(),
3080            ]
3081        });
3082        &SCHEMA[..]
3083    }
3084    fn eval<'a, 'b, 'c>(
3085        &self,
3086        args: &'c [ArgumentHandle<'a, 'b>],
3087        _ctx: &dyn FunctionContext<'b>,
3088    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3089        let number_f = coerce_num(&scalar_like_value(&args[0])?)?.trunc() as i64; // number of failures
3090        let number_s = coerce_num(&scalar_like_value(&args[1])?)?.trunc() as i64; // number of successes
3091        let prob_s = coerce_num(&scalar_like_value(&args[2])?)?; // probability of success
3092        let cumulative = coerce_num(&scalar_like_value(&args[3])?)? != 0.0;
3093
3094        if number_f < 0 || number_s < 1 || prob_s <= 0.0 || prob_s >= 1.0 {
3095            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3096                ExcelError::new_num(),
3097            )));
3098        }
3099
3100        let result = if cumulative {
3101            // CDF: sum from i=0 to number_f of P(X=i)
3102            // This is equivalent to I_{prob_s}(number_s, number_f + 1) using regularized beta
3103            beta_i(prob_s, number_s as f64, (number_f + 1) as f64)
3104        } else {
3105            // PMF: C(number_f + number_s - 1, number_s - 1) * prob_s^number_s * (1-prob_s)^number_f
3106            // = C(k + r - 1, r - 1) * p^r * (1-p)^k where k = number_f, r = number_s
3107            let ln_prob = ln_binom(number_f + number_s - 1, number_s - 1)
3108                + (number_s as f64) * prob_s.ln()
3109                + (number_f as f64) * (1.0 - prob_s).ln();
3110            ln_prob.exp()
3111        };
3112
3113        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3114            result,
3115        )))
3116    }
3117}
3118
3119/// HYPGEOM.DIST(sample_s, number_sample, population_s, number_pop, cumulative) - Hypergeometric distribution
3120/// sample_s = number of successes in sample
3121/// number_sample = sample size
3122/// population_s = number of successes in population
3123/// number_pop = population size
3124#[derive(Debug)]
3125pub struct HypgeomDistFn;
3126impl Function for HypgeomDistFn {
3127    func_caps!(PURE);
3128    fn name(&self) -> &'static str {
3129        "HYPGEOM.DIST"
3130    }
3131    fn min_args(&self) -> usize {
3132        5
3133    }
3134    fn arg_schema(&self) -> &'static [ArgSchema] {
3135        use std::sync::LazyLock;
3136        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
3137            vec![
3138                ArgSchema::number_lenient_scalar(),
3139                ArgSchema::number_lenient_scalar(),
3140                ArgSchema::number_lenient_scalar(),
3141                ArgSchema::number_lenient_scalar(),
3142                ArgSchema::number_lenient_scalar(),
3143            ]
3144        });
3145        &SCHEMA[..]
3146    }
3147    fn eval<'a, 'b, 'c>(
3148        &self,
3149        args: &'c [ArgumentHandle<'a, 'b>],
3150        _ctx: &dyn FunctionContext<'b>,
3151    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3152        let sample_s = coerce_num(&scalar_like_value(&args[0])?)?.trunc() as i64; // successes in sample
3153        let number_sample = coerce_num(&scalar_like_value(&args[1])?)?.trunc() as i64; // sample size
3154        let population_s = coerce_num(&scalar_like_value(&args[2])?)?.trunc() as i64; // successes in population
3155        let number_pop = coerce_num(&scalar_like_value(&args[3])?)?.trunc() as i64; // population size
3156        let cumulative = coerce_num(&scalar_like_value(&args[4])?)? != 0.0;
3157
3158        // Validation
3159        if number_pop <= 0
3160            || population_s < 0
3161            || population_s > number_pop
3162            || number_sample < 0
3163            || number_sample > number_pop
3164            || sample_s < 0
3165        {
3166            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3167                ExcelError::new_num(),
3168            )));
3169        }
3170
3171        // sample_s must be at least max(0, number_sample - (number_pop - population_s))
3172        // and at most min(number_sample, population_s)
3173        let min_successes = 0.max(number_sample - (number_pop - population_s));
3174        let max_successes = number_sample.min(population_s);
3175
3176        if sample_s < min_successes || sample_s > max_successes {
3177            // Return 0 for PMF, or appropriate CDF value
3178            if cumulative {
3179                if sample_s < min_successes {
3180                    return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0)));
3181                } else {
3182                    return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(1.0)));
3183                }
3184            } else {
3185                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0)));
3186            }
3187        }
3188
3189        let result = if cumulative {
3190            // CDF: sum from i=min_successes to sample_s of P(X=i)
3191            let mut sum = 0.0;
3192            for i in min_successes..=sample_s {
3193                sum += hypgeom_pmf(i, number_sample, population_s, number_pop);
3194            }
3195            sum
3196        } else {
3197            // PMF: C(population_s, sample_s) * C(number_pop - population_s, number_sample - sample_s) / C(number_pop, number_sample)
3198            hypgeom_pmf(sample_s, number_sample, population_s, number_pop)
3199        };
3200
3201        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3202            result,
3203        )))
3204    }
3205}
3206
3207/// Helper: Hypergeometric PMF
3208fn hypgeom_pmf(k: i64, n: i64, k_pop: i64, n_pop: i64) -> f64 {
3209    // P(X=k) = C(K, k) * C(N-K, n-k) / C(N, n)
3210    // Using logs to avoid overflow
3211    let ln_prob = ln_binom(k_pop, k) + ln_binom(n_pop - k_pop, n - k) - ln_binom(n_pop, n);
3212    ln_prob.exp()
3213}
3214
3215/* ═══════════════════════════════════════════════════════════════════════════
3216COVARIANCE AND CORRELATION FUNCTIONS
3217═══════════════════════════════════════════════════════════════════════════ */
3218
3219/// COVARIANCE.P(array1, array2) - Population covariance
3220#[derive(Debug)]
3221pub struct CovariancePFn;
3222impl Function for CovariancePFn {
3223    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3224    fn name(&self) -> &'static str {
3225        "COVARIANCE.P"
3226    }
3227    fn aliases(&self) -> &'static [&'static str] {
3228        &["COVAR"]
3229    }
3230    fn min_args(&self) -> usize {
3231        2
3232    }
3233    fn arg_schema(&self) -> &'static [ArgSchema] {
3234        &ARG_RANGE_NUM_LENIENT_ONE[..]
3235    }
3236    fn eval<'a, 'b, 'c>(
3237        &self,
3238        args: &'c [ArgumentHandle<'a, 'b>],
3239        _ctx: &dyn FunctionContext<'b>,
3240    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3241        let (y, x) = match collect_paired_arrays(args) {
3242            Ok(v) => v,
3243            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
3244        };
3245
3246        let n = x.len() as f64;
3247        let mean_x = x.iter().sum::<f64>() / n;
3248        let mean_y = y.iter().sum::<f64>() / n;
3249
3250        let mut sum_xy = 0.0;
3251        for i in 0..x.len() {
3252            let dx = x[i] - mean_x;
3253            let dy = y[i] - mean_y;
3254            sum_xy += dx * dy;
3255        }
3256
3257        // Population covariance divides by n
3258        let covar = sum_xy / n;
3259        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3260            covar,
3261        )))
3262    }
3263}
3264
3265/// COVARIANCE.S(array1, array2) - Sample covariance
3266#[derive(Debug)]
3267pub struct CovarianceSFn;
3268impl Function for CovarianceSFn {
3269    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3270    fn name(&self) -> &'static str {
3271        "COVARIANCE.S"
3272    }
3273    fn min_args(&self) -> usize {
3274        2
3275    }
3276    fn arg_schema(&self) -> &'static [ArgSchema] {
3277        &ARG_RANGE_NUM_LENIENT_ONE[..]
3278    }
3279    fn eval<'a, 'b, 'c>(
3280        &self,
3281        args: &'c [ArgumentHandle<'a, 'b>],
3282        _ctx: &dyn FunctionContext<'b>,
3283    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3284        let (y, x) = match collect_paired_arrays(args) {
3285            Ok(v) => v,
3286            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
3287        };
3288
3289        let n = x.len();
3290        if n < 2 {
3291            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3292                ExcelError::new_div(),
3293            )));
3294        }
3295
3296        let mean_x = x.iter().sum::<f64>() / n as f64;
3297        let mean_y = y.iter().sum::<f64>() / n as f64;
3298
3299        let mut sum_xy = 0.0;
3300        for i in 0..n {
3301            let dx = x[i] - mean_x;
3302            let dy = y[i] - mean_y;
3303            sum_xy += dx * dy;
3304        }
3305
3306        // Sample covariance divides by (n - 1)
3307        let covar = sum_xy / (n - 1) as f64;
3308        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3309            covar,
3310        )))
3311    }
3312}
3313
3314/// PEARSON(array1, array2) - Pearson correlation coefficient (same as CORREL)
3315#[derive(Debug)]
3316pub struct PearsonFn;
3317impl Function for PearsonFn {
3318    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3319    fn name(&self) -> &'static str {
3320        "PEARSON"
3321    }
3322    fn min_args(&self) -> usize {
3323        2
3324    }
3325    fn arg_schema(&self) -> &'static [ArgSchema] {
3326        &ARG_RANGE_NUM_LENIENT_ONE[..]
3327    }
3328    fn eval<'a, 'b, 'c>(
3329        &self,
3330        args: &'c [ArgumentHandle<'a, 'b>],
3331        _ctx: &dyn FunctionContext<'b>,
3332    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3333        let (y, x) = match collect_paired_arrays(args) {
3334            Ok(v) => v,
3335            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
3336        };
3337
3338        let n = x.len() as f64;
3339        let mean_x = x.iter().sum::<f64>() / n;
3340        let mean_y = y.iter().sum::<f64>() / n;
3341
3342        let mut sum_xy = 0.0;
3343        let mut sum_x2 = 0.0;
3344        let mut sum_y2 = 0.0;
3345
3346        for i in 0..x.len() {
3347            let dx = x[i] - mean_x;
3348            let dy = y[i] - mean_y;
3349            sum_xy += dx * dy;
3350            sum_x2 += dx * dx;
3351            sum_y2 += dy * dy;
3352        }
3353
3354        let denom = (sum_x2 * sum_y2).sqrt();
3355        if denom == 0.0 {
3356            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3357                ExcelError::new_div(),
3358            )));
3359        }
3360
3361        let correl = sum_xy / denom;
3362        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3363            correl,
3364        )))
3365    }
3366}
3367
3368/// RSQ(known_y's, known_x's) - R-squared value (square of correlation)
3369#[derive(Debug)]
3370pub struct RsqFn;
3371impl Function for RsqFn {
3372    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3373    fn name(&self) -> &'static str {
3374        "RSQ"
3375    }
3376    fn min_args(&self) -> usize {
3377        2
3378    }
3379    fn arg_schema(&self) -> &'static [ArgSchema] {
3380        &ARG_RANGE_NUM_LENIENT_ONE[..]
3381    }
3382    fn eval<'a, 'b, 'c>(
3383        &self,
3384        args: &'c [ArgumentHandle<'a, 'b>],
3385        _ctx: &dyn FunctionContext<'b>,
3386    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3387        let (y, x) = match collect_paired_arrays(args) {
3388            Ok(v) => v,
3389            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
3390        };
3391
3392        let n = x.len() as f64;
3393        let mean_x = x.iter().sum::<f64>() / n;
3394        let mean_y = y.iter().sum::<f64>() / n;
3395
3396        let mut sum_xy = 0.0;
3397        let mut sum_x2 = 0.0;
3398        let mut sum_y2 = 0.0;
3399
3400        for i in 0..x.len() {
3401            let dx = x[i] - mean_x;
3402            let dy = y[i] - mean_y;
3403            sum_xy += dx * dy;
3404            sum_x2 += dx * dx;
3405            sum_y2 += dy * dy;
3406        }
3407
3408        let denom = sum_x2 * sum_y2;
3409        if denom == 0.0 {
3410            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3411                ExcelError::new_div(),
3412            )));
3413        }
3414
3415        // R-squared = r^2 = (sum_xy)^2 / (sum_x2 * sum_y2)
3416        let rsq = (sum_xy * sum_xy) / denom;
3417        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(rsq)))
3418    }
3419}
3420
3421/// STEYX(known_y's, known_x's) - Standard error of the predicted y-value
3422#[derive(Debug)]
3423pub struct SteyxFn;
3424impl Function for SteyxFn {
3425    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3426    fn name(&self) -> &'static str {
3427        "STEYX"
3428    }
3429    fn min_args(&self) -> usize {
3430        2
3431    }
3432    fn arg_schema(&self) -> &'static [ArgSchema] {
3433        &ARG_RANGE_NUM_LENIENT_ONE[..]
3434    }
3435    fn eval<'a, 'b, 'c>(
3436        &self,
3437        args: &'c [ArgumentHandle<'a, 'b>],
3438        _ctx: &dyn FunctionContext<'b>,
3439    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3440        let (y, x) = match collect_paired_arrays(args) {
3441            Ok(v) => v,
3442            Err(e) => return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(e))),
3443        };
3444
3445        let n = x.len();
3446        if n < 3 {
3447            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3448                ExcelError::new_div(),
3449            )));
3450        }
3451
3452        let n_f = n as f64;
3453        let mean_x = x.iter().sum::<f64>() / n_f;
3454        let mean_y = y.iter().sum::<f64>() / n_f;
3455
3456        let mut sum_xy = 0.0;
3457        let mut sum_x2 = 0.0;
3458        let mut sum_y2 = 0.0;
3459
3460        for i in 0..n {
3461            let dx = x[i] - mean_x;
3462            let dy = y[i] - mean_y;
3463            sum_xy += dx * dy;
3464            sum_x2 += dx * dx;
3465            sum_y2 += dy * dy;
3466        }
3467
3468        if sum_x2 == 0.0 {
3469            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3470                ExcelError::new_div(),
3471            )));
3472        }
3473
3474        // STEYX = sqrt((sum_y2 - (sum_xy)^2 / sum_x2) / (n - 2))
3475        let sse = sum_y2 - (sum_xy * sum_xy) / sum_x2;
3476        if sse < 0.0 {
3477            // This can happen due to floating point errors; return 0 in such case
3478            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(0.0)));
3479        }
3480        let steyx = (sse / (n_f - 2.0)).sqrt();
3481        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3482            steyx,
3483        )))
3484    }
3485}
3486
3487/* ─────────────────────────── SKEW ──────────────────────────── */
3488
3489/// SKEW(number1, [number2], ...) - Skewness of a distribution
3490#[derive(Debug)]
3491pub struct SkewFn;
3492impl Function for SkewFn {
3493    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3494    fn name(&self) -> &'static str {
3495        "SKEW"
3496    }
3497    fn min_args(&self) -> usize {
3498        1
3499    }
3500    fn variadic(&self) -> bool {
3501        true
3502    }
3503    fn arg_schema(&self) -> &'static [ArgSchema] {
3504        &ARG_RANGE_NUM_LENIENT_ONE[..]
3505    }
3506    fn eval<'a, 'b, 'c>(
3507        &self,
3508        args: &'c [ArgumentHandle<'a, 'b>],
3509        _ctx: &dyn FunctionContext<'b>,
3510    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3511        let nums = collect_numeric_stats(args)?;
3512        let n = nums.len();
3513
3514        // SKEW requires at least 3 data points
3515        if n < 3 {
3516            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3517                ExcelError::new_div(),
3518            )));
3519        }
3520
3521        let n_f = n as f64;
3522        let mean = nums.iter().sum::<f64>() / n_f;
3523
3524        // Calculate sample standard deviation
3525        let mut sum_sq = 0.0;
3526        for &v in &nums {
3527            let d = v - mean;
3528            sum_sq += d * d;
3529        }
3530        let stdev = (sum_sq / (n_f - 1.0)).sqrt();
3531
3532        if stdev == 0.0 {
3533            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3534                ExcelError::new_div(),
3535            )));
3536        }
3537
3538        // Calculate sum of cubed deviations normalized by stdev
3539        let mut sum_cubed = 0.0;
3540        for &v in &nums {
3541            let d = (v - mean) / stdev;
3542            sum_cubed += d * d * d;
3543        }
3544
3545        // Excel SKEW formula: n / ((n-1)*(n-2)) * sum((xi - mean)/stdev)^3
3546        let skew = (n_f / ((n_f - 1.0) * (n_f - 2.0))) * sum_cubed;
3547        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(skew)))
3548    }
3549}
3550
3551/* ─────────────────────────── KURT ──────────────────────────── */
3552
3553/// KURT(number1, [number2], ...) - Kurtosis of a distribution
3554#[derive(Debug)]
3555pub struct KurtFn;
3556impl Function for KurtFn {
3557    func_caps!(PURE, NUMERIC_ONLY, REDUCTION);
3558    fn name(&self) -> &'static str {
3559        "KURT"
3560    }
3561    fn min_args(&self) -> usize {
3562        1
3563    }
3564    fn variadic(&self) -> bool {
3565        true
3566    }
3567    fn arg_schema(&self) -> &'static [ArgSchema] {
3568        &ARG_RANGE_NUM_LENIENT_ONE[..]
3569    }
3570    fn eval<'a, 'b, 'c>(
3571        &self,
3572        args: &'c [ArgumentHandle<'a, 'b>],
3573        _ctx: &dyn FunctionContext<'b>,
3574    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3575        let nums = collect_numeric_stats(args)?;
3576        let n = nums.len();
3577
3578        // KURT requires at least 4 data points
3579        if n < 4 {
3580            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3581                ExcelError::new_div(),
3582            )));
3583        }
3584
3585        let n_f = n as f64;
3586        let mean = nums.iter().sum::<f64>() / n_f;
3587
3588        // Calculate sample standard deviation
3589        let mut sum_sq = 0.0;
3590        for &v in &nums {
3591            let d = v - mean;
3592            sum_sq += d * d;
3593        }
3594        let stdev = (sum_sq / (n_f - 1.0)).sqrt();
3595
3596        if stdev == 0.0 {
3597            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3598                ExcelError::new_div(),
3599            )));
3600        }
3601
3602        // Calculate sum of fourth powers of deviations normalized by stdev
3603        let mut sum_fourth = 0.0;
3604        for &v in &nums {
3605            let d = (v - mean) / stdev;
3606            sum_fourth += d * d * d * d;
3607        }
3608
3609        // Excel KURT formula (excess kurtosis):
3610        // n*(n+1) / ((n-1)*(n-2)*(n-3)) * sum((xi - mean)/stdev)^4 - 3*(n-1)^2 / ((n-2)*(n-3))
3611        let term1 = (n_f * (n_f + 1.0)) / ((n_f - 1.0) * (n_f - 2.0) * (n_f - 3.0)) * sum_fourth;
3612        let term2 = (3.0 * (n_f - 1.0) * (n_f - 1.0)) / ((n_f - 2.0) * (n_f - 3.0));
3613        let kurt = term1 - term2;
3614        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(kurt)))
3615    }
3616}
3617
3618/* ─────────────────────────── FISHER ──────────────────────────── */
3619
3620/// FISHER(x) - Fisher transformation
3621#[derive(Debug)]
3622pub struct FisherFn;
3623impl Function for FisherFn {
3624    func_caps!(PURE, NUMERIC_ONLY);
3625    fn name(&self) -> &'static str {
3626        "FISHER"
3627    }
3628    fn min_args(&self) -> usize {
3629        1
3630    }
3631    fn arg_schema(&self) -> &'static [ArgSchema] {
3632        &ARG_RANGE_NUM_LENIENT_ONE[..]
3633    }
3634    fn eval<'a, 'b, 'c>(
3635        &self,
3636        args: &'c [ArgumentHandle<'a, 'b>],
3637        _ctx: &dyn FunctionContext<'b>,
3638    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3639        let x = coerce_num(&scalar_like_value(&args[0])?)?;
3640
3641        // FISHER requires -1 < x < 1
3642        if x <= -1.0 || x >= 1.0 {
3643            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3644                ExcelError::new_num(),
3645            )));
3646        }
3647
3648        // Fisher transformation: 0.5 * ln((1 + x) / (1 - x))
3649        let fisher = 0.5 * ((1.0 + x) / (1.0 - x)).ln();
3650        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3651            fisher,
3652        )))
3653    }
3654}
3655
3656/* ─────────────────────────── FISHERINV ──────────────────────────── */
3657
3658/// FISHERINV(y) - Inverse Fisher transformation
3659#[derive(Debug)]
3660pub struct FisherInvFn;
3661impl Function for FisherInvFn {
3662    func_caps!(PURE, NUMERIC_ONLY);
3663    fn name(&self) -> &'static str {
3664        "FISHERINV"
3665    }
3666    fn min_args(&self) -> usize {
3667        1
3668    }
3669    fn arg_schema(&self) -> &'static [ArgSchema] {
3670        &ARG_RANGE_NUM_LENIENT_ONE[..]
3671    }
3672    fn eval<'a, 'b, 'c>(
3673        &self,
3674        args: &'c [ArgumentHandle<'a, 'b>],
3675        _ctx: &dyn FunctionContext<'b>,
3676    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3677        let y = coerce_num(&scalar_like_value(&args[0])?)?;
3678
3679        // Inverse Fisher transformation: (e^(2y) - 1) / (e^(2y) + 1)
3680        let e2y = (2.0 * y).exp();
3681        let fisherinv = (e2y - 1.0) / (e2y + 1.0);
3682        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3683            fisherinv,
3684        )))
3685    }
3686}
3687
3688/* ─────────────────────────── FORECAST.LINEAR ──────────────────────────── */
3689
3690/// FORECAST.LINEAR(x, known_y's, known_x's) - Returns predicted y value for x using linear regression
3691/// The formula is: y = intercept + slope * x
3692/// where slope = sum((xi - mean_x)(yi - mean_y)) / sum((xi - mean_x)^2)
3693/// and intercept = mean_y - slope * mean_x
3694#[derive(Debug)]
3695pub struct ForecastLinearFn;
3696impl Function for ForecastLinearFn {
3697    func_caps!(PURE, NUMERIC_ONLY);
3698    fn name(&self) -> &'static str {
3699        "FORECAST.LINEAR"
3700    }
3701    fn aliases(&self) -> &'static [&'static str] {
3702        &["FORECAST"]
3703    }
3704    fn min_args(&self) -> usize {
3705        3
3706    }
3707    fn arg_schema(&self) -> &'static [ArgSchema] {
3708        &ARG_RANGE_NUM_LENIENT_ONE[..]
3709    }
3710    fn eval<'a, 'b, 'c>(
3711        &self,
3712        args: &'c [ArgumentHandle<'a, 'b>],
3713        _ctx: &dyn FunctionContext<'b>,
3714    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3715        // args[0] = x value to forecast
3716        // args[1] = known_y's
3717        // args[2] = known_x's
3718        let x = match coerce_num(&scalar_like_value(&args[0])?) {
3719            Ok(n) => n,
3720            Err(_) => {
3721                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3722                    ExcelError::new_value(),
3723                )));
3724            }
3725        };
3726
3727        let y_vals = collect_numeric_stats(&args[1..2])?;
3728        let x_vals = collect_numeric_stats(&args[2..3])?;
3729
3730        // Arrays must have same length
3731        if y_vals.len() != x_vals.len() {
3732            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3733                ExcelError::new_na(),
3734            )));
3735        }
3736
3737        if y_vals.is_empty() {
3738            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3739                ExcelError::new_na(),
3740            )));
3741        }
3742
3743        let n = x_vals.len() as f64;
3744        let mean_x = x_vals.iter().sum::<f64>() / n;
3745        let mean_y = y_vals.iter().sum::<f64>() / n;
3746
3747        let mut sum_xy = 0.0;
3748        let mut sum_x2 = 0.0;
3749
3750        for i in 0..x_vals.len() {
3751            let dx = x_vals[i] - mean_x;
3752            let dy = y_vals[i] - mean_y;
3753            sum_xy += dx * dy;
3754            sum_x2 += dx * dx;
3755        }
3756
3757        if sum_x2 == 0.0 {
3758            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3759                ExcelError::new_div(),
3760            )));
3761        }
3762
3763        let slope = sum_xy / sum_x2;
3764        let intercept = mean_y - slope * mean_x;
3765        let forecast = intercept + slope * x;
3766
3767        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
3768            forecast,
3769        )))
3770    }
3771}
3772
3773/* ─────────────────────────── LINEST ──────────────────────────── */
3774
3775/// LINEST(known_y's, known_x's, [const], [stats]) - Returns statistics describing the linear trend
3776/// With stats=FALSE (default): returns 1x2 array [slope, intercept]
3777/// With stats=TRUE: returns 5x2 array with regression statistics
3778#[derive(Debug)]
3779pub struct LinestFn;
3780impl Function for LinestFn {
3781    func_caps!(PURE, NUMERIC_ONLY);
3782    fn name(&self) -> &'static str {
3783        "LINEST"
3784    }
3785    fn min_args(&self) -> usize {
3786        1
3787    }
3788    fn variadic(&self) -> bool {
3789        true
3790    }
3791    fn arg_schema(&self) -> &'static [ArgSchema] {
3792        &ARG_RANGE_NUM_LENIENT_ONE[..]
3793    }
3794    fn eval<'a, 'b, 'c>(
3795        &self,
3796        args: &'c [ArgumentHandle<'a, 'b>],
3797        _ctx: &dyn FunctionContext<'b>,
3798    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
3799        // args[0] = known_y's (required)
3800        // args[1] = known_x's (optional, defaults to {1,2,3,...})
3801        // args[2] = const (optional, default TRUE - whether to compute intercept)
3802        // args[3] = stats (optional, default FALSE - whether to return additional statistics)
3803
3804        let y_vals = collect_numeric_stats(&args[0..1])?;
3805
3806        if y_vals.is_empty() {
3807            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3808                ExcelError::new_na(),
3809            )));
3810        }
3811
3812        // Get known_x's or generate default {1, 2, 3, ...}
3813        let x_vals = if args.len() >= 2 {
3814            collect_numeric_stats(&args[1..2])?
3815        } else {
3816            (1..=y_vals.len()).map(|i| i as f64).collect()
3817        };
3818
3819        // Arrays must have same length
3820        if y_vals.len() != x_vals.len() {
3821            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3822                ExcelError::new_ref(),
3823            )));
3824        }
3825
3826        // Parse const argument (default TRUE)
3827        let use_const = if args.len() >= 3 {
3828            match scalar_like_value(&args[2])? {
3829                LiteralValue::Boolean(b) => b,
3830                LiteralValue::Number(n) => n != 0.0,
3831                LiteralValue::Int(i) => i != 0,
3832                _ => true,
3833            }
3834        } else {
3835            true
3836        };
3837
3838        // Parse stats argument (default FALSE)
3839        let return_stats = if args.len() >= 4 {
3840            match scalar_like_value(&args[3])? {
3841                LiteralValue::Boolean(b) => b,
3842                LiteralValue::Number(n) => n != 0.0,
3843                LiteralValue::Int(i) => i != 0,
3844                _ => false,
3845            }
3846        } else {
3847            false
3848        };
3849
3850        let n = x_vals.len() as f64;
3851
3852        // Calculate regression coefficients
3853        let (slope, intercept) = if use_const {
3854            // Normal linear regression with intercept
3855            let mean_x = x_vals.iter().sum::<f64>() / n;
3856            let mean_y = y_vals.iter().sum::<f64>() / n;
3857
3858            let mut sum_xy = 0.0;
3859            let mut sum_x2 = 0.0;
3860
3861            for i in 0..x_vals.len() {
3862                let dx = x_vals[i] - mean_x;
3863                let dy = y_vals[i] - mean_y;
3864                sum_xy += dx * dy;
3865                sum_x2 += dx * dx;
3866            }
3867
3868            if sum_x2 == 0.0 {
3869                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3870                    ExcelError::new_div(),
3871                )));
3872            }
3873
3874            let slope = sum_xy / sum_x2;
3875            let intercept = mean_y - slope * mean_x;
3876            (slope, intercept)
3877        } else {
3878            // Regression through origin (intercept = 0)
3879            let mut sum_xy = 0.0;
3880            let mut sum_x2 = 0.0;
3881
3882            for i in 0..x_vals.len() {
3883                sum_xy += x_vals[i] * y_vals[i];
3884                sum_x2 += x_vals[i] * x_vals[i];
3885            }
3886
3887            if sum_x2 == 0.0 {
3888                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
3889                    ExcelError::new_div(),
3890                )));
3891            }
3892
3893            let slope = sum_xy / sum_x2;
3894            (slope, 0.0)
3895        };
3896
3897        if !return_stats {
3898            // Return just slope and intercept as 1x2 array: [[slope, intercept]]
3899            let row = vec![LiteralValue::Number(slope), LiteralValue::Number(intercept)];
3900            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(vec![
3901                row,
3902            ])));
3903        }
3904
3905        // Calculate additional statistics for stats=TRUE
3906        // Row 1: [slope, intercept]
3907        // Row 2: [se_slope, se_intercept]
3908        // Row 3: [r_squared, se_y]
3909        // Row 4: [F_statistic, df]
3910        // Row 5: [ss_reg, ss_resid]
3911
3912        let mean_y = y_vals.iter().sum::<f64>() / n;
3913
3914        // Calculate residuals and sums of squares
3915        let mut ss_resid = 0.0; // Sum of squared residuals
3916        let mut ss_tot = 0.0; // Total sum of squares
3917
3918        for i in 0..x_vals.len() {
3919            let y_pred = slope * x_vals[i] + intercept;
3920            let residual = y_vals[i] - y_pred;
3921            ss_resid += residual * residual;
3922            let dy_tot = y_vals[i] - mean_y;
3923            ss_tot += dy_tot * dy_tot;
3924        }
3925
3926        let ss_reg = ss_tot - ss_resid; // Regression sum of squares
3927
3928        // R-squared
3929        let r_squared = if ss_tot == 0.0 {
3930            1.0 // Perfect fit or all y values are the same
3931        } else {
3932            1.0 - (ss_resid / ss_tot)
3933        };
3934
3935        // Degrees of freedom
3936        let df = if use_const {
3937            (n as i64 - 2).max(1) as f64 // n - k - 1 where k=1 (one predictor)
3938        } else {
3939            (n as i64 - 1).max(1) as f64 // n - k when no intercept
3940        };
3941
3942        // Standard error of y estimate
3943        let se_y = if df > 0.0 {
3944            (ss_resid / df).sqrt()
3945        } else {
3946            0.0
3947        };
3948
3949        // Standard errors of coefficients
3950        let mean_x = x_vals.iter().sum::<f64>() / n;
3951        let mut sum_x2_centered = 0.0;
3952        let mut sum_x2_raw = 0.0;
3953        for &xi in &x_vals {
3954            sum_x2_centered += (xi - mean_x).powi(2);
3955            sum_x2_raw += xi * xi;
3956        }
3957
3958        let se_slope = if sum_x2_centered > 0.0 && df > 0.0 {
3959            se_y / sum_x2_centered.sqrt()
3960        } else {
3961            f64::NAN
3962        };
3963
3964        let se_intercept = if use_const && sum_x2_centered > 0.0 && df > 0.0 {
3965            se_y * (sum_x2_raw / (n * sum_x2_centered)).sqrt()
3966        } else {
3967            f64::NAN
3968        };
3969
3970        // F-statistic
3971        let f_stat = if ss_resid > 0.0 && df > 0.0 {
3972            (ss_reg / 1.0) / (ss_resid / df) // MSR / MSE
3973        } else if ss_resid == 0.0 {
3974            f64::INFINITY // Perfect fit
3975        } else {
3976            f64::NAN
3977        };
3978
3979        // Build 5x2 result array
3980        let rows = vec![
3981            vec![LiteralValue::Number(slope), LiteralValue::Number(intercept)],
3982            vec![
3983                LiteralValue::Number(se_slope),
3984                LiteralValue::Number(se_intercept),
3985            ],
3986            vec![LiteralValue::Number(r_squared), LiteralValue::Number(se_y)],
3987            vec![LiteralValue::Number(f_stat), LiteralValue::Number(df)],
3988            vec![LiteralValue::Number(ss_reg), LiteralValue::Number(ss_resid)],
3989        ];
3990
3991        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(rows)))
3992    }
3993}
3994
3995/* ─────────────────────────── CONFIDENCE.NORM ──────────────────────────── */
3996
3997/// CONFIDENCE.NORM(alpha, standard_dev, size) - Returns the confidence interval for a population mean
3998/// using a normal distribution.
3999/// Formula: z_crit * standard_dev / sqrt(size), where z_crit = NORM.S.INV(1 - alpha/2)
4000#[derive(Debug)]
4001pub struct ConfidenceNormFn;
4002impl Function for ConfidenceNormFn {
4003    func_caps!(PURE);
4004    fn name(&self) -> &'static str {
4005        "CONFIDENCE.NORM"
4006    }
4007    fn aliases(&self) -> &'static [&'static str] {
4008        &["CONFIDENCE"]
4009    }
4010    fn min_args(&self) -> usize {
4011        3
4012    }
4013    fn arg_schema(&self) -> &'static [ArgSchema] {
4014        use std::sync::LazyLock;
4015        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
4016            vec![
4017                ArgSchema::number_lenient_scalar(),
4018                ArgSchema::number_lenient_scalar(),
4019                ArgSchema::number_lenient_scalar(),
4020            ]
4021        });
4022        &SCHEMA[..]
4023    }
4024    fn eval<'a, 'b, 'c>(
4025        &self,
4026        args: &'c [ArgumentHandle<'a, 'b>],
4027        _ctx: &dyn FunctionContext<'b>,
4028    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4029        let alpha = coerce_num(&scalar_like_value(&args[0])?)?;
4030        let std_dev = coerce_num(&scalar_like_value(&args[1])?)?;
4031        let size = coerce_num(&scalar_like_value(&args[2])?)?;
4032
4033        // Validate inputs
4034        if alpha <= 0.0 || alpha >= 1.0 {
4035            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4036                ExcelError::new_num(),
4037            )));
4038        }
4039        if std_dev <= 0.0 {
4040            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4041                ExcelError::new_num(),
4042            )));
4043        }
4044        if size < 1.0 {
4045            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4046                ExcelError::new_num(),
4047            )));
4048        }
4049
4050        // z_crit = NORM.S.INV(1 - alpha/2)
4051        let z_crit = match std_norm_inv(1.0 - alpha / 2.0) {
4052            Some(z) => z,
4053            None => {
4054                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4055                    ExcelError::new_num(),
4056                )));
4057            }
4058        };
4059
4060        let result = z_crit * std_dev / size.sqrt();
4061        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
4062            result,
4063        )))
4064    }
4065}
4066
4067/* ─────────────────────────── CONFIDENCE.T ──────────────────────────── */
4068
4069/// CONFIDENCE.T(alpha, standard_dev, size) - Returns the confidence interval for a population mean
4070/// using a Student's t-distribution.
4071/// Formula: t_crit * standard_dev / sqrt(size), where t_crit = T.INV(1 - alpha/2, size - 1)
4072#[derive(Debug)]
4073pub struct ConfidenceTFn;
4074impl Function for ConfidenceTFn {
4075    func_caps!(PURE);
4076    fn name(&self) -> &'static str {
4077        "CONFIDENCE.T"
4078    }
4079    fn min_args(&self) -> usize {
4080        3
4081    }
4082    fn arg_schema(&self) -> &'static [ArgSchema] {
4083        use std::sync::LazyLock;
4084        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
4085            vec![
4086                ArgSchema::number_lenient_scalar(),
4087                ArgSchema::number_lenient_scalar(),
4088                ArgSchema::number_lenient_scalar(),
4089            ]
4090        });
4091        &SCHEMA[..]
4092    }
4093    fn eval<'a, 'b, 'c>(
4094        &self,
4095        args: &'c [ArgumentHandle<'a, 'b>],
4096        _ctx: &dyn FunctionContext<'b>,
4097    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4098        let alpha = coerce_num(&scalar_like_value(&args[0])?)?;
4099        let std_dev = coerce_num(&scalar_like_value(&args[1])?)?;
4100        let size = coerce_num(&scalar_like_value(&args[2])?)?;
4101
4102        // Validate inputs - size must be >= 2 for t-distribution (df = size - 1 >= 1)
4103        if alpha <= 0.0 || alpha >= 1.0 {
4104            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4105                ExcelError::new_num(),
4106            )));
4107        }
4108        if std_dev <= 0.0 {
4109            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4110                ExcelError::new_num(),
4111            )));
4112        }
4113        if size < 2.0 {
4114            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4115                ExcelError::new_num(),
4116            )));
4117        }
4118
4119        let df = size - 1.0;
4120
4121        // t_crit = T.INV(1 - alpha/2, df)
4122        let t_crit = match t_inv(1.0 - alpha / 2.0, df) {
4123            Some(t) => t,
4124            None => {
4125                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4126                    ExcelError::new_num(),
4127                )));
4128            }
4129        };
4130
4131        let result = t_crit * std_dev / size.sqrt();
4132        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
4133            result,
4134        )))
4135    }
4136}
4137
4138/* ─────────────────────────── Z.TEST ──────────────────────────── */
4139
4140/// Z.TEST(array, x, [sigma]) - Returns the one-tailed P-value of a z-test.
4141/// z = (mean(array) - x) / (sigma / sqrt(n))
4142/// Returns 1 - NORM.S.DIST(z, TRUE)
4143/// If sigma is omitted, uses the population standard deviation of the array.
4144#[derive(Debug)]
4145pub struct ZTestFn;
4146impl Function for ZTestFn {
4147    func_caps!(PURE);
4148    fn name(&self) -> &'static str {
4149        "Z.TEST"
4150    }
4151    fn aliases(&self) -> &'static [&'static str] {
4152        &["ZTEST"]
4153    }
4154    fn min_args(&self) -> usize {
4155        2
4156    }
4157    fn variadic(&self) -> bool {
4158        true
4159    }
4160    fn arg_schema(&self) -> &'static [ArgSchema] {
4161        use std::sync::LazyLock;
4162        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
4163            vec![
4164                {
4165                    let mut s = ArgSchema::number_lenient_scalar();
4166                    s.shape = crate::args::ShapeKind::Range;
4167                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
4168                    s
4169                },
4170                ArgSchema::number_lenient_scalar(),
4171                ArgSchema::number_lenient_scalar(), // optional sigma
4172            ]
4173        });
4174        &SCHEMA[..]
4175    }
4176    fn eval<'a, 'b, 'c>(
4177        &self,
4178        args: &'c [ArgumentHandle<'a, 'b>],
4179        _ctx: &dyn FunctionContext<'b>,
4180    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4181        // Collect numeric values from the array argument
4182        let data = collect_numeric_stats(&args[0..1])?;
4183
4184        if data.is_empty() {
4185            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4186                ExcelError::new_na(),
4187            )));
4188        }
4189
4190        let x = coerce_num(&scalar_like_value(&args[1])?)?;
4191
4192        let n = data.len() as f64;
4193        let mean: f64 = data.iter().sum::<f64>() / n;
4194
4195        // Calculate sigma: use provided value or compute population std dev
4196        let sigma = if args.len() > 2 {
4197            let s = coerce_num(&scalar_like_value(&args[2])?)?;
4198            if s <= 0.0 {
4199                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4200                    ExcelError::new_num(),
4201                )));
4202            }
4203            s
4204        } else {
4205            // Population standard deviation
4206            let variance: f64 = data.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n;
4207            let std_dev = variance.sqrt();
4208            if std_dev == 0.0 {
4209                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4210                    ExcelError::new_div(),
4211                )));
4212            }
4213            std_dev
4214        };
4215
4216        // z = (mean - x) / (sigma / sqrt(n))
4217        let z = (mean - x) / (sigma / n.sqrt());
4218
4219        // P-value = 1 - NORM.S.DIST(z, TRUE)
4220        let p_value = 1.0 - std_norm_cdf(z);
4221
4222        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
4223            p_value,
4224        )))
4225    }
4226}
4227
4228/* ─────────────────────────── TREND ──────────────────────────── */
4229
4230/// TREND(known_y's, [known_x's], [new_x's], [const]) - Returns y values along a linear trend
4231/// Uses linear regression y = mx + b
4232/// - If new_x's provided, calculates trend values for those x's
4233/// - If new_x's omitted, uses known_x's
4234/// - const=TRUE (default): calculate intercept normally
4235/// - const=FALSE: force intercept through origin
4236#[derive(Debug)]
4237pub struct TrendFn;
4238impl Function for TrendFn {
4239    func_caps!(PURE, NUMERIC_ONLY);
4240    fn name(&self) -> &'static str {
4241        "TREND"
4242    }
4243    fn min_args(&self) -> usize {
4244        1
4245    }
4246    fn variadic(&self) -> bool {
4247        true
4248    }
4249    fn arg_schema(&self) -> &'static [ArgSchema] {
4250        &ARG_RANGE_NUM_LENIENT_ONE[..]
4251    }
4252    fn eval<'a, 'b, 'c>(
4253        &self,
4254        args: &'c [ArgumentHandle<'a, 'b>],
4255        _ctx: &dyn FunctionContext<'b>,
4256    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4257        // TREND: args[0] = known_y's (required)
4258        // args[1] = known_x's (optional, defaults to {1,2,3,...})
4259        // args[2] = new_x's (optional, defaults to known_x's)
4260        // args[3] = const (optional, default TRUE - whether to compute intercept)
4261
4262        let y_vals = collect_numeric_stats(&args[0..1])?;
4263
4264        if y_vals.is_empty() {
4265            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4266                ExcelError::new_na(),
4267            )));
4268        }
4269
4270        // Helper to check if argument is empty/omitted
4271        // Note: Empty arguments are represented as empty text strings by the parser
4272        fn is_arg_empty(arg: &ArgumentHandle) -> bool {
4273            match scalar_like_value(arg) {
4274                Ok(LiteralValue::Empty) => true,
4275                Ok(LiteralValue::Text(s)) if s.is_empty() => true,
4276                _ => false,
4277            }
4278        }
4279
4280        // Get known_x's or generate default {1, 2, 3, ...}
4281        let x_vals = if args.len() >= 2 && !is_arg_empty(&args[1]) {
4282            collect_numeric_stats(&args[1..2])?
4283        } else {
4284            (1..=y_vals.len()).map(|i| i as f64).collect()
4285        };
4286
4287        // Arrays must have same length
4288        if y_vals.len() != x_vals.len() {
4289            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4290                ExcelError::new_ref(),
4291            )));
4292        }
4293
4294        // Get new_x's or use known_x's - check if argument is empty/omitted
4295        let new_x_vals = if args.len() >= 3 && !is_arg_empty(&args[2]) {
4296            collect_numeric_stats(&args[2..3])?
4297        } else {
4298            x_vals.clone()
4299        };
4300
4301        if new_x_vals.is_empty() {
4302            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4303                ExcelError::new_na(),
4304            )));
4305        }
4306
4307        // Parse const argument (default TRUE)
4308        let use_const = if args.len() >= 4 {
4309            match scalar_like_value(&args[3])? {
4310                LiteralValue::Boolean(b) => b,
4311                LiteralValue::Number(n) => n != 0.0,
4312                LiteralValue::Int(i) => i != 0,
4313                LiteralValue::Empty => true, // empty defaults to TRUE
4314                _ => true,
4315            }
4316        } else {
4317            true
4318        };
4319
4320        let n = x_vals.len() as f64;
4321
4322        // Calculate regression coefficients
4323        let (slope, intercept) = if use_const {
4324            // Normal linear regression with intercept
4325            let mean_x = x_vals.iter().sum::<f64>() / n;
4326            let mean_y = y_vals.iter().sum::<f64>() / n;
4327
4328            let mut sum_xy = 0.0;
4329            let mut sum_x2 = 0.0;
4330
4331            for i in 0..x_vals.len() {
4332                let dx = x_vals[i] - mean_x;
4333                let dy = y_vals[i] - mean_y;
4334                sum_xy += dx * dy;
4335                sum_x2 += dx * dx;
4336            }
4337
4338            if sum_x2 == 0.0 {
4339                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4340                    ExcelError::new_div(),
4341                )));
4342            }
4343
4344            let slope = sum_xy / sum_x2;
4345            let intercept = mean_y - slope * mean_x;
4346            (slope, intercept)
4347        } else {
4348            // Regression through origin (intercept = 0)
4349            let mut sum_xy = 0.0;
4350            let mut sum_x2 = 0.0;
4351
4352            for i in 0..x_vals.len() {
4353                sum_xy += x_vals[i] * y_vals[i];
4354                sum_x2 += x_vals[i] * x_vals[i];
4355            }
4356
4357            if sum_x2 == 0.0 {
4358                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4359                    ExcelError::new_div(),
4360                )));
4361            }
4362
4363            let slope = sum_xy / sum_x2;
4364            (slope, 0.0)
4365        };
4366
4367        // Calculate predicted y values for new_x's
4368        let predicted: Vec<LiteralValue> = new_x_vals
4369            .iter()
4370            .map(|&x| LiteralValue::Number(slope * x + intercept))
4371            .collect();
4372
4373        // Return as 1xN array (row vector)
4374        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(vec![
4375            predicted,
4376        ])))
4377    }
4378}
4379
4380/* ─────────────────────────── GROWTH ──────────────────────────── */
4381
4382/// GROWTH(known_y's, [known_x's], [new_x's], [const]) - Returns values along exponential growth trend
4383/// Uses exponential regression y = b * m^x
4384/// - Similar parameters to TREND but for exponential growth
4385/// - const=TRUE: calculate b normally
4386/// - const=FALSE: force b = 1
4387#[derive(Debug)]
4388pub struct GrowthFn;
4389impl Function for GrowthFn {
4390    func_caps!(PURE, NUMERIC_ONLY);
4391    fn name(&self) -> &'static str {
4392        "GROWTH"
4393    }
4394    fn min_args(&self) -> usize {
4395        1
4396    }
4397    fn variadic(&self) -> bool {
4398        true
4399    }
4400    fn arg_schema(&self) -> &'static [ArgSchema] {
4401        &ARG_RANGE_NUM_LENIENT_ONE[..]
4402    }
4403    fn eval<'a, 'b, 'c>(
4404        &self,
4405        args: &'c [ArgumentHandle<'a, 'b>],
4406        _ctx: &dyn FunctionContext<'b>,
4407    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4408        // GROWTH: args[0] = known_y's (required)
4409        // args[1] = known_x's (optional, defaults to {1,2,3,...})
4410        // args[2] = new_x's (optional, defaults to known_x's)
4411        // args[3] = const (optional, default TRUE - whether to compute intercept)
4412
4413        let y_vals = collect_numeric_stats(&args[0..1])?;
4414
4415        if y_vals.is_empty() {
4416            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4417                ExcelError::new_na(),
4418            )));
4419        }
4420
4421        // Check that all y values are positive (required for log transformation)
4422        for &y in &y_vals {
4423            if y <= 0.0 {
4424                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4425                    ExcelError::new_num(),
4426                )));
4427            }
4428        }
4429
4430        // Helper to check if argument is empty/omitted
4431        // Note: Empty arguments are represented as empty text strings by the parser
4432        fn is_arg_empty(arg: &ArgumentHandle) -> bool {
4433            match scalar_like_value(arg) {
4434                Ok(LiteralValue::Empty) => true,
4435                Ok(LiteralValue::Text(s)) if s.is_empty() => true,
4436                _ => false,
4437            }
4438        }
4439
4440        // Get known_x's or generate default {1, 2, 3, ...}
4441        let x_vals = if args.len() >= 2 && !is_arg_empty(&args[1]) {
4442            collect_numeric_stats(&args[1..2])?
4443        } else {
4444            (1..=y_vals.len()).map(|i| i as f64).collect()
4445        };
4446
4447        // Arrays must have same length
4448        if y_vals.len() != x_vals.len() {
4449            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4450                ExcelError::new_ref(),
4451            )));
4452        }
4453
4454        // Get new_x's or use known_x's - check if argument is empty/omitted
4455        let new_x_vals = if args.len() >= 3 && !is_arg_empty(&args[2]) {
4456            collect_numeric_stats(&args[2..3])?
4457        } else {
4458            x_vals.clone()
4459        };
4460
4461        if new_x_vals.is_empty() {
4462            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4463                ExcelError::new_na(),
4464            )));
4465        }
4466
4467        // Parse const argument (default TRUE)
4468        let use_const = if args.len() >= 4 {
4469            match scalar_like_value(&args[3])? {
4470                LiteralValue::Boolean(b) => b,
4471                LiteralValue::Number(n) => n != 0.0,
4472                LiteralValue::Int(i) => i != 0,
4473                LiteralValue::Empty => true, // empty defaults to TRUE
4474                _ => true,
4475            }
4476        } else {
4477            true
4478        };
4479
4480        // Transform to log space: ln(y) = ln(b) + x*ln(m)
4481        // This is linear regression on log-transformed y values
4482        let ln_y_vals: Vec<f64> = y_vals.iter().map(|&y| y.ln()).collect();
4483
4484        let n = x_vals.len() as f64;
4485
4486        // Calculate regression coefficients in log space
4487        let (ln_m, ln_b) = if use_const {
4488            // Normal linear regression with intercept
4489            let mean_x = x_vals.iter().sum::<f64>() / n;
4490            let mean_ln_y = ln_y_vals.iter().sum::<f64>() / n;
4491
4492            let mut sum_xy = 0.0;
4493            let mut sum_x2 = 0.0;
4494
4495            for i in 0..x_vals.len() {
4496                let dx = x_vals[i] - mean_x;
4497                let dy = ln_y_vals[i] - mean_ln_y;
4498                sum_xy += dx * dy;
4499                sum_x2 += dx * dx;
4500            }
4501
4502            if sum_x2 == 0.0 {
4503                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4504                    ExcelError::new_div(),
4505                )));
4506            }
4507
4508            let ln_m = sum_xy / sum_x2;
4509            let ln_b = mean_ln_y - ln_m * mean_x;
4510            (ln_m, ln_b)
4511        } else {
4512            // Regression through origin in log space (ln_b = 0, so b = 1)
4513            let mut sum_xy = 0.0;
4514            let mut sum_x2 = 0.0;
4515
4516            for i in 0..x_vals.len() {
4517                sum_xy += x_vals[i] * ln_y_vals[i];
4518                sum_x2 += x_vals[i] * x_vals[i];
4519            }
4520
4521            if sum_x2 == 0.0 {
4522                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4523                    ExcelError::new_div(),
4524                )));
4525            }
4526
4527            let ln_m = sum_xy / sum_x2;
4528            (ln_m, 0.0)
4529        };
4530
4531        // Convert back from log space: m = e^ln_m, b = e^ln_b
4532        let m = ln_m.exp();
4533        let b = ln_b.exp();
4534
4535        // Calculate predicted y values: y = b * m^x
4536        let predicted: Vec<LiteralValue> = new_x_vals
4537            .iter()
4538            .map(|&x| LiteralValue::Number(b * m.powf(x)))
4539            .collect();
4540
4541        // Return as 1xN array (row vector)
4542        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(vec![
4543            predicted,
4544        ])))
4545    }
4546}
4547
4548/* ─────────────────────────── LOGEST ──────────────────────────── */
4549
4550/// LOGEST(known_y's, [known_x's], [const], [stats]) - Returns parameters of exponential curve
4551/// Returns array: [[m, b]] when stats=FALSE
4552/// Returns 5x2 array with statistics when stats=TRUE (like LINEST)
4553/// The exponential curve is y = b * m^x
4554#[derive(Debug)]
4555pub struct LogestFn;
4556impl Function for LogestFn {
4557    func_caps!(PURE, NUMERIC_ONLY);
4558    fn name(&self) -> &'static str {
4559        "LOGEST"
4560    }
4561    fn min_args(&self) -> usize {
4562        1
4563    }
4564    fn variadic(&self) -> bool {
4565        true
4566    }
4567    fn arg_schema(&self) -> &'static [ArgSchema] {
4568        &ARG_RANGE_NUM_LENIENT_ONE[..]
4569    }
4570    fn eval<'a, 'b, 'c>(
4571        &self,
4572        args: &'c [ArgumentHandle<'a, 'b>],
4573        _ctx: &dyn FunctionContext<'b>,
4574    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4575        // args[0] = known_y's (required)
4576        // args[1] = known_x's (optional, defaults to {1,2,3,...})
4577        // args[2] = const (optional, default TRUE - whether to compute b)
4578        // args[3] = stats (optional, default FALSE - whether to return additional statistics)
4579
4580        let y_vals = collect_numeric_stats(&args[0..1])?;
4581
4582        if y_vals.is_empty() {
4583            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4584                ExcelError::new_na(),
4585            )));
4586        }
4587
4588        // Check that all y values are positive (required for log transformation)
4589        for &y in &y_vals {
4590            if y <= 0.0 {
4591                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4592                    ExcelError::new_num(),
4593                )));
4594            }
4595        }
4596
4597        // Get known_x's or generate default {1, 2, 3, ...}
4598        let x_vals = if args.len() >= 2 {
4599            collect_numeric_stats(&args[1..2])?
4600        } else {
4601            (1..=y_vals.len()).map(|i| i as f64).collect()
4602        };
4603
4604        // Arrays must have same length
4605        if y_vals.len() != x_vals.len() {
4606            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4607                ExcelError::new_ref(),
4608            )));
4609        }
4610
4611        // Parse const argument (default TRUE)
4612        let use_const = if args.len() >= 3 {
4613            match scalar_like_value(&args[2])? {
4614                LiteralValue::Boolean(b) => b,
4615                LiteralValue::Number(n) => n != 0.0,
4616                LiteralValue::Int(i) => i != 0,
4617                _ => true,
4618            }
4619        } else {
4620            true
4621        };
4622
4623        // Parse stats argument (default FALSE)
4624        let return_stats = if args.len() >= 4 {
4625            match scalar_like_value(&args[3])? {
4626                LiteralValue::Boolean(b) => b,
4627                LiteralValue::Number(n) => n != 0.0,
4628                LiteralValue::Int(i) => i != 0,
4629                _ => false,
4630            }
4631        } else {
4632            false
4633        };
4634
4635        // Transform to log space: ln(y) = ln(b) + x*ln(m)
4636        let ln_y_vals: Vec<f64> = y_vals.iter().map(|&y| y.ln()).collect();
4637
4638        let n = x_vals.len() as f64;
4639
4640        // Calculate regression coefficients in log space
4641        let (ln_m, ln_b) = if use_const {
4642            // Normal linear regression with intercept
4643            let mean_x = x_vals.iter().sum::<f64>() / n;
4644            let mean_ln_y = ln_y_vals.iter().sum::<f64>() / n;
4645
4646            let mut sum_xy = 0.0;
4647            let mut sum_x2 = 0.0;
4648
4649            for i in 0..x_vals.len() {
4650                let dx = x_vals[i] - mean_x;
4651                let dy = ln_y_vals[i] - mean_ln_y;
4652                sum_xy += dx * dy;
4653                sum_x2 += dx * dx;
4654            }
4655
4656            if sum_x2 == 0.0 {
4657                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4658                    ExcelError::new_div(),
4659                )));
4660            }
4661
4662            let ln_m = sum_xy / sum_x2;
4663            let ln_b = mean_ln_y - ln_m * mean_x;
4664            (ln_m, ln_b)
4665        } else {
4666            // Regression through origin in log space (ln_b = 0, so b = 1)
4667            let mut sum_xy = 0.0;
4668            let mut sum_x2 = 0.0;
4669
4670            for i in 0..x_vals.len() {
4671                sum_xy += x_vals[i] * ln_y_vals[i];
4672                sum_x2 += x_vals[i] * x_vals[i];
4673            }
4674
4675            if sum_x2 == 0.0 {
4676                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4677                    ExcelError::new_div(),
4678                )));
4679            }
4680
4681            let ln_m = sum_xy / sum_x2;
4682            (ln_m, 0.0)
4683        };
4684
4685        // Convert from log space to get m and b
4686        let m = ln_m.exp();
4687        let b = ln_b.exp();
4688
4689        if !return_stats {
4690            // Return just m and b as 1x2 array: [[m, b]]
4691            let row = vec![LiteralValue::Number(m), LiteralValue::Number(b)];
4692            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(vec![
4693                row,
4694            ])));
4695        }
4696
4697        // Calculate additional statistics for stats=TRUE
4698        // Statistics are computed in log space, then converted
4699        // Row 1: [m, b]
4700        // Row 2: [se_m, se_b] - standard errors (converted from log space)
4701        // Row 3: [r_squared, se_y] - R-squared and standard error of y estimate
4702        // Row 4: [F_statistic, df] - F-statistic and degrees of freedom
4703        // Row 5: [ss_reg, ss_resid] - regression sum of squares and residual sum of squares
4704
4705        let mean_ln_y = ln_y_vals.iter().sum::<f64>() / n;
4706
4707        // Calculate residuals and sums of squares in log space
4708        let mut ss_resid = 0.0;
4709        let mut ss_tot = 0.0;
4710
4711        for i in 0..x_vals.len() {
4712            let ln_y_pred = ln_m * x_vals[i] + ln_b;
4713            let residual = ln_y_vals[i] - ln_y_pred;
4714            ss_resid += residual * residual;
4715            let dy_tot = ln_y_vals[i] - mean_ln_y;
4716            ss_tot += dy_tot * dy_tot;
4717        }
4718
4719        let ss_reg = ss_tot - ss_resid;
4720
4721        // R-squared (same in both spaces for transformed regression)
4722        let r_squared = if ss_tot == 0.0 {
4723            1.0
4724        } else {
4725            1.0 - (ss_resid / ss_tot)
4726        };
4727
4728        // Degrees of freedom
4729        let df = if use_const {
4730            (n as i64 - 2).max(1) as f64
4731        } else {
4732            (n as i64 - 1).max(1) as f64
4733        };
4734
4735        // Standard error of y estimate (in log space)
4736        let se_ln_y = if df > 0.0 {
4737            (ss_resid / df).sqrt()
4738        } else {
4739            0.0
4740        };
4741
4742        // Standard errors of coefficients in log space
4743        let mean_x = x_vals.iter().sum::<f64>() / n;
4744        let mut sum_x2_centered = 0.0;
4745        let mut sum_x2_raw = 0.0;
4746        for &xi in &x_vals {
4747            sum_x2_centered += (xi - mean_x).powi(2);
4748            sum_x2_raw += xi * xi;
4749        }
4750
4751        let se_ln_m = if sum_x2_centered > 0.0 && df > 0.0 {
4752            se_ln_y / sum_x2_centered.sqrt()
4753        } else {
4754            f64::NAN
4755        };
4756
4757        let se_ln_b = if use_const && sum_x2_centered > 0.0 && df > 0.0 {
4758            se_ln_y * (sum_x2_raw / (n * sum_x2_centered)).sqrt()
4759        } else {
4760            f64::NAN
4761        };
4762
4763        // Convert standard errors: se_m = m * se_ln_m (delta method approximation)
4764        let se_m = m * se_ln_m;
4765        let se_b = b * se_ln_b;
4766
4767        // Standard error of y estimate - convert from log space
4768        // This is an approximation; for exponential models, se_y in original space varies with x
4769        let se_y = se_ln_y;
4770
4771        // F-statistic
4772        let f_stat = if ss_resid > 0.0 && df > 0.0 {
4773            (ss_reg / 1.0) / (ss_resid / df)
4774        } else if ss_resid == 0.0 {
4775            f64::INFINITY
4776        } else {
4777            f64::NAN
4778        };
4779
4780        // Build 5x2 result array
4781        let rows = vec![
4782            vec![LiteralValue::Number(m), LiteralValue::Number(b)],
4783            vec![LiteralValue::Number(se_m), LiteralValue::Number(se_b)],
4784            vec![LiteralValue::Number(r_squared), LiteralValue::Number(se_y)],
4785            vec![LiteralValue::Number(f_stat), LiteralValue::Number(df)],
4786            vec![LiteralValue::Number(ss_reg), LiteralValue::Number(ss_resid)],
4787        ];
4788
4789        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(rows)))
4790    }
4791}
4792
4793/* ─────────────────────────── PERCENTRANK ──────────────────────────── */
4794
4795/// PERCENTRANK.INC(array, x, [significance]) - Returns percentage rank (inclusive)
4796/// Returns rank of x in array as percentage (0 to 1 inclusive)
4797/// Uses interpolation for values between data points
4798/// significance: number of significant digits (default 3)
4799#[derive(Debug)]
4800pub struct PercentRankIncFn;
4801impl Function for PercentRankIncFn {
4802    func_caps!(PURE, NUMERIC_ONLY);
4803    fn name(&self) -> &'static str {
4804        "PERCENTRANK.INC"
4805    }
4806    fn aliases(&self) -> &'static [&'static str] {
4807        &["PERCENTRANK"]
4808    }
4809    fn min_args(&self) -> usize {
4810        2
4811    }
4812    fn variadic(&self) -> bool {
4813        true
4814    }
4815    fn arg_schema(&self) -> &'static [ArgSchema] {
4816        &ARG_RANGE_NUM_LENIENT_ONE[..]
4817    }
4818    fn eval<'a, 'b, 'c>(
4819        &self,
4820        args: &'c [ArgumentHandle<'a, 'b>],
4821        _ctx: &dyn FunctionContext<'b>,
4822    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4823        if args.len() < 2 {
4824            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4825                ExcelError::new_num(),
4826            )));
4827        }
4828
4829        // Get x value (the value to find the rank of)
4830        let x = match coerce_num(&scalar_like_value(&args[1])?) {
4831            Ok(n) => n,
4832            Err(_) => {
4833                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4834                    ExcelError::new_num(),
4835                )));
4836            }
4837        };
4838
4839        // Get optional significance (default 3)
4840        let significance = if args.len() > 2 {
4841            match coerce_num(&scalar_like_value(&args[2])?) {
4842                Ok(n) => {
4843                    let s = n as i32;
4844                    if s < 1 {
4845                        return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4846                            ExcelError::new_num(),
4847                        )));
4848                    }
4849                    s as u32
4850                }
4851                Err(_) => 3,
4852            }
4853        } else {
4854            3
4855        };
4856
4857        // Collect and sort the data array
4858        let mut nums = collect_numeric_stats(&args[0..1])?;
4859        if nums.is_empty() {
4860            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4861                ExcelError::new_num(),
4862            )));
4863        }
4864        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
4865
4866        let n = nums.len();
4867
4868        // Check if x is outside the range
4869        if x < nums[0] || x > nums[n - 1] {
4870            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4871                ExcelError::new_na(),
4872            )));
4873        }
4874
4875        // Find the rank using linear interpolation
4876        // For PERCENTRANK.INC, the formula is: rank = (position) / (n-1)
4877        // where position is 0-based and uses linear interpolation
4878        let rank = if n == 1 {
4879            // Single element - rank is 0 (or 1.0 if we want, but Excel returns 0)
4880            0.0
4881        } else {
4882            let mut rank_val = 0.0;
4883            for i in 0..n - 1 {
4884                if (nums[i] - x).abs() < 1e-12 {
4885                    // Exact match at position i
4886                    rank_val = (i as f64) / ((n - 1) as f64);
4887                    break;
4888                } else if nums[i] < x && x < nums[i + 1] {
4889                    // Interpolate between positions i and i+1
4890                    let frac = (x - nums[i]) / (nums[i + 1] - nums[i]);
4891                    rank_val = ((i as f64) + frac) / ((n - 1) as f64);
4892                    break;
4893                } else if i == n - 2 && (nums[n - 1] - x).abs() < 1e-12 {
4894                    // Exact match at last position
4895                    rank_val = 1.0;
4896                }
4897            }
4898            rank_val
4899        };
4900
4901        // Truncate to significance decimal places
4902        let multiplier = 10_f64.powi(significance as i32);
4903        let truncated = (rank * multiplier).trunc() / multiplier;
4904
4905        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
4906            truncated,
4907        )))
4908    }
4909}
4910
4911/// PERCENTRANK.EXC(array, x, [significance]) - Returns percentage rank (exclusive)
4912/// Same as PERCENTRANK.INC but excludes 0 and 1 from range
4913/// Range is 1/(n+1) to n/(n+1)
4914#[derive(Debug)]
4915pub struct PercentRankExcFn;
4916impl Function for PercentRankExcFn {
4917    func_caps!(PURE, NUMERIC_ONLY);
4918    fn name(&self) -> &'static str {
4919        "PERCENTRANK.EXC"
4920    }
4921    fn min_args(&self) -> usize {
4922        2
4923    }
4924    fn variadic(&self) -> bool {
4925        true
4926    }
4927    fn arg_schema(&self) -> &'static [ArgSchema] {
4928        &ARG_RANGE_NUM_LENIENT_ONE[..]
4929    }
4930    fn eval<'a, 'b, 'c>(
4931        &self,
4932        args: &'c [ArgumentHandle<'a, 'b>],
4933        _ctx: &dyn FunctionContext<'b>,
4934    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
4935        if args.len() < 2 {
4936            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4937                ExcelError::new_num(),
4938            )));
4939        }
4940
4941        // Get x value (the value to find the rank of)
4942        let x = match coerce_num(&scalar_like_value(&args[1])?) {
4943            Ok(n) => n,
4944            Err(_) => {
4945                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4946                    ExcelError::new_num(),
4947                )));
4948            }
4949        };
4950
4951        // Get optional significance (default 3)
4952        let significance = if args.len() > 2 {
4953            match coerce_num(&scalar_like_value(&args[2])?) {
4954                Ok(n) => {
4955                    let s = n as i32;
4956                    if s < 1 {
4957                        return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4958                            ExcelError::new_num(),
4959                        )));
4960                    }
4961                    s as u32
4962                }
4963                Err(_) => 3,
4964            }
4965        } else {
4966            3
4967        };
4968
4969        // Collect and sort the data array
4970        let mut nums = collect_numeric_stats(&args[0..1])?;
4971        if nums.is_empty() {
4972            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4973                ExcelError::new_num(),
4974            )));
4975        }
4976        nums.sort_by(|a, b| a.partial_cmp(b).unwrap());
4977
4978        let n = nums.len();
4979
4980        // Check if x is outside the range
4981        if x < nums[0] || x > nums[n - 1] {
4982            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
4983                ExcelError::new_na(),
4984            )));
4985        }
4986
4987        // For PERCENTRANK.EXC, the formula is: rank = position / (n+1)
4988        // where position is 1-based and uses linear interpolation
4989        let rank = {
4990            let mut rank_val = 0.0;
4991            for i in 0..n {
4992                if (nums[i] - x).abs() < 1e-12 {
4993                    // Exact match at position i (1-based: i+1)
4994                    rank_val = ((i + 1) as f64) / ((n + 1) as f64);
4995                    break;
4996                } else if i < n - 1 && nums[i] < x && x < nums[i + 1] {
4997                    // Interpolate between positions i and i+1 (1-based: i+1 and i+2)
4998                    let frac = (x - nums[i]) / (nums[i + 1] - nums[i]);
4999                    let position = ((i + 1) as f64) + frac;
5000                    rank_val = position / ((n + 1) as f64);
5001                    break;
5002                }
5003            }
5004            rank_val
5005        };
5006
5007        // Truncate to significance decimal places
5008        let multiplier = 10_f64.powi(significance as i32);
5009        let truncated = (rank * multiplier).trunc() / multiplier;
5010
5011        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
5012            truncated,
5013        )))
5014    }
5015}
5016
5017/* ─────────────────────────── FREQUENCY ──────────────────────────── */
5018
5019/// FREQUENCY(data_array, bins_array) - Returns frequency distribution
5020/// Returns vertical array of frequencies
5021/// Counts values in each bin: <= bin[0], (bin[0], bin[1]], ..., > bin[n-1]
5022/// Returns array with one more element than bins_array
5023#[derive(Debug)]
5024pub struct FrequencyFn;
5025impl Function for FrequencyFn {
5026    func_caps!(PURE, NUMERIC_ONLY);
5027    fn name(&self) -> &'static str {
5028        "FREQUENCY"
5029    }
5030    fn min_args(&self) -> usize {
5031        2
5032    }
5033    fn variadic(&self) -> bool {
5034        false
5035    }
5036    fn arg_schema(&self) -> &'static [ArgSchema] {
5037        &ARG_RANGE_NUM_LENIENT_ONE[..]
5038    }
5039    fn eval<'a, 'b, 'c>(
5040        &self,
5041        args: &'c [ArgumentHandle<'a, 'b>],
5042        _ctx: &dyn FunctionContext<'b>,
5043    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5044        if args.len() < 2 {
5045            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5046                ExcelError::new_num(),
5047            )));
5048        }
5049
5050        // Collect data array
5051        let data = collect_numeric_stats(&args[0..1])?;
5052
5053        // Collect bins array
5054        let mut bins = collect_numeric_stats(&args[1..2])?;
5055
5056        // Handle empty bins - return single count of all data
5057        if bins.is_empty() {
5058            let rows = vec![vec![LiteralValue::Number(data.len() as f64)]];
5059            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(rows)));
5060        }
5061
5062        // Sort bins
5063        bins.sort_by(|a, b| a.partial_cmp(b).unwrap());
5064
5065        // Calculate frequencies
5066        // Result has bins.len() + 1 elements
5067        let mut frequencies = vec![0usize; bins.len() + 1];
5068
5069        for &value in &data {
5070            // Find which bin the value belongs to
5071            let mut found = false;
5072            for (i, &bin) in bins.iter().enumerate() {
5073                if i == 0 {
5074                    // First bin: count values <= bins[0]
5075                    if value <= bin {
5076                        frequencies[0] += 1;
5077                        found = true;
5078                        break;
5079                    }
5080                } else {
5081                    // Intermediate bins: count values > bins[i-1] AND <= bins[i]
5082                    if value <= bin {
5083                        frequencies[i] += 1;
5084                        found = true;
5085                        break;
5086                    }
5087                }
5088            }
5089            // Last bin: values > bins[last]
5090            if !found {
5091                frequencies[bins.len()] += 1;
5092            }
5093        }
5094
5095        // Return as vertical array (column vector)
5096        let rows: Vec<Vec<LiteralValue>> = frequencies
5097            .into_iter()
5098            .map(|f| vec![LiteralValue::Number(f as f64)])
5099            .collect();
5100
5101        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Array(rows)))
5102    }
5103}
5104
5105/* ─────────────────────────── T.DIST.2T ──────────────────────────── */
5106
5107/// T.DIST.2T(x, deg_freedom) - Returns the two-tailed Student's t-distribution
5108/// Returns P(|T| > x) = 2 * (1 - t_cdf(|x|, df))
5109#[derive(Debug)]
5110pub struct TDist2TFn;
5111impl Function for TDist2TFn {
5112    func_caps!(PURE);
5113    fn name(&self) -> &'static str {
5114        "T.DIST.2T"
5115    }
5116    fn min_args(&self) -> usize {
5117        2
5118    }
5119    fn arg_schema(&self) -> &'static [ArgSchema] {
5120        use std::sync::LazyLock;
5121        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
5122            vec![
5123                ArgSchema::number_lenient_scalar(),
5124                ArgSchema::number_lenient_scalar(),
5125            ]
5126        });
5127        &SCHEMA[..]
5128    }
5129    fn eval<'a, 'b, 'c>(
5130        &self,
5131        args: &'c [ArgumentHandle<'a, 'b>],
5132        _ctx: &dyn FunctionContext<'b>,
5133    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5134        let x = coerce_num(&scalar_like_value(&args[0])?)?;
5135        let df = coerce_num(&scalar_like_value(&args[1])?)?;
5136
5137        // x must be non-negative for T.DIST.2T, df must be >= 1
5138        if x < 0.0 || df < 1.0 {
5139            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5140                ExcelError::new_num(),
5141            )));
5142        }
5143
5144        // Two-tailed: P(|T| > x) = 2 * (1 - t_cdf(x, df))
5145        let p = 2.0 * (1.0 - t_cdf(x, df));
5146        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(p)))
5147    }
5148}
5149
5150/* ─────────────────────────── T.INV.2T ──────────────────────────── */
5151
5152/// T.INV.2T(probability, deg_freedom) - Returns the two-tailed inverse of Student's t-distribution
5153/// Returns the value t such that P(|T| > t) = probability
5154/// This is equivalent to t_inv(1 - probability/2, df)
5155#[derive(Debug)]
5156pub struct TInv2TFn;
5157impl Function for TInv2TFn {
5158    func_caps!(PURE);
5159    fn name(&self) -> &'static str {
5160        "T.INV.2T"
5161    }
5162    fn aliases(&self) -> &'static [&'static str] {
5163        &["TINV"]
5164    }
5165    fn min_args(&self) -> usize {
5166        2
5167    }
5168    fn arg_schema(&self) -> &'static [ArgSchema] {
5169        use std::sync::LazyLock;
5170        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
5171            vec![
5172                ArgSchema::number_lenient_scalar(),
5173                ArgSchema::number_lenient_scalar(),
5174            ]
5175        });
5176        &SCHEMA[..]
5177    }
5178    fn eval<'a, 'b, 'c>(
5179        &self,
5180        args: &'c [ArgumentHandle<'a, 'b>],
5181        _ctx: &dyn FunctionContext<'b>,
5182    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5183        let p = coerce_num(&scalar_like_value(&args[0])?)?;
5184        let df = coerce_num(&scalar_like_value(&args[1])?)?;
5185
5186        // probability must be in (0, 1], df >= 1
5187        if p <= 0.0 || p > 1.0 || df < 1.0 {
5188            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5189                ExcelError::new_num(),
5190            )));
5191        }
5192
5193        // For two-tailed: we want t such that P(|T| > t) = p
5194        // P(|T| > t) = 2 * (1 - F(t)) where F is CDF
5195        // So 1 - F(t) = p/2, meaning F(t) = 1 - p/2
5196        // Thus t = t_inv(1 - p/2, df)
5197        match t_inv(1.0 - p / 2.0, df) {
5198            Some(result) => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(
5199                result,
5200            ))),
5201            None => Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5202                ExcelError::new_num(),
5203            ))),
5204        }
5205    }
5206}
5207
5208/* ─────────────────────────── T.TEST ──────────────────────────── */
5209
5210/// T.TEST(array1, array2, tails, type) - Returns probability for Student's t-test
5211/// tails: 1 for one-tailed, 2 for two-tailed
5212/// type: 1=paired, 2=two-sample equal variance, 3=two-sample unequal variance (Welch's)
5213#[derive(Debug)]
5214pub struct TTestFn;
5215impl Function for TTestFn {
5216    func_caps!(PURE);
5217    fn name(&self) -> &'static str {
5218        "T.TEST"
5219    }
5220    fn aliases(&self) -> &'static [&'static str] {
5221        &["TTEST"]
5222    }
5223    fn min_args(&self) -> usize {
5224        4
5225    }
5226    fn arg_schema(&self) -> &'static [ArgSchema] {
5227        use std::sync::LazyLock;
5228        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
5229            vec![
5230                {
5231                    let mut s = ArgSchema::number_lenient_scalar();
5232                    s.shape = crate::args::ShapeKind::Range;
5233                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5234                    s
5235                },
5236                {
5237                    let mut s = ArgSchema::number_lenient_scalar();
5238                    s.shape = crate::args::ShapeKind::Range;
5239                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5240                    s
5241                },
5242                ArgSchema::number_lenient_scalar(), // tails
5243                ArgSchema::number_lenient_scalar(), // type
5244            ]
5245        });
5246        &SCHEMA[..]
5247    }
5248    fn eval<'a, 'b, 'c>(
5249        &self,
5250        args: &'c [ArgumentHandle<'a, 'b>],
5251        _ctx: &dyn FunctionContext<'b>,
5252    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5253        let array1 = collect_numeric_stats(&args[0..1])?;
5254        let array2 = collect_numeric_stats(&args[1..2])?;
5255        let tails = coerce_num(&scalar_like_value(&args[2])?)? as i32;
5256        let test_type = coerce_num(&scalar_like_value(&args[3])?)? as i32;
5257
5258        // Validate tails (1 or 2) and type (1, 2, or 3)
5259        if !(1..=2).contains(&tails) || !(1..=3).contains(&test_type) {
5260            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5261                ExcelError::new_num(),
5262            )));
5263        }
5264
5265        let n1 = array1.len();
5266        let n2 = array2.len();
5267
5268        // For paired test, arrays must have same length
5269        if test_type == 1 && n1 != n2 {
5270            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5271                ExcelError::new_na(),
5272            )));
5273        }
5274
5275        // Need at least 2 data points for meaningful t-test
5276        if n1 < 2 || n2 < 2 {
5277            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5278                ExcelError::new_num(),
5279            )));
5280        }
5281
5282        let (t_stat, df) = match test_type {
5283            1 => {
5284                // Paired t-test
5285                let n = n1 as f64;
5286                let diffs: Vec<f64> = array1
5287                    .iter()
5288                    .zip(array2.iter())
5289                    .map(|(a, b)| a - b)
5290                    .collect();
5291                let mean_diff = diffs.iter().sum::<f64>() / n;
5292                let var_diff =
5293                    diffs.iter().map(|d| (d - mean_diff).powi(2)).sum::<f64>() / (n - 1.0);
5294                if var_diff == 0.0 {
5295                    return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5296                        ExcelError::new_div(),
5297                    )));
5298                }
5299                let se = (var_diff / n).sqrt();
5300                (mean_diff / se, n - 1.0)
5301            }
5302            2 => {
5303                // Two-sample equal variance (pooled)
5304                let n1f = n1 as f64;
5305                let n2f = n2 as f64;
5306                let mean1 = array1.iter().sum::<f64>() / n1f;
5307                let mean2 = array2.iter().sum::<f64>() / n2f;
5308                let var1 = array1.iter().map(|x| (x - mean1).powi(2)).sum::<f64>() / (n1f - 1.0);
5309                let var2 = array2.iter().map(|x| (x - mean2).powi(2)).sum::<f64>() / (n2f - 1.0);
5310
5311                // Pooled variance
5312                let sp2 = ((n1f - 1.0) * var1 + (n2f - 1.0) * var2) / (n1f + n2f - 2.0);
5313                if sp2 == 0.0 {
5314                    return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5315                        ExcelError::new_div(),
5316                    )));
5317                }
5318                let se = (sp2 * (1.0 / n1f + 1.0 / n2f)).sqrt();
5319                ((mean1 - mean2) / se, n1f + n2f - 2.0)
5320            }
5321            3 => {
5322                // Welch's t-test (unequal variance)
5323                let n1f = n1 as f64;
5324                let n2f = n2 as f64;
5325                let mean1 = array1.iter().sum::<f64>() / n1f;
5326                let mean2 = array2.iter().sum::<f64>() / n2f;
5327                let var1 = array1.iter().map(|x| (x - mean1).powi(2)).sum::<f64>() / (n1f - 1.0);
5328                let var2 = array2.iter().map(|x| (x - mean2).powi(2)).sum::<f64>() / (n2f - 1.0);
5329
5330                let s1_n = var1 / n1f;
5331                let s2_n = var2 / n2f;
5332                let se = (s1_n + s2_n).sqrt();
5333                if se == 0.0 {
5334                    return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5335                        ExcelError::new_div(),
5336                    )));
5337                }
5338
5339                // Welch-Satterthwaite degrees of freedom
5340                let df_num = (s1_n + s2_n).powi(2);
5341                let df_denom = s1_n.powi(2) / (n1f - 1.0) + s2_n.powi(2) / (n2f - 1.0);
5342                let df = if df_denom == 0.0 {
5343                    1.0
5344                } else {
5345                    df_num / df_denom
5346                };
5347                ((mean1 - mean2) / se, df)
5348            }
5349            _ => unreachable!(),
5350        };
5351
5352        // Calculate p-value based on tails
5353        let t_abs = t_stat.abs();
5354        let p = if tails == 1 {
5355            1.0 - t_cdf(t_abs, df)
5356        } else {
5357            2.0 * (1.0 - t_cdf(t_abs, df))
5358        };
5359
5360        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(p)))
5361    }
5362}
5363
5364/* ─────────────────────────── F.TEST ──────────────────────────── */
5365
5366/// F.TEST(array1, array2) - Returns result of F-test for comparing variances
5367/// Returns the two-tailed probability that variances are not significantly different
5368/// F = larger_variance / smaller_variance
5369#[derive(Debug)]
5370pub struct FTestFn;
5371impl Function for FTestFn {
5372    func_caps!(PURE);
5373    fn name(&self) -> &'static str {
5374        "F.TEST"
5375    }
5376    fn aliases(&self) -> &'static [&'static str] {
5377        &["FTEST"]
5378    }
5379    fn min_args(&self) -> usize {
5380        2
5381    }
5382    fn arg_schema(&self) -> &'static [ArgSchema] {
5383        use std::sync::LazyLock;
5384        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
5385            vec![
5386                {
5387                    let mut s = ArgSchema::number_lenient_scalar();
5388                    s.shape = crate::args::ShapeKind::Range;
5389                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5390                    s
5391                },
5392                {
5393                    let mut s = ArgSchema::number_lenient_scalar();
5394                    s.shape = crate::args::ShapeKind::Range;
5395                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5396                    s
5397                },
5398            ]
5399        });
5400        &SCHEMA[..]
5401    }
5402    fn eval<'a, 'b, 'c>(
5403        &self,
5404        args: &'c [ArgumentHandle<'a, 'b>],
5405        _ctx: &dyn FunctionContext<'b>,
5406    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5407        let array1 = collect_numeric_stats(&args[0..1])?;
5408        let array2 = collect_numeric_stats(&args[1..2])?;
5409
5410        let n1 = array1.len();
5411        let n2 = array2.len();
5412
5413        // Need at least 2 points in each array
5414        if n1 < 2 || n2 < 2 {
5415            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5416                ExcelError::new_div(),
5417            )));
5418        }
5419
5420        let n1f = n1 as f64;
5421        let n2f = n2 as f64;
5422
5423        let mean1 = array1.iter().sum::<f64>() / n1f;
5424        let mean2 = array2.iter().sum::<f64>() / n2f;
5425
5426        let var1 = array1.iter().map(|x| (x - mean1).powi(2)).sum::<f64>() / (n1f - 1.0);
5427        let var2 = array2.iter().map(|x| (x - mean2).powi(2)).sum::<f64>() / (n2f - 1.0);
5428
5429        // Handle zero variance
5430        if var1 == 0.0 || var2 == 0.0 {
5431            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5432                ExcelError::new_div(),
5433            )));
5434        }
5435
5436        // F-statistic: Excel's F.TEST uses var1/var2 (order matters for degrees of freedom)
5437        // and returns two-tailed p-value
5438        let f = var1 / var2;
5439        let df1 = n1f - 1.0;
5440        let df2 = n2f - 1.0;
5441
5442        // Two-tailed p-value: min(F.DIST(f), 1-F.DIST(f)) * 2
5443        let p_lower = f_cdf(f, df1, df2);
5444        let p_upper = 1.0 - p_lower;
5445        let p = 2.0 * p_lower.min(p_upper);
5446
5447        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(p)))
5448    }
5449}
5450
5451/* ─────────────────────────── CHISQ.TEST ──────────────────────────── */
5452
5453/// CHISQ.TEST(actual_range, expected_range) - Returns chi-squared test for independence
5454/// Returns p-value from chi-squared distribution
5455#[derive(Debug)]
5456pub struct ChisqTestFn;
5457impl Function for ChisqTestFn {
5458    func_caps!(PURE);
5459    fn name(&self) -> &'static str {
5460        "CHISQ.TEST"
5461    }
5462    fn aliases(&self) -> &'static [&'static str] {
5463        &["CHITEST"]
5464    }
5465    fn min_args(&self) -> usize {
5466        2
5467    }
5468    fn arg_schema(&self) -> &'static [ArgSchema] {
5469        use std::sync::LazyLock;
5470        static SCHEMA: LazyLock<Vec<ArgSchema>> = LazyLock::new(|| {
5471            vec![
5472                {
5473                    let mut s = ArgSchema::number_lenient_scalar();
5474                    s.shape = crate::args::ShapeKind::Range;
5475                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5476                    s
5477                },
5478                {
5479                    let mut s = ArgSchema::number_lenient_scalar();
5480                    s.shape = crate::args::ShapeKind::Range;
5481                    s.coercion = formualizer_common::CoercionPolicy::NumberLenientText;
5482                    s
5483                },
5484            ]
5485        });
5486        &SCHEMA[..]
5487    }
5488    fn eval<'a, 'b, 'c>(
5489        &self,
5490        args: &'c [ArgumentHandle<'a, 'b>],
5491        _ctx: &dyn FunctionContext<'b>,
5492    ) -> Result<crate::traits::CalcValue<'b>, ExcelError> {
5493        let actual = collect_numeric_stats(&args[0..1])?;
5494        let expected = collect_numeric_stats(&args[1..2])?;
5495
5496        // Arrays must have same length
5497        if actual.len() != expected.len() {
5498            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5499                ExcelError::new_na(),
5500            )));
5501        }
5502
5503        if actual.is_empty() {
5504            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5505                ExcelError::new_na(),
5506            )));
5507        }
5508
5509        // Calculate chi-squared statistic: sum((observed - expected)^2 / expected)
5510        let mut chi_sq = 0.0;
5511        for (obs, exp) in actual.iter().zip(expected.iter()) {
5512            if *exp <= 0.0 {
5513                return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5514                    ExcelError::new_num(),
5515                )));
5516            }
5517            chi_sq += (obs - exp).powi(2) / exp;
5518        }
5519
5520        // Degrees of freedom = number of categories - 1
5521        let df = (actual.len() - 1) as f64;
5522
5523        if df < 1.0 {
5524            return Ok(crate::traits::CalcValue::Scalar(LiteralValue::Error(
5525                ExcelError::new_num(),
5526            )));
5527        }
5528
5529        // P-value = 1 - CHISQ.DIST(chi_sq, df, TRUE) = right-tail probability
5530        let p = 1.0 - chisq_cdf(chi_sq, df);
5531
5532        Ok(crate::traits::CalcValue::Scalar(LiteralValue::Number(p)))
5533    }
5534}
5535
5536pub fn register_builtins() {
5537    use std::sync::Arc;
5538    crate::function_registry::register_function(Arc::new(ForecastLinearFn));
5539    crate::function_registry::register_function(Arc::new(LinestFn));
5540    crate::function_registry::register_function(Arc::new(LARGE));
5541    crate::function_registry::register_function(Arc::new(SMALL));
5542    crate::function_registry::register_function(Arc::new(MEDIAN));
5543    crate::function_registry::register_function(Arc::new(StdevSample));
5544    crate::function_registry::register_function(Arc::new(StdevPop));
5545    crate::function_registry::register_function(Arc::new(VarSample));
5546    crate::function_registry::register_function(Arc::new(VarPop));
5547    crate::function_registry::register_function(Arc::new(PercentileInc));
5548    crate::function_registry::register_function(Arc::new(PercentileExc));
5549    crate::function_registry::register_function(Arc::new(QuartileInc));
5550    crate::function_registry::register_function(Arc::new(QuartileExc));
5551    crate::function_registry::register_function(Arc::new(RankEqFn));
5552    crate::function_registry::register_function(Arc::new(RankAvgFn));
5553    crate::function_registry::register_function(Arc::new(ModeSingleFn));
5554    crate::function_registry::register_function(Arc::new(ModeMultiFn));
5555    crate::function_registry::register_function(Arc::new(ProductFn));
5556    crate::function_registry::register_function(Arc::new(GeomeanFn));
5557    crate::function_registry::register_function(Arc::new(HarmeanFn));
5558    crate::function_registry::register_function(Arc::new(AvedevFn));
5559    crate::function_registry::register_function(Arc::new(DevsqFn));
5560    crate::function_registry::register_function(Arc::new(MaxIfsFn));
5561    crate::function_registry::register_function(Arc::new(MinIfsFn));
5562    crate::function_registry::register_function(Arc::new(TrimmeanFn));
5563    crate::function_registry::register_function(Arc::new(CorrelFn));
5564    crate::function_registry::register_function(Arc::new(SlopeFn));
5565    crate::function_registry::register_function(Arc::new(InterceptFn));
5566    // Covariance and correlation functions
5567    crate::function_registry::register_function(Arc::new(CovariancePFn));
5568    crate::function_registry::register_function(Arc::new(CovarianceSFn));
5569    crate::function_registry::register_function(Arc::new(PearsonFn));
5570    crate::function_registry::register_function(Arc::new(RsqFn));
5571    crate::function_registry::register_function(Arc::new(SteyxFn));
5572    crate::function_registry::register_function(Arc::new(SkewFn));
5573    crate::function_registry::register_function(Arc::new(KurtFn));
5574    crate::function_registry::register_function(Arc::new(FisherFn));
5575    crate::function_registry::register_function(Arc::new(FisherInvFn));
5576    // Statistical distributions
5577    crate::function_registry::register_function(Arc::new(NormSDistFn));
5578    crate::function_registry::register_function(Arc::new(NormSInvFn));
5579    crate::function_registry::register_function(Arc::new(NormDistFn));
5580    crate::function_registry::register_function(Arc::new(NormInvFn));
5581    crate::function_registry::register_function(Arc::new(LognormDistFn));
5582    crate::function_registry::register_function(Arc::new(LognormInvFn));
5583    crate::function_registry::register_function(Arc::new(PhiFn));
5584    crate::function_registry::register_function(Arc::new(GaussFn));
5585    crate::function_registry::register_function(Arc::new(StandardizeFn));
5586    crate::function_registry::register_function(Arc::new(TDistFn));
5587    crate::function_registry::register_function(Arc::new(TInvFn));
5588    crate::function_registry::register_function(Arc::new(ChisqDistFn));
5589    crate::function_registry::register_function(Arc::new(ChisqInvFn));
5590    crate::function_registry::register_function(Arc::new(FDistFn));
5591    crate::function_registry::register_function(Arc::new(FInvFn));
5592    // Discrete distributions
5593    crate::function_registry::register_function(Arc::new(BinomDistFn));
5594    crate::function_registry::register_function(Arc::new(PoissonDistFn));
5595    crate::function_registry::register_function(Arc::new(ExponDistFn));
5596    crate::function_registry::register_function(Arc::new(GammaDistFn));
5597    // Additional distributions
5598    crate::function_registry::register_function(Arc::new(WeibullDistFn));
5599    crate::function_registry::register_function(Arc::new(BetaDistFn));
5600    crate::function_registry::register_function(Arc::new(NegbinomDistFn));
5601    crate::function_registry::register_function(Arc::new(HypgeomDistFn));
5602    // Confidence intervals and hypothesis testing
5603    crate::function_registry::register_function(Arc::new(ConfidenceNormFn));
5604    crate::function_registry::register_function(Arc::new(ConfidenceTFn));
5605    crate::function_registry::register_function(Arc::new(ZTestFn));
5606    // Regression and trend functions
5607    crate::function_registry::register_function(Arc::new(TrendFn));
5608    crate::function_registry::register_function(Arc::new(GrowthFn));
5609    crate::function_registry::register_function(Arc::new(LogestFn));
5610    // Percent rank and frequency functions
5611    crate::function_registry::register_function(Arc::new(PercentRankIncFn));
5612    crate::function_registry::register_function(Arc::new(PercentRankExcFn));
5613    crate::function_registry::register_function(Arc::new(FrequencyFn));
5614    // Hypothesis testing functions
5615    crate::function_registry::register_function(Arc::new(TDist2TFn));
5616    crate::function_registry::register_function(Arc::new(TInv2TFn));
5617    crate::function_registry::register_function(Arc::new(TTestFn));
5618    crate::function_registry::register_function(Arc::new(FTestFn));
5619    crate::function_registry::register_function(Arc::new(ChisqTestFn));
5620}
5621
5622#[cfg(test)]
5623mod tests_basic_stats {
5624    use super::*;
5625    use crate::test_workbook::TestWorkbook;
5626    use crate::traits::ArgumentHandle;
5627    use formualizer_common::LiteralValue;
5628    use formualizer_parse::parser::{ASTNode, ASTNodeType};
5629    fn interp(wb: &TestWorkbook) -> crate::interpreter::Interpreter<'_> {
5630        wb.interpreter()
5631    }
5632    fn arr(vals: Vec<f64>) -> ASTNode {
5633        ASTNode::new(
5634            ASTNodeType::Literal(LiteralValue::Array(vec![
5635                vals.into_iter().map(LiteralValue::Number).collect(),
5636            ])),
5637            None,
5638        )
5639    }
5640    #[test]
5641    fn median_even() {
5642        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(MEDIAN));
5643        let ctx = interp(&wb);
5644        let node = arr(vec![1.0, 3.0, 5.0, 7.0]);
5645        let f = ctx.context.get_function("", "MEDIAN").unwrap();
5646        let out = f
5647            .dispatch(
5648                &[ArgumentHandle::new(&node, &ctx)],
5649                &ctx.function_context(None),
5650            )
5651            .unwrap();
5652        assert_eq!(out, LiteralValue::Number(4.0));
5653    }
5654    #[test]
5655    fn median_odd() {
5656        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(MEDIAN));
5657        let ctx = interp(&wb);
5658        let node = arr(vec![1.0, 9.0, 5.0]);
5659        let f = ctx.context.get_function("", "MEDIAN").unwrap();
5660        let out = f
5661            .dispatch(
5662                &[ArgumentHandle::new(&node, &ctx)],
5663                &ctx.function_context(None),
5664            )
5665            .unwrap();
5666        assert_eq!(out, LiteralValue::Number(5.0));
5667    }
5668    #[test]
5669    fn large_basic() {
5670        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(LARGE));
5671        let ctx = interp(&wb);
5672        let nums = arr(vec![10.0, 20.0, 30.0]);
5673        let k = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(2.0)), None);
5674        let f = ctx.context.get_function("", "LARGE").unwrap();
5675        let out = f
5676            .dispatch(
5677                &[
5678                    ArgumentHandle::new(&nums, &ctx),
5679                    ArgumentHandle::new(&k, &ctx),
5680                ],
5681                &ctx.function_context(None),
5682            )
5683            .unwrap();
5684        assert_eq!(out, LiteralValue::Number(20.0));
5685    }
5686    #[test]
5687    fn small_basic() {
5688        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(SMALL));
5689        let ctx = interp(&wb);
5690        let nums = arr(vec![10.0, 20.0, 30.0]);
5691        let k = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(2.0)), None);
5692        let f = ctx.context.get_function("", "SMALL").unwrap();
5693        let out = f
5694            .dispatch(
5695                &[
5696                    ArgumentHandle::new(&nums, &ctx),
5697                    ArgumentHandle::new(&k, &ctx),
5698                ],
5699                &ctx.function_context(None),
5700            )
5701            .unwrap();
5702        assert_eq!(out, LiteralValue::Number(20.0));
5703    }
5704    #[test]
5705    fn percentile_inc_quarter() {
5706        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(PercentileInc));
5707        let ctx = interp(&wb);
5708        let nums = arr(vec![1.0, 2.0, 3.0, 4.0]);
5709        let p = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(0.25)), None);
5710        let f = ctx.context.get_function("", "PERCENTILE.INC").unwrap();
5711        match f
5712            .dispatch(
5713                &[
5714                    ArgumentHandle::new(&nums, &ctx),
5715                    ArgumentHandle::new(&p, &ctx),
5716                ],
5717                &ctx.function_context(None),
5718            )
5719            .unwrap()
5720            .into_literal()
5721        {
5722            LiteralValue::Number(v) => assert!((v - 1.75).abs() < 1e-9),
5723            other => panic!("unexpected {other:?}"),
5724        }
5725    }
5726    #[test]
5727    fn rank_eq_descending() {
5728        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(RankEqFn));
5729        let ctx = interp(&wb);
5730        // target 5 among {10,5,1} descending => ranks 1,2,3 => expect 2
5731        let target = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(5.0)), None);
5732        let arr_node = arr(vec![10.0, 5.0, 1.0]);
5733        let f = ctx.context.get_function("", "RANK.EQ").unwrap();
5734        let out = f
5735            .dispatch(
5736                &[
5737                    ArgumentHandle::new(&target, &ctx),
5738                    ArgumentHandle::new(&arr_node, &ctx),
5739                ],
5740                &ctx.function_context(None),
5741            )
5742            .unwrap();
5743        assert_eq!(out, LiteralValue::Number(2.0));
5744    }
5745    #[test]
5746    fn rank_eq_ascending_order_arg() {
5747        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(RankEqFn));
5748        let ctx = interp(&wb);
5749        // ascending order=1: array {1,5,10}; target 5 => rank 2
5750        let target = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(5.0)), None);
5751        let arr_node = arr(vec![1.0, 5.0, 10.0]);
5752        let order = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(1.0)), None);
5753        let f = ctx.context.get_function("", "RANK.EQ").unwrap();
5754        let out = f
5755            .dispatch(
5756                &[
5757                    ArgumentHandle::new(&target, &ctx),
5758                    ArgumentHandle::new(&arr_node, &ctx),
5759                    ArgumentHandle::new(&order, &ctx),
5760                ],
5761                &ctx.function_context(None),
5762            )
5763            .unwrap();
5764        assert_eq!(out, LiteralValue::Number(2.0));
5765    }
5766    #[test]
5767    fn rank_avg_ties() {
5768        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(RankAvgFn));
5769        let ctx = interp(&wb);
5770        // descending array {5,5,1} target 5 positions 1 and 2 avg -> 1.5
5771        let target = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(5.0)), None);
5772        let arr_node = arr(vec![5.0, 5.0, 1.0]);
5773        let f = ctx.context.get_function("", "RANK.AVG").unwrap();
5774        let out = f
5775            .dispatch(
5776                &[
5777                    ArgumentHandle::new(&target, &ctx),
5778                    ArgumentHandle::new(&arr_node, &ctx),
5779                ],
5780                &ctx.function_context(None),
5781            )
5782            .unwrap()
5783            .into_literal();
5784        match out {
5785            LiteralValue::Number(v) => assert!((v - 1.5).abs() < 1e-12),
5786            other => panic!("expected number got {other:?}"),
5787        }
5788    }
5789    #[test]
5790    fn stdev_var_sample_population() {
5791        let wb = TestWorkbook::new()
5792            .with_function(std::sync::Arc::new(StdevSample))
5793            .with_function(std::sync::Arc::new(StdevPop))
5794            .with_function(std::sync::Arc::new(VarSample))
5795            .with_function(std::sync::Arc::new(VarPop));
5796        let ctx = interp(&wb);
5797        let arr_node = arr(vec![2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0]); // variance population = 4, sample = 4.571428...
5798        let stdev_p = ctx.context.get_function("", "STDEV.P").unwrap();
5799        let stdev_s = ctx.context.get_function("", "STDEV.S").unwrap();
5800        let var_p = ctx.context.get_function("", "VAR.P").unwrap();
5801        let var_s = ctx.context.get_function("", "VAR.S").unwrap();
5802        let args = [ArgumentHandle::new(&arr_node, &ctx)];
5803        match var_p
5804            .dispatch(&args, &ctx.function_context(None))
5805            .unwrap()
5806            .into_literal()
5807        {
5808            LiteralValue::Number(v) => assert!((v - 4.0).abs() < 1e-12),
5809            other => panic!("unexpected {other:?}"),
5810        }
5811        match var_s
5812            .dispatch(&args, &ctx.function_context(None))
5813            .unwrap()
5814            .into_literal()
5815        {
5816            LiteralValue::Number(v) => assert!((v - 4.571428571428571).abs() < 1e-9),
5817            other => panic!("unexpected {other:?}"),
5818        }
5819        match stdev_p
5820            .dispatch(&args, &ctx.function_context(None))
5821            .unwrap()
5822            .into_literal()
5823        {
5824            LiteralValue::Number(v) => assert!((v - 2.0).abs() < 1e-12),
5825            other => panic!("unexpected {other:?}"),
5826        }
5827        match stdev_s
5828            .dispatch(&args, &ctx.function_context(None))
5829            .unwrap()
5830            .into_literal()
5831        {
5832            LiteralValue::Number(v) => assert!((v - 2.138089935).abs() < 1e-9),
5833            other => panic!("unexpected {other:?}"),
5834        }
5835    }
5836    #[test]
5837    fn quartile_inc_exc() {
5838        let wb = TestWorkbook::new()
5839            .with_function(std::sync::Arc::new(QuartileInc))
5840            .with_function(std::sync::Arc::new(QuartileExc));
5841        let ctx = interp(&wb);
5842        let arr_node = arr(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
5843        let q1 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(1.0)), None);
5844        let q2 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(2.0)), None);
5845        let f_inc = ctx.context.get_function("", "QUARTILE.INC").unwrap();
5846        let f_exc = ctx.context.get_function("", "QUARTILE.EXC").unwrap();
5847        let arg_inc_q1 = [
5848            ArgumentHandle::new(&arr_node, &ctx),
5849            ArgumentHandle::new(&q1, &ctx),
5850        ];
5851        let arg_inc_q2 = [
5852            ArgumentHandle::new(&arr_node, &ctx),
5853            ArgumentHandle::new(&q2, &ctx),
5854        ];
5855        match f_inc
5856            .dispatch(&arg_inc_q1, &ctx.function_context(None))
5857            .unwrap()
5858            .into_literal()
5859        {
5860            LiteralValue::Number(v) => assert!((v - 2.0).abs() < 1e-12),
5861            other => panic!("unexpected {other:?}"),
5862        }
5863        match f_inc
5864            .dispatch(&arg_inc_q2, &ctx.function_context(None))
5865            .unwrap()
5866            .into_literal()
5867        {
5868            LiteralValue::Number(v) => assert!((v - 3.0).abs() < 1e-12),
5869            other => panic!("unexpected {other:?}"),
5870        }
5871        // QUARTILE.EXC Q1 for 5-point set uses exclusive percentile => 1.5
5872        match f_exc
5873            .dispatch(&arg_inc_q1, &ctx.function_context(None))
5874            .unwrap()
5875            .into_literal()
5876        {
5877            LiteralValue::Number(v) => assert!((v - 1.5).abs() < 1e-12),
5878            other => panic!("unexpected {other:?}"),
5879        }
5880        match f_exc
5881            .dispatch(&arg_inc_q2, &ctx.function_context(None))
5882            .unwrap()
5883            .into_literal()
5884        {
5885            LiteralValue::Number(v) => assert!((v - 3.0).abs() < 1e-12),
5886            other => panic!("unexpected {other:?}"),
5887        }
5888    }
5889
5890    // --- eval()/dispatch equivalence tests for variance / stdev ---
5891    #[test]
5892    fn fold_equivalence_var_stdev() {
5893        use crate::function::Function as _; // trait import
5894        let wb = TestWorkbook::new()
5895            .with_function(std::sync::Arc::new(VarSample))
5896            .with_function(std::sync::Arc::new(VarPop))
5897            .with_function(std::sync::Arc::new(StdevSample))
5898            .with_function(std::sync::Arc::new(StdevPop));
5899        let ctx = interp(&wb);
5900        let arr_node = arr(vec![1.0, 2.0, 5.0, 5.0, 9.0]);
5901        let args = [ArgumentHandle::new(&arr_node, &ctx)];
5902
5903        let var_s_fn = VarSample; // concrete instance to call eval()
5904        let var_p_fn = VarPop;
5905        let stdev_s_fn = StdevSample;
5906        let stdev_p_fn = StdevPop;
5907
5908        let fctx = ctx.function_context(None);
5909        // Dispatch results (will use fold path)
5910        let disp_var_s = ctx
5911            .context
5912            .get_function("", "VAR.S")
5913            .unwrap()
5914            .dispatch(&args, &fctx)
5915            .unwrap()
5916            .into_literal();
5917        let disp_var_p = ctx
5918            .context
5919            .get_function("", "VAR.P")
5920            .unwrap()
5921            .dispatch(&args, &fctx)
5922            .unwrap()
5923            .into_literal();
5924        let disp_stdev_s = ctx
5925            .context
5926            .get_function("", "STDEV.S")
5927            .unwrap()
5928            .dispatch(&args, &fctx)
5929            .unwrap()
5930            .into_literal();
5931        let disp_stdev_p = ctx
5932            .context
5933            .get_function("", "STDEV.P")
5934            .unwrap()
5935            .dispatch(&args, &fctx)
5936            .unwrap()
5937            .into_literal();
5938
5939        // Scalar path results
5940        let scalar_var_s = var_s_fn.eval(&args, &fctx).unwrap().into_literal();
5941        let scalar_var_p = var_p_fn.eval(&args, &fctx).unwrap().into_literal();
5942        let scalar_stdev_s = stdev_s_fn.eval(&args, &fctx).unwrap().into_literal();
5943        let scalar_stdev_p = stdev_p_fn.eval(&args, &fctx).unwrap().into_literal();
5944
5945        fn assert_close(a: &LiteralValue, b: &LiteralValue) {
5946            match (a, b) {
5947                (LiteralValue::Number(x), LiteralValue::Number(y)) => {
5948                    assert!((x - y).abs() < 1e-12, "mismatch {x} vs {y}")
5949                }
5950                _ => assert_eq!(a, b),
5951            }
5952        }
5953        assert_close(&disp_var_s, &scalar_var_s);
5954        assert_close(&disp_var_p, &scalar_var_p);
5955        assert_close(&disp_stdev_s, &scalar_stdev_s);
5956        assert_close(&disp_stdev_p, &scalar_stdev_p);
5957    }
5958
5959    #[test]
5960    fn fold_equivalence_edge_cases() {
5961        use crate::function::Function as _;
5962        let wb = TestWorkbook::new()
5963            .with_function(std::sync::Arc::new(VarSample))
5964            .with_function(std::sync::Arc::new(VarPop))
5965            .with_function(std::sync::Arc::new(StdevSample))
5966            .with_function(std::sync::Arc::new(StdevPop));
5967        let ctx = interp(&wb);
5968        // Single numeric element -> sample variance/div0, population variance 0
5969        let single = arr(vec![42.0]);
5970        let args_single = [ArgumentHandle::new(&single, &ctx)];
5971        let fctx = ctx.function_context(None);
5972        let disp_var_s = ctx
5973            .context
5974            .get_function("", "VAR.S")
5975            .unwrap()
5976            .dispatch(&args_single, &fctx)
5977            .unwrap();
5978        let scalar_var_s = VarSample.eval(&args_single, &fctx).unwrap().into_literal();
5979        assert_eq!(disp_var_s, scalar_var_s);
5980        let disp_var_p = ctx
5981            .context
5982            .get_function("", "VAR.P")
5983            .unwrap()
5984            .dispatch(&args_single, &fctx)
5985            .unwrap();
5986        let scalar_var_p = VarPop.eval(&args_single, &fctx).unwrap().into_literal();
5987        assert_eq!(disp_var_p, scalar_var_p);
5988        let disp_stdev_p = ctx
5989            .context
5990            .get_function("", "STDEV.P")
5991            .unwrap()
5992            .dispatch(&args_single, &fctx)
5993            .unwrap();
5994        let scalar_stdev_p = StdevPop.eval(&args_single, &fctx).unwrap().into_literal();
5995        assert_eq!(disp_stdev_p, scalar_stdev_p);
5996        let disp_stdev_s = ctx
5997            .context
5998            .get_function("", "STDEV.S")
5999            .unwrap()
6000            .dispatch(&args_single, &fctx)
6001            .unwrap();
6002        let scalar_stdev_s = StdevSample
6003            .eval(&args_single, &fctx)
6004            .unwrap()
6005            .into_literal();
6006        assert_eq!(disp_stdev_s, scalar_stdev_s);
6007    }
6008
6009    #[test]
6010    fn legacy_aliases_match_modern() {
6011        let wb = TestWorkbook::new()
6012            .with_function(std::sync::Arc::new(PercentileInc))
6013            .with_function(std::sync::Arc::new(QuartileInc))
6014            .with_function(std::sync::Arc::new(RankEqFn));
6015        let ctx = interp(&wb);
6016        let arr_node = arr(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
6017        let p = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(0.4)), None);
6018        let q2 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(2.0)), None);
6019        let target = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(4.0)), None);
6020        let args_p = [
6021            ArgumentHandle::new(&arr_node, &ctx),
6022            ArgumentHandle::new(&p, &ctx),
6023        ];
6024        let args_q = [
6025            ArgumentHandle::new(&arr_node, &ctx),
6026            ArgumentHandle::new(&q2, &ctx),
6027        ];
6028        let args_rank = [
6029            ArgumentHandle::new(&target, &ctx),
6030            ArgumentHandle::new(&arr_node, &ctx),
6031        ];
6032        let modern_p = ctx
6033            .context
6034            .get_function("", "PERCENTILE.INC")
6035            .unwrap()
6036            .dispatch(&args_p, &ctx.function_context(None))
6037            .unwrap()
6038            .into_literal();
6039        let legacy_p = ctx
6040            .context
6041            .get_function("", "PERCENTILE")
6042            .unwrap()
6043            .dispatch(&args_p, &ctx.function_context(None))
6044            .unwrap()
6045            .into_literal();
6046        assert_eq!(modern_p, legacy_p);
6047        let modern_q = ctx
6048            .context
6049            .get_function("", "QUARTILE.INC")
6050            .unwrap()
6051            .dispatch(&args_q, &ctx.function_context(None))
6052            .unwrap()
6053            .into_literal();
6054        let legacy_q = ctx
6055            .context
6056            .get_function("", "QUARTILE")
6057            .unwrap()
6058            .dispatch(&args_q, &ctx.function_context(None))
6059            .unwrap()
6060            .into_literal();
6061        assert_eq!(modern_q, legacy_q);
6062        let modern_rank = ctx
6063            .context
6064            .get_function("", "RANK.EQ")
6065            .unwrap()
6066            .dispatch(&args_rank, &ctx.function_context(None))
6067            .unwrap()
6068            .into_literal();
6069        let legacy_rank = ctx
6070            .context
6071            .get_function("", "RANK")
6072            .unwrap()
6073            .dispatch(&args_rank, &ctx.function_context(None))
6074            .unwrap()
6075            .into_literal();
6076        assert_eq!(modern_rank, legacy_rank);
6077    }
6078
6079    #[test]
6080    fn mode_single_basic_and_alias() {
6081        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(ModeSingleFn));
6082        let ctx = interp(&wb);
6083        let arr_node = arr(vec![5.0, 2.0, 2.0, 3.0, 3.0, 3.0]);
6084        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6085        let mode_sngl = ctx
6086            .context
6087            .get_function("", "MODE.SNGL")
6088            .unwrap()
6089            .dispatch(&args, &ctx.function_context(None))
6090            .unwrap()
6091            .into_literal();
6092        assert_eq!(mode_sngl, LiteralValue::Number(3.0));
6093        let mode_alias = ctx
6094            .context
6095            .get_function("", "MODE")
6096            .unwrap()
6097            .dispatch(&args, &ctx.function_context(None))
6098            .unwrap()
6099            .into_literal();
6100        assert_eq!(mode_alias, mode_sngl);
6101    }
6102
6103    #[test]
6104    fn mode_single_no_duplicates() {
6105        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(ModeSingleFn));
6106        let ctx = interp(&wb);
6107        let arr_node = arr(vec![1.0, 2.0, 3.0]);
6108        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6109        let out = ctx
6110            .context
6111            .get_function("", "MODE.SNGL")
6112            .unwrap()
6113            .dispatch(&args, &ctx.function_context(None))
6114            .unwrap()
6115            .into_literal();
6116        match out {
6117            LiteralValue::Error(e) => assert!(e.to_string().contains("#N/A")),
6118            _ => panic!("expected #N/A"),
6119        }
6120    }
6121
6122    #[test]
6123    fn mode_multi_basic() {
6124        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(ModeMultiFn));
6125        let ctx = interp(&wb);
6126        let arr_node = arr(vec![2.0, 3.0, 2.0, 4.0, 3.0, 5.0, 2.0, 3.0]);
6127        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6128        let out = ctx
6129            .context
6130            .get_function("", "MODE.MULT")
6131            .unwrap()
6132            .dispatch(&args, &ctx.function_context(None))
6133            .unwrap()
6134            .into_literal();
6135        let expected = LiteralValue::Array(vec![
6136            vec![LiteralValue::Number(2.0)],
6137            vec![LiteralValue::Number(3.0)],
6138        ]);
6139        assert_eq!(out, expected);
6140    }
6141
6142    #[test]
6143    fn large_small_fold_vs_scalar() {
6144        let wb = TestWorkbook::new()
6145            .with_function(std::sync::Arc::new(LARGE))
6146            .with_function(std::sync::Arc::new(SMALL));
6147        let ctx = interp(&wb);
6148        let arr_node = arr(vec![10.0, 5.0, 7.0, 12.0, 9.0]);
6149        let k_node = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(2.0)), None);
6150        let args = [
6151            ArgumentHandle::new(&arr_node, &ctx),
6152            ArgumentHandle::new(&k_node, &ctx),
6153        ];
6154        let f_large = ctx.context.get_function("", "LARGE").unwrap();
6155        let disp_large = f_large
6156            .dispatch(&args, &ctx.function_context(None))
6157            .unwrap()
6158            .into_literal();
6159        let scalar_large = LARGE
6160            .eval(&args, &ctx.function_context(None))
6161            .unwrap()
6162            .into_literal();
6163        assert_eq!(disp_large, scalar_large);
6164        let f_small = ctx.context.get_function("", "SMALL").unwrap();
6165        let disp_small = f_small
6166            .dispatch(&args, &ctx.function_context(None))
6167            .unwrap()
6168            .into_literal();
6169        let scalar_small = SMALL
6170            .eval(&args, &ctx.function_context(None))
6171            .unwrap()
6172            .into_literal();
6173        assert_eq!(disp_small, scalar_small);
6174    }
6175
6176    #[test]
6177    fn mode_fold_vs_scalar() {
6178        let wb = TestWorkbook::new()
6179            .with_function(std::sync::Arc::new(ModeSingleFn))
6180            .with_function(std::sync::Arc::new(ModeMultiFn));
6181        let ctx = interp(&wb);
6182        let arr_node = arr(vec![2.0, 3.0, 2.0, 4.0, 3.0, 3.0, 2.0]);
6183        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6184        let f_single = ctx.context.get_function("", "MODE.SNGL").unwrap();
6185        let disp_single = f_single
6186            .dispatch(&args, &ctx.function_context(None))
6187            .unwrap()
6188            .into_literal();
6189        let scalar_single = ModeSingleFn
6190            .eval(&args, &ctx.function_context(None))
6191            .unwrap()
6192            .into_literal();
6193        assert_eq!(disp_single, scalar_single);
6194        let f_multi = ctx.context.get_function("", "MODE.MULT").unwrap();
6195        let disp_multi = f_multi
6196            .dispatch(&args, &ctx.function_context(None))
6197            .unwrap()
6198            .into_literal();
6199        let scalar_multi = ModeMultiFn
6200            .eval(&args, &ctx.function_context(None))
6201            .unwrap()
6202            .into_literal();
6203        assert_eq!(disp_multi, scalar_multi);
6204    }
6205
6206    #[test]
6207    fn median_fold_vs_scalar_even() {
6208        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(MEDIAN));
6209        let ctx = interp(&wb);
6210        let arr_node = arr(vec![7.0, 1.0, 9.0, 5.0]); // sorted: 1,5,7,9 median=(5+7)/2=6
6211        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6212        let f_med = ctx.context.get_function("", "MEDIAN").unwrap();
6213        let disp = f_med
6214            .dispatch(&args, &ctx.function_context(None))
6215            .unwrap()
6216            .into_literal();
6217        let scalar = MEDIAN
6218            .eval(&args, &ctx.function_context(None))
6219            .unwrap()
6220            .into_literal();
6221        assert_eq!(disp, scalar);
6222        assert_eq!(disp, LiteralValue::Number(6.0));
6223    }
6224
6225    #[test]
6226    fn median_fold_vs_scalar_odd() {
6227        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(MEDIAN));
6228        let ctx = interp(&wb);
6229        let arr_node = arr(vec![9.0, 2.0, 5.0]); // sorted 2,5,9 median=5
6230        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6231        let f_med = ctx.context.get_function("", "MEDIAN").unwrap();
6232        let disp = f_med
6233            .dispatch(&args, &ctx.function_context(None))
6234            .unwrap()
6235            .into_literal();
6236        let scalar = MEDIAN
6237            .eval(&args, &ctx.function_context(None))
6238            .unwrap()
6239            .into_literal();
6240        assert_eq!(disp, scalar);
6241        assert_eq!(disp, LiteralValue::Number(5.0));
6242    }
6243
6244    // Newly added edge case tests for statistical semantics.
6245    #[test]
6246    fn percentile_inc_edges() {
6247        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(PercentileInc));
6248        let ctx = interp(&wb);
6249        let arr_node = arr(vec![10.0, 20.0, 30.0, 40.0]);
6250        let p0 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(0.0)), None);
6251        let p1 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(1.0)), None);
6252        let f = ctx.context.get_function("", "PERCENTILE.INC").unwrap();
6253        let args0 = [
6254            ArgumentHandle::new(&arr_node, &ctx),
6255            ArgumentHandle::new(&p0, &ctx),
6256        ];
6257        let args1 = [
6258            ArgumentHandle::new(&arr_node, &ctx),
6259            ArgumentHandle::new(&p1, &ctx),
6260        ];
6261        assert_eq!(
6262            f.dispatch(&args0, &ctx.function_context(None))
6263                .unwrap()
6264                .into_literal(),
6265            LiteralValue::Number(10.0)
6266        );
6267        assert_eq!(
6268            f.dispatch(&args1, &ctx.function_context(None))
6269                .unwrap()
6270                .into_literal(),
6271            LiteralValue::Number(40.0)
6272        );
6273    }
6274
6275    #[test]
6276    fn percentile_exc_invalid() {
6277        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(PercentileExc));
6278        let ctx = interp(&wb);
6279        let arr_node = arr(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
6280        let p_bad0 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(0.0)), None);
6281        let p_bad1 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(1.0)), None);
6282        let f = ctx.context.get_function("", "PERCENTILE.EXC").unwrap();
6283        for bad in [&p_bad0, &p_bad1] {
6284            let args = [
6285                ArgumentHandle::new(&arr_node, &ctx),
6286                ArgumentHandle::new(bad, &ctx),
6287            ];
6288            match f
6289                .dispatch(&args, &ctx.function_context(None))
6290                .unwrap()
6291                .into_literal()
6292            {
6293                LiteralValue::Error(e) => assert!(e.to_string().contains("#NUM!")),
6294                other => panic!("expected #NUM! got {other:?}"),
6295            }
6296        }
6297    }
6298
6299    #[test]
6300    fn quartile_invalids() {
6301        let wb = TestWorkbook::new()
6302            .with_function(std::sync::Arc::new(QuartileInc))
6303            .with_function(std::sync::Arc::new(QuartileExc));
6304        let ctx = interp(&wb);
6305        let arr_node = arr(vec![1.0, 2.0, 3.0]);
6306        // QUARTILE.INC invalid q=5
6307        let q5 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(5.0)), None);
6308        let args_bad_inc = [
6309            ArgumentHandle::new(&arr_node, &ctx),
6310            ArgumentHandle::new(&q5, &ctx),
6311        ];
6312        match ctx
6313            .context
6314            .get_function("", "QUARTILE.INC")
6315            .unwrap()
6316            .dispatch(&args_bad_inc, &ctx.function_context(None))
6317            .unwrap()
6318            .into_literal()
6319        {
6320            LiteralValue::Error(e) => assert!(e.to_string().contains("#NUM!")),
6321            other => panic!("expected #NUM! {other:?}"),
6322        }
6323        // QUARTILE.EXC invalid q=0
6324        let q0 = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(0.0)), None);
6325        let args_bad_exc = [
6326            ArgumentHandle::new(&arr_node, &ctx),
6327            ArgumentHandle::new(&q0, &ctx),
6328        ];
6329        match ctx
6330            .context
6331            .get_function("", "QUARTILE.EXC")
6332            .unwrap()
6333            .dispatch(&args_bad_exc, &ctx.function_context(None))
6334            .unwrap()
6335            .into_literal()
6336        {
6337            LiteralValue::Error(e) => assert!(e.to_string().contains("#NUM!")),
6338            other => panic!("expected #NUM! {other:?}"),
6339        }
6340    }
6341
6342    #[test]
6343    fn rank_target_not_found() {
6344        let wb = TestWorkbook::new()
6345            .with_function(std::sync::Arc::new(RankEqFn))
6346            .with_function(std::sync::Arc::new(RankAvgFn));
6347        let ctx = interp(&wb);
6348        let arr_node = arr(vec![1.0, 2.0, 3.0]);
6349        let target = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(4.0)), None);
6350        let args = [
6351            ArgumentHandle::new(&target, &ctx),
6352            ArgumentHandle::new(&arr_node, &ctx),
6353        ];
6354        for name in ["RANK.EQ", "RANK.AVG"] {
6355            match ctx
6356                .context
6357                .get_function("", name)
6358                .unwrap()
6359                .dispatch(&args, &ctx.function_context(None))
6360                .unwrap()
6361                .into_literal()
6362            {
6363                LiteralValue::Error(e) => assert!(e.to_string().contains("#N/A")),
6364                other => panic!("expected #N/A {other:?}"),
6365            }
6366        }
6367    }
6368
6369    #[test]
6370    fn mode_mult_ordering() {
6371        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(ModeMultiFn));
6372        let ctx = interp(&wb);
6373        // Two modes with same frequency; ensure ascending ordering in array result
6374        let arr_node = arr(vec![5.0, 2.0, 2.0, 5.0, 3.0, 7.0, 5.0, 2.0]); // 2 and 5 appear 4 times each
6375        let args = [ArgumentHandle::new(&arr_node, &ctx)];
6376        let out = ctx
6377            .context
6378            .get_function("", "MODE.MULT")
6379            .unwrap()
6380            .dispatch(&args, &ctx.function_context(None))
6381            .unwrap()
6382            .into_literal();
6383        match out {
6384            LiteralValue::Array(rows) => {
6385                let vals: Vec<f64> = rows
6386                    .into_iter()
6387                    .map(|r| {
6388                        if let LiteralValue::Number(n) = r[0] {
6389                            n
6390                        } else {
6391                            panic!("expected number")
6392                        }
6393                    })
6394                    .collect();
6395                assert_eq!(vals, vec![2.0, 5.0]);
6396            }
6397            other => panic!("expected array {other:?}"),
6398        }
6399    }
6400
6401    #[test]
6402    fn boolean_and_text_in_range_are_ignored() {
6403        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(StdevPop));
6404        let ctx = interp(&wb);
6405        let mixed = ASTNode::new(
6406            ASTNodeType::Literal(LiteralValue::Array(vec![vec![
6407                LiteralValue::Number(1.0),
6408                LiteralValue::Text("ABC".into()),
6409                LiteralValue::Boolean(true),
6410                LiteralValue::Number(4.0),
6411            ]])),
6412            None,
6413        );
6414        let f = ctx.context.get_function("", "STDEV.P").unwrap();
6415        let out = f
6416            .dispatch(
6417                &[ArgumentHandle::new(&mixed, &ctx)],
6418                &ctx.function_context(None),
6419            )
6420            .unwrap()
6421            .into_literal();
6422        // NOTE: Inline array literal is treated as a direct scalar argument (not a range reference),
6423        // so boolean TRUE is coerced to 1. Dataset becomes {1,1,4}; population stdev = sqrt(6/3)=sqrt(2).
6424        match out {
6425            LiteralValue::Number(v) => {
6426                assert!((v - 2f64.sqrt()).abs() < 1e-12, "expected sqrt(2) got {v}")
6427            }
6428            other => panic!("unexpected {other:?}"),
6429        }
6430    }
6431
6432    #[test]
6433    fn boolean_direct_arg_coerces() {
6434        let wb = TestWorkbook::new().with_function(std::sync::Arc::new(StdevPop));
6435        let ctx = interp(&wb);
6436        let one = ASTNode::new(ASTNodeType::Literal(LiteralValue::Number(1.0)), None);
6437        let t = ASTNode::new(ASTNodeType::Literal(LiteralValue::Boolean(true)), None);
6438        let f = ctx.context.get_function("", "STDEV.P").unwrap();
6439        let args = [
6440            ArgumentHandle::new(&one, &ctx),
6441            ArgumentHandle::new(&t, &ctx),
6442        ];
6443        let out = f
6444            .dispatch(&args, &ctx.function_context(None))
6445            .unwrap()
6446            .into_literal();
6447        assert_eq!(out, LiteralValue::Number(0.0));
6448    }
6449}
formualizer_eval/builtins/stats/mod.rs

formualizer_eval/builtins/stats/
mod.rs
