kalc_lib/math.rs
1use crate::cas::isolate;
2#[cfg(feature = "fastrand")]
3use crate::complex::{rand_gamma, rand_norm};
4use crate::{
5 complex::{
6 LimSide::{Both, Left, Right},
7 NumStr,
8 NumStr::{
9 And, Comma, Converse, Conversion, Division, Equal, Exponent, Func, Greater,
10 GreaterEqual, Implies, InternalMultiplication, LeftBracket, LeftCurlyBracket, Lesser,
11 LesserEqual, Matrix, Minus, Modulo, Multiplication, Nand, NearEqual, Nor, Not,
12 NotEqual, Num, Or, Plus, PlusMinus, Range, RightBracket, RightCurlyBracket, Root,
13 ShiftLeft, ShiftRight, Tetration, Vector, Xor,
14 },
15 about_eq, add, and, area, atan, binomial, change_basis, cofactor, coordinate, cube, cubic,
16 determinant, digamma, div, eigenvalues, eigenvectors, eq, erf, erfc, eta, euleriannumbers,
17 euleriannumbersint, extrema, gamma, gcd, ge, generalized_eigenvectors, gt, hsv2rgb,
18 identity, implies, incomplete_beta, incomplete_gamma, inverse, iter, jcf, kernel, lambertw,
19 length, limit, lower_incomplete_gamma, minors, mul_units, mvec, nand, ne, nor, not,
20 nth_prime, or, pow_nth, prime_factors, quadratic, quartic, range, rcf, recursion,
21 regularized_incomplete_beta, rem, root, rref, shl, shr, slog, slope, solve, sort, sort_mat,
22 sqr, sub, subfactorial, sum, surface_area, taylor, tetration, to, to_cyl, to_polar, trace,
23 transpose, unity, variance, xor, zeta,
24 },
25 fraction::{c_to_rational, rationalize},
26 misc::do_math_with_var,
27 units::{AngleType, Number, Options, Units},
28};
29use rug::{
30 Complex, Float, Integer,
31 float::{
32 Constant::Pi,
33 Special::{Infinity, Nan},
34 },
35 integer::IsPrime,
36 ops::Pow,
37};
38#[cfg(feature = "fastrand")]
39use std::ops::Rem;
40use std::{cmp::Ordering, time::SystemTime};
41pub fn do_math(
42 mut function: Vec<NumStr>,
43 options: Options,
44 mut func_vars: Vec<(String, Vec<NumStr>)>,
45) -> Result<NumStr, &'static str> {
46 if function.is_empty() {
47 return Err(" ");
48 }
49 compute_funcvars(&mut function, options, &mut func_vars);
50 let mut i = 0;
51 while i < function.len() {
52 match &function[i] {
53 LeftCurlyBracket => {
54 let mut j = i + 1;
55 let mut count = 1;
56 while count > 0 {
57 if j >= function.len() {
58 return Err("curly bracket err");
59 }
60 match &function[j] {
61 LeftCurlyBracket => count += 1,
62 RightCurlyBracket => count -= 1,
63 _ => {}
64 }
65 j += 1;
66 }
67 if i + 1 == j - 1 {
68 return Err("no interior vector");
69 }
70 let mut single = 0;
71 let v = &function[i + 1..j - 1];
72 let mut vec = Vec::new();
73 let mut mat = Vec::<Vec<Number>>::new();
74 for (f, n) in v.iter().enumerate() {
75 match &n {
76 Comma if count == 0 => {
77 let z = do_math(v[single..f].to_vec(), options, func_vars.clone())?;
78 match z {
79 Num(n) => vec.push(*n),
80 Vector(n) => mat.push(n),
81 _ => return Err("broken matrix"),
82 }
83 single = f + 1;
84 }
85 LeftBracket | LeftCurlyBracket => count += 1,
86 RightBracket | RightCurlyBracket => count -= 1,
87 _ => {}
88 }
89 }
90 if single != v.len() {
91 let z = do_math(v[single..].to_vec(), options, func_vars.clone())?;
92 match z {
93 Num(n) => vec.push(*n),
94 Vector(n) => mat.push(n),
95 _ => return Err("broken matrix"),
96 }
97 }
98 function.drain(i..j);
99 if !mat.is_empty() {
100 if vec.is_empty() {
101 function.insert(i, Matrix(mat));
102 } else {
103 return Err("vector err");
104 }
105 } else {
106 function.insert(i, Vector(vec));
107 }
108 }
109 LeftBracket => {
110 let mut j = i + 1;
111 let mut count = 1;
112 while count > 0 {
113 if j >= function.len() {
114 return Err("round bracket err");
115 }
116 match &function[j] {
117 LeftBracket => count += 1,
118 RightBracket => count -= 1,
119 _ => {}
120 }
121 j += 1;
122 }
123 if i + 1 == j - 1 {
124 return Err("no interior bracket");
125 }
126 if i != 0 {
127 if let Func(k) = &function[i - 1] {
128 if matches!(
129 k.as_str(),
130 "next"
131 | "log"
132 | "exp"
133 | "ζ"
134 | "polygamma"
135 | "digamma"
136 | "inter"
137 | "interpolate"
138 | "lobf"
139 | "plane"
140 | "lineofbestfit"
141 | "ψ"
142 | "rotate"
143 | "multinomial"
144 | "gcd"
145 | "gcf"
146 | "lcm"
147 | "ssrt"
148 | "W"
149 | "unity"
150 | "productlog"
151 | "lambertw"
152 | "ceil"
153 | "floor"
154 | "round"
155 | "int"
156 | "trunc"
157 | "frac"
158 | "fract"
159 | "slog"
160 | "root"
161 | "atan"
162 | "arctan"
163 | "atan2"
164 | "hypot"
165 | "normP"
166 | "normD"
167 | "betaP"
168 | "betaC"
169 | "bi"
170 | "binomial"
171 | "angle"
172 | "cross"
173 | "dot"
174 | "part"
175 | "proj"
176 | "project"
177 | "oproj"
178 | "oproject"
179 | "C"
180 | "P"
181 | "Ap"
182 | "An"
183 | "gamma"
184 | "Γ"
185 | "γ"
186 | "lower_gamma"
187 | "ph"
188 | "pochhammer"
189 | "Β"
190 | "B"
191 | "beta"
192 | "I"
193 | "quad"
194 | "quadratic"
195 | "cubic"
196 | "quartic"
197 | "percentilerank"
198 | "percentile"
199 | "eigenvalues"
200 | "eigenvectors"
201 | "generalized_eigenvectors"
202 | "change_basis"
203 | "coordinate"
204 | "mod"
205 | "covariance"
206 | "cov"
207 | "rand_norm"
208 | "rand_uniform"
209 | "rand_int"
210 | "rand_gamma"
211 | "rand_beta"
212 | "gamma_pdf"
213 | "gamma_cdf"
214 | "beta_cdf"
215 | "beta_pdf"
216 | "norm_cdf"
217 | "norm_pdf"
218 | "lognorm_cdf"
219 | "binomial_cdf"
220 | "geometric_cdf"
221 | "lognorm_pdf"
222 | "binomial_pmf"
223 | "geometric_pmf"
224 | "rand_lognorm"
225 | "rand_binomial"
226 | "poisson_pmf"
227 | "poisson_cdf"
228 | "rand_neg_binomial"
229 | "neg_binomial_cdf"
230 | "neg_binomial_pmf"
231 | "hypergeometric_pmf"
232 | "hypergeometric_cdf"
233 | "rand_hypergeometric"
234 | "neg_hypergeometric_pmf"
235 | "neg_hypergeometric_cdf"
236 | "rand_neg_hypergeometric"
237 | "union"
238 | "poly"
239 | "polynomial"
240 | "intersection"
241 | "set_difference"
242 | "symmetric_difference"
243 | "cartesian_product"
244 | "remove"
245 | "extend"
246 | "link"
247 | "subset"
248 | "element"
249 ) {
250 function.remove(j - 1);
251 function.remove(i);
252 let v = function.drain(i..j - 2).collect::<Vec<NumStr>>();
253 count = 0;
254 let mut place = Vec::new();
255 for (f, n) in v.iter().enumerate() {
256 match n {
257 Comma if count == 0 => place.push(f),
258 LeftBracket | LeftCurlyBracket => count += 1,
259 RightBracket | RightCurlyBracket => count -= 1,
260 _ => {}
261 }
262 }
263 if !place.is_empty() {
264 let mut func = vec![function[i - 1].clone()];
265 func.push(do_math(
266 v[..place[0]].to_vec(),
267 options,
268 func_vars.clone(),
269 )?);
270 for (k, l) in place.iter().enumerate() {
271 func.push(do_math(
272 v[l + 1..if k + 1 != place.len() {
273 place[k + 1]
274 } else {
275 v.len()
276 }]
277 .to_vec(),
278 options,
279 func_vars.clone(),
280 )?);
281 }
282 function[i - 1] = do_math(func, options, func_vars.clone())?;
283 } else {
284 let v = vec![
285 function[i - 1].clone(),
286 do_math(v, options, func_vars.clone())?,
287 ];
288 function[i - 1] = do_math(v, options, func_vars.clone())?;
289 }
290 continue;
291 } else if matches!(
292 k.as_str(),
293 "sum"
294 | "area"
295 | "surfacearea"
296 | "sarea"
297 | "solve"
298 | "∫"
299 | "length"
300 | "slope"
301 | "taylor"
302 | "iter"
303 | "extrema"
304 | "summation"
305 | "prod"
306 | "product"
307 | "Σ"
308 | "Π"
309 | "vec"
310 | "mat"
311 | "piecewise"
312 | "pw"
313 | "D"
314 | "integrate"
315 | "arclength"
316 | "lim"
317 | "limit"
318 | "set"
319 | "isolate"
320 ) {
321 i = j - 1;
322 continue;
323 } else {
324 function.remove(j - 1);
325 function.remove(i);
326 function[i - 1] = do_math(
327 vec![
328 function[i - 1].clone(),
329 do_math(
330 function.drain(i..j - 2).collect::<Vec<NumStr>>(),
331 options,
332 func_vars.clone(),
333 )?,
334 ],
335 options,
336 func_vars.clone(),
337 )?;
338 continue;
339 }
340 }
341 }
342 function.remove(j - 1);
343 function[i] = do_math(
344 function.drain(i + 1..j - 1).collect::<Vec<NumStr>>(),
345 options,
346 func_vars.clone(),
347 )?;
348 }
349 _ => {}
350 }
351 i += 1;
352 }
353 if function.len() == 1 {
354 if let Func(s) = &function[0] {
355 if !matches!(s.as_str(), "rnd" | "rand" | "epoch") {
356 return Ok(function[0].clone());
357 }
358 } else {
359 return Ok(function[0].clone());
360 }
361 }
362 i = 0;
363 let to_deg = match options.angle {
364 AngleType::Degrees => 180 / Complex::with_val(options.prec, Pi),
365 AngleType::Radians => Complex::with_val(options.prec, 1),
366 AngleType::Gradians => 200 / Complex::with_val(options.prec, Pi),
367 };
368 while i < function.len().saturating_sub(1) {
369 if let Func(s) = &function[i].clone() {
370 if !matches!(s.as_str(), "epoch" | "rnd" | "rand") {
371 if matches!(
372 s.as_str(),
373 "sum"
374 | "area"
375 | "sarea"
376 | "surfacearea"
377 | "solve"
378 | "∫"
379 | "length"
380 | "slope"
381 | "taylor"
382 | "iter"
383 | "extrema"
384 | "product"
385 | "prod"
386 | "summation"
387 | "Σ"
388 | "Π"
389 | "vec"
390 | "mat"
391 | "piecewise"
392 | "pw"
393 | "D"
394 | "integrate"
395 | "arclength"
396 | "lim"
397 | "limit"
398 | "set"
399 | "isolate"
400 ) {
401 let mut place = Vec::new();
402 let mut count = 0;
403 let mut count2 = 0;
404 for (f, n) in function[i + 2..].iter().enumerate() {
405 if s == "piecewise" || s == "pw" {
406 match n {
407 Comma if (count == 0 || count == 1) && count2 == 0 => {
408 place.push(f + i + 2);
409 }
410 LeftBracket => {
411 count += 1;
412 count2 += 1;
413 }
414 LeftCurlyBracket => count += 1,
415 RightCurlyBracket => {
416 if count == 0 {
417 place.push(f + i + 2);
418 break;
419 }
420 count -= 1;
421 }
422 RightBracket => {
423 if count == 0 {
424 place.push(f + i + 2);
425 break;
426 }
427 count -= 1;
428 count2 -= 1;
429 }
430 _ => {}
431 }
432 } else {
433 match n {
434 Comma if count == 0 => {
435 place.push(f + i + 2);
436 }
437 LeftBracket | LeftCurlyBracket => {
438 count += 1;
439 }
440 RightBracket | RightCurlyBracket => {
441 if count == 0 {
442 place.push(f + i + 2);
443 break;
444 }
445 count -= 1;
446 }
447 _ => {}
448 }
449 }
450 }
451 match (
452 s.as_str(),
453 if place.is_empty() || place[0] - 1 == i + 1 {
454 Func(String::new())
455 } else {
456 function[place[0] - 1].clone()
457 },
458 ) {
459 ("iter", Func(var)) if place.len() == 4 || place.len() == 5 => {
460 function[i] = iter(
461 function[place[0] + 1..place[1]].to_vec(),
462 func_vars.clone(),
463 options,
464 var.to_string(),
465 do_math(
466 function[place[1] + 1..place[2]].to_vec(),
467 options,
468 func_vars.clone(),
469 )?,
470 do_math(
471 function[place[2] + 1..place[3]].to_vec(),
472 options,
473 func_vars.clone(),
474 )?
475 .num()?
476 .number
477 .real()
478 .clone(),
479 place.len() == 5,
480 )?;
481 function.drain(i + 1..=*place.last().unwrap());
482 }
483 ("extrema", Func(var)) if place.len() == 2 || place.len() == 3 => {
484 function[i] = extrema(
485 function[place[0] + 1..place[1]].to_vec(),
486 func_vars.clone(),
487 options,
488 var.to_string(),
489 if place.len() == 3 {
490 do_math(
491 function[place[1] + 1..place[2]].to_vec(),
492 options,
493 func_vars.clone(),
494 )?
495 .num()?
496 } else {
497 Number::from(Complex::new(options.prec), None)
498 },
499 )?;
500 function.drain(i + 1..=*place.last().unwrap());
501 }
502 ("solve", Func(var)) if place.len() == 2 || place.len() == 3 => {
503 function[i] = solve(
504 function[place[0] + 1..place[1]].to_vec(),
505 func_vars.clone(),
506 options,
507 var.to_string(),
508 if place.len() == 3 {
509 do_math(
510 function[place[1] + 1..place[2]].to_vec(),
511 options,
512 func_vars.clone(),
513 )?
514 .num()?
515 } else {
516 Number::from(Complex::new(options.prec), None)
517 },
518 )?;
519 function.drain(i + 1..=*place.last().unwrap());
520 }
521 ("set", Func(var)) if place.len() == 3 => {
522 function[i] = do_math_with_var(
523 function[place[0] + 1..place[1]].to_vec(),
524 options,
525 func_vars.clone(),
526 &var,
527 do_math(
528 function[place[1] + 1..place[2]].to_vec(),
529 options,
530 func_vars.clone(),
531 )?,
532 )?;
533 function.drain(i + 1..=*place.last().unwrap());
534 }
535 ("lim" | "limit", Func(var)) if place.len() == 3 || place.len() == 4 => {
536 function[i] = limit(
537 function[place[0] + 1..place[1]].to_vec(),
538 func_vars.clone(),
539 options,
540 var.to_string(),
541 do_math(
542 function[place[1] + 1..place[2]].to_vec(),
543 options,
544 func_vars.clone(),
545 )?
546 .num()?,
547 if place.len() == 4 {
548 match do_math(
549 function[place[2] + 1..place[3]].to_vec(),
550 options,
551 func_vars.clone(),
552 )?
553 .num()?
554 .number
555 .real()
556 .cmp0()
557 {
558 Some(Ordering::Less) => Left,
559 Some(Ordering::Greater) => Right,
560 _ => Both,
561 }
562 } else {
563 Both
564 },
565 )?;
566 function.drain(i + 1..=*place.last().unwrap());
567 }
568 ("surfacearea" | "sarea", Func(var)) if place.len() == 7 => {
569 if let Func(var2) = &function[place[0] + 1] {
570 function[i] = NumStr::new(surface_area(
571 function[place[1] + 1..place[2]].to_vec(),
572 func_vars.clone(),
573 options,
574 var.to_string(),
575 do_math(
576 function[place[4] + 1..place[5]].to_vec(),
577 options,
578 func_vars.clone(),
579 )?
580 .num()?
581 .number,
582 do_math(
583 function[place[5] + 1..place[6]].to_vec(),
584 options,
585 func_vars.clone(),
586 )?
587 .num()?,
588 var2.to_string(),
589 function[place[2] + 1..place[3]].to_vec(),
590 function[place[3] + 1..place[4]].to_vec(),
591 )?);
592 function.drain(i + 1..=*place.last().unwrap());
593 } else {
594 return Err("bad var");
595 }
596 }
597 ("length" | "arclength", Func(var)) if place.len() == 4 => {
598 function[i] = NumStr::new(length(
599 function[place[0] + 1..place[1]].to_vec(),
600 func_vars.clone(),
601 options,
602 var.to_string(),
603 do_math(
604 function[place[1] + 1..place[2]].to_vec(),
605 options,
606 func_vars.clone(),
607 )?
608 .num()?
609 .number,
610 do_math(
611 function[place[2] + 1..place[3]].to_vec(),
612 options,
613 func_vars.clone(),
614 )?
615 .num()?,
616 )?);
617 function.drain(i + 1..=*place.last().unwrap());
618 }
619 ("∫" | "area" | "integrate", Func(var))
620 if place.len() == 4 || place.len() == 5 || place.len() == 6 =>
621 {
622 function[i] = area(
623 function[place[0] + 1..place[1]].to_vec(),
624 func_vars.clone(),
625 options,
626 var.to_string(),
627 do_math(
628 function[place[1] + 1..place[2]].to_vec(),
629 options,
630 func_vars.clone(),
631 )?
632 .num()?
633 .number,
634 do_math(
635 function[place[2] + 1..place[3]].to_vec(),
636 options,
637 func_vars.clone(),
638 )?
639 .num()?,
640 if place.len() == 5 {
641 do_math(
642 function[place[3] + 1..place[4]].to_vec(),
643 options,
644 func_vars.clone(),
645 )?
646 .num()?
647 .number
648 } else {
649 Complex::with_val(options.prec, 1)
650 },
651 place.len() != 6,
652 )?;
653 function.drain(i + 1..=*place.last().unwrap());
654 }
655 ("isolate", Func(var)) if place.len() == 2 => {
656 function[i] = isolate(
657 &function[place[0] + 1..place[1]],
658 func_vars.clone(),
659 options,
660 var.to_string(),
661 )?;
662 function.drain(i + 1..=*place.last().unwrap());
663 }
664 ("slope" | "D", Func(var))
665 if place.len() == 3 || place.len() == 4 || place.len() == 5 =>
666 {
667 function[i] = slope(
668 function[place[0] + 1..place[1]].to_vec(),
669 func_vars.clone(),
670 options,
671 var.to_string(),
672 do_math(
673 function[place[1] + 1..place[2]].to_vec(),
674 options,
675 func_vars.clone(),
676 )?
677 .num()?,
678 place.len() != 5,
679 if place.len() >= 4 {
680 do_math(
681 function[place[2] + 1..place[3]].to_vec(),
682 options,
683 func_vars.clone(),
684 )?
685 .num()?
686 .number
687 .real()
688 .to_integer()
689 .unwrap_or_default()
690 .to_u32()
691 .unwrap_or_default()
692 } else {
693 1
694 },
695 )?;
696 function.drain(i + 1..=*place.last().unwrap());
697 }
698 ("taylor", Func(var)) if place.len() == 4 || place.len() == 5 => {
699 function[i] = taylor(
700 function[place[0] + 1..place[1]].to_vec(),
701 func_vars.clone(),
702 options,
703 var.to_string(),
704 if place.len() == 5 {
705 Some(
706 do_math(
707 function[place[3] + 1..place[4]].to_vec(),
708 options,
709 func_vars.clone(),
710 )?
711 .num()?,
712 )
713 } else {
714 None
715 },
716 do_math(
717 function[place[1] + 1..place[2]].to_vec(),
718 options,
719 func_vars.clone(),
720 )?
721 .num()?,
722 do_math(
723 function[place[2] + 1..place[3]].to_vec(),
724 options,
725 func_vars.clone(),
726 )?
727 .num()?
728 .number
729 .real()
730 .to_integer()
731 .unwrap_or_default()
732 .to_usize()
733 .unwrap_or_default(),
734 )?;
735 function.drain(i + 1..=*place.last().unwrap());
736 }
737 ("pw" | "piecewise", _) if !place.is_empty() => {
738 let mut ans = None;
739 let mut start = i + 3;
740 for (i, end) in place[0..if place.len() % 2 == 1 {
741 place.len()
742 } else {
743 place.len().saturating_sub(1)
744 }]
745 .iter()
746 .enumerate()
747 {
748 if i + 1 == place.len()
749 || (i % 2 == 0
750 && do_math(
751 function[*end + 1..place[i + 1] - 1].to_vec(),
752 options,
753 func_vars.clone(),
754 )?
755 .num()?
756 .number
757 .real()
758 == &1.0)
759 {
760 ans = Some(recursion(
761 func_vars.clone(),
762 function[if i + 1 == place.len() {
763 start.saturating_sub(1)
764 } else {
765 start
766 }..*end]
767 .to_vec(),
768 options,
769 )?);
770 break;
771 } else {
772 start = end + 2;
773 }
774 }
775 function[i] = if let Some(n) = ans {
776 n
777 } else {
778 NumStr::new(Number::from(
779 Complex::with_val(options.prec, Nan),
780 None,
781 ))
782 };
783 function.drain(i + 1..=*place.last().unwrap());
784 }
785 (
786 "sum" | "product" | "prod" | "summation" | "Σ" | "Π" | "vec" | "mat",
787 Func(var),
788 ) if place.len() == 4 => {
789 let start = do_math(
790 function[place[1] + 1..place[2]].to_vec(),
791 options,
792 func_vars.clone(),
793 )?
794 .num()?
795 .number;
796 let end = do_math(
797 function[place[2] + 1..place[3]].to_vec(),
798 options,
799 func_vars.clone(),
800 )?
801 .num()?
802 .number;
803 if !start.imag().is_zero() || !end.imag().is_zero() {
804 return Err("imag start/end");
805 }
806 if !start.real().clone().fract().is_zero()
807 || !end.real().clone().fract().is_zero()
808 {
809 return Err("fractional start/end");
810 }
811 let start = start.real();
812 let end = end.real();
813 function[i] = match s.as_str() {
814 "vec" | "mat" => mvec(
815 function[place[0] + 1..place[1]].to_vec(),
816 func_vars.clone(),
817 &var,
818 start
819 .to_integer()
820 .unwrap_or_default()
821 .to_isize()
822 .unwrap_or_default(),
823 end.to_integer()
824 .unwrap_or_default()
825 .to_isize()
826 .unwrap_or_default(),
827 s == "vec",
828 options,
829 )?,
830 _ => sum(
831 function[place[0] + 1..place[1]].to_vec(),
832 func_vars.clone(),
833 &var,
834 start.clone(),
835 end.clone(),
836 !(s == "sum" || s == "summation" || s == "Σ"),
837 options,
838 )?,
839 };
840 function.drain(i + 1..=place[3]);
841 }
842 ("sum" | "summation" | "Σ", _) if place.len() <= 1 => {
843 function[i] = match if place.is_empty() {
844 Ok(function.remove(i + 1).clone())
845 } else {
846 do_math(
847 function.drain(i + 1..=place[0]).collect::<Vec<NumStr>>(),
848 options,
849 func_vars.clone(),
850 )
851 } {
852 Ok(Num(a)) => Num(a),
853 Ok(Vector(a)) => NumStr::new(Number::from(
854 a.iter().fold(Complex::new(options.prec), |sum, val| {
855 sum + val.number.clone()
856 }),
857 None,
858 )),
859 Ok(Matrix(a)) => NumStr::new(Number::from(
860 a.iter()
861 .flatten()
862 .fold(Complex::new(options.prec), |sum, val| {
863 sum + val.number.clone()
864 }),
865 None,
866 )),
867 _ => return Err("sum err"),
868 };
869 }
870 ("product" | "prod" | "Π", _) if place.len() <= 1 => {
871 function[i] = match if place.is_empty() {
872 Ok(function.remove(i + 1).clone())
873 } else {
874 do_math(
875 function.drain(i + 1..=place[0]).collect::<Vec<NumStr>>(),
876 options,
877 func_vars.clone(),
878 )
879 } {
880 Ok(Num(a)) => Num(a),
881 Ok(Vector(a)) => NumStr::new(Number::from(
882 a.iter()
883 .fold(Complex::with_val(options.prec, 1), |sum, val| {
884 sum * val.number.clone()
885 }),
886 None,
887 )),
888 Ok(Matrix(a)) => NumStr::new(Number::from(
889 a.iter()
890 .flatten()
891 .fold(Complex::with_val(options.prec, 1), |sum, val| {
892 sum * val.number.clone()
893 }),
894 None,
895 )),
896 _ => return Err("prod err"),
897 };
898 }
899 (_, _) => return Err("arg/var err with sum/prod/vec/slope or similar"),
900 }
901 } else {
902 let arg = function.remove(i + 1);
903 function[i] = match arg.clone() {
904 Matrix(a) => match s.as_str() {
905 "plane" => {
906 if a.len() != 3 || a.iter().any(|a| a.len() != 3) {
907 return Err("dimensions too high");
908 }
909 let x1 = &a[0][0].number;
910 let y1 = &a[0][1].number;
911 let z1 = &a[0][2].number;
912 let x2 = &a[1][0].number;
913 let y2 = &a[1][1].number;
914 let z2 = &a[1][2].number;
915 let x3 = &a[2][0].number;
916 let y3 = &a[2][1].number;
917 let z3 = &a[2][2].number;
918 let t1 = z2 - z3.clone() * y2 / y3;
919 let t2 = x2 - x3.clone() * y2 / y3;
920 let t3 = t1 / t2.clone();
921 let c1 = z1 - x1 * t3.clone() - y1.clone() / y3 * (z3 - x3 * t3);
922 let c2 = 1 - x1 / t2.clone() + x1.clone() * y2 / y3 / t2.clone()
923 - y1.clone() / y3
924 + x3.clone() * y1 / y3 / t2.clone()
925 - x3.clone() * y1 * y2 / y3 / y3 / t2.clone();
926 let c: Complex = c1 / c2;
927 let b = (z3
928 - c.clone()
929 - x3 * (z2 - c.clone() - y2.clone() / y3 * (z3 - c.clone()))
930 / t2)
931 / y3;
932 let a = (z2 - c.clone() - b.clone() * y2) / x2;
933 if function.len() > i + 2 {
934 let x = function.remove(i + 1).num()?.number;
935 let y = function.remove(i + 1).num()?.number;
936 NumStr::new(Number::from(a * x + b * y + c, None))
937 } else {
938 Vector(vec![
939 Number::from(a, None),
940 Number::from(b, None),
941 Number::from(c, None),
942 ])
943 }
944 }
945 "lobf" | "lineofbestfit" => {
946 if a.is_empty() || a.iter().any(|a| a.len() != 2) {
947 return Err("dimensions too high");
948 }
949 let mut xsum = Complex::new(options.prec);
950 let mut ysum = Complex::new(options.prec);
951 let mut xxsum = Complex::new(options.prec);
952 let mut xysum = Complex::new(options.prec);
953 for row in &a {
954 let x = row[0].number.clone();
955 let y = row[1].number.clone();
956 xsum += x.clone();
957 ysum += y.clone();
958 xxsum += sqr(x.clone());
959 xysum += x * y;
960 }
961 let m: Complex = (a.len() * xysum - xsum.clone() * ysum.clone())
962 / (a.len() * xxsum - sqr(xsum.clone()));
963 let b = (ysum - m.clone() * xsum) / a.len();
964 if function.len() > i + 1 {
965 let x = function.remove(i + 1).num()?.number;
966 NumStr::new(Number::from(m * x + b, a[0][1].units))
967 } else {
968 Vector(vec![Number::from(m, None), Number::from(b, None)])
969 }
970 }
971 "inter" | "interpolate" => {
972 if function.len() > i + 1 {
973 if !a.is_empty() && a.iter().all(|a| a.len() == 2) {
974 let x = function.remove(i + 1).num()?.number;
975 let mut sum = Complex::new(options.prec);
976 for i in 0..a.len() {
977 let mut prod = Complex::with_val(options.prec, 1);
978 for j in 0..a.len() {
979 if j != i {
980 prod *= (x.clone() - a[j][0].number.clone())
981 / (a[i][0].number.clone()
982 - a[j][0].number.clone())
983 }
984 }
985 sum += prod * a[i][1].number.clone()
986 }
987 NumStr::new(Number::from(sum, a[0][1].units))
988 } else {
989 return Err("dimensions too high");
990 }
991 } else {
992 return Err("no x value given");
993 }
994 }
995 "sort" => Matrix(sort_mat(a, options.prec)),
996 "max" => {
997 let mut vec = Vec::with_capacity(a.len().max(1));
998 for j in a {
999 let mut max = j[0].clone();
1000 for i in j {
1001 if i.number.real() > max.number.real() {
1002 max = i
1003 }
1004 }
1005 vec.push(max)
1006 }
1007 Vector(vec)
1008 }
1009 "min" => {
1010 let mut vec = Vec::with_capacity(a.len().max(1));
1011 for j in a {
1012 let mut min = j[0].clone();
1013 for i in j {
1014 if i.number.real() < min.number.real() {
1015 min = i
1016 }
1017 }
1018 vec.push(min)
1019 }
1020 Vector(vec)
1021 }
1022 "flatten" => Vector(a.into_iter().flatten().collect::<Vec<Number>>()),
1023 "cofactor" | "cofactors" | "cof" => Matrix(cofactor(&a)?),
1024 "minor" | "minors" => Matrix(minors(&a)?),
1025 "adjugate" | "adj" => Matrix(transpose(&cofactor(&a)?)),
1026 "inverse" | "inv" => Matrix(inverse(&a)?),
1027 "transpose" | "trans" => Matrix(transpose(&a)),
1028 "len" => NumStr::new(Number::from(
1029 Complex::with_val(options.prec, a.len()),
1030 None,
1031 )),
1032 "wid" | "width" => NumStr::new(Number::from(
1033 Complex::with_val(options.prec, a[0].len()),
1034 None,
1035 )),
1036 "tr" | "trace" => NumStr::new(trace(&a)),
1037 "det" | "norm" | "determinant" => NumStr::new(determinant(&a)?),
1038 "part" => {
1039 if function.len() > i + 2 {
1040 match (function.remove(i + 1), function.remove(i + 1)) {
1041 (Num(b), Num(c)) => {
1042 let b = b.number;
1043 let c = c.number;
1044 let n1 =
1045 b.clone().real().to_integer().unwrap_or_default();
1046 let getcol = n1 == -1;
1047 let n1 = n1.to_usize().unwrap_or_default();
1048 let n2 = c
1049 .clone()
1050 .real()
1051 .to_integer()
1052 .unwrap_or_default()
1053 .to_usize()
1054 .unwrap_or_default();
1055 if getcol {
1056 if a.iter().all(|a| n2 < a.len()) {
1057 Vector(
1058 a.iter()
1059 .map(|a| a[n2].clone())
1060 .collect::<Vec<Number>>(),
1061 )
1062 } else {
1063 return Err("out of range");
1064 }
1065 } else if n1 < a.len() && n2 < a[n1].len() {
1066 NumStr::new(a[n1][n2].clone())
1067 } else {
1068 return Err("not in matrix");
1069 }
1070 }
1071 (Num(b), Vector(c)) => {
1072 let b = b.number;
1073 let n1 =
1074 b.clone().real().to_integer().unwrap_or_default();
1075 let getcol = n1 == -1;
1076 let n1 = n1.to_usize().unwrap_or_default();
1077 if getcol {
1078 let mut mat = Vec::with_capacity(c.len().max(1));
1079 for n in c {
1080 let n = n
1081 .number
1082 .clone()
1083 .real()
1084 .to_integer()
1085 .unwrap_or_default()
1086 .to_usize()
1087 .unwrap_or_default();
1088 if a.iter().all(|a| n < a.len()) {
1089 mat.push(
1090 a.iter()
1091 .map(|a| a[n].clone())
1092 .collect::<Vec<Number>>(),
1093 )
1094 } else {
1095 return Err("out of range");
1096 }
1097 }
1098 Matrix(transpose(&mat))
1099 } else {
1100 let mut vec = Vec::with_capacity(c.len().max(1));
1101 for n in c {
1102 let n2 = n
1103 .number
1104 .clone()
1105 .real()
1106 .to_integer()
1107 .unwrap_or_default()
1108 .to_usize()
1109 .unwrap_or_default();
1110 if n1 < a.len() && n2 < a[n1].len() {
1111 vec.push(a[n1][n2].clone())
1112 } else {
1113 return Err("not in matrix");
1114 }
1115 }
1116 Vector(vec)
1117 }
1118 }
1119 (Vector(b), Num(c)) => {
1120 let c = c.number;
1121 let n2 = c
1122 .clone()
1123 .real()
1124 .to_integer()
1125 .unwrap_or_default()
1126 .to_usize()
1127 .unwrap_or_default();
1128 let mut vec = Vec::with_capacity(b.len().max(1));
1129 for n in b {
1130 let n1 = n
1131 .number
1132 .clone()
1133 .real()
1134 .to_integer()
1135 .unwrap_or_default()
1136 .to_usize()
1137 .unwrap_or_default();
1138 if n1 < a.len() && n2 < a[n1].len() {
1139 vec.push(a[n1][n2].clone())
1140 } else {
1141 return Err("not in matrix");
1142 }
1143 }
1144 Vector(vec)
1145 }
1146 (Vector(b), Vector(c)) => {
1147 let mut mat = Vec::with_capacity(b.len().max(1));
1148 for g in b {
1149 let mut vec = Vec::with_capacity(c.len().max(1));
1150 let n1 = g
1151 .number
1152 .clone()
1153 .real()
1154 .to_integer()
1155 .unwrap_or_default()
1156 .to_usize()
1157 .unwrap_or_default();
1158 for n in c.clone() {
1159 let n2 = n
1160 .number
1161 .clone()
1162 .real()
1163 .to_integer()
1164 .unwrap_or_default()
1165 .to_usize()
1166 .unwrap_or_default();
1167 if n1 < a.len() && n2 < a[n1].len() {
1168 vec.push(a[n1][n2].clone())
1169 } else {
1170 return Err("not in matrix");
1171 }
1172 }
1173 mat.push(vec);
1174 }
1175 Matrix(mat)
1176 }
1177 _ => return Err("wrong part num"),
1178 }
1179 } else if function.len() > i + 1 {
1180 match function.remove(i + 1) {
1181 Num(b) => {
1182 let b = b.number;
1183 let n = b
1184 .clone()
1185 .real()
1186 .to_integer()
1187 .unwrap_or_default()
1188 .to_usize()
1189 .unwrap_or_default();
1190 if n < a.len() {
1191 Vector(a[n].clone())
1192 } else {
1193 return Err("out of range");
1194 }
1195 }
1196 Vector(b) => {
1197 let mut vec = Vec::with_capacity(b.len().max(1));
1198 for i in b {
1199 let n = i
1200 .number
1201 .clone()
1202 .real()
1203 .to_integer()
1204 .unwrap_or_default()
1205 .to_usize()
1206 .unwrap_or_default();
1207 if n < a.len() {
1208 vec.push(a[n].clone());
1209 } else {
1210 return Err("out of range");
1211 }
1212 }
1213 Matrix(vec)
1214 }
1215 _ => return Err("non num/vec"),
1216 }
1217 } else {
1218 return Err("no arg");
1219 }
1220 }
1221 "weighted_mean" => {
1222 if a.iter().any(|a| a.len() != 2) {
1223 return Err("bad data");
1224 }
1225 NumStr::new(Number::from(
1226 a.iter().fold(Complex::new(options.prec), |sum, val| {
1227 sum + val[0].number.clone() * val[1].number.clone()
1228 }) / a.iter().fold(Complex::new(options.prec), |sum, val| {
1229 sum + val[1].number.clone()
1230 }),
1231 a[0][0].units,
1232 ))
1233 }
1234 "mean" | "μ" => NumStr::new(Number::from(
1235 a.iter()
1236 .flatten()
1237 .fold(Complex::new(options.prec), |sum, val| {
1238 sum + val.number.clone()
1239 })
1240 / a.iter().fold(0, |sum, a| sum + a.len()),
1241 a[0][0].units,
1242 )),
1243 "mode" => {
1244 let mut most = (Vec::with_capacity(a.len().max(1)), 0);
1245 for i in a.iter().flatten() {
1246 let mut count = 0;
1247 for j in a.iter().flatten() {
1248 if i == j {
1249 count += 1;
1250 }
1251 }
1252 if count > most.1 {
1253 most.0.clear();
1254 most.0.push(i.clone());
1255 most.1 = count;
1256 }
1257 if count == most.1 && !most.0.iter().any(|j| i == j) {
1258 most.0.push(i.clone())
1259 }
1260 }
1261 if most.0.len() == 1 {
1262 NumStr::new(most.0[0].clone())
1263 } else {
1264 Vector(most.0)
1265 }
1266 }
1267 "median" => {
1268 let a = sort(a.iter().flatten().cloned().collect::<Vec<Number>>());
1269 if a.len() % 2 == 0 {
1270 Vector(vec![a[a.len() / 2 - 1].clone(), a[a.len() / 2].clone()])
1271 } else {
1272 NumStr::new(a[a.len() / 2].clone())
1273 }
1274 }
1275 "all" => {
1276 let mut res = true;
1277 for a in a.iter().flatten() {
1278 if !(a.number.imag().is_zero() && a.number.real() == &1) {
1279 res = false
1280 }
1281 }
1282 NumStr::new(Number::from(
1283 Complex::with_val(options.prec, res as u8),
1284 None,
1285 ))
1286 }
1287 "any" => {
1288 let mut res = false;
1289 for a in a.iter().flatten() {
1290 if a.number.imag().is_zero() && a.number.real() == &1 {
1291 res = true
1292 }
1293 }
1294 NumStr::new(Number::from(
1295 Complex::with_val(options.prec, res as u8),
1296 None,
1297 ))
1298 }
1299 "eigenvalues" => {
1300 if function.len() > i + 1 && !matches!(&function[i + 1], Func(_)) {
1301 function.remove(i + 1);
1302 eigenvalues(&a, true)?
1303 } else {
1304 eigenvalues(&a, false)?
1305 }
1306 }
1307 "eigenvectors" => {
1308 if function.len() > i + 1 && !matches!(&function[i + 1], Func(_)) {
1309 function.remove(i + 1);
1310 eigenvectors(&a, true)?
1311 } else {
1312 eigenvectors(&a, false)?
1313 }
1314 }
1315 "rref" => Matrix(rref(a)?),
1316 "ker" => Matrix(kernel(a)?),
1317 "ran" => Matrix(range(a)?),
1318 "generalized_eigenvectors" => {
1319 if function.len() > i + 1 && !matches!(&function[i + 1], Func(_)) {
1320 function.remove(i + 1);
1321 generalized_eigenvectors(&a, true)?
1322 } else {
1323 generalized_eigenvectors(&a, false)?
1324 }
1325 }
1326 "rcf" => rcf(a)?,
1327 "jcf" => jcf(a)?,
1328 "change_basis" => {
1329 if function.len() > i + 1 && !matches!(&function[i + 1], Func(_)) {
1330 let beta = function.remove(i + 1).mat()?;
1331 if function.len() > i + 1
1332 && !matches!(&function[i + 1], Func(_))
1333 {
1334 change_basis(a, &beta, &function.remove(i + 1).mat()?)?
1335 } else {
1336 let i = identity(a.len(), a[0][0].number.prec().0);
1337 change_basis(a, &i, &beta)?
1338 }
1339 } else {
1340 return Err("missing arg");
1341 }
1342 }
1343 "coordinate" => {
1344 if function.len() > i + 1 && !matches!(&function[i + 1], Func(_)) {
1345 let v = function.remove(i + 1).vec()?;
1346 coordinate(v, a)?
1347 } else {
1348 return Err("missing arg");
1349 }
1350 }
1351 "rank" => NumStr::new(Number::from(
1352 Complex::with_val(a[0][0].number.prec(), range(a)?.len()),
1353 None,
1354 )),
1355 "null" => NumStr::new(Number::from(
1356 Complex::with_val(a[0][0].number.prec(), kernel(a)?.len()),
1357 None,
1358 )),
1359 "to_list" => {
1360 let mut vec = Vec::new();
1361 for a in a {
1362 if a.len() != 2 {
1363 return Err("bad list");
1364 }
1365 for _ in 0..a[1]
1366 .number
1367 .real()
1368 .to_integer()
1369 .unwrap_or_default()
1370 .to_usize()
1371 .unwrap_or_default()
1372 {
1373 vec.push(a[0].clone())
1374 }
1375 }
1376 if vec.is_empty() {
1377 return Err("bad list");
1378 }
1379 Vector(sort(vec))
1380 }
1381 "to_freq" => {
1382 if a.is_empty() {
1383 return Err("bad list");
1384 }
1385 let mut a = sort(a.into_iter().flatten().collect::<Vec<Number>>());
1386 let mut last = a[0].clone();
1387 let mut count = 1;
1388 a.remove(0);
1389 let mut mat = Vec::new();
1390 for a in a {
1391 if a != last {
1392 mat.push(vec![
1393 last.clone(),
1394 Number::from(
1395 Complex::with_val(options.prec, count),
1396 None,
1397 ),
1398 ]);
1399 last = a;
1400 count = 0;
1401 }
1402 count += 1;
1403 }
1404 mat.push(vec![
1405 last.clone(),
1406 Number::from(Complex::with_val(options.prec, count), None),
1407 ]);
1408 Matrix(mat)
1409 }
1410 #[cfg(feature = "fastrand")]
1411 "rand_weighted" => {
1412 if a.iter().any(|a| a.len() != 2) {
1413 return Err("bad data");
1414 }
1415 let mut sum = Integer::new();
1416 for i in &a {
1417 sum += i[1].number.real().to_integer().unwrap_or_default();
1418 }
1419 let n = sum.to_u128().unwrap_or_default();
1420 if n == 0 {
1421 return Err("bad data");
1422 }
1423 let max = u128::MAX - u128::MAX.rem(n);
1424 let mut rnd = u128::MAX;
1425 while rnd >= max {
1426 rnd = fastrand::u128(..);
1427 }
1428 rnd = rnd.rem(n) + 1;
1429 let mut num =
1430 Number::from(Complex::with_val(options.prec, Nan), None);
1431 for i in &a {
1432 rnd = rnd.saturating_sub(
1433 i[1].number
1434 .real()
1435 .to_integer()
1436 .unwrap_or_default()
1437 .to_u128()
1438 .unwrap_or_default(),
1439 );
1440 if rnd == 0 {
1441 num = i[0].clone();
1442 break;
1443 }
1444 }
1445 NumStr::new(num)
1446 }
1447 #[cfg(feature = "fastrand")]
1448 "roll" => {
1449 let mut sum: Integer = Integer::new();
1450 for i in a {
1451 if i.len() != 2 {
1452 return Err("bad dice data");
1453 }
1454 let a = i[0].number.real().to_integer().unwrap_or_default();
1455 if a > u128::MAX || a == 0 {
1456 return Err("dice too large or bad dice data");
1457 }
1458 let n = a.to_u128().unwrap_or_default();
1459 let max = u128::MAX - u128::MAX.rem(n);
1460 let end = i[1]
1461 .number
1462 .real()
1463 .to_integer()
1464 .unwrap_or_default()
1465 .to_u128()
1466 .unwrap_or_default();
1467 let mut i = 0;
1468 while i < end {
1469 let rnd = fastrand::u128(..);
1470 if rnd < max {
1471 sum += rnd.rem(n) + 1;
1472 i += 1;
1473 }
1474 }
1475 }
1476 NumStr::new(Number::from(
1477 Complex::with_val(options.prec, sum),
1478 None,
1479 ))
1480 }
1481 "dice" => {
1482 let mut faces = Vec::new();
1483 for a in a {
1484 if a.len() != 2 {
1485 return Err("bad list");
1486 }
1487 for _ in 0..a[1]
1488 .number
1489 .real()
1490 .to_integer()
1491 .unwrap_or_default()
1492 .to_usize()
1493 .unwrap_or_default()
1494 {
1495 faces.push(
1496 a[0].number
1497 .real()
1498 .to_integer()
1499 .unwrap_or_default()
1500 .to_usize()
1501 .unwrap_or_default(),
1502 )
1503 }
1504 }
1505 if faces.iter().any(|c| c == &0) {
1506 return Err("bad face value");
1507 }
1508 if faces.is_empty() {
1509 Vector(vec![Number::from(
1510 Complex::with_val(options.prec, 1),
1511 None,
1512 )])
1513 } else if faces.len() == 1 {
1514 Vector(vec![
1515 Number::from(
1516 Complex::with_val(options.prec, 1),
1517 None
1518 );
1519 faces[0]
1520 ])
1521 } else {
1522 let mut last = vec![Integer::from(1); faces[0]];
1523 let mut current = last.clone();
1524 for i in 1..faces.len() {
1525 current = Vec::new();
1526 for p in 0..=faces[0..=i].iter().sum::<usize>() - i - 1 {
1527 let value = last[(p + 1).saturating_sub(faces[i])
1528 ..=p.min(faces[0..i].iter().sum::<usize>() - i)]
1529 .iter()
1530 .sum::<Integer>();
1531 current.push(value)
1532 }
1533 last.clone_from(¤t)
1534 }
1535 Vector(
1536 current
1537 .iter()
1538 .map(|a| {
1539 Number::from(
1540 Complex::with_val(options.prec, a),
1541 None,
1542 )
1543 })
1544 .collect::<Vec<Number>>(),
1545 )
1546 }
1547 }
1548 "poly" | "polynomial" => {
1549 if i + 1 < function.len() {
1550 if a.is_empty() {
1551 NumStr::new(Number::from(Complex::new(options.prec), None))
1552 } else {
1553 let x = function.remove(i + 1).num()?;
1554 let mut sum =
1555 vec![
1556 Number::from(Complex::new(options.prec), None);
1557 a.len()
1558 ];
1559 for (i, v) in transpose(&a).iter().rev().enumerate() {
1560 for (s, a) in sum.iter_mut().zip(v.iter().map(|a| {
1561 Number::from(
1562 a.number.clone() * x.number.clone().pow(i),
1563 mul_units(
1564 a.units,
1565 Some(
1566 x.units
1567 .unwrap_or_default()
1568 .pow(i as f64),
1569 ),
1570 ),
1571 )
1572 })) {
1573 *s = add(s, &a)
1574 }
1575 }
1576 Vector(sum)
1577 }
1578 } else {
1579 return Err("not enough args");
1580 }
1581 }
1582 "norm_combine" => {
1583 if a.is_empty() {
1584 return Err("empty vec");
1585 }
1586 let mut mu = Complex::new(options.prec);
1587 let mut var = Complex::new(options.prec);
1588 let mut ws = Complex::new(options.prec);
1589 for d in &a {
1590 if d.len() < 2 {
1591 return Err("not enough data");
1592 }
1593 if d.len() == 2 {
1594 mu += d[0].number.clone();
1595 var += d[1].number.clone().pow(2);
1596 ws += 1;
1597 } else {
1598 mu += d[0].number.clone() * d[2].number.clone();
1599 var += d[1].number.clone().pow(2) * d[2].number.clone();
1600 ws += d[2].number.clone();
1601 }
1602 }
1603 mu /= ws.clone();
1604 for d in &a {
1605 if d.len() == 2 {
1606 var += (mu.clone() - d[0].number.clone()).pow(2);
1607 } else {
1608 var += (mu.clone() - d[0].number.clone()).pow(2)
1609 * d[2].number.clone();
1610 }
1611 }
1612 var /= ws;
1613 Vector(vec![Number::from(mu, None), Number::from(var.sqrt(), None)])
1614 }
1615 _ => do_functions(arg, options, &mut function, i, &to_deg, s)?,
1616 },
1617 Vector(a) => match s.as_str() {
1618 "transpose" | "trans" => Matrix(transpose(&[a])),
1619 "to_list" => {
1620 let mut vec = Vec::new();
1621 for (i, a) in a.iter().enumerate() {
1622 let num =
1623 Number::from(Complex::with_val(options.prec, i + 1), None);
1624 for _ in 1..=a
1625 .number
1626 .real()
1627 .to_integer()
1628 .unwrap_or_default()
1629 .to_usize()
1630 .unwrap_or_default()
1631 {
1632 vec.push(num.clone())
1633 }
1634 }
1635 if vec.is_empty() {
1636 return Err("bad list");
1637 }
1638 Vector(sort(vec))
1639 }
1640 #[cfg(feature = "fastrand")]
1641 "roll" => {
1642 let mut sum: Integer = Integer::new();
1643 let mut i = 0;
1644 while i < a.len() {
1645 let a = a[i].number.real().to_integer().unwrap_or_default();
1646 if a > u128::MAX as f64 || a == 0 {
1647 return Err("dice too large or bad dice data");
1648 }
1649 let n = a.to_u128().unwrap_or_default();
1650 let max = u128::MAX - u128::MAX.rem(n);
1651 let rnd = fastrand::u128(..);
1652 if rnd < max {
1653 sum += rnd.rem(n) + 1;
1654 i += 1;
1655 }
1656 }
1657 NumStr::new(Number::from(
1658 Complex::with_val(options.prec, sum),
1659 None,
1660 ))
1661 }
1662 "dice" => {
1663 let faces = a
1664 .iter()
1665 .map(|c| {
1666 c.number
1667 .real()
1668 .to_integer()
1669 .unwrap_or_default()
1670 .to_usize()
1671 .unwrap_or_default()
1672 })
1673 .collect::<Vec<usize>>();
1674 if faces.iter().any(|c| c == &0) {
1675 return Err("bad face value");
1676 }
1677 let mut last = vec![Integer::from(1); faces[0]];
1678 if a.is_empty() {
1679 Vector(vec![Number::from(
1680 Complex::with_val(options.prec, 1),
1681 None,
1682 )])
1683 } else if faces.len() == 1 {
1684 Vector(
1685 last.iter()
1686 .map(|a| {
1687 Number::from(
1688 Complex::with_val(options.prec, a),
1689 None,
1690 )
1691 })
1692 .collect::<Vec<Number>>(),
1693 )
1694 } else {
1695 let mut current = Vec::new();
1696 for i in 1..faces.len() {
1697 current = Vec::new();
1698 for p in 0..=faces[0..=i].iter().sum::<usize>() - i - 1 {
1699 let value = last[(p + 1).saturating_sub(faces[i])
1700 ..=p.min(faces[0..i].iter().sum::<usize>() - i)]
1701 .iter()
1702 .sum::<Integer>();
1703 current.push(value)
1704 }
1705 last.clone_from(¤t);
1706 }
1707 Vector(
1708 current
1709 .iter()
1710 .map(|a| {
1711 Number::from(
1712 Complex::with_val(options.prec, a),
1713 None,
1714 )
1715 })
1716 .collect::<Vec<Number>>(),
1717 )
1718 }
1719 }
1720 "quartiles" => {
1721 if a.len() < 2 {
1722 return Err("not enough data");
1723 }
1724 let units = a[0].units;
1725 let a = sort(a);
1726 let half1 = &a[0..a.len() / 2];
1727 let half2 = if a.len() % 2 == 0 {
1728 &a[a.len() / 2..a.len()]
1729 } else {
1730 &a[a.len() / 2 + 1..a.len()]
1731 };
1732 if half1.len() % 2 == 0 {
1733 Vector(vec![
1734 Number::from(
1735 (half1[half1.len() / 2 - 1].number.clone()
1736 + half1[half1.len() / 2].number.clone())
1737 / 2,
1738 units,
1739 ),
1740 Number::from(
1741 if a.len() % 2 == 0 {
1742 (half1[half1.len().saturating_sub(1)]
1743 .number
1744 .clone()
1745 + half2[0].number.clone())
1746 / 2
1747 } else {
1748 a[a.len() / 2].number.clone()
1749 },
1750 units,
1751 ),
1752 Number::from(
1753 (half2[half2.len() / 2 - 1].number.clone()
1754 + half2[half2.len() / 2].number.clone())
1755 / 2,
1756 units,
1757 ),
1758 ])
1759 } else {
1760 Vector(vec![
1761 Number::from(half1[half1.len() / 2].number.clone(), units),
1762 Number::from(
1763 if a.len() % 2 == 0 {
1764 (half1[half1.len().saturating_sub(1)]
1765 .number
1766 .clone()
1767 + half2[0].number.clone())
1768 / 2
1769 } else {
1770 a[a.len() / 2].number.clone()
1771 },
1772 units,
1773 ),
1774 Number::from(half2[half2.len() / 2].number.clone(), units),
1775 ])
1776 }
1777 }
1778 "percentile" => {
1779 if function.len() <= i + 1 {
1780 return Err("not enough input");
1781 }
1782 let b = function.remove(i + 1).num()?.number;
1783 let r: Float = (b.real().clone() / 100) * a.len();
1784 let r = r
1785 .ceil()
1786 .to_integer()
1787 .unwrap_or_default()
1788 .to_usize()
1789 .unwrap_or_default();
1790 if r > a.len() {
1791 return Err("bad input");
1792 }
1793 NumStr::new(sort(a)[r.saturating_sub(1)].clone())
1794 }
1795 "percentilerank" => {
1796 if function.len() <= i + 1 {
1797 return Err("not enough input");
1798 }
1799 let mut cf = 0;
1800 let mut f = 0;
1801 let b = function.remove(i + 1).num()?.number;
1802 for a in sort(a.clone()) {
1803 if a.number.real() < b.real() {
1804 cf += 1;
1805 } else if a.number == b {
1806 f += 1;
1807 } else {
1808 break;
1809 }
1810 }
1811 NumStr::new(Number::from(
1812 100 * (cf + Complex::with_val(options.prec, f) / 2) / a.len(),
1813 None,
1814 ))
1815 }
1816 "to_freq" => {
1817 if a.is_empty() {
1818 return Err("bad list");
1819 }
1820 let mut a = sort(a);
1821 let mut last = a[0].clone();
1822 let mut count = 1;
1823 a.remove(0);
1824 let mut mat = Vec::new();
1825 for a in a {
1826 if a != last {
1827 mat.push(vec![
1828 last.clone(),
1829 Number::from(
1830 Complex::with_val(options.prec, count),
1831 None,
1832 ),
1833 ]);
1834 last = a;
1835 count = 0;
1836 }
1837 count += 1;
1838 }
1839 mat.push(vec![
1840 last.clone(),
1841 Number::from(Complex::with_val(options.prec, count), None),
1842 ]);
1843 Matrix(mat)
1844 }
1845 "kurtosis" => {
1846 fn quar(z: Complex) -> Complex {
1847 if z.imag().is_zero() {
1848 z.pow(4)
1849 } else {
1850 z.clone() * &z * &z * z
1851 }
1852 }
1853 let mean = a.iter().fold(Complex::new(options.prec), |sum, val| {
1854 sum + val.number.clone()
1855 }) / a.len();
1856 NumStr::new(Number::from(
1857 a.iter().fold(Complex::new(options.prec), |sum, val| {
1858 sum + quar(val.number.clone() - mean.clone())
1859 }) / a.len()
1860 / sqr(variance(&a, Some(mean.clone()), options.prec).number)
1861 - 3,
1862 None,
1863 ))
1864 }
1865 "skew" | "skewness" => {
1866 let mean = a.iter().fold(Complex::new(options.prec), |sum, val| {
1867 sum + val.number.clone()
1868 }) / a.len();
1869 NumStr::new(Number::from(
1870 a.iter().fold(Complex::new(options.prec), |sum, val| {
1871 sum + cube(val.number.clone() - mean.clone())
1872 }) / a.len()
1873 / cube(
1874 variance(&a, Some(mean.clone()), options.prec)
1875 .number
1876 .sqrt(),
1877 ),
1878 None,
1879 ))
1880 }
1881 "sd" | "standarddeviation" | "σ" => NumStr::new(Number::from(
1882 variance(&a, None, options.prec).number.sqrt(),
1883 a[0].units,
1884 )),
1885 "variance" | "var" => NumStr::new(variance(&a, None, options.prec)),
1886 "covariance" | "cov" => {
1887 let mut sum = Complex::new(options.prec);
1888 if function.len() <= i + 1 {
1889 return Err("not enough input");
1890 }
1891 let b = function.remove(i + 1).vec()?;
1892 if a.len() != b.len() {
1893 return Err("different sized data sets");
1894 }
1895 let ma = a.iter().fold(Complex::new(options.prec), |sum, val| {
1896 sum + val.number.clone()
1897 }) / a.len();
1898 let mb = b.iter().fold(Complex::new(options.prec), |sum, val| {
1899 sum + val.number.clone()
1900 }) / b.len();
1901 for (a, b) in a.iter().zip(b.iter()) {
1902 sum += (a.number.clone() - ma.clone())
1903 * (b.number.clone() - mb.clone());
1904 }
1905 NumStr::new(Number::from(
1906 sum / (a.len() - 1),
1907 mul_units(a[0].units, b[0].units),
1908 ))
1909 }
1910 "all" => {
1911 let mut res = true;
1912 for a in a {
1913 if !(a.number.imag().is_zero() && a.number.real() == &1) {
1914 res = false
1915 }
1916 }
1917 NumStr::new(Number::from(
1918 Complex::with_val(options.prec, res as u8),
1919 None,
1920 ))
1921 }
1922 "any" => {
1923 let mut res = false;
1924 for a in a {
1925 if a.number.imag().is_zero() && a.number.real() == &1 {
1926 res = true
1927 }
1928 }
1929 NumStr::new(Number::from(
1930 Complex::with_val(options.prec, res as u8),
1931 None,
1932 ))
1933 }
1934 "sort" => Vector(sort(a)),
1935 "hsv_to_rgb" => {
1936 if a.len() == 3 {
1937 Vector(hsv2rgb(
1938 a[0].number.real().clone(),
1939 a[1].number.real().clone(),
1940 a[2].number.real().clone(),
1941 ))
1942 } else {
1943 return Err("not 3 length");
1944 }
1945 }
1946 "mean" | "μ" => NumStr::new(Number::from(
1947 a.iter().fold(Complex::new(options.prec), |sum, val| {
1948 sum + val.number.clone()
1949 }) / a.len(),
1950 a[0].units,
1951 )),
1952 "geo_mean" => NumStr::new(Number::from(
1953 pow_nth(
1954 a.iter()
1955 .fold(Complex::with_val(options.prec, 1), |sum, val| {
1956 sum * val.number.clone()
1957 }),
1958 Complex::with_val(options.prec, 1) / a.len(),
1959 ),
1960 a[0].units,
1961 )),
1962 "median" => {
1963 let a = sort(a);
1964 if a.len() % 2 == 0 {
1965 NumStr::new(Number::from(
1966 (a[a.len() / 2 - 1].number.clone()
1967 + a[a.len() / 2].number.clone())
1968 / 2,
1969 a[0].units,
1970 ))
1971 } else {
1972 NumStr::new(Number::from(
1973 a[a.len() / 2].number.clone(),
1974 a[0].units,
1975 ))
1976 }
1977 }
1978 "mode" => {
1979 let mut most = (Vec::new(), 0);
1980 for i in &a {
1981 let mut count = 0;
1982 for j in &a {
1983 if i == j {
1984 count += 1;
1985 }
1986 }
1987 if count > most.1 {
1988 most = (vec![i.clone()], count);
1989 }
1990 if count == most.1 && !most.0.iter().any(|j| i == j) {
1991 most.0.push(i.clone())
1992 }
1993 }
1994 if most.0.len() == 1 {
1995 NumStr::new(most.0[0].clone())
1996 } else {
1997 Vector(most.0)
1998 }
1999 }
2000 "max" => {
2001 let mut max = a[0].clone();
2002 for i in a {
2003 if i.number.real() > max.number.real() {
2004 max = i
2005 }
2006 }
2007 NumStr::new(max)
2008 }
2009 "min" => {
2010 let mut min = a[0].clone();
2011 for i in a {
2012 if i.number.real() < min.number.real() {
2013 min = i
2014 }
2015 }
2016 NumStr::new(min)
2017 }
2018 "reverse" => Vector(a.iter().rev().cloned().collect()),
2019 "len" => NumStr::new(Number::from(
2020 Complex::with_val(options.prec, a.len()),
2021 None,
2022 )),
2023 "norm" => {
2024 let units = a[0].units;
2025 let mut n = Complex::new(options.prec);
2026 for i in a {
2027 n += sqr(i.number.abs());
2028 }
2029 NumStr::new(Number::from(n.sqrt(), units))
2030 }
2031 "normalize" => {
2032 let mut n = Complex::new(options.prec);
2033 for i in a.clone() {
2034 n += sqr(i.number);
2035 }
2036 Vector(
2037 a.iter()
2038 .map(|x| {
2039 Number::from(x.number.clone() / n.clone().sqrt(), None)
2040 })
2041 .collect(),
2042 )
2043 }
2044 "car" | "cartesian" => {
2045 if a.len() == 2 {
2046 let t = a[1].number.clone() / to_deg.clone();
2047 Vector(vec![
2048 Number::from(
2049 a[0].number.clone() * t.clone().cos(),
2050 a[0].units,
2051 ),
2052 Number::from(
2053 a[0].number.clone() * t.clone().sin(),
2054 a[0].units,
2055 ),
2056 ])
2057 } else if a.len() == 3 {
2058 let t1 = a[1].number.clone() / to_deg.clone();
2059 let t2 = a[2].number.clone() / to_deg.clone();
2060 Vector(vec![
2061 Number::from(
2062 a[0].number.clone()
2063 * t1.clone().sin()
2064 * t2.clone().cos(),
2065 a[0].units,
2066 ),
2067 Number::from(
2068 a[0].number.clone()
2069 * t1.clone().sin()
2070 * t2.clone().sin(),
2071 a[0].units,
2072 ),
2073 Number::from(
2074 a[0].number.clone() * t1.clone().cos(),
2075 a[0].units,
2076 ),
2077 ])
2078 } else {
2079 return Err("incorrect polar form");
2080 }
2081 }
2082 "polar" | "pol" => Vector(to_polar(a.clone(), to_deg.clone())),
2083 "cyl" | "cylinder" => Vector(to_cyl(a.clone(), to_deg.clone())),
2084 "angle" => {
2085 if function.len() > i + 1 {
2086 let b = function.remove(i + 1).vec()?;
2087 if a.len() == 3 && b.len() == 3 {
2088 let c: Complex = sqr(a[0].number.clone())
2089 + sqr(a[1].number.clone())
2090 + sqr(a[2].number.clone());
2091 let d: Complex = sqr(b[0].number.clone())
2092 + sqr(b[1].number.clone())
2093 + sqr(b[2].number.clone());
2094 NumStr::new(Number::from(
2095 ((a[0].number.clone() * b[0].number.clone()
2096 + a[1].number.clone() * b[1].number.clone()
2097 + a[2].number.clone() * b[2].number.clone())
2098 / (c.sqrt() * d.sqrt()))
2099 .acos()
2100 * to_deg.clone(),
2101 Some(Units {
2102 angle: 1.0,
2103 ..Units::default()
2104 }),
2105 ))
2106 } else if a.len() == 2 && b.len() == 2 {
2107 let c: Complex =
2108 sqr(a[0].number.clone()) + sqr(a[1].number.clone());
2109 let d: Complex =
2110 sqr(b[0].number.clone()) + sqr(b[1].number.clone());
2111 NumStr::new(Number::from(
2112 ((a[0].number.clone() * b[0].number.clone()
2113 + a[1].number.clone() * b[1].number.clone())
2114 / (c.sqrt() * d.sqrt()))
2115 .acos()
2116 * to_deg.clone(),
2117 Some(Units {
2118 angle: 1.0,
2119 ..Units::default()
2120 }),
2121 ))
2122 } else {
2123 return Err("cant decern angles");
2124 }
2125 } else {
2126 let vec = to_polar(a, to_deg.clone());
2127 if vec.len() == 3 {
2128 Vector(vec[1..=2].to_vec())
2129 } else if vec.len() == 2 {
2130 NumStr::new(vec[1].clone())
2131 } else {
2132 return Err("cant decern angles");
2133 }
2134 }
2135 }
2136 "cross" => {
2137 if function.len() > i + 1 {
2138 let b = function.remove(i + 1).vec()?;
2139 let units = mul_units(a[0].units, b[0].units);
2140 if a.len() == 3 && b.len() == 3 {
2141 Vector(vec![
2142 Number::from(
2143 a[1].number.clone() * &b[2].number
2144 - a[2].number.clone() * &b[1].number,
2145 units,
2146 ),
2147 Number::from(
2148 a[2].number.clone() * &b[0].number
2149 - a[0].number.clone() * &b[2].number,
2150 units,
2151 ),
2152 Number::from(
2153 a[0].number.clone() * &b[1].number
2154 - a[1].number.clone() * &b[0].number,
2155 units,
2156 ),
2157 ])
2158 } else if a.len() == 2 && b.len() == 2 {
2159 NumStr::new(Number::from(
2160 a[0].number.clone() * &b[1].number
2161 - a[1].number.clone() * &b[0].number,
2162 units,
2163 ))
2164 } else {
2165 return Err("cant cross");
2166 }
2167 } else {
2168 return Err("no args");
2169 }
2170 }
2171 "project" | "proj" => {
2172 if function.len() > i + 1 {
2173 let b = function.remove(i + 1).clone().vec()?;
2174 if b.len() == a.len() {
2175 let mut dot = Complex::new(options.prec);
2176 for i in a
2177 .iter()
2178 .zip(b.iter())
2179 .map(|(a, b)| a.number.clone() * b.number.clone())
2180 {
2181 dot += i;
2182 }
2183 let mut norm = Complex::new(options.prec);
2184 for i in &b {
2185 norm += sqr(i.number.clone().abs());
2186 }
2187 Vector(
2188 b.iter()
2189 .map(|n| {
2190 Number::from(
2191 dot.clone() / norm.clone()
2192 * n.number.clone(),
2193 n.units,
2194 )
2195 })
2196 .collect::<Vec<Number>>(),
2197 )
2198 } else {
2199 return Err("cant project");
2200 }
2201 } else {
2202 return Err("no args");
2203 }
2204 }
2205 "oproject" | "oproj" => {
2206 if function.len() > i + 1 {
2207 let b = function.remove(i + 1).clone().vec()?;
2208 if b.len() == a.len() {
2209 let mut dot = Complex::new(options.prec);
2210 for i in a
2211 .iter()
2212 .zip(b.iter())
2213 .map(|(a, b)| a.number.clone() * b.number.clone())
2214 {
2215 dot += i;
2216 }
2217 let mut norm = Complex::new(options.prec);
2218 for i in &b {
2219 norm += sqr(i.number.clone().abs());
2220 }
2221 Vector(
2222 b.iter()
2223 .zip(a.iter())
2224 .map(|(n, a)| {
2225 Number::from(
2226 a.number.clone()
2227 - dot.clone() / norm.clone()
2228 * n.number.clone(),
2229 a.units,
2230 )
2231 })
2232 .collect::<Vec<Number>>(),
2233 )
2234 } else {
2235 return Err("cant project");
2236 }
2237 } else {
2238 return Err("no args");
2239 }
2240 }
2241 "dot" => {
2242 if function.len() > i + 1 {
2243 let mut n = Complex::new(options.prec);
2244 let b = function.remove(i + 1).vec()?;
2245 for i in a
2246 .iter()
2247 .zip(b.iter())
2248 .map(|(a, b)| a.number.clone() * b.number.clone())
2249 {
2250 n += i;
2251 }
2252 NumStr::new(Number::from(n, mul_units(a[0].units, b[0].units)))
2253 } else {
2254 return Err("no args");
2255 }
2256 }
2257 "part" => {
2258 if function.len() > i + 1 {
2259 match function.remove(i + 1) {
2260 Num(b) => {
2261 let b = b.number;
2262 let n = b
2263 .clone()
2264 .real()
2265 .to_integer()
2266 .unwrap_or_default()
2267 .to_usize()
2268 .unwrap_or_default();
2269 if n < a.len() {
2270 NumStr::new(a[n].clone())
2271 } else {
2272 return Err("out of range");
2273 }
2274 }
2275 Vector(b) => {
2276 let mut vec = Vec::new();
2277 for i in b {
2278 let n = i
2279 .number
2280 .clone()
2281 .real()
2282 .to_integer()
2283 .unwrap_or_default()
2284 .to_usize()
2285 .unwrap_or_default();
2286 if n < a.len() {
2287 vec.push(a[n].clone());
2288 } else {
2289 return Err("out of range");
2290 }
2291 }
2292 Vector(vec)
2293 }
2294 _ => return Err("non num/vec"),
2295 }
2296 } else {
2297 return Err("no args");
2298 }
2299 }
2300 "split" => Matrix(
2301 a.iter()
2302 .map(|a| {
2303 vec![
2304 Number::from(a.number.real().clone().into(), None),
2305 Number::from(a.number.imag().clone().into(), None),
2306 ]
2307 })
2308 .collect::<Vec<Vec<Number>>>(),
2309 ),
2310 "uniq" => {
2311 let mut a = a;
2312 a.dedup();
2313 Vector(a)
2314 }
2315 "factors" | "factor" => {
2316 let mut fail = false;
2317 let mut mat = Vec::new();
2318 for num in a {
2319 let num = num.number;
2320 if num.imag().clone().is_zero() {
2321 if num.real().clone().fract().is_zero() {
2322 let mut vec = Vec::new();
2323 let n = num
2324 .real()
2325 .to_integer()
2326 .unwrap_or_default()
2327 .to_usize()
2328 .unwrap_or_default();
2329 for i in 1..=n {
2330 if n % i == 0 {
2331 vec.push(Number::from(
2332 Complex::with_val(options.prec, i),
2333 None,
2334 ));
2335 }
2336 }
2337 mat.push(vec);
2338 } else {
2339 fail = true;
2340 break;
2341 }
2342 } else {
2343 fail = true;
2344 break;
2345 }
2346 }
2347 if fail {
2348 NumStr::new(Number::from(
2349 Complex::with_val(options.prec, Nan),
2350 None,
2351 ))
2352 } else {
2353 Matrix(mat)
2354 }
2355 }
2356 "union" => {
2357 if function.len() > i + 1 {
2358 if let Vector(b) = function.remove(i + 1) {
2359 let mut a = a;
2360 a.extend(b);
2361 a = sort(a);
2362 a.dedup();
2363 Vector(a)
2364 } else {
2365 return Err("arg not vector");
2366 }
2367 } else {
2368 return Err("arg not vector");
2369 }
2370 }
2371 "intersection" => {
2372 if function.len() > i + 1 {
2373 if let Vector(b) = function.remove(i + 1) {
2374 let mut v = Vec::new();
2375 'main: for n1 in a {
2376 for n2 in &b {
2377 if &n1 == n2 {
2378 v.push(n2.clone());
2379 continue 'main;
2380 }
2381 }
2382 }
2383 Vector(v)
2384 } else {
2385 return Err("arg not vector");
2386 }
2387 } else {
2388 return Err("arg not vector");
2389 }
2390 }
2391 "set_difference" => {
2392 if function.len() > i + 1 {
2393 if let Vector(b) = function.remove(i + 1) {
2394 let mut v = Vec::new();
2395 'main: for n1 in a {
2396 for n2 in &b {
2397 if &n1 == n2 {
2398 continue 'main;
2399 }
2400 }
2401 v.push(n1);
2402 }
2403 Vector(v)
2404 } else {
2405 return Err("arg not vector");
2406 }
2407 } else {
2408 return Err("arg not vector");
2409 }
2410 }
2411 "symmetric_difference" => {
2412 if function.len() > i + 1 {
2413 if let Vector(b) = function.remove(i + 1) {
2414 let mut a1 = Vec::new();
2415 'main: for n1 in a.clone() {
2416 for n2 in &b {
2417 if &n1 == n2 {
2418 continue 'main;
2419 }
2420 }
2421 a1.push(n1);
2422 }
2423 let mut a2 = Vec::new();
2424 'main: for n1 in b.clone() {
2425 for n2 in &a {
2426 if &n1 == n2 {
2427 continue 'main;
2428 }
2429 }
2430 a2.push(n1);
2431 }
2432 a2.extend(a1);
2433 let mut u = sort(a2);
2434 u.dedup();
2435 Vector(u)
2436 } else {
2437 return Err("arg not vector");
2438 }
2439 } else {
2440 return Err("arg not vector");
2441 }
2442 }
2443 "cartesian_product" => {
2444 if function.len() > i + 1 {
2445 if let Vector(b) = function.remove(i + 1) {
2446 let mut m = Vec::new();
2447 for n1 in a {
2448 for n2 in b.clone() {
2449 m.push(vec![n1.clone(), n2])
2450 }
2451 }
2452 Matrix(m)
2453 } else {
2454 return Err("arg not vector");
2455 }
2456 } else {
2457 return Err("arg not vector");
2458 }
2459 }
2460 "power_set" => {
2461 let mut m = Vec::new();
2462 for i in 0..1 << a.len() {
2463 let mut v = Vec::new();
2464 for (j, n) in a.iter().enumerate() {
2465 if (i >> j) & 1 == 1 {
2466 v.push(n.clone())
2467 }
2468 }
2469 m.push(v);
2470 }
2471 Matrix(m)
2472 }
2473 "set_fix" => {
2474 let mut a = sort(a);
2475 a.dedup();
2476 Vector(a)
2477 }
2478 "subset" => {
2479 if function.len() > i + 1 {
2480 if let Vector(b) = function.remove(i + 1) {
2481 let mut ainb = true;
2482 'main: for n1 in b {
2483 for n2 in &a {
2484 if &n1 == n2 {
2485 continue 'main;
2486 }
2487 }
2488 ainb = false;
2489 }
2490 NumStr::new(Number::from(
2491 Complex::with_val(options.prec, ainb as u8),
2492 None,
2493 ))
2494 } else {
2495 return Err("arg not vector");
2496 }
2497 } else {
2498 return Err("arg not vector");
2499 }
2500 }
2501 "element" => {
2502 if function.len() > i + 1 {
2503 if let Num(b) = function.remove(i + 1) {
2504 let mut ainb = false;
2505 for n2 in &a {
2506 if &*b == n2 {
2507 ainb = true;
2508 }
2509 }
2510 NumStr::new(Number::from(
2511 Complex::with_val(options.prec, ainb as u8),
2512 None,
2513 ))
2514 } else {
2515 return Err("arg not vector");
2516 }
2517 } else {
2518 return Err("arg not vector");
2519 }
2520 }
2521 "extend" | "link" => {
2522 if function.len() > i + 1 {
2523 match function.remove(i + 1) {
2524 Num(n) => {
2525 let mut a = a;
2526 a.push(*n);
2527 Vector(a)
2528 }
2529 Vector(v) => {
2530 let mut a = a;
2531 a.extend(v);
2532 Vector(a)
2533 }
2534 _ => return Err("bad arg"),
2535 }
2536 } else {
2537 return Err("bad arg");
2538 }
2539 }
2540 "remove" => {
2541 if function.len() > i + 1 {
2542 match function.remove(i + 1) {
2543 Num(n) => {
2544 let mut a = a;
2545 let n = n
2546 .number
2547 .real()
2548 .to_integer()
2549 .unwrap_or_default()
2550 .to_usize()
2551 .unwrap_or_default();
2552 if n >= a.len() {
2553 return Err("bad range");
2554 }
2555 a.remove(n);
2556 Vector(a)
2557 }
2558 Vector(v) => {
2559 let mut a = a;
2560 for n in v {
2561 let n = n
2562 .number
2563 .real()
2564 .to_integer()
2565 .unwrap_or_default()
2566 .to_usize()
2567 .unwrap_or_default();
2568 if n >= a.len() {
2569 return Err("bad range");
2570 }
2571 a.remove(n);
2572 }
2573 Vector(a)
2574 }
2575 _ => return Err("bad range"),
2576 }
2577 } else {
2578 return Err("bad range");
2579 }
2580 }
2581 "rationalize" => Matrix(
2582 a.iter()
2583 .map(|c| c_to_rational(c.number.clone(), options))
2584 .collect::<Vec<Vec<Number>>>(),
2585 ),
2586 "prime_factors" => {
2587 if a.len() != 2 {
2588 return Err("expected vector length 2");
2589 }
2590 let mut p1 = prime_factors(
2591 a[0].number.real().to_integer().unwrap_or_default(),
2592 );
2593 let p2 = prime_factors(
2594 a[1].number.real().to_integer().unwrap_or_default(),
2595 );
2596 for p in p1.iter_mut() {
2597 for m in &p2 {
2598 if p.0 == m.0 {
2599 p.1 -= m.1
2600 }
2601 }
2602 }
2603 for m in p2.iter() {
2604 if !p1.iter().any(|p| p.0 == m.0) {
2605 p1.push((m.0.clone(), -m.1));
2606 }
2607 }
2608 Matrix(
2609 p1.iter()
2610 .map(|p| {
2611 vec![
2612 Number::from(
2613 Complex::with_val(options.prec, p.0.clone()),
2614 None,
2615 ),
2616 Number::from(
2617 Complex::with_val(options.prec, p.1),
2618 None,
2619 ),
2620 ]
2621 })
2622 .collect::<Vec<Vec<Number>>>(),
2623 )
2624 }
2625 "poly" | "polynomial" => {
2626 if i + 1 < function.len() {
2627 let x = function.remove(i + 1).num()?;
2628 let mut sum = Number::from(Complex::new(options.prec), None);
2629 for (i, a) in a.iter().rev().enumerate() {
2630 let n = Number::from(
2631 a.number.clone() * x.number.clone().pow(i),
2632 mul_units(
2633 a.units,
2634 Some(x.units.unwrap_or_default().pow(i as f64)),
2635 ),
2636 );
2637 if i == 0 { sum = n } else { sum = add(&sum, &n) }
2638 }
2639 NumStr::new(sum)
2640 } else {
2641 return Err("not enough args");
2642 }
2643 }
2644 _ => do_functions(arg, options, &mut function, i, &to_deg, s)?,
2645 },
2646 _ => match s.as_str() {
2647 "rationalize" => Vector(c_to_rational(arg.num()?.number, options)),
2648 "domain_coloring_rgb" => {
2649 let pi = Float::with_val(options.prec, Pi);
2650 let num = arg.num()?.number;
2651 let hue: Float = 1 + (-num.clone()).arg().real().clone() / π
2652 let sat: Float = (1 + num.clone().abs().real().clone().fract()) / 2;
2653 let val: Float = {
2654 let (r, i) = (num * &pi).into_real_imag();
2655 let t1: Float = r.sin();
2656 let t2: Float = i.sin();
2657 (t1 * t2).abs().pow(0.125)
2658 };
2659 Vector(hsv2rgb(3 * hue, sat, val))
2660 }
2661 "multinomial" => {
2662 let mut numerator: Complex = arg.num()?.number + 1;
2663 let mut divisor = gamma(numerator.clone());
2664 while i + 1 < function.len() && !matches!(&function[i + 1], Func(_))
2665 {
2666 let temp = function.remove(i + 1).num()?.number;
2667 numerator += temp.clone();
2668 let temp = temp.clone() + 1;
2669 divisor *= gamma(temp);
2670 }
2671 NumStr::new(Number::from(gamma(numerator) / divisor, None))
2672 }
2673 "Β" | "B" | "beta" => {
2674 if i + 1 < function.len() {
2675 let a = arg.num()?.number;
2676 let b = function.remove(i + 1).num()?.number;
2677 if i + 1 < function.len() {
2678 let x = function.remove(i + 1).num()?.number;
2679 NumStr::new(Number::from(incomplete_beta(a, b, x), None))
2680 } else if a.imag().is_zero() && b.imag().is_zero() {
2681 NumStr::new(Number::from(
2682 gamma(a.clone()) * gamma(b.clone())
2683 / gamma(a + b.clone()),
2684 None,
2685 ))
2686 } else {
2687 NumStr::new(Number::from(
2688 incomplete_beta(
2689 Complex::with_val(options.prec, 1),
2690 a,
2691 b,
2692 ),
2693 None,
2694 ))
2695 }
2696 } else {
2697 return Err("not enough args");
2698 }
2699 }
2700 "gamma_pdf" => {
2701 if i + 2 < function.len() {
2702 let x = arg.num()?.number;
2703 let a = function.remove(i + 1).num()?.number;
2704 let b = function.remove(i + 1).num()?.number;
2705 NumStr::new(Number::from(
2706 pow_nth(b.clone(), -a.clone()) / gamma(a.clone())
2707 * pow_nth(x.clone(), a - 1)
2708 * (-x / b).exp(),
2709 None,
2710 ))
2711 } else {
2712 return Err("not enough args");
2713 }
2714 }
2715 "gamma_cdf" => {
2716 if i + 2 < function.len() {
2717 let x = arg.num()?.number;
2718 let a = function.remove(i + 1).num()?.number;
2719 let b = function.remove(i + 1).num()?.number;
2720 NumStr::new(Number::from(
2721 1 - incomplete_gamma(a.clone(), x / b) / gamma(a),
2722 None,
2723 ))
2724 } else {
2725 return Err("not enough args");
2726 }
2727 }
2728 "I" | "betaC" | "beta_cdf" => {
2729 if i + 2 < function.len() {
2730 let x = arg.num()?.number;
2731 let a = function.remove(i + 1).num()?.number;
2732 let b = function.remove(i + 1).num()?.number;
2733 if i + 1 < function.len() {
2734 let z = function.remove(i + 1).num()?.number;
2735 NumStr::new(Number::from(
2736 regularized_incomplete_beta(a, b.clone(), z.clone())
2737 - regularized_incomplete_beta(x, b, z),
2738 None,
2739 ))
2740 } else {
2741 NumStr::new(Number::from(
2742 regularized_incomplete_beta(x, a, b),
2743 None,
2744 ))
2745 }
2746 } else {
2747 return Err("not enough args");
2748 }
2749 }
2750 "betaP" | "beta_pdf" => {
2751 if i + 2 < function.len() {
2752 let x = arg.num()?.number;
2753 let alpha = function.remove(i + 1).num()?.number;
2754 let beta = function.remove(i + 1).num()?.number;
2755 let c: Complex = 1 - x.clone();
2756 NumStr::new(Number::from(
2757 gamma(alpha.clone() + beta.clone())
2758 * pow_nth(x, alpha.clone() - 1)
2759 * pow_nth(c, beta.clone() - 1)
2760 / (gamma(alpha) * gamma(beta)),
2761 None,
2762 ))
2763 } else {
2764 return Err("not enough args");
2765 }
2766 }
2767 "normP" | "norm_pdf" => {
2768 if i + 2 < function.len() {
2769 let x = arg.num()?.number;
2770 let mu = function.remove(i + 1).num()?.number;
2771 let sigma = function.remove(i + 1).num()?.number;
2772 let n: Complex = sqr(x - mu);
2773 let n: Complex = -n / (2 * sqr(sigma.clone()));
2774 let tau: Complex = 2 * Complex::with_val(options.prec, Pi);
2775 NumStr::new(Number::from(n.exp() / (sigma * tau.sqrt()), None))
2776 } else {
2777 return Err("not enough args");
2778 }
2779 }
2780 "normD" | "norm_cdf" => {
2781 let mut a = arg.num()?.number;
2782 if i + 2 < function.len() {
2783 a -= function.remove(i + 1).num()?.number;
2784 a /= function.remove(i + 1).num()?.number;
2785 }
2786 if a.imag().is_zero() {
2787 let two = Float::with_val(options.prec, 2);
2788 NumStr::new(Number::from(
2789 ((-a / two.clone().sqrt()).real().clone().erfc() / two)
2790 .into(),
2791 None,
2792 ))
2793 } else {
2794 let two = Float::with_val(options.prec, 2);
2795 NumStr::new(Number::from(
2796 erf(-a / two.clone().sqrt()) / two,
2797 None,
2798 ))
2799 }
2800 }
2801 "lognorm_cdf" => {
2802 let mut a = arg.num()?.number.ln();
2803 if i + 2 < function.len() {
2804 a -= function.remove(i + 1).num()?.number;
2805 a /= function.remove(i + 1).num()?.number;
2806 }
2807 if a.imag().is_zero() {
2808 let two = Float::with_val(options.prec, 2);
2809 NumStr::new(Number::from(
2810 ((-a / two.clone().sqrt()).real().clone().erfc() / two)
2811 .into(),
2812 None,
2813 ))
2814 } else {
2815 let two = Float::with_val(options.prec, 2);
2816 NumStr::new(Number::from(
2817 erf(-a / two.clone().sqrt()) / two,
2818 None,
2819 ))
2820 }
2821 }
2822 "lognorm_pdf" => {
2823 if i + 2 < function.len() {
2824 let x = arg.num()?.number;
2825 let mu = function.remove(i + 1).num()?.number;
2826 let sigma = function.remove(i + 1).num()?.number;
2827 let n: Complex = sqr(x.clone().ln() - mu);
2828 let n: Complex = -n / (2 * sqr(sigma.clone()));
2829 let tau: Complex = 2 * Complex::with_val(options.prec, Pi);
2830 NumStr::new(Number::from(
2831 n.exp() / (sigma * tau.sqrt() * x),
2832 None,
2833 ))
2834 } else {
2835 return Err("not enough args");
2836 }
2837 }
2838 "poisson_pmf" => {
2839 if i + 1 < function.len() {
2840 let k = arg.num()?.number;
2841 let l = function.remove(i + 1).num()?.number;
2842 NumStr::new(Number::from(
2843 pow_nth(l.clone(), k.clone()) * (-l).exp() / gamma(k + 1),
2844 None,
2845 ))
2846 } else {
2847 return Err("not enough args");
2848 }
2849 }
2850 "poisson_cdf" => {
2851 if i + 1 < function.len() {
2852 let k = arg.num()?.number;
2853 let l = function.remove(i + 1).num()?.number;
2854 NumStr::new(Number::from(
2855 incomplete_gamma(k.clone() + 1, l) / gamma(k + 1),
2856 None,
2857 ))
2858 } else {
2859 return Err("not enough args");
2860 }
2861 }
2862 "binomial_cdf" => {
2863 if i + 2 < function.len() {
2864 let k = arg.num()?.number;
2865 let n = function.remove(i + 1).num()?.number;
2866 let p = function.remove(i + 1).num()?.number;
2867 let q: Complex = 1 - p.clone();
2868 NumStr::new(Number::from(
2869 regularized_incomplete_beta(q, n - k.clone(), 1 + k),
2870 None,
2871 ))
2872 } else {
2873 return Err("not enough args");
2874 }
2875 }
2876 "binomial_pmf" => {
2877 if i + 2 < function.len() {
2878 let k = arg.num()?.number;
2879 let n = function.remove(i + 1).num()?.number;
2880 let p = function.remove(i + 1).num()?.number;
2881 let q: Complex = 1 - p.clone();
2882 NumStr::new(Number::from(
2883 binomial(n.clone(), k.clone())
2884 * pow_nth(p, k.clone())
2885 * pow_nth(q, n - k),
2886 None,
2887 ))
2888 } else {
2889 return Err("not enough args");
2890 }
2891 }
2892 "neg_binomial_cdf" => {
2893 if i + 2 < function.len() {
2894 let k = arg.num()?.number;
2895 let r = function.remove(i + 1).num()?.number;
2896 let p = function.remove(i + 1).num()?.number;
2897 NumStr::new(Number::from(
2898 regularized_incomplete_beta(p, r, k + 1),
2899 None,
2900 ))
2901 } else {
2902 return Err("not enough args");
2903 }
2904 }
2905 "neg_binomial_pmf" => {
2906 if i + 2 < function.len() {
2907 let k = arg.num()?.number;
2908 let r = function.remove(i + 1).num()?.number;
2909 let p = function.remove(i + 1).num()?.number;
2910 let q: Complex = 1 - p.clone();
2911 NumStr::new(Number::from(
2912 binomial(k.clone() + r.clone() - 1, k.clone())
2913 * pow_nth(p, r)
2914 * pow_nth(q, k),
2915 None,
2916 ))
2917 } else {
2918 return Err("not enough args");
2919 }
2920 }
2921 "hypergeometric_cdf" => {
2922 if i + 3 < function.len() {
2923 let mut k = arg.num()?.number.real().clone().floor();
2924 let pop = function.remove(i + 1).num()?.number;
2925 let success = function.remove(i + 1).num()?.number;
2926 let draws = function.remove(i + 1).num()?.number;
2927 let mut sum = Complex::new(options.prec);
2928 while k >= 0 {
2929 sum += binomial(success.clone(), k.clone().into())
2930 * binomial(
2931 pop.clone() - success.clone(),
2932 draws.clone() - k.clone(),
2933 )
2934 / binomial(pop.clone(), draws.clone());
2935 k -= 1
2936 }
2937 NumStr::new(Number::from(
2938 Complex::with_val(options.prec, sum),
2939 None,
2940 ))
2941 } else {
2942 return Err("not enough args");
2943 }
2944 }
2945 "hypergeometric_pmf" => {
2946 if i + 3 < function.len() {
2947 let k = arg.num()?.number;
2948 let pop = function.remove(i + 1).num()?.number;
2949 let success = function.remove(i + 1).num()?.number;
2950 let draws = function.remove(i + 1).num()?.number;
2951 NumStr::new(Number::from(
2952 binomial(success.clone(), k.clone())
2953 * binomial(
2954 pop.clone() - success.clone(),
2955 draws.clone() - k,
2956 )
2957 / binomial(pop, draws),
2958 None,
2959 ))
2960 } else {
2961 return Err("not enough args");
2962 }
2963 }
2964 #[cfg(feature = "fastrand")]
2965 "rand_hypergeometric" => {
2966 if i + 2 < function.len() {
2967 let mut pop = arg.num()?.number.real().clone();
2968 let mut success =
2969 function.remove(i + 1).num()?.number.real().clone();
2970 let mut draws =
2971 function.remove(i + 1).num()?.number.real().clone();
2972 let mut sum = Integer::new();
2973 while draws > 0 {
2974 if success.clone() / pop.clone()
2975 > Float::with_val(options.prec, fastrand::u128(..))
2976 / u128::MAX
2977 {
2978 sum += 1;
2979 success -= 1;
2980 }
2981 pop -= 1;
2982 draws -= 1
2983 }
2984 NumStr::new(Number::from(
2985 Complex::with_val(options.prec, sum),
2986 None,
2987 ))
2988 } else {
2989 return Err("not enough args");
2990 }
2991 }
2992 "neg_hypergeometric_cdf" => {
2993 if i + 3 < function.len() {
2994 let mut k = arg.num()?.number.real().clone().floor();
2995 let pop = function.remove(i + 1).num()?.number;
2996 let success = function.remove(i + 1).num()?.number;
2997 let fails = function.remove(i + 1).num()?.number;
2998 let mut sum = Complex::new(options.prec);
2999 while k >= 0 {
3000 sum += binomial(
3001 k.clone() + fails.clone() - 1,
3002 k.clone().into(),
3003 ) * binomial(
3004 pop.clone() - k.clone() - fails.clone(),
3005 success.clone() - k.clone(),
3006 ) / binomial(pop.clone(), success.clone());
3007 k -= 1
3008 }
3009 NumStr::new(Number::from(
3010 Complex::with_val(options.prec, sum),
3011 None,
3012 ))
3013 } else {
3014 return Err("not enough args");
3015 }
3016 }
3017 "neg_hypergeometric_pmf" => {
3018 if i + 3 < function.len() {
3019 let k = arg.num()?.number;
3020 let pop = function.remove(i + 1).num()?.number;
3021 let success = function.remove(i + 1).num()?.number;
3022 let fails = function.remove(i + 1).num()?.number;
3023 NumStr::new(Number::from(
3024 binomial(k.clone() + fails.clone() - 1, k.clone())
3025 * binomial(
3026 pop.clone() - k.clone() - fails,
3027 success.clone() - k,
3028 )
3029 / binomial(pop, success),
3030 None,
3031 ))
3032 } else {
3033 return Err("not enough args");
3034 }
3035 }
3036 #[cfg(feature = "fastrand")]
3037 "rand_neg_hypergeometric" => {
3038 if i + 2 < function.len() {
3039 let mut pop = arg.num()?.number.real().clone();
3040 let mut success =
3041 function.remove(i + 1).num()?.number.real().clone();
3042 let mut fails =
3043 function.remove(i + 1).num()?.number.real().clone();
3044 let mut sum = Integer::new();
3045 while fails > 0 {
3046 if success.clone() / pop.clone()
3047 > Float::with_val(options.prec, fastrand::u128(..))
3048 / u128::MAX
3049 {
3050 sum += 1;
3051 success -= 1;
3052 } else {
3053 fails -= 1
3054 }
3055 pop -= 1;
3056 }
3057 NumStr::new(Number::from(
3058 Complex::with_val(options.prec, sum),
3059 None,
3060 ))
3061 } else {
3062 return Err("not enough args");
3063 }
3064 }
3065 "geometric_cdf" => {
3066 if i + 1 < function.len() {
3067 let k = arg.num()?.number;
3068 let p = function.remove(i + 1).num()?.number;
3069 let q: Complex = 1 - p.clone();
3070 NumStr::new(Number::from(1 - pow_nth(q, k), None))
3071 } else {
3072 return Err("not enough args");
3073 }
3074 }
3075 "geometric_pmf" => {
3076 if i + 1 < function.len() {
3077 let k = arg.num()?.number;
3078 let p = function.remove(i + 1).num()?.number;
3079 let q: Complex = 1 - p.clone();
3080 NumStr::new(Number::from(pow_nth(q, k - 1) * p, None))
3081 } else {
3082 return Err("not enough args");
3083 }
3084 }
3085 "quartic" => {
3086 if i + 5 < function.len() {
3087 let a = arg.num()?;
3088 let b = function.remove(i + 1).num()?;
3089 let c = function.remove(i + 1).num()?;
3090 let d = function.remove(i + 1).num()?;
3091 let e = function.remove(i + 1).num()?;
3092 function.remove(i + 1);
3093 let n = quartic(a, b, c, d, e, true);
3094 if n.len() == 1 {
3095 NumStr::new(n[0].clone())
3096 } else {
3097 Vector(n)
3098 }
3099 } else if i + 4 < function.len() {
3100 let a = arg.num()?;
3101 let b = function.remove(i + 1).num()?;
3102 let c = function.remove(i + 1).num()?;
3103 let d = function.remove(i + 1).num()?;
3104 let e = function.remove(i + 1).num()?;
3105 let n = quartic(a, b, c, d, e, false);
3106 if n.len() == 1 {
3107 NumStr::new(n[0].clone())
3108 } else {
3109 Vector(n)
3110 }
3111 } else if i + 3 < function.len() {
3112 let b = arg.num()?;
3113 let c = function.remove(i + 1).num()?;
3114 let d = function.remove(i + 1).num()?;
3115 let e = function.remove(i + 1).num()?;
3116 let n = quartic(
3117 Number::from(Complex::with_val(options.prec, 1), None),
3118 b,
3119 c,
3120 d,
3121 e,
3122 false,
3123 );
3124 if n.len() == 1 {
3125 NumStr::new(n[0].clone())
3126 } else {
3127 Vector(n)
3128 }
3129 } else {
3130 return Err("not enough args");
3131 }
3132 }
3133 "cubic" => {
3134 if i + 4 < function.len() {
3135 let a = arg.num()?;
3136 let b = function.remove(i + 1).num()?;
3137 let c = function.remove(i + 1).num()?;
3138 let d = function.remove(i + 1).num()?;
3139 function.remove(i + 1);
3140 let n = cubic(a, b, c, d, true);
3141 if n.len() == 1 {
3142 NumStr::new(n[0].clone())
3143 } else {
3144 Vector(n)
3145 }
3146 } else if i + 3 < function.len() {
3147 let a = arg.num()?;
3148 let b = function.remove(i + 1).num()?;
3149 let c = function.remove(i + 1).num()?;
3150 let d = function.remove(i + 1).num()?;
3151 let n = cubic(a, b, c, d, false);
3152 if n.len() == 1 {
3153 NumStr::new(n[0].clone())
3154 } else {
3155 Vector(n)
3156 }
3157 } else if i + 2 < function.len() {
3158 let b = arg.num()?;
3159 let c = function.remove(i + 1).num()?;
3160 let d = function.remove(i + 1).num()?;
3161 let n = cubic(
3162 Number::from(Complex::with_val(options.prec, 1), None),
3163 b,
3164 c,
3165 d,
3166 false,
3167 );
3168 if n.len() == 1 {
3169 NumStr::new(n[0].clone())
3170 } else {
3171 Vector(n)
3172 }
3173 } else {
3174 return Err("not enough args");
3175 }
3176 }
3177 "quad" | "quadratic" => {
3178 if i + 3 < function.len() {
3179 let a = arg.num()?;
3180 let b = function.remove(i + 1).num()?;
3181 let c = function.remove(i + 1).num()?;
3182 function.remove(i + 1);
3183 let n = quadratic(a, b, c, true);
3184 if n.is_empty() {
3185 NumStr::new(Number::from(
3186 Complex::with_val(options.prec, Nan),
3187 None,
3188 ))
3189 } else if n.len() == 1 {
3190 NumStr::new(n[0].clone())
3191 } else {
3192 Vector(n)
3193 }
3194 } else if i + 2 < function.len() {
3195 let a = arg.num()?;
3196 let b = function.remove(i + 1).num()?;
3197 let c = function.remove(i + 1).num()?;
3198 let n = quadratic(a, b, c, false);
3199 if n.len() == 1 {
3200 NumStr::new(n[0].clone())
3201 } else {
3202 Vector(n)
3203 }
3204 } else if i + 1 < function.len() {
3205 let b = arg.num()?;
3206 let c = function.remove(i + 1).num()?;
3207 let n = quadratic(
3208 Number::from(Complex::with_val(options.prec, 1), None),
3209 b,
3210 c,
3211 false,
3212 );
3213 if n.len() == 1 {
3214 NumStr::new(n[0].clone())
3215 } else {
3216 Vector(n)
3217 }
3218 } else {
3219 return Err("not enough args");
3220 }
3221 }
3222 "cossin" => {
3223 let (a, b) = arg.num()?.number.sin_cos(Complex::new(options.prec));
3224 Vector(vec![Number::from(b, None), Number::from(a, None)])
3225 }
3226 "sincos" => {
3227 let (a, b) = arg.num()?.number.sin_cos(Complex::new(options.prec));
3228 Vector(vec![Number::from(a, None), Number::from(b, None)])
3229 }
3230 "split" => {
3231 let a = arg.num()?.number;
3232 Vector(vec![
3233 Number::from(a.real().clone().into(), None),
3234 Number::from(a.imag().clone().into(), None),
3235 ])
3236 }
3237 "ceil" => {
3238 let a = arg.num()?;
3239 let units = a.units;
3240 let m = if i + 1 < function.len() {
3241 Complex::with_val(a.number.prec(), options.base.1)
3242 .pow(function.remove(i + 1).num()?.number)
3243 } else {
3244 Complex::with_val(a.number.prec(), 1)
3245 };
3246 let a = a.number * m.clone();
3247 NumStr::new(Number::from(
3248 Complex::with_val(
3249 options.prec,
3250 (a.real().clone().ceil(), a.imag().clone().ceil()),
3251 ) / m,
3252 units,
3253 ))
3254 }
3255 "floor" => {
3256 let a = arg.num()?;
3257 let units = a.units;
3258 let m = if i + 1 < function.len() {
3259 Complex::with_val(a.number.prec(), options.base.1)
3260 .pow(function.remove(i + 1).num()?.number)
3261 } else {
3262 Complex::with_val(a.number.prec(), 1)
3263 };
3264 let a = a.number * m.clone();
3265 NumStr::new(Number::from(
3266 Complex::with_val(
3267 options.prec,
3268 (a.real().clone().floor(), a.imag().clone().floor()),
3269 ) / m,
3270 units,
3271 ))
3272 }
3273 "round" => {
3274 let a = arg.num()?;
3275 let units = a.units;
3276 let m = if i + 1 < function.len() {
3277 Complex::with_val(a.number.prec(), options.base.1)
3278 .pow(function.remove(i + 1).num()?.number)
3279 } else {
3280 Complex::with_val(a.number.prec(), 1)
3281 };
3282 let a = a.number * m.clone();
3283 NumStr::new(Number::from(
3284 Complex::with_val(
3285 options.prec,
3286 (a.real().clone().round(), a.imag().clone().round()),
3287 ) / m,
3288 units,
3289 ))
3290 }
3291 "int" | "trunc" => {
3292 let a = arg.num()?;
3293 let units = a.units;
3294 let m = if i + 1 < function.len() {
3295 Complex::with_val(a.number.prec(), options.base.1)
3296 .pow(function.remove(i + 1).num()?.number)
3297 } else {
3298 Complex::with_val(a.number.prec(), 1)
3299 };
3300 let a = a.number * m.clone();
3301 NumStr::new(Number::from(
3302 Complex::with_val(
3303 options.prec,
3304 (a.real().clone().trunc(), a.imag().clone().trunc()),
3305 ) / m,
3306 units,
3307 ))
3308 }
3309 "frac" | "fract" => {
3310 let a = arg.num()?;
3311 let units = a.units;
3312 let m = if i + 1 < function.len() {
3313 Complex::with_val(a.number.prec(), options.base.1)
3314 .pow(function.remove(i + 1).num()?.number)
3315 } else {
3316 Complex::with_val(a.number.prec(), 1)
3317 };
3318 let a = a.number * m.clone();
3319 NumStr::new(Number::from(
3320 Complex::with_val(
3321 options.prec,
3322 (a.real().clone().fract(), a.imag().clone().fract()),
3323 ) / m,
3324 units,
3325 ))
3326 }
3327 "iden" | "identity" => Matrix(identity(
3328 arg.num()?
3329 .number
3330 .real()
3331 .to_integer()
3332 .unwrap_or_default()
3333 .to_usize()
3334 .unwrap_or_default(),
3335 options.prec,
3336 )),
3337 "rotate" => {
3338 if i + 2 < function.len() {
3339 let (sina, cosa) = (arg.num()?.number / to_deg.clone())
3340 .sin_cos(Complex::new(options.prec));
3341 let (sinb, cosb) = (function.remove(i + 1).num()?.number
3342 / to_deg.clone())
3343 .sin_cos(Complex::new(options.prec));
3344 let (sinc, cosc) = (function.remove(i + 1).num()?.number
3345 / to_deg.clone())
3346 .sin_cos(Complex::new(options.prec));
3347 Matrix(vec![
3348 vec![
3349 Number::from(cosa.clone() * cosb.clone(), None),
3350 Number::from(
3351 cosa.clone() * sinb.clone() * sinc.clone()
3352 - sina.clone() * cosc.clone(),
3353 None,
3354 ),
3355 Number::from(
3356 cosa.clone() * sinb.clone() * cosc.clone()
3357 + sina.clone() * sinc.clone(),
3358 None,
3359 ),
3360 ],
3361 vec![
3362 Number::from(sina.clone() * cosb.clone(), None),
3363 Number::from(
3364 sina.clone() * sinb.clone() * sinc.clone()
3365 + cosa.clone() * cosc.clone(),
3366 None,
3367 ),
3368 Number::from(
3369 sina.clone() * sinb.clone() * cosc.clone()
3370 - cosa.clone() * sinc.clone(),
3371 None,
3372 ),
3373 ],
3374 vec![
3375 Number::from(-sinb.clone(), None),
3376 Number::from(cosb.clone() * sinc.clone(), None),
3377 Number::from(cosb.clone() * cosc.clone(), None),
3378 ],
3379 ])
3380 } else {
3381 let (sin, cos) = (arg.num()?.number / to_deg.clone())
3382 .sin_cos(Complex::new(options.prec));
3383 Matrix(vec![
3384 vec![
3385 Number::from(cos.clone(), None),
3386 Number::from(-sin.clone(), None),
3387 ],
3388 vec![Number::from(sin, None), Number::from(cos, None)],
3389 ])
3390 }
3391 }
3392 "prime_factors" => {
3393 let a = arg.num()?.number;
3394 if a.imag().is_zero() {
3395 if a.real().clone().fract().is_zero() {
3396 let n = a.real().to_integer().unwrap_or_default();
3397 let m = prime_factors(n);
3398 if m.is_empty() {
3399 NumStr::new(Number::from(
3400 Complex::with_val(options.prec, Nan),
3401 None,
3402 ))
3403 } else {
3404 Matrix(
3405 m.iter()
3406 .map(|(a, b)| {
3407 vec![
3408 Number::from(
3409 Complex::with_val(options.prec, a),
3410 None,
3411 ),
3412 Number::from(
3413 Complex::with_val(options.prec, b),
3414 None,
3415 ),
3416 ]
3417 })
3418 .collect::<Vec<Vec<Number>>>(),
3419 )
3420 }
3421 } else if let Some(a) = rationalize(a.real().clone(), options) {
3422 let mut p1 = prime_factors(a.0);
3423 let p2 = prime_factors(a.1);
3424 for p in p1.iter_mut() {
3425 for m in &p2 {
3426 if p.0 == m.0 {
3427 p.1 -= m.1
3428 }
3429 }
3430 }
3431 for m in p2.iter() {
3432 if !p1.iter().any(|p| p.0 == m.0) {
3433 p1.push((m.0.clone(), -m.1));
3434 }
3435 }
3436 Matrix(
3437 p1.iter()
3438 .map(|p| {
3439 vec![
3440 Number::from(
3441 Complex::with_val(
3442 options.prec,
3443 p.0.clone(),
3444 ),
3445 None,
3446 ),
3447 Number::from(
3448 Complex::with_val(options.prec, p.1),
3449 None,
3450 ),
3451 ]
3452 })
3453 .collect::<Vec<Vec<Number>>>(),
3454 )
3455 } else {
3456 NumStr::new(Number::from(
3457 Complex::with_val(options.prec, Nan),
3458 None,
3459 ))
3460 }
3461 } else {
3462 NumStr::new(Number::from(
3463 Complex::with_val(options.prec, Nan),
3464 None,
3465 ))
3466 }
3467 }
3468 "factors" | "factor" => {
3469 let a = arg.num()?.number;
3470 if a.imag().clone().is_zero() {
3471 if a.real().clone().fract().is_zero() {
3472 let mut vec = Vec::new();
3473 let n = a
3474 .real()
3475 .to_integer()
3476 .unwrap_or_default()
3477 .to_usize()
3478 .unwrap_or_default();
3479 for i in 1..=n {
3480 if n % i == 0 {
3481 vec.push(Number::from(
3482 Complex::with_val(options.prec, i),
3483 None,
3484 ));
3485 }
3486 }
3487 Vector(vec)
3488 } else {
3489 NumStr::new(Number::from(
3490 Complex::with_val(options.prec, Nan),
3491 None,
3492 ))
3493 }
3494 } else {
3495 NumStr::new(Number::from(
3496 Complex::with_val(options.prec, Nan),
3497 None,
3498 ))
3499 }
3500 }
3501 "unity" => {
3502 let vec = if i + 1 < function.len()
3503 && !matches!(&function[i + 1], Func(_))
3504 {
3505 unity(arg.num()?.number, function.remove(i + 1).num()?.number)
3506 } else {
3507 unity(Complex::with_val(options.prec, 1), arg.num()?.number)
3508 };
3509 if vec.is_empty() {
3510 Vector(vec![Number::from(
3511 Complex::with_val(options.prec, Nan),
3512 None,
3513 )])
3514 } else {
3515 Vector(vec)
3516 }
3517 }
3518 _ => do_functions(arg, options, &mut function, i, &to_deg, s)?,
3519 },
3520 }
3521 }
3522 }
3523 }
3524 i += 1;
3525 }
3526 i = 0;
3527 while i < function.len() {
3528 if let Func(s) = &function[i] {
3529 function[i] = match s.as_str() {
3530 #[cfg(feature = "fastrand")]
3531 "rnd" | "rand" => NumStr::new(Number::from(
3532 Complex::with_val(options.prec, fastrand::u128(..)) / u128::MAX,
3533 None,
3534 )),
3535 "epoch" => NumStr::new(Number::from(
3536 Complex::with_val(
3537 options.prec,
3538 match SystemTime::now().duration_since(SystemTime::UNIX_EPOCH) {
3539 Ok(n) => n.as_nanos(),
3540 _ => return Err("epoch fail"),
3541 },
3542 ) / 1000000000,
3543 Some(Units {
3544 second: 1.0,
3545 ..Units::default()
3546 }),
3547 )),
3548 _ => {
3549 i += 1;
3550 continue;
3551 }
3552 }
3553 } else {
3554 i += 1;
3555 }
3556 }
3557 i = 1;
3558 while i < function.len().saturating_sub(1) {
3559 function[i] = match &function[i] {
3560 Modulo => function[i - 1].func(&function[i + 1], rem)?,
3561 Range => to(&function[i - 1], &function[i + 1])?,
3562 _ => {
3563 i += 1;
3564 continue;
3565 }
3566 };
3567 function.remove(i + 1);
3568 function.remove(i - 1);
3569 }
3570 i = function.len().saturating_sub(2);
3571 while i != 0 {
3572 function[i] = match &function[i] {
3573 Exponent => function[i - 1].pow(&function[i + 1])?,
3574 Tetration => function[i - 1].func(&function[i + 1], tetration)?,
3575 Root => function[i - 1].func(&function[i + 1], root)?,
3576 _ => {
3577 i -= 1;
3578 continue;
3579 }
3580 };
3581 function.remove(i + 1);
3582 function.remove(i - 1);
3583 i = i.saturating_sub(2);
3584 }
3585 i = 1;
3586 while i < function.len().saturating_sub(1) {
3587 function[i] = match &function[i] {
3588 InternalMultiplication => function[i - 1].mul(&function[i + 1])?,
3589 _ => {
3590 i += 1;
3591 continue;
3592 }
3593 };
3594 function.remove(i + 1);
3595 function.remove(i - 1);
3596 }
3597 i = 1;
3598 while i < function.len().saturating_sub(1) {
3599 function[i] = match &function[i] {
3600 Multiplication => function[i - 1].mul(&function[i + 1])?,
3601 Division => function[i - 1].func(&function[i + 1], div)?,
3602 _ => {
3603 i += 1;
3604 continue;
3605 }
3606 };
3607 function.remove(i + 1);
3608 function.remove(i - 1);
3609 }
3610 i = 1;
3611 while i < function.len().saturating_sub(1) {
3612 function[i] = match &function[i] {
3613 PlusMinus => function[i - 1].pm(&function[i + 1])?,
3614 Plus => function[i - 1].func(&function[i + 1], add)?,
3615 Minus => function[i - 1].func(&function[i + 1], sub)?,
3616 _ => {
3617 i += 1;
3618 continue;
3619 }
3620 };
3621 function.remove(i + 1);
3622 function.remove(i - 1);
3623 }
3624 if options.units {
3625 i = 1;
3626 while i < function.len().saturating_sub(1) {
3627 function[i] = match &function[i] {
3628 Conversion => function[i - 1].func(&function[i + 1], div)?,
3629 _ => {
3630 i += 1;
3631 continue;
3632 }
3633 };
3634 function.remove(i + 1);
3635 function.remove(i - 1);
3636 }
3637 }
3638 i = 1;
3639 while i < function.len().saturating_sub(1) {
3640 function[i] = match &function[i] {
3641 Lesser => {
3642 if i + 3 < function.len()
3643 && matches!(
3644 &function[i + 2],
3645 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3646 )
3647 {
3648 function[i - 1] = function[i + 1].func(&function[i - 1], gt)?;
3649 function[i] = And;
3650 continue;
3651 }
3652 function[i + 1].func(&function[i - 1], gt)?
3653 }
3654 LesserEqual => {
3655 if i + 3 < function.len()
3656 && matches!(
3657 &function[i + 2],
3658 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3659 )
3660 {
3661 function[i - 1] = function[i + 1].func(&function[i - 1], ge)?;
3662 function[i] = And;
3663 continue;
3664 }
3665 function[i + 1].func(&function[i - 1], ge)?
3666 }
3667 Greater => {
3668 if i + 3 < function.len()
3669 && matches!(
3670 &function[i + 2],
3671 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3672 )
3673 {
3674 function[i - 1] = function[i - 1].func(&function[i + 1], gt)?;
3675 function[i] = And;
3676 continue;
3677 }
3678 function[i - 1].func(&function[i + 1], gt)?
3679 }
3680 GreaterEqual => {
3681 if i + 3 < function.len()
3682 && matches!(
3683 &function[i + 2],
3684 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3685 )
3686 {
3687 function[i - 1] = function[i - 1].func(&function[i + 1], ge)?;
3688 function[i] = And;
3689 continue;
3690 }
3691 function[i - 1].func(&function[i + 1], ge)?
3692 }
3693 Equal => {
3694 if i + 3 < function.len()
3695 && matches!(
3696 &function[i + 2],
3697 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3698 )
3699 {
3700 function[i - 1] = function[i - 1].func(&function[i + 1], eq)?;
3701 function[i] = And;
3702 continue;
3703 }
3704 function[i - 1].func(&function[i + 1], eq)?
3705 }
3706 NotEqual => {
3707 if i + 3 < function.len()
3708 && matches!(
3709 &function[i + 2],
3710 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3711 )
3712 {
3713 function[i - 1] = function[i - 1].func(&function[i + 1], ne)?;
3714 function[i] = And;
3715 continue;
3716 }
3717 function[i - 1].func(&function[i + 1], ne)?
3718 }
3719 NearEqual => {
3720 if i + 3 < function.len()
3721 && matches!(
3722 &function[i + 2],
3723 Lesser | Greater | Equal | NearEqual | LesserEqual | GreaterEqual
3724 )
3725 {
3726 function[i - 1] = function[i - 1].func(&function[i + 1], about_eq)?;
3727 function[i] = And;
3728 continue;
3729 }
3730 function[i - 1].func(&function[i + 1], about_eq)?
3731 }
3732 ShiftRight => function[i - 1].func(&function[i + 1], shr)?,
3733 ShiftLeft => function[i - 1].func(&function[i + 1], shl)?,
3734 _ => {
3735 i += 1;
3736 continue;
3737 }
3738 };
3739 function.remove(i + 1);
3740 function.remove(i - 1);
3741 }
3742 if !function.is_empty() && function[0] == Not {
3743 function[0] = not(&function.remove(1))?;
3744 }
3745 i = 1;
3746 while i < function.len().saturating_sub(1) {
3747 function[i] = match &function[i] {
3748 And => function[i - 1].func(&function[i + 1], and)?,
3749 Or => function[i - 1].func(&function[i + 1], or)?,
3750 Xor => function[i - 1].func(&function[i + 1], xor)?,
3751 Implies => function[i - 1].func(&function[i + 1], implies)?,
3752 Nand => function[i - 1].func(&function[i + 1], nand)?,
3753 Nor => function[i - 1].func(&function[i + 1], nor)?,
3754 Converse => function[i + 1].func(&function[i - 1], implies)?,
3755 Not => {
3756 function[i] = not(&function.remove(i + 1))?;
3757 continue;
3758 }
3759 _ => {
3760 i += 1;
3761 continue;
3762 }
3763 };
3764 function.remove(i + 1);
3765 function.remove(i - 1);
3766 }
3767 if function.len() == 1 {
3768 Ok(function[0].clone())
3769 } else {
3770 Err("failed to compute")
3771 }
3772}
3773fn do_functions(
3774 a: NumStr,
3775 options: Options,
3776 function: &mut Vec<NumStr>,
3777 k: usize,
3778 to_deg: &Complex,
3779 s: &str,
3780) -> Result<NumStr, &'static str> {
3781 if function.len() > k + 1 {
3782 match (a.clone(), function[k + 1].clone()) {
3783 (Num(a), Num(b)) => {
3784 function.remove(k + 1);
3785 return Ok(NumStr::new(functions(
3786 *a,
3787 Some(*b),
3788 to_deg.clone(),
3789 s,
3790 options,
3791 )?));
3792 }
3793 (Vector(a), Vector(b)) => {
3794 function.remove(k + 1);
3795 let mut mat = Vec::new();
3796 for a in a {
3797 let mut vec = Vec::new();
3798 for b in &b {
3799 vec.push(functions(
3800 a.clone(),
3801 Some(b.clone()),
3802 to_deg.clone(),
3803 s,
3804 options,
3805 )?)
3806 }
3807 mat.push(vec);
3808 }
3809 return Ok(Matrix(mat));
3810 }
3811 (Matrix(a), Matrix(b))
3812 if a.len() == b.len() && (0..a.len()).all(|i| a[i].len() == b[i].len()) =>
3813 {
3814 function.remove(k + 1);
3815 let mut mat = Vec::new();
3816 for i in 0..a.len() {
3817 let mut vec = Vec::new();
3818 for j in 0..a[i].len() {
3819 vec.push(functions(
3820 a[i][j].clone(),
3821 Some(b[i][j].clone()),
3822 to_deg.clone(),
3823 s,
3824 options,
3825 )?)
3826 }
3827 mat.push(vec.clone());
3828 }
3829 return Ok(Matrix(mat));
3830 }
3831 (Num(a), Vector(b)) => {
3832 function.remove(k + 1);
3833 let mut vec = Vec::new();
3834 for i in b {
3835 vec.push(functions(*a.clone(), Some(i), to_deg.clone(), s, options)?)
3836 }
3837 return Ok(Vector(vec));
3838 }
3839 (Vector(a), Num(b)) => {
3840 function.remove(k + 1);
3841 let mut vec = Vec::new();
3842 for i in a {
3843 vec.push(functions(i, Some(*b.clone()), to_deg.clone(), s, options)?)
3844 }
3845 return Ok(Vector(vec));
3846 }
3847 (Num(a), Matrix(b)) => {
3848 function.remove(k + 1);
3849 let mut mat = Vec::new();
3850 for i in b {
3851 let mut vec = Vec::new();
3852 for j in i {
3853 vec.push(functions(*a.clone(), Some(j), to_deg.clone(), s, options)?)
3854 }
3855 mat.push(vec.clone());
3856 }
3857 return Ok(Matrix(mat));
3858 }
3859 (Matrix(a), Num(b)) => {
3860 function.remove(k + 1);
3861 let mut mat = Vec::new();
3862 for i in a {
3863 let mut vec = Vec::new();
3864 for j in i {
3865 vec.push(functions(j, Some(*b.clone()), to_deg.clone(), s, options)?)
3866 }
3867 mat.push(vec.clone());
3868 }
3869 return Ok(Matrix(mat));
3870 }
3871 (Matrix(a), Vector(b)) if a.len() == b.len() => {
3872 function.remove(k + 1);
3873 let mut mat = Vec::new();
3874 for i in 0..a.len() {
3875 let mut vec = Vec::new();
3876 for j in 0..a[i].len() {
3877 vec.push(functions(
3878 a[i][j].clone(),
3879 Some(b[i].clone()),
3880 to_deg.clone(),
3881 s,
3882 options,
3883 )?)
3884 }
3885 mat.push(vec.clone());
3886 }
3887 return Ok(Matrix(mat));
3888 }
3889 (Vector(a), Matrix(b)) if a.len() == b.len() => {
3890 function.remove(k + 1);
3891 let mut mat = Vec::new();
3892 for i in 0..b.len() {
3893 let mut vec = Vec::new();
3894 for j in 0..b[i].len() {
3895 vec.push(functions(
3896 a[i].clone(),
3897 Some(b[i][j].clone()),
3898 to_deg.clone(),
3899 s,
3900 options,
3901 )?)
3902 }
3903 mat.push(vec.clone());
3904 }
3905 return Ok(Matrix(mat));
3906 }
3907 _ => {}
3908 }
3909 }
3910 match a {
3911 Matrix(a) => {
3912 let mut mat = Vec::new();
3913 for i in a {
3914 let mut vec = Vec::new();
3915 for j in i {
3916 vec.push(functions(j, None, to_deg.clone(), s, options)?)
3917 }
3918 mat.push(vec.clone());
3919 }
3920 Ok(Matrix(mat))
3921 }
3922 Vector(a) => {
3923 let mut vec = Vec::new();
3924 for i in a {
3925 vec.push(functions(i, None, to_deg.clone(), s, options)?)
3926 }
3927 Ok(Vector(vec))
3928 }
3929 Num(a) => Ok(NumStr::new(functions(
3930 *a,
3931 None,
3932 to_deg.clone(),
3933 s,
3934 options,
3935 )?)),
3936 _ => Err("str err1"),
3937 }
3938}
3939fn functions(
3940 mut a: Number,
3941 mut c: Option<Number>,
3942 to_deg: Complex,
3943 s: &str,
3944 options: Options,
3945) -> Result<Number, &'static str> {
3946 if a.number.imag().is_zero() && !a.number.imag().is_sign_positive() {
3947 a.number = Complex::with_val(a.number.prec(), a.number.real())
3948 }
3949 let n = if matches!(
3950 s,
3951 "root"
3952 | "sqrt"
3953 | "asquare"
3954 | "exp"
3955 | "aln"
3956 | "square"
3957 | "asqrt"
3958 | "cube"
3959 | "acbrt"
3960 | "asin"
3961 | "arcsin"
3962 | "acsc"
3963 | "arccsc"
3964 | "acos"
3965 | "arccos"
3966 | "asec"
3967 | "arcsec"
3968 | "atan2"
3969 | "hypot"
3970 | "atan"
3971 | "arctan"
3972 | "acot"
3973 | "arccot"
3974 | "ceil"
3975 | "floor"
3976 | "round"
3977 | "frac"
3978 | "fract"
3979 | "cbrt"
3980 | "acube"
3981 | "units"
3982 | "int"
3983 | "trunc"
3984 | "recip"
3985 | "abs"
3986 | "norm"
3987 | "onlyreal"
3988 | "onlyre"
3989 | "ore"
3990 | "onlyimag"
3991 | "onlyim"
3992 | "oim"
3993 | "re"
3994 | "real"
3995 | "im"
3996 | "imag"
3997 | "next"
3998 | "rand_norm"
3999 | "rand_uniform"
4000 | "rand_int"
4001 ) {
4002 if let Some(ref b) = c {
4003 if b.number.imag().is_zero() && !b.number.imag().is_sign_positive() {
4004 c = Some(Number::from(
4005 Complex::with_val(b.number.prec(), b.number.real()),
4006 b.units,
4007 ))
4008 }
4009 }
4010 match s {
4011 "sqrt" | "asquare" => Number::from(a.number.sqrt(), a.units.map(|a| a.root(2.0))),
4012 "root" => {
4013 if let Some(b) = c {
4014 root(&a, &b)
4015 } else {
4016 Number::from(a.number.sqrt(), a.units.map(|a| a.root(2.0)))
4017 }
4018 }
4019 "exp" | "aln" => {
4020 if let Some(b) = c {
4021 Number::from(
4022 pow_nth(a.number, b.number.clone()),
4023 a.units.map(|a| a.pow(b.number.real().to_f64())),
4024 )
4025 } else {
4026 Number::from(a.number.exp(), None)
4027 }
4028 }
4029 "square" | "asqrt" => Number::from(sqr(a.number), a.units.map(|a| a.pow(2.0))),
4030 "cube" | "acbrt" => Number::from(cube(a.number), a.units.map(|a| a.pow(3.0))),
4031 "asin" | "arcsin" => Number::from(
4032 a.number.clone().asin() * to_deg,
4033 Some(Units {
4034 angle: 1.0,
4035 ..Units::default()
4036 }),
4037 ),
4038 "acsc" | "arccsc" => Number::from(
4039 a.number.clone().recip().asin() * to_deg,
4040 Some(Units {
4041 angle: 1.0,
4042 ..Units::default()
4043 }),
4044 ),
4045 "acos" | "arccos" => Number::from(
4046 a.number.clone().acos() * to_deg,
4047 Some(Units {
4048 angle: 1.0,
4049 ..Units::default()
4050 }),
4051 ),
4052 "asec" | "arcsec" => Number::from(
4053 a.number.clone().recip().acos() * to_deg,
4054 Some(Units {
4055 angle: 1.0,
4056 ..Units::default()
4057 }),
4058 ),
4059 "atan2" => {
4060 if let Some(b) = c {
4061 Number::from(
4062 atan(b.number, a.number) * to_deg,
4063 Some(Units {
4064 angle: 1.0,
4065 ..Units::default()
4066 }),
4067 )
4068 } else {
4069 return Err("not enough args");
4070 }
4071 }
4072 "hypot" => {
4073 if let Some(b) = c {
4074 Number::from(
4075 (b.number.clone() * b.number + a.number.clone() * a.number).sqrt(),
4076 a.units,
4077 )
4078 } else {
4079 return Err("not enough args");
4080 }
4081 }
4082 "atan" | "arctan" => {
4083 if let Some(b) = c {
4084 Number::from(
4085 atan(a.number, b.number) * to_deg,
4086 Some(Units {
4087 angle: 1.0,
4088 ..Units::default()
4089 }),
4090 )
4091 } else {
4092 Number::from(
4093 a.number.atan() * to_deg,
4094 Some(Units {
4095 angle: 1.0,
4096 ..Units::default()
4097 }),
4098 )
4099 }
4100 }
4101 "acot" | "arccot" => Number::from(
4102 a.number.recip().atan() * to_deg,
4103 Some(Units {
4104 angle: 1.0,
4105 ..Units::default()
4106 }),
4107 ),
4108 "cbrt" | "acube" => Number::from(
4109 if a.number.imag().is_zero() {
4110 if a.number.real().is_zero() {
4111 Complex::new(options.prec)
4112 } else if a.number.real().is_sign_positive() {
4113 pow_nth(a.number, Complex::with_val(options.prec, 3).recip())
4114 } else {
4115 -pow_nth(-a.number, Complex::with_val(options.prec, 3).recip())
4116 }
4117 } else {
4118 pow_nth(a.number, Complex::with_val(options.prec, 3).recip())
4119 },
4120 a.units.map(|a| a.root(3.0)),
4121 ),
4122 "abs" | "norm" => Number::from(a.number.abs(), a.units),
4123 "recip" => Number::from(a.number.recip(), a.units.map(|a| a.pow(-1.0))),
4124 "units" => Number::from(Complex::with_val(options.prec, 1), a.units),
4125 "onlyreal" | "onlyre" | "ore" => {
4126 if -a.number.imag().clone().abs().log10() > a.number.prec().0 / 4 {
4127 Number::from(a.number.real().clone().into(), a.units)
4128 } else {
4129 Number::from(Complex::with_val(options.prec, Nan), None)
4130 }
4131 }
4132 "onlyimag" | "onlyim" | "oim" => {
4133 if -a.number.real().clone().abs().log10() > a.number.prec().0 / 4 {
4134 Number::from(a.number.imag().clone().into(), a.units)
4135 } else {
4136 Number::from(Complex::with_val(options.prec, Nan), None)
4137 }
4138 }
4139 "re" | "real" => Number::from(a.number.real().clone().into(), a.units),
4140 "im" | "imag" => Number::from(a.number.imag().clone().into(), a.units),
4141 "next" => {
4142 if let Some(b) = c {
4143 let mut real: Float = a.number.real().clone();
4144 let imag: Float = a.number.imag().clone();
4145 if b.number.real().is_infinite() {
4146 if b.number.real().is_sign_positive() {
4147 real.next_up()
4148 } else {
4149 real.next_down()
4150 }
4151 } else {
4152 real.next_toward(b.number.real());
4153 }
4154 Number::from(Complex::with_val(options.prec, (real, imag)), a.units)
4155 } else {
4156 let mut real: Float = a.number.real().clone();
4157 let imag: Float = a.number.imag().clone();
4158 real.next_up();
4159 Number::from(Complex::with_val(options.prec, (real, imag)), a.units)
4160 }
4161 }
4162 #[cfg(feature = "fastrand")]
4163 "rand_norm" => {
4164 if let Some(b) = c {
4165 Number::from(rand_norm(a.number, b.number), a.units)
4166 } else {
4167 return Err("not enough args");
4168 }
4169 }
4170 #[cfg(feature = "fastrand")]
4171 "rand_uniform" => {
4172 if let Some(b) = c {
4173 let units = a.units;
4174 let a = a.number;
4175 let b = b.number;
4176 Number::from(
4177 (b.clone() - a.clone()) * Float::with_val(options.prec, fastrand::u128(..))
4178 / u128::MAX
4179 + if a.real() < b.real() { a } else { b },
4180 units,
4181 )
4182 } else {
4183 return Err("not enough args");
4184 }
4185 }
4186 #[cfg(feature = "fastrand")]
4187 "rand_int" => {
4188 if let Some(b) = c {
4189 let units = a.units;
4190 let ar = a
4191 .number
4192 .real()
4193 .to_integer()
4194 .unwrap_or_default()
4195 .to_i128()
4196 .unwrap_or_default();
4197 let br = b
4198 .number
4199 .real()
4200 .to_integer()
4201 .unwrap_or_default()
4202 .to_i128()
4203 .unwrap_or_default();
4204 let ai = a
4205 .number
4206 .imag()
4207 .to_integer()
4208 .unwrap_or_default()
4209 .to_i128()
4210 .unwrap_or_default();
4211 let bi = b
4212 .number
4213 .imag()
4214 .to_integer()
4215 .unwrap_or_default()
4216 .to_i128()
4217 .unwrap_or_default();
4218 Number::from(
4219 Complex::with_val(
4220 options.prec,
4221 (
4222 fastrand::i128(ar.min(br)..=br.max(ar)),
4223 fastrand::i128(ai.min(bi)..=bi.max(ai)),
4224 ),
4225 ),
4226 units,
4227 )
4228 } else {
4229 return Err("not enough args");
4230 }
4231 }
4232 _ => return Err("unreachable"),
4233 }
4234 } else {
4235 let a = a.number;
4236 let mut d = None;
4237 if let Some(ref b) = c {
4238 if b.number.imag().is_zero() && !b.number.imag().is_sign_positive() {
4239 d = Some(Complex::with_val(b.number.prec(), b.number.real()))
4240 } else {
4241 d = Some(b.number.clone())
4242 }
4243 }
4244 Number::from(
4245 match s {
4246 "sin" => (a / to_deg).sin(),
4247 "csc" => (a / to_deg).sin().recip(),
4248 "cos" => (a / to_deg).cos(),
4249 "sec" => (a / to_deg).cos().recip(),
4250 "tan" => (a / to_deg).tan(),
4251 "cot" => (a / to_deg).tan().recip(),
4252 "sinh" => a.sinh(),
4253 "csch" => a.sinh().recip(),
4254 "cosh" => a.cosh(),
4255 "sech" => a.cosh().recip(),
4256 "tanh" => a.tanh(),
4257 "coth" => a.tanh().recip(),
4258 "asinh" | "arcsinh" => a.asinh(),
4259 "acsch" | "arccsch" => a.recip().asinh(),
4260 "acosh" | "arccosh" => a.acosh(),
4261 "asech" | "arcsech" => a.clone().recip().acosh(),
4262 "atanh" | "arctanh" => a.clone().atanh(),
4263 "acoth" | "arccoth" => a.clone().recip().atanh(),
4264 "cis" => {
4265 let b = a / to_deg.clone();
4266 b.clone().cos() + b.sin() * Complex::with_val(options.prec, (0.0, 1.0))
4267 }
4268 "ln" | "aexp" => a.ln(),
4269 "W" | "productlog" | "lambertw" => {
4270 if let Some(b) = d {
4271 lambertw(b, a.real().to_integer().unwrap_or_default())
4272 } else {
4273 lambertw(a, Integer::new())
4274 }
4275 }
4276 "log" => {
4277 let a = a.ln();
4278 if let Some(b) = d {
4279 let b = b.ln();
4280 if a.is_zero() {
4281 if b.is_zero() {
4282 Complex::with_val(options.prec, Nan)
4283 } else {
4284 Complex::with_val(options.prec, Infinity)
4285 }
4286 } else if b.real().is_infinite() {
4287 -Complex::with_val(options.prec, Infinity)
4288 } else {
4289 b / a
4290 }
4291 } else {
4292 a
4293 }
4294 }
4295 "ssrt" => {
4296 if let Some(b) = d {
4297 let b = b.ln();
4298 b.clone() / lambertw(b, a.real().to_integer().unwrap_or_default())
4299 } else {
4300 let a = a.ln();
4301 a.clone() / lambertw(a, Integer::new())
4302 }
4303 }
4304 "slog" => {
4305 if let Some(b) = d {
4306 if a.real() > &1 {
4307 slog(&a, &b)
4308 } else {
4309 Complex::with_val(options.prec, Nan)
4310 }
4311 } else {
4312 return Err("no args");
4313 }
4314 }
4315 "Ap" => {
4316 if let Some(b) = d {
4317 if a.real().is_sign_negative() {
4318 Complex::with_val(options.prec, Nan)
4319 } else {
4320 pub fn pow_n(z: Complex, n: usize) -> Complex {
4321 if !z.imag().is_zero() && n <= 256 {
4322 let mut p = z.clone();
4323 for _ in 1..n {
4324 p *= &z;
4325 }
4326 p
4327 } else {
4328 z.pow(n)
4329 }
4330 }
4331 let mut sum = Complex::new(options.prec);
4332 let n = a
4333 .real()
4334 .to_integer()
4335 .unwrap_or_default()
4336 .to_u32()
4337 .unwrap_or_default();
4338 for k in 0..=n {
4339 sum += pow_n(b.clone(), k as usize) * euleriannumbersint(n, k)
4340 }
4341 sum
4342 }
4343 } else {
4344 return Err("no args");
4345 }
4346 }
4347 "An" => {
4348 if let Some(b) = d {
4349 euleriannumbers(
4350 a,
4351 b.real()
4352 .to_integer()
4353 .unwrap_or_default()
4354 .to_i32()
4355 .unwrap_or_default(),
4356 )
4357 } else {
4358 return Err("no args");
4359 }
4360 }
4361 "P" => {
4362 if let Some(b) = d {
4363 if a.real().is_sign_negative()
4364 && a.real().clone().fract().is_zero()
4365 && a.imag().is_zero()
4366 && b.imag().is_zero()
4367 {
4368 let a = a + Complex::with_val(
4369 options.prec,
4370 (0, Float::with_val(options.prec, 0.5).pow(options.prec / 2)),
4371 );
4372 (gamma(a.clone() + 1) / gamma(a.clone() - b + 1))
4373 .real()
4374 .clone()
4375 .into()
4376 } else {
4377 gamma(a.clone() + 1) / gamma(a.clone() - b + 1)
4378 }
4379 } else {
4380 return Err("no args");
4381 }
4382 }
4383 "C" | "bi" | "binomial" => {
4384 if let Some(b) = d {
4385 binomial(a, b)
4386 } else {
4387 return Err("no args");
4388 }
4389 }
4390 "pochhammer" | "ph" => {
4391 if let Some(b) = d {
4392 if !a.real().is_sign_positive() && a.imag().is_zero() && b.imag().is_zero()
4393 {
4394 pow_nth(Complex::with_val(options.prec, -1), b.clone())
4395 * gamma(1 - a.clone())
4396 / gamma(1 - a - b)
4397 } else {
4398 gamma(b.clone() + a.clone()) / gamma(a.clone())
4399 }
4400 } else {
4401 return Err("not enough args");
4402 }
4403 }
4404 "lower_gamma" | "γ" => {
4405 if let Some(b) = d {
4406 lower_incomplete_gamma(a, b)
4407 } else {
4408 return Err("not enough args");
4409 }
4410 }
4411 "gamma" | "Γ" => {
4412 if let Some(b) = d {
4413 incomplete_gamma(a, b)
4414 } else {
4415 gamma(a)
4416 }
4417 }
4418 "sgn" | "sign" => {
4419 if a.is_zero() {
4420 Complex::new(options.prec)
4421 } else {
4422 a.clone() / a.abs()
4423 }
4424 }
4425 "arg" => a.arg(),
4426 "doublefact" | "doublefactorial" => {
4427 let two = Complex::with_val(options.prec, 2);
4428 let pi = Complex::with_val(options.prec, Pi);
4429 two.pow(a.clone() / 2 + (1 - (pi.clone() * a.clone()).cos()) / 4)
4430 * pi.clone().pow(((pi * a.clone()).cos() - 1) / 4)
4431 * gamma(a.clone() / 2 + 1)
4432 }
4433 "fact" | "factorial" => gamma(a.clone() + 1),
4434 "subfact" | "subfactorial" => {
4435 let near = a.real().clone().fract().is_zero();
4436 let neg = a.real().is_sign_negative();
4437 let neari = a.imag().is_zero();
4438 if !neari || neg || !near {
4439 if neg && near && neari {
4440 let a = a + Complex::with_val(
4441 options.prec,
4442 (0, Float::with_val(options.prec, 0.5).pow(options.prec / 8)),
4443 );
4444 subfactorial(a)
4445 } else {
4446 subfactorial(a)
4447 }
4448 } else if a.real().is_zero() {
4449 Complex::with_val(options.prec, 1)
4450 } else {
4451 (gamma(a.clone() + 1) / Float::with_val(options.prec, 1).exp())
4452 .real()
4453 .clone()
4454 .round()
4455 .into()
4456 }
4457 }
4458 "sinc" => a.clone().sin() / a,
4459 "conj" | "conjugate" => a.conj(),
4460 "erf" => {
4461 if a.imag().is_zero() {
4462 a.real().clone().erf().into()
4463 } else {
4464 erf(a)
4465 }
4466 }
4467 "erfc" => {
4468 if a.imag().is_zero() {
4469 a.real().clone().erfc().into()
4470 } else {
4471 erfc(a)
4472 }
4473 }
4474 "erfi" => {
4475 let i = Complex::with_val(options.prec, (0, 1));
4476 -i.clone() * erf(i * a)
4477 }
4478 "ai" => {
4479 if a.imag().is_zero() {
4480 a.real().clone().ai().into()
4481 } else {
4482 Complex::with_val(options.prec, Nan)
4483 }
4484 }
4485 "trigamma" => digamma(a, 1),
4486 "digamma" | "polygamma" | "ψ" => {
4487 if let Some(b) = d {
4488 digamma(
4489 b,
4490 a.real()
4491 .to_integer()
4492 .unwrap_or_default()
4493 .to_u32()
4494 .unwrap_or_default(),
4495 )
4496 } else if a.imag().is_zero() {
4497 if a.real().is_sign_negative() && a.real().clone().fract().is_zero() {
4498 Complex::with_val(options.prec, Infinity)
4499 } else {
4500 a.real().clone().digamma().into()
4501 }
4502 } else {
4503 digamma(a, 0)
4504 }
4505 }
4506 "eta" | "η" => eta(a),
4507 "zeta" | "ζ" => {
4508 if a.imag().is_zero() {
4509 a.real().clone().zeta().into()
4510 } else {
4511 zeta(a)
4512 }
4513 }
4514 "nth_prime" | "prime" => {
4515 if a.imag().is_zero() && a.real().clone().fract().is_zero() {
4516 Complex::with_val(
4517 options.prec,
4518 nth_prime(a.real().to_integer().unwrap_or_default()),
4519 )
4520 } else {
4521 Complex::with_val(options.prec, Nan)
4522 }
4523 }
4524 "mod" => {
4525 if let Some(b) = d {
4526 let c = a.clone() / b.clone();
4527 let c = Complex::with_val(
4528 options.prec,
4529 (c.real().clone().floor(), c.imag().clone().floor()),
4530 );
4531 a - b * c
4532 } else {
4533 return Err("not enough args");
4534 }
4535 }
4536 "lcm" => {
4537 if let Some(b) = d {
4538 if !a.real().is_finite()
4539 || !b.real().is_finite()
4540 || b.real().is_zero()
4541 || a.real().is_zero()
4542 {
4543 Complex::with_val(options.prec, Nan)
4544 } else {
4545 let a = a.real().clone().abs().to_integer().unwrap_or_default();
4546 let b = b.real().clone().abs().to_integer().unwrap_or_default();
4547 Complex::with_val(options.prec, a.clone() * b.clone() / gcd(a, b))
4548 }
4549 } else {
4550 return Err("not enough args");
4551 }
4552 }
4553 "gcd" | "gcf" => {
4554 if let Some(b) = d {
4555 if !a.real().is_finite()
4556 || !b.real().is_finite()
4557 || b.real().is_zero()
4558 || a.real().is_zero()
4559 {
4560 Complex::with_val(options.prec, Nan)
4561 } else {
4562 Complex::with_val(
4563 options.prec,
4564 gcd(
4565 a.real().clone().abs().to_integer().unwrap_or_default(),
4566 b.real().clone().abs().to_integer().unwrap_or_default(),
4567 ),
4568 )
4569 }
4570 } else {
4571 return Err("not enough args");
4572 }
4573 }
4574 "is_nan" => Complex::with_val(options.prec, a.real().is_nan() as u8),
4575 "is_inf" => Complex::with_val(options.prec, a.real().is_infinite() as u8),
4576 "is_fin" | "is_finite" => {
4577 Complex::with_val(options.prec, a.real().is_finite() as u8)
4578 }
4579 "isprime" | "is_prime" => {
4580 if a.imag().is_zero()
4581 && a.real().clone().fract().is_zero()
4582 && a.real().is_finite()
4583 {
4584 Complex::with_val(
4585 options.prec,
4586 (a.real()
4587 .to_integer()
4588 .unwrap_or_default()
4589 .is_probably_prime(100)
4590 != IsPrime::No) as u8,
4591 )
4592 } else {
4593 Complex::new(options.prec)
4594 }
4595 }
4596 #[cfg(feature = "fastrand")]
4597 "rand_gamma" => {
4598 if let Some(b) = d {
4599 rand_gamma(a.real().clone(), b.real().clone()).into()
4600 } else {
4601 return Err("not enough args");
4602 }
4603 }
4604 #[cfg(feature = "fastrand")]
4605 "rand_beta" => {
4606 if let Some(b) = d {
4607 let x = rand_gamma(a.real().clone(), Float::with_val(options.prec, 1));
4608 let y = rand_gamma(b.real().clone(), Float::with_val(options.prec, 1));
4609 (x.clone() / (x + y)).into()
4610 } else {
4611 return Err("not enough args");
4612 }
4613 }
4614 #[cfg(feature = "fastrand")]
4615 "rand_lognorm" => {
4616 if let Some(b) = d {
4617 rand_norm(a, b).exp()
4618 } else {
4619 return Err("not enough args");
4620 }
4621 }
4622 #[cfg(feature = "fastrand")]
4623 "rand_bernoulli" => {
4624 if *a.real() > Float::with_val(options.prec, fastrand::u128(..)) / u128::MAX {
4625 Complex::with_val(options.prec, 1)
4626 } else {
4627 Complex::new(options.prec)
4628 }
4629 }
4630 #[cfg(feature = "fastrand")]
4631 "rand_binomial" => {
4632 if let Some(b) = d {
4633 let p = b.real();
4634 let mut n = a.real().to_integer().unwrap_or_default();
4635 let mut sum = Integer::new();
4636 while n > 0 {
4637 if *p > Float::with_val(options.prec, fastrand::u128(..)) / u128::MAX {
4638 sum += 1;
4639 }
4640 n -= 1
4641 }
4642 Complex::with_val(options.prec, sum)
4643 } else {
4644 return Err("not enough args");
4645 }
4646 }
4647 #[cfg(feature = "fastrand")]
4648 "rand_neg_binomial" => {
4649 if let Some(b) = d {
4650 let p = b.real();
4651 if *p <= 0 {
4652 return Err("p must be greater then 0");
4653 }
4654 let mut r = a.real().to_integer().unwrap_or_default();
4655 let mut sum = Integer::new();
4656 while r > 0 {
4657 if *p > Float::with_val(options.prec, fastrand::u128(..)) / u128::MAX {
4658 r -= 1
4659 } else {
4660 sum += 1;
4661 }
4662 }
4663 Complex::with_val(options.prec, sum)
4664 } else {
4665 return Err("not enough args");
4666 }
4667 }
4668 #[cfg(feature = "fastrand")]
4669 "rand_geometric" => {
4670 let q: Float = 1 - a.real().clone();
4671 let n: Float = (Float::with_val(options.prec, fastrand::u128(..)) / u128::MAX)
4672 .ln()
4673 / q.ln();
4674 n.ceil().into()
4675 }
4676 #[cfg(feature = "fastrand")]
4677 "rand_poisson" => {
4678 let mut prod = Complex::with_val(options.prec, 1);
4679 let lim = (-a).exp();
4680 let mut n = Integer::new();
4681 while lim.real() < prod.real() {
4682 prod *= Float::with_val(options.prec, fastrand::u128(..)) / u128::MAX;
4683 n += 1;
4684 }
4685 Complex::with_val(options.prec, n - 1)
4686 }
4687 _ => {
4688 return Err("wrong input type");
4689 }
4690 },
4691 None,
4692 )
4693 };
4694 if n.number.imag().is_zero() && !n.number.imag().is_sign_positive() {
4695 Ok(Number::from(n.number.real().clone().into(), n.units))
4696 } else {
4697 Ok(n)
4698 }
4699}
4700pub fn compute_funcvars(
4701 function: &mut [NumStr],
4702 options: Options,
4703 func_vars: &mut Vec<(String, Vec<NumStr>)>,
4704) {
4705 let mut i = 0;
4706 while i < func_vars.len() {
4707 let v = func_vars[i].clone();
4708 let mut cont = false;
4709 for f in function.iter_mut() {
4710 if let Func(s) = &f {
4711 if *s == v.0 {
4712 cont = true;
4713 break;
4714 }
4715 }
4716 }
4717 if !cont && i + 1 < func_vars.len() {
4718 for fv in func_vars[i + 1..].iter_mut() {
4719 for f in fv.1.iter_mut() {
4720 if let Func(s) = &f {
4721 if *s == v.0 {
4722 cont = true;
4723 break;
4724 }
4725 }
4726 }
4727 }
4728 }
4729 if cont
4730 && (v.1.len() != 1
4731 || (if let Func(s) = &v.1[0] {
4732 matches!(s.as_str(), "rnd" | "rand" | "epoch")
4733 } else {
4734 false
4735 }))
4736 && !v.0.ends_with(')')
4737 {
4738 if let Ok(n) = do_math(v.1.clone(), options, func_vars[..i].to_vec()) {
4739 for f in function.iter_mut() {
4740 if let Func(s) = &f {
4741 if *s == v.0 {
4742 *f = n.clone();
4743 }
4744 }
4745 }
4746 if i + 1 < func_vars.len() {
4747 for fv in func_vars[i + 1..].iter_mut() {
4748 for f in fv.1.iter_mut() {
4749 if let Func(s) = &f {
4750 if *s == v.0 {
4751 *f = n.clone();
4752 }
4753 }
4754 }
4755 }
4756 }
4757 func_vars.remove(i);
4758 continue;
4759 }
4760 }
4761 i += 1;
4762 }
4763}