kalc_lib/
complex.rs

1use crate::{
2    complex::NumStr::{
3        Comma, Division, Exponent, Func, LeftBracket, LeftCurlyBracket, Matrix, Minus,
4        Multiplication, Num, Plus, RightBracket, RightCurlyBracket, Vector,
5    },
6    math::do_math,
7    misc::{do_math_with_var, place_funcvar, place_var},
8    parse::simplify,
9    units::{Number, Options, Units},
10};
11#[cfg(feature = "rayon")]
12use rayon::iter::{IntoParallelIterator, ParallelIterator};
13use rug::{
14    Complex, Float, Integer,
15    float::{
16        Constant::Pi,
17        Special::{Infinity, Nan},
18    },
19    integer::IsPrime,
20    ops::Pow,
21};
22#[cfg(feature = "serde")]
23use serde::{Deserialize, Serialize};
24use std::cmp::Ordering;
25#[derive(Clone, PartialEq, Debug)]
26#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
27pub enum NumStr {
28    Num(Box<Number>),
29    Func(String),
30    Vector(Vec<Number>),
31    Matrix(Vec<Vec<Number>>),
32    LeftBracket,
33    RightBracket,
34    LeftCurlyBracket,
35    RightCurlyBracket,
36    Comma,
37    Plus,
38    Minus,
39    PlusMinus,
40    Multiplication,
41    Division,
42    InternalMultiplication,
43    Tetration,
44    Root,
45    Exponent,
46    Equal,
47    NotEqual,
48    Greater,
49    GreaterEqual,
50    Lesser,
51    LesserEqual,
52    Modulo,
53    Range,
54    Conversion,
55    NearEqual,
56    ShiftLeft,
57    ShiftRight,
58    And,
59    Or,
60    Not,
61    Xor,
62    Nand,
63    Implies,
64    Nor,
65    Converse,
66}
67impl NumStr {
68    pub fn new(n: Number) -> Self {
69        Num(Box::new(n))
70    }
71}
72impl Number {
73    pub fn from(number: Complex, units: Option<Units>) -> Number {
74        Self { number, units }
75    }
76    pub fn set_prec(&mut self, prec: u32) {
77        self.number.set_prec(prec)
78    }
79    pub fn from_f64(f: f64, opt: &Options) -> Number {
80        let number = Complex::with_val(opt.prec, f);
81        let units = None;
82        Self { number, units }
83    }
84    pub fn new(opt: &Options) -> Number {
85        let number = Complex::new(opt.prec);
86        let units = None;
87        Self { number, units }
88    }
89    pub fn real(&self) -> &Float {
90        self.number.real()
91    }
92    pub fn imag(&self) -> &Float {
93        self.number.imag()
94    }
95}
96pub fn add(a: &Number, b: &Number) -> Number {
97    Number::from(
98        a.number.clone() + b.number.clone(),
99        if a.units == b.units { a.units } else { None },
100    )
101}
102pub fn sub(a: &Number, b: &Number) -> Number {
103    Number::from(
104        a.number.clone() - b.number.clone(),
105        if a.units == b.units { a.units } else { None },
106    )
107}
108pub fn set_prec(function: &mut [NumStr], func_vars: &mut [(String, Vec<NumStr>)], prec: u32) {
109    function.iter_mut().for_each(|n| n.set_prec(prec));
110    func_vars
111        .iter_mut()
112        .for_each(|(_, f)| f.iter_mut().for_each(|n| n.set_prec(prec)));
113}
114impl NumStr {
115    pub fn set_prec(&mut self, prec: u32) {
116        match self {
117            Num(n) => n.set_prec(prec),
118            Vector(v) => v.iter_mut().for_each(|n| n.set_prec(prec)),
119            Matrix(m) => m
120                .iter_mut()
121                .for_each(|v| v.iter_mut().for_each(|n| n.set_prec(prec))),
122            _ => {}
123        }
124    }
125    pub fn mul(&self, b: &Self) -> Result<Self, &'static str> {
126        fn m(a: &Number, b: &Number) -> Number {
127            Number::from(
128                {
129                    let a = a.number.clone();
130                    let b = b.number.clone();
131                    if a.real().is_infinite() || b.real().is_infinite() {
132                        if (a.real().is_infinite() && b.is_zero())
133                            || (b.real().is_infinite() && a.is_zero())
134                        {
135                            Complex::with_val(a.prec(), Nan)
136                        } else {
137                            match (a.real().is_sign_positive(), b.real().is_sign_positive()) {
138                                (true, true) | (false, false) => {
139                                    Complex::with_val(a.prec(), Infinity)
140                                }
141                                (false, true) | (true, false) => {
142                                    -Complex::with_val(a.prec(), Infinity)
143                                }
144                            }
145                        }
146                    } else {
147                        a * b.clone()
148                    }
149                },
150                match (a.units, b.units) {
151                    (Some(a), Some(b)) => Some(a.mul(&b)),
152                    (Some(a), None) | (None, Some(a)) => Some(a),
153                    (None, None) => None,
154                },
155            )
156        }
157        Ok(match (self, b) {
158            (Num(a), Num(b)) => NumStr::new(m(a, b)),
159            (Num(b), Vector(a)) | (Vector(a), Num(b)) => {
160                Vector(a.iter().map(|a| m(a, b)).collect())
161            }
162            (Vector(a), Vector(b)) if a.len() == b.len() => {
163                Vector(a.iter().zip(b.iter()).map(|(a, b)| m(a, b)).collect())
164            }
165            (Num(b), Matrix(a)) | (Matrix(a), Num(b)) => Matrix(
166                a.iter()
167                    .map(|a| a.iter().map(|a| m(a, b)).collect())
168                    .collect(),
169            ),
170            (Matrix(a), Vector(b)) if a.iter().all(|a| a.len() == b.len()) => Vector(
171                a.iter()
172                    .map(|a| {
173                        let mut iter = a.iter().zip(b.iter()).map(|(a, b)| m(a, b));
174                        let mut sum = iter.next().unwrap();
175                        for val in iter {
176                            sum = add(&sum, &val)
177                        }
178                        sum
179                    })
180                    .collect::<Vec<Number>>(),
181            ),
182            (Vector(a), Matrix(b)) if a.len() == b.len() => Vector(
183                transpose(b)
184                    .iter()
185                    .map(|b| {
186                        let mut iter = b.iter().zip(a.iter()).map(|(a, b)| m(a, b));
187                        let mut sum = iter.next().unwrap();
188                        for val in iter {
189                            sum = add(&sum, &val)
190                        }
191                        sum
192                    })
193                    .collect::<Vec<Number>>(),
194            ),
195            (Matrix(a), Matrix(b)) if a.iter().all(|a| a.len() == b.len()) => Matrix(
196                a.iter()
197                    .map(|a| {
198                        transpose(b)
199                            .iter()
200                            .map(|b| {
201                                let mut iter = a.iter().zip(b.iter()).map(|(a, b)| m(a, b));
202                                let mut sum = iter.next().unwrap();
203                                for val in iter {
204                                    sum = add(&sum, &val)
205                                }
206                                sum
207                            })
208                            .collect::<Vec<Number>>()
209                    })
210                    .collect(),
211            ),
212            _ => return Err("mul err"),
213        })
214    }
215    pub fn pm(&self, b: &Self) -> Result<Self, &'static str> {
216        Ok(match (self, b) {
217            (Num(a), Num(b)) => Vector(vec![add(a, b), sub(a, b)]),
218            (Num(a), Vector(b)) => Vector(
219                b.iter()
220                    .map(|b| add(a, b))
221                    .chain(b.iter().map(|b| sub(a, b)))
222                    .collect(),
223            ),
224            (Vector(b), Num(a)) => Vector(
225                b.iter()
226                    .map(|b| add(b, a))
227                    .chain(b.iter().map(|b| sub(b, a)))
228                    .collect(),
229            ),
230            (Vector(a), Vector(b)) if a.len() == b.len() => Vector(
231                a.iter()
232                    .zip(b.iter())
233                    .map(|(a, b)| add(a, b))
234                    .chain(a.iter().zip(b.iter()).map(|(a, b)| sub(a, b)))
235                    .collect(),
236            ),
237            (Matrix(a), Num(b)) => Vector(
238                a.iter()
239                    .flat_map(|a| {
240                        a.iter()
241                            .map(|a| add(a, b))
242                            .chain(a.iter().map(|a| sub(a, b)))
243                            .collect::<Vec<Number>>()
244                    })
245                    .collect::<Vec<Number>>(),
246            ),
247            (Num(b), Matrix(a)) => Vector(
248                a.iter()
249                    .flat_map(|a| {
250                        a.iter()
251                            .map(|a| add(b, a))
252                            .chain(a.iter().map(|a| sub(b, a)))
253                            .collect::<Vec<Number>>()
254                    })
255                    .collect::<Vec<Number>>(),
256            ),
257            _ => return Err("plus-minus unsupported"),
258        })
259    }
260    pub fn pow(&self, b: &Self) -> Result<Self, &'static str> {
261        fn p(a: &Number, b: &Number) -> Number {
262            Number::from(
263                {
264                    let a = a.number.clone();
265                    let b = b.number.clone();
266                    if a.real().is_infinite() {
267                        if b.is_zero() {
268                            Complex::with_val(a.prec(), Nan)
269                        } else if b.real().is_sign_positive() {
270                            Complex::with_val(a.prec(), Infinity)
271                        } else {
272                            Complex::new(a.prec())
273                        }
274                    } else if b.real().is_infinite() {
275                        if a.clone().abs() == 1 {
276                            Complex::with_val(a.prec(), Nan)
277                        } else if b.real().is_sign_positive() == a.real().clone().trunc().is_zero()
278                        {
279                            Complex::new(a.prec())
280                        } else {
281                            Complex::with_val(a.prec(), Infinity)
282                        }
283                    } else if a.is_zero() && b.real().is_zero() {
284                        Complex::with_val(a.prec(), Nan)
285                    } else {
286                        pow_nth(a, b)
287                    }
288                },
289                match (a.units, b.units) {
290                    (Some(a), None) => Some(a.pow(b.number.real().to_f64())),
291                    _ => None,
292                },
293            )
294        }
295        Ok(match (self, b) {
296            (Num(a), Num(b)) => NumStr::new(p(a, b)),
297            (Num(a), Vector(b)) => Vector(b.iter().map(|b| p(a, b)).collect()),
298            (Vector(a), Num(b)) => Vector(a.iter().map(|a| p(a, b)).collect()),
299            (Vector(a), Vector(b)) if a.len() == b.len() => {
300                Vector(a.iter().zip(b.iter()).map(|(a, b)| p(a, b)).collect())
301            }
302            (Num(a), Matrix(b)) => Matrix(
303                b.iter()
304                    .map(|b| b.iter().map(|b| p(a, b)).collect())
305                    .collect(),
306            ),
307            (Matrix(a), Num(b)) if a.iter().all(|c| a.len() == c.len()) => {
308                let b = b.number.clone();
309                if b.imag().is_zero() && b.real().clone().fract().is_zero() {
310                    if b.real().is_zero() {
311                        let mut mat = Vec::with_capacity(a.len().max(1));
312                        for i in 0..a.len() {
313                            let mut vec = Vec::with_capacity(a.len().max(1));
314                            for j in 0..a.len() {
315                                vec.push(if i == j {
316                                    Number::from(Complex::with_val(a[0][0].number.prec(), 1), None)
317                                } else {
318                                    Number::from(Complex::new(a[0][0].number.prec()), None)
319                                })
320                            }
321                            mat.push(vec);
322                        }
323                        Matrix(mat)
324                    } else {
325                        let mut mat = Matrix(a.clone());
326                        let c = b
327                            .real()
328                            .clone()
329                            .abs()
330                            .to_integer()
331                            .unwrap_or_default()
332                            .to_usize()
333                            .unwrap_or_default();
334                        for _ in 1..c {
335                            mat = mat.mul(&Matrix(a.clone()))?;
336                        }
337                        if b.real().is_sign_positive() {
338                            mat
339                        } else {
340                            Matrix(inverse(&mat.mat()?)?)
341                        }
342                    }
343                } else {
344                    return Err("no imag/fractional support for matrix powers");
345                }
346            }
347            (Vector(b), Matrix(a)) if b.len() == a.len() => Matrix(
348                a.iter()
349                    .zip(b.iter())
350                    .map(|(a, b)| a.iter().map(|a| p(b, a)).collect())
351                    .collect(),
352            ),
353            (Matrix(a), Vector(b)) if a.len() == b.len() => Matrix(
354                a.iter()
355                    .zip(b.iter())
356                    .map(|(a, b)| a.iter().map(|a| p(a, b)).collect())
357                    .collect(),
358            ),
359            (Matrix(a), Vector(b)) if a.iter().all(|a| a.len() == b.len()) => Matrix(
360                a.iter()
361                    .map(|a| a.iter().zip(b.iter()).map(|(a, b)| p(a, b)).collect())
362                    .collect(),
363            ),
364            (Matrix(a), Matrix(b)) if a[0].len() == b[0].len() && a.len() == b.len() => Matrix(
365                a.iter()
366                    .zip(b.iter())
367                    .map(|(a, b)| {
368                        a.iter()
369                            .zip(b.iter())
370                            .map(|(a, b)| p(a, b))
371                            .collect::<Vec<Number>>()
372                    })
373                    .collect(),
374            ),
375            _ => return Err("pow err"),
376        })
377    }
378    pub fn func<F>(&self, b: &Self, func: F) -> Result<Self, &'static str>
379    where
380        F: Fn(&Number, &Number) -> Number,
381    {
382        Ok(match (self, b) {
383            (Num(a), Num(b)) => NumStr::new(func(a, b)),
384            (Num(a), Vector(b)) => Vector(b.iter().map(|b| func(a, b)).collect()),
385            (Vector(a), Num(b)) => Vector(a.iter().map(|a| func(a, b)).collect()),
386            (Num(a), Matrix(b)) => Matrix(
387                b.iter()
388                    .map(|b| b.iter().map(|b| func(a, b)).collect())
389                    .collect(),
390            ),
391            (Matrix(a), Num(b)) => Matrix(
392                a.iter()
393                    .map(|a| a.iter().map(|a| func(a, b)).collect())
394                    .collect(),
395            ),
396            (Vector(a), Vector(b)) if a.len() == b.len() => {
397                Vector(a.iter().zip(b.iter()).map(|(a, b)| func(a, b)).collect())
398            }
399            (Vector(b), Matrix(a)) if b.len() == a.len() => Matrix(
400                a.iter()
401                    .zip(b.iter())
402                    .map(|(a, b)| a.iter().map(|a| func(b, a)).collect())
403                    .collect(),
404            ),
405            (Matrix(a), Vector(b)) if a.len() == b.len() => Matrix(
406                a.iter()
407                    .zip(b.iter())
408                    .map(|(a, b)| a.iter().map(|a| func(a, b)).collect())
409                    .collect(),
410            ),
411            (Matrix(a), Vector(b)) if a.iter().all(|a| a.len() == b.len()) => Matrix(
412                a.iter()
413                    .map(|a| a.iter().zip(b.iter()).map(|(a, b)| func(a, b)).collect())
414                    .collect(),
415            ),
416            (Matrix(a), Matrix(b)) if a[0].len() == b[0].len() && a.len() == b.len() => Matrix(
417                a.iter()
418                    .zip(b.iter())
419                    .map(|(a, b)| {
420                        a.iter()
421                            .zip(b.iter())
422                            .map(|(a, b)| func(a, b))
423                            .collect::<Vec<Number>>()
424                    })
425                    .collect(),
426            ),
427            _ => return Err("operation err"),
428        })
429    }
430    pub fn str_is(&self, s: &str) -> bool {
431        match self {
432            Func(s2) => s == s2,
433            _ => false,
434        }
435    }
436    pub fn num(&self) -> Result<Number, &'static str> {
437        match self {
438            Num(n) => Ok(*n.clone()),
439            _ => Err("failed to get number"),
440        }
441    }
442    pub fn vec(&self) -> Result<Vec<Number>, &'static str> {
443        match self {
444            Vector(v) => Ok(v.clone()),
445            _ => Err("failed to get vector"),
446        }
447    }
448    pub fn mat(&self) -> Result<Vec<Vec<Number>>, &'static str> {
449        match self {
450            Matrix(m) => Ok(m.clone()),
451            _ => Err("failed to get matrix"),
452        }
453    }
454}
455pub fn and(a: &Number, b: &Number) -> Number {
456    let a = &a.number;
457    let b = &b.number;
458    Number::from(
459        Complex::with_val(
460            a.prec(),
461            (a.imag().is_zero() && b.imag().is_zero() && a.real() == &1 && b.real() == &1) as u8,
462        ),
463        None,
464    )
465}
466pub fn or(a: &Number, b: &Number) -> Number {
467    let a = &a.number;
468    let b = &b.number;
469    Number::from(
470        Complex::with_val(
471            a.prec(),
472            (a.imag().is_zero() && b.imag().is_zero() && (a.real() == &1 || b.real() == &1)) as u8,
473        ),
474        None,
475    )
476}
477pub fn xor(a: &Number, b: &Number) -> Number {
478    let a = &a.number;
479    let b = &b.number;
480    Number::from(
481        Complex::with_val(
482            a.prec(),
483            (a.imag().is_zero() && b.imag().is_zero() && ((a.real() == &1) ^ (b.real() == &1)))
484                as u8,
485        ),
486        None,
487    )
488}
489pub fn nand(a: &Number, b: &Number) -> Number {
490    let a = &a.number;
491    let b = &b.number;
492    Number::from(
493        Complex::with_val(
494            a.prec(),
495            (a.imag().is_zero() && b.imag().is_zero() && (a.real() != &1 || b.real() != &1)) as u8,
496        ),
497        None,
498    )
499}
500pub fn nor(a: &Number, b: &Number) -> Number {
501    let a = &a.number;
502    let b = &b.number;
503    Number::from(
504        Complex::with_val(
505            a.prec(),
506            (a.imag().is_zero() && b.imag().is_zero() && (a.real() != &1 && b.real() != &1)) as u8,
507        ),
508        None,
509    )
510}
511pub fn implies(a: &Number, b: &Number) -> Number {
512    let a = &a.number;
513    let b = &b.number;
514    Number::from(
515        Complex::with_val(
516            a.prec(),
517            (a.imag().is_zero() && b.imag().is_zero() && (a.real() != &1 || b.real() == &1)) as u8,
518        ),
519        None,
520    )
521}
522pub fn not(a: &NumStr) -> Result<NumStr, &'static str> {
523    match a {
524        Num(a) => {
525            let a = &a.number;
526            Ok(NumStr::new(Number::from(
527                Complex::with_val(a.prec(), (a.imag().is_zero() && a.real() != &1) as u8),
528                None,
529            )))
530        }
531        Vector(v) => {
532            let mut o = Vec::with_capacity(v.len().max(1));
533            for a in v {
534                o.push(Number::from(
535                    Complex::with_val(a.number.prec(), a.number.is_zero() && a.number.real() != &1),
536                    None,
537                ));
538            }
539            Ok(Vector(o))
540        }
541        Matrix(m) => {
542            let mut k = Vec::with_capacity(m.len().max(1));
543            for v in m {
544                let mut o = Vec::with_capacity(v.len().max(1));
545                for a in v {
546                    o.push(Number::from(
547                        Complex::with_val(
548                            a.number.prec(),
549                            a.number.is_zero() && a.number.real() != &1,
550                        ),
551                        None,
552                    ));
553                }
554                k.push(o)
555            }
556            Ok(Matrix(k))
557        }
558        _ => Err("bad not input"),
559    }
560}
561pub fn div(a: &Number, b: &Number) -> Number {
562    Number::from(
563        {
564            let a = a.number.clone();
565            let b = b.number.clone();
566            if b.is_zero() || a.real().is_infinite() {
567                if a.is_zero() || b.real().is_infinite() {
568                    Complex::with_val(a.prec(), Nan)
569                } else if a.real().is_sign_positive() == b.real().is_sign_positive() {
570                    Complex::with_val(a.prec(), Infinity)
571                } else {
572                    -Complex::with_val(a.prec(), Infinity)
573                }
574            } else {
575                a / b.clone()
576            }
577        },
578        match (a.units, b.units) {
579            (Some(a), Some(b)) => Some(a.div(&b)),
580            (Some(a), None) => Some(a),
581            (None, Some(b)) => Some(Units::default().div(&b)),
582            (None, None) => None,
583        },
584    )
585}
586pub fn root(a: &Number, b: &Number) -> Number {
587    Number::from(
588        {
589            let a = a.number.clone();
590            let b = b.number.clone();
591            if b == 2 {
592                a.sqrt()
593            } else if b.imag().is_zero()
594                && b.real().clone() % 2 == 0
595                && b.real().clone().fract().is_zero()
596                && a.imag().is_zero()
597                && a.real().is_sign_negative()
598            {
599                -pow_nth(-a, b.recip())
600            } else {
601                pow_nth(a, b.recip())
602            }
603        },
604        match (a.units, b.units) {
605            (Some(a), None) => Some(a.root(b.number.real().to_f64())),
606            _ => None,
607        },
608    )
609}
610pub fn unity(y: Complex, x: Complex) -> Vec<Number> {
611    if y.is_zero() {
612        let n = x
613            .into_real_imag()
614            .0
615            .trunc()
616            .to_integer()
617            .unwrap_or_default()
618            .to_i128()
619            .unwrap_or_default() as usize;
620        return vec![Number::from(Complex::new(y.prec()), None); n];
621    }
622    let y = y.ln();
623    let mut vec: Vec<Number> = Vec::new();
624    let taui: Complex = 2 * Complex::with_val(x.prec(), (0, Pi));
625    let r: Float = x.imag().clone().pow(2) / 2;
626    let r: Float = x.real().clone() / 2 + r / x.real().clone();
627    let n = (x.imag().clone() * y.real() / x.real() - y.imag()) / taui.imag().clone();
628    let left: Float = n.clone() - r.clone();
629    let right: Float = r + n;
630    for k in if left < right {
631        left.clone()
632            .trunc()
633            .to_integer()
634            .unwrap_or_default()
635            .to_i128()
636            .unwrap_or_default()
637            ..=if left.clone().fract().is_zero() && right.clone().fract().is_zero() {
638                right
639                    .clone()
640                    .trunc()
641                    .to_integer()
642                    .unwrap_or_default()
643                    .to_i128()
644                    .unwrap_or_default()
645                    - 1
646            } else {
647                right
648                    .clone()
649                    .trunc()
650                    .to_integer()
651                    .unwrap_or_default()
652                    .to_i128()
653                    .unwrap_or_default()
654            }
655    } else {
656        right
657            .clone()
658            .trunc()
659            .to_integer()
660            .unwrap_or_default()
661            .to_i128()
662            .unwrap_or_default()
663            ..=if left.clone().fract().is_zero() && right.clone().fract().is_zero() {
664                left.clone()
665                    .trunc()
666                    .to_integer()
667                    .unwrap_or_default()
668                    .to_i128()
669                    .unwrap_or_default()
670                    - 1
671            } else {
672                left.clone()
673                    .trunc()
674                    .to_integer()
675                    .unwrap_or_default()
676                    .to_i128()
677                    .unwrap_or_default()
678            }
679    } {
680        let r: Complex = (y.clone() + k * taui.clone()) / x.clone();
681        let r: Complex = r.exp();
682        vec.push(Number::from(r, None))
683    }
684    vec
685}
686pub fn shl(a: &Number, b: &Number) -> Number {
687    Number::from(
688        a.number.clone() * pow_nth(Complex::with_val(a.number.prec(), 2), b.number.clone()),
689        None,
690    )
691}
692pub fn shr(a: &Number, b: &Number) -> Number {
693    Number::from(
694        a.number.clone() * pow_nth(Complex::with_val(a.number.prec(), 2), -b.number.clone()),
695        None,
696    )
697}
698pub fn ne(a: &Number, b: &Number) -> Number {
699    let ua = a.units;
700    let ub = b.units;
701    let a = a.number.clone();
702    let b = b.number.clone();
703    let c: Complex = a.clone() - b.clone();
704    let int = Integer::from(10).pow(a.prec().0 / 4);
705    let re: Float = c.real().clone() * int.clone();
706    let re: Float = re.round() / int.clone();
707    let im: Float = c.imag().clone() * int.clone();
708    let im: Float = im.round() / int;
709    Number::from(
710        Complex::with_val(
711            a.prec(),
712            ((!(re.is_zero()
713                || (a.real().is_infinite()
714                    && b.real().is_infinite()
715                    && a.real().is_sign_positive() == b.real().is_sign_positive()))
716                || !(im.is_zero()
717                    || (a.imag().is_infinite()
718                        && b.imag().is_infinite()
719                        && a.imag().is_sign_positive() == b.imag().is_sign_positive())))
720                || match (ua, ub) {
721                    (Some(a), Some(b)) => a != b,
722                    (Some(a), None) | (None, Some(a)) => !a.is_none(),
723                    (None, None) => false,
724                }) as u8,
725        ),
726        None,
727    )
728}
729pub fn eq(a: &Number, b: &Number) -> Number {
730    let ua = a.units;
731    let ub = b.units;
732    let a = a.number.clone();
733    let b = b.number.clone();
734    let c: Complex = a.clone() - b.clone();
735    let int = Integer::from(10).pow(a.prec().0 / 4);
736    let re: Float = c.real().clone() * int.clone();
737    let re: Float = re.round() / int.clone();
738    let im: Float = c.imag().clone() * int.clone();
739    let im: Float = im.round() / int;
740    Number::from(
741        Complex::with_val(
742            a.prec(),
743            ((re.is_zero()
744                || (a.real().is_infinite()
745                    && b.real().is_infinite()
746                    && a.real().is_sign_positive() == b.real().is_sign_positive())
747                    && im.is_zero()
748                || (a.imag().is_infinite()
749                    && b.imag().is_infinite()
750                    && a.imag().is_sign_positive() == b.imag().is_sign_positive()))
751                && match (ua, ub) {
752                    (Some(a), Some(b)) => a == b,
753                    (Some(a), None) | (None, Some(a)) => a.is_none(),
754                    (None, None) => true,
755                }) as u8,
756        ),
757        None,
758    )
759}
760pub fn about_eq(a: &Number, b: &Number) -> Number {
761    let ua = a.units;
762    let ub = b.units;
763    let a = a.number.clone();
764    let b = b.number.clone();
765    Number::from(
766        Complex::with_val(
767            a.prec(),
768            ((a.real().is_sign_positive() == b.real().is_sign_positive()
769                && a.imag().is_sign_positive() == b.imag().is_sign_positive()
770                && a.real().clone().abs().log10().floor()
771                    == b.real().clone().abs().log10().floor()
772                && a.imag().clone().abs().log10().floor()
773                    == b.imag().clone().abs().log10().floor())
774                && match (ua, ub) {
775                    (Some(a), Some(b)) => a == b,
776                    (Some(a), None) | (None, Some(a)) => a.is_none(),
777                    (None, None) => true,
778                }) as u8,
779        ),
780        None,
781    )
782}
783pub fn ge(a: &Number, b: &Number) -> Number {
784    Number::from(
785        Complex::with_val(a.number.prec(), (a.number.real() >= b.number.real()) as u8),
786        None,
787    )
788}
789pub fn gt(a: &Number, b: &Number) -> Number {
790    Number::from(
791        Complex::with_val(a.number.prec(), (a.number.real() > b.number.real()) as u8),
792        None,
793    )
794}
795pub fn rem(a: &Number, b: &Number) -> Number {
796    let a = &a.number;
797    let b = &b.number;
798    let c = a.clone() / b.clone();
799    let c = Complex::with_val(
800        a.prec(),
801        (c.real().clone().floor(), c.imag().clone().floor()),
802    );
803    Number::from(a - b * c, None)
804}
805pub fn digamma(mut z: Complex, mut n: u32) -> Complex {
806    n += 1;
807    let oop = z.prec();
808    let op = oop.0 / 2;
809    let prec = n * op;
810    z.set_prec(prec);
811    let h: Float = Float::with_val(prec, 0.5).pow(op / 2);
812    let num = Integer::from(n);
813    let mut sum = Complex::new(prec);
814    for k in 0..=n {
815        if k % 2 == 0 {
816            sum += num.clone().binomial(k) * gamma(z.clone() + h.clone() * (n - k)).ln()
817        } else {
818            sum -= num.clone().binomial(k) * gamma(z.clone() + h.clone() * (n - k)).ln()
819        }
820    }
821    let mut n = sum * Float::with_val(prec, 2).pow(op / 2 * n);
822    n.set_prec(oop);
823    n
824}
825pub fn gamma(a: Complex) -> Complex {
826    if !a.imag().is_zero() {
827        gamma0(a)
828    } else if a.real().is_sign_negative() && a.real().clone().fract().is_zero() {
829        Complex::with_val(a.prec(), Infinity)
830    } else {
831        a.real().clone().gamma().into()
832    }
833}
834pub fn tetration(a: &Number, b: &Number) -> Number {
835    let a = a.number.clone();
836    let b = b.number.clone();
837    Number::from(
838        if b.real().clone().fract().is_zero() {
839            if b == 0 {
840                Complex::with_val(a.prec(), 1)
841            } else if b.real().is_sign_positive() {
842                (1..b
843                    .real()
844                    .to_integer()
845                    .unwrap_or_default()
846                    .to_usize()
847                    .unwrap_or_default())
848                    .fold(a.clone(), |tetration, _| pow_nth(a.clone(), tetration))
849            } else if b == -1 {
850                Complex::new(a.prec())
851            } else {
852                Complex::with_val(a.prec(), (Infinity, Nan))
853            }
854        } else {
855            tetration_recursion(a.clone(), b.clone())
856        },
857        None,
858    )
859}
860fn tetration_recursion(a: Complex, b: Complex) -> Complex {
861    if b.real().is_sign_positive() {
862        pow_nth(a.clone(), tetration_recursion(a, b - 1))
863    } else if b.real() <= &-1 {
864        tetration_recursion(a.clone(), b + 1).ln() / a.ln()
865    } else {
866        let aln = a.abs().clone().ln();
867        1 + b.clone() * (aln.clone() * (2 + b.clone()) - b) / (1 + aln)
868    }
869}
870pub fn slog(a: &Complex, b: &Complex) -> Complex {
871    if b.real() <= &0 {
872        let z = &pow_nth(a.clone(), b.clone());
873        if z.real() <= b.real() {
874            Complex::with_val(a.prec(), Nan)
875        } else {
876            slog(a, z) - 1
877        }
878    } else if b.real() > &1 {
879        let z = &(b.clone().ln() / a.clone().ln());
880        if z.real() >= b.real() {
881            Complex::with_val(a.prec(), Nan)
882        } else {
883            slog(a, z) + 1
884        }
885    } else {
886        let a = a.clone().ln();
887        b.clone() * (2 * a.clone() + b * (1 - a.clone())) / (1 + a) - 1
888    }
889}
890pub fn atan(a: Complex, b: Complex) -> Complex {
891    if a.imag().is_zero() && b.imag().is_zero() {
892        Complex::with_val(a.prec(), (a.real(), b.real())).arg()
893    } else {
894        let i = Complex::with_val(a.prec(), (0, 1));
895        let abs: Complex = sqr(a.clone()) + sqr(b.clone());
896        -i.clone() * ((a + b * i) / abs.sqrt()).ln()
897    }
898}
899pub fn to_polar(mut a: Vec<Number>, to_deg: Complex) -> Vec<Number> {
900    if a.len() == 1 {
901        a.push(Number::from(Complex::new(a[0].number.prec()), None));
902    }
903    if a.len() != 2 && a.len() != 3 {
904        a
905    } else if a.len() == 2 {
906        if a[1].number.is_zero() {
907            if a[0].number.is_zero() {
908                vec![
909                    Number::from(Complex::new(a[0].number.prec()), None),
910                    Number::from(
911                        Complex::new(a[0].number.prec()),
912                        Some(Units {
913                            angle: 1.0,
914                            ..Units::default()
915                        }),
916                    ),
917                ]
918            } else {
919                vec![
920                    Number::from(a[0].number.clone().abs(), a[0].units),
921                    Number::from(
922                        if a[0].number.real().is_sign_positive() {
923                            Complex::new(a[0].number.prec())
924                        } else {
925                            to_deg * Float::with_val(a[0].number.prec().0, Pi)
926                        },
927                        Some(Units {
928                            angle: 1.0,
929                            ..Units::default()
930                        }),
931                    ),
932                ]
933            }
934        } else {
935            let mut n: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
936            n = n.sqrt();
937            vec![
938                Number::from(n.clone(), a[0].units),
939                Number::from(
940                    atan(a[0].number.clone(), a[1].number.clone()) * to_deg,
941                    Some(Units {
942                        angle: 1.0,
943                        ..Units::default()
944                    }),
945                ),
946            ]
947        }
948    } else if a[1].number.is_zero() {
949        if a[0].number.is_zero() {
950            if a[2].number.is_zero() {
951                vec![
952                    Number::from(Complex::new(a[0].number.prec()), None),
953                    Number::from(
954                        Complex::new(a[0].number.prec()),
955                        Some(Units {
956                            angle: 1.0,
957                            ..Units::default()
958                        }),
959                    ),
960                    Number::from(
961                        Complex::new(a[0].number.prec()),
962                        Some(Units {
963                            angle: 1.0,
964                            ..Units::default()
965                        }),
966                    ),
967                ]
968            } else {
969                vec![
970                    Number::from(a[2].number.clone().abs(), a[2].units),
971                    Number::from(
972                        Complex::new(a[0].number.prec()),
973                        Some(Units {
974                            angle: 1.0,
975                            ..Units::default()
976                        }),
977                    ),
978                    Number::from(
979                        Complex::new(a[0].number.prec()),
980                        Some(Units {
981                            angle: 1.0,
982                            ..Units::default()
983                        }),
984                    ),
985                ]
986            }
987        } else {
988            let nxy: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
989            let mut n: Complex = nxy.clone() + sqr(a[2].number.clone());
990            n = n.sqrt();
991            vec![
992                Number::from(n.clone(), a[0].units),
993                Number::from(
994                    atan(a[2].number.clone(), nxy.sqrt()) * to_deg.clone(),
995                    Some(Units {
996                        angle: 1.0,
997                        ..Units::default()
998                    }),
999                ),
1000                Number::from(
1001                    Complex::new(a[0].number.prec()),
1002                    Some(Units {
1003                        angle: 1.0,
1004                        ..Units::default()
1005                    }),
1006                ),
1007            ]
1008        }
1009    } else {
1010        let nxy: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
1011        let mut n: Complex = nxy.clone() + sqr(a[2].number.clone());
1012        n = n.sqrt();
1013        vec![
1014            Number::from(n.clone(), a[0].units),
1015            Number::from(
1016                atan(a[2].number.clone(), nxy.sqrt()) * to_deg.clone(),
1017                Some(Units {
1018                    angle: 1.0,
1019                    ..Units::default()
1020                }),
1021            ),
1022            Number::from(
1023                atan(a[0].number.clone(), a[1].number.clone()) * to_deg.clone(),
1024                Some(Units {
1025                    angle: 1.0,
1026                    ..Units::default()
1027                }),
1028            ),
1029        ]
1030    }
1031}
1032pub fn to_cyl(mut a: Vec<Number>, to_deg: Complex) -> Vec<Number> {
1033    if a.len() == 1 {
1034        a.push(Number::from(Complex::new(a[0].number.prec()), None));
1035    }
1036    if a.len() != 2 && a.len() != 3 {
1037        a
1038    } else if a.len() == 2 {
1039        if a[1].number.is_zero() {
1040            if a[0].number.is_zero() {
1041                vec![
1042                    Number::from(Complex::new(a[0].number.prec()), None),
1043                    Number::from(
1044                        Complex::new(a[0].number.prec()),
1045                        Some(Units {
1046                            angle: 1.0,
1047                            ..Units::default()
1048                        }),
1049                    ),
1050                ]
1051            } else {
1052                vec![
1053                    Number::from(a[0].number.clone().abs(), a[0].units),
1054                    Number::from(
1055                        if a[0].number.real().is_sign_positive() {
1056                            Complex::new(a[0].number.prec())
1057                        } else {
1058                            to_deg * Float::with_val(a[0].number.prec().0, Pi)
1059                        },
1060                        Some(Units {
1061                            angle: 1.0,
1062                            ..Units::default()
1063                        }),
1064                    ),
1065                ]
1066            }
1067        } else {
1068            let mut n: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
1069            n = n.sqrt();
1070            vec![
1071                Number::from(n.clone(), a[0].units),
1072                Number::from(
1073                    atan(a[0].number.clone(), a[1].number.clone()) * to_deg,
1074                    Some(Units {
1075                        angle: 1.0,
1076                        ..Units::default()
1077                    }),
1078                ),
1079            ]
1080        }
1081    } else if a[1].number.is_zero() {
1082        if a[0].number.is_zero() {
1083            if a[2].number.is_zero() {
1084                vec![
1085                    Number::from(Complex::new(a[0].number.prec()), None),
1086                    Number::from(
1087                        Complex::new(a[0].number.prec()),
1088                        Some(Units {
1089                            angle: 1.0,
1090                            ..Units::default()
1091                        }),
1092                    ),
1093                    Number::from(Complex::new(a[0].number.prec()), None),
1094                ]
1095            } else {
1096                vec![
1097                    Number::from(Complex::new(a[0].number.prec()), None),
1098                    Number::from(
1099                        Complex::new(a[0].number.prec()),
1100                        Some(Units {
1101                            angle: 1.0,
1102                            ..Units::default()
1103                        }),
1104                    ),
1105                    Number::from(a[2].number.clone().abs(), a[2].units),
1106                ]
1107            }
1108        } else {
1109            let nxy: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
1110            vec![
1111                Number::from(nxy.sqrt(), a[0].units),
1112                Number::from(
1113                    atan(a[0].number.clone(), a[1].number.clone()) * to_deg.clone(),
1114                    Some(Units {
1115                        angle: 1.0,
1116                        ..Units::default()
1117                    }),
1118                ),
1119                a[2].clone(),
1120            ]
1121        }
1122    } else {
1123        let nxy: Complex = sqr(a[0].number.clone()) + sqr(a[1].number.clone());
1124        vec![
1125            Number::from(nxy.sqrt(), a[0].units),
1126            Number::from(
1127                atan(a[0].number.clone(), a[1].number.clone()) * to_deg.clone(),
1128                Some(Units {
1129                    angle: 1.0,
1130                    ..Units::default()
1131                }),
1132            ),
1133            a[2].clone(),
1134        ]
1135    }
1136}
1137pub fn to(a: &NumStr, b: &NumStr) -> Result<NumStr, &'static str> {
1138    Ok(match (a, b) {
1139        (Num(a), Num(b)) => {
1140            let prec = a.number.prec();
1141            let a = a
1142                .number
1143                .real()
1144                .to_integer()
1145                .unwrap_or_default()
1146                .to_isize()
1147                .unwrap_or_default();
1148            let b = b
1149                .number
1150                .real()
1151                .to_integer()
1152                .unwrap_or_default()
1153                .to_isize()
1154                .unwrap_or_default();
1155            let vec: Vec<Number> = if a < b {
1156                (a..=b)
1157                    .map(|a| Number::from(Complex::with_val(prec, a), None))
1158                    .collect()
1159            } else {
1160                (b..=a)
1161                    .rev()
1162                    .map(|a| Number::from(Complex::with_val(prec, a), None))
1163                    .collect()
1164            };
1165            if vec.is_empty() {
1166                return Err("start range greater then end range");
1167            }
1168            Vector(vec)
1169        }
1170        (Vector(a), Num(b)) => {
1171            let prec = b.number.prec();
1172            let b = b
1173                .number
1174                .real()
1175                .to_integer()
1176                .unwrap_or_default()
1177                .to_isize()
1178                .unwrap_or_default();
1179            let mat: Vec<Vec<Number>> = a
1180                .iter()
1181                .map(|a| {
1182                    let a = a
1183                        .number
1184                        .real()
1185                        .to_integer()
1186                        .unwrap_or_default()
1187                        .to_isize()
1188                        .unwrap_or_default();
1189                    if a < b {
1190                        (a..=b)
1191                            .map(|a| Number::from(Complex::with_val(prec, a), None))
1192                            .collect()
1193                    } else {
1194                        (b..=a)
1195                            .rev()
1196                            .map(|a| Number::from(Complex::with_val(prec, a), None))
1197                            .collect()
1198                    }
1199                })
1200                .collect();
1201            if mat.is_empty() || mat.iter().any(|vec| vec.is_empty()) {
1202                return Err("start range greater then end range");
1203            }
1204            Matrix(mat)
1205        }
1206        (Num(a), Vector(b)) => {
1207            let prec = a.number.prec();
1208            let a = a
1209                .number
1210                .real()
1211                .to_integer()
1212                .unwrap_or_default()
1213                .to_isize()
1214                .unwrap_or_default();
1215            let mat: Vec<Vec<Number>> = b
1216                .iter()
1217                .map(|b| {
1218                    let b = b
1219                        .number
1220                        .real()
1221                        .to_integer()
1222                        .unwrap_or_default()
1223                        .to_isize()
1224                        .unwrap_or_default();
1225                    if a < b {
1226                        (a..=b)
1227                            .map(|a| Number::from(Complex::with_val(prec, a), None))
1228                            .collect()
1229                    } else {
1230                        (b..=a)
1231                            .rev()
1232                            .map(|a| Number::from(Complex::with_val(prec, a), None))
1233                            .collect()
1234                    }
1235                })
1236                .collect();
1237            if mat.is_empty() || mat.iter().any(|vec| vec.is_empty()) {
1238                return Err("start range greater then end range");
1239            }
1240            Matrix(mat)
1241        }
1242        _ => return Err(".. err"),
1243    })
1244}
1245pub fn mvec(
1246    function: Vec<NumStr>,
1247    func_vars: Vec<(String, Vec<NumStr>)>,
1248    var: &str,
1249    start: isize,
1250    end: isize,
1251    mvec: bool,
1252    options: Options,
1253) -> Result<NumStr, &'static str> {
1254    let mut vec = Vec::new();
1255    let mut mat = Vec::new();
1256    let mut test = true;
1257    let mut body = |z: isize| -> Result<(), &'static str> {
1258        match do_math_with_var(
1259            function.clone(),
1260            options,
1261            func_vars.clone(),
1262            var,
1263            NumStr::new(Number::from(Complex::with_val(options.prec, z), None)),
1264        )? {
1265            Num(n) => {
1266                if test && vec.is_empty() {
1267                    test = false;
1268                    vec = Vec::with_capacity((end - start + 1).max(1) as usize)
1269                }
1270                vec.push(*n)
1271            }
1272            Vector(v) if mvec => {
1273                if test && vec.is_empty() {
1274                    test = false;
1275                    vec = Vec::with_capacity((end - start + 1).max(1) as usize * v.len())
1276                }
1277                vec.extend(v)
1278            }
1279            Vector(v) => {
1280                if test && mat.is_empty() {
1281                    test = false;
1282                    mat = Vec::with_capacity((end - start + 1).max(1) as usize)
1283                }
1284                mat.push(v)
1285            }
1286            Matrix(m) if !mvec => {
1287                if test && mat.is_empty() {
1288                    test = false;
1289                    mat = Vec::with_capacity((end - start + 1).max(1) as usize * m.len().max(1))
1290                }
1291                mat.extend(m)
1292            }
1293            _ => return Err("cant create 3d matrix"),
1294        }
1295        Ok(())
1296    };
1297    if start < end {
1298        for z in start..=end {
1299            body(z)?
1300        }
1301    } else {
1302        for z in (end..=start).rev() {
1303            body(z)?
1304        }
1305    }
1306    if mat.is_empty() {
1307        if vec.is_empty() {
1308            Err("start>end")
1309        } else {
1310            Ok(Vector(vec))
1311        }
1312    } else {
1313        Ok(Matrix(mat))
1314    }
1315}
1316pub fn sum(
1317    function: Vec<NumStr>,
1318    func_vars: Vec<(String, Vec<NumStr>)>,
1319    var: &str,
1320    mut start: Float,
1321    mut end: Float,
1322    product: bool,
1323    options: Options,
1324) -> Result<NumStr, &'static str> {
1325    if start > end {
1326        (start, end) = (end, start)
1327    }
1328    if end.is_infinite() {
1329        let mut j = 0;
1330        let mut n = start
1331            .to_integer()
1332            .unwrap_or_default()
1333            .to_isize()
1334            .unwrap_or_default();
1335        let k = if end.is_sign_positive() { 1 } else { -1 };
1336        let mut last = NumStr::new(Number::from(Complex::with_val(options.prec, Nan), None));
1337        let mut value = do_math_with_var(
1338            function.clone(),
1339            options,
1340            func_vars.clone(),
1341            var,
1342            NumStr::new(Number::from(Complex::with_val(options.prec, n), None)),
1343        )?;
1344        while last != value {
1345            if j > 10000 {
1346                return Ok(NumStr::new(Number::from(
1347                    Complex::with_val(options.prec, Nan),
1348                    None,
1349                )));
1350            }
1351            j += 1;
1352            n += k;
1353            last = value;
1354            let math = do_math_with_var(
1355                function.clone(),
1356                options,
1357                func_vars.clone(),
1358                var,
1359                NumStr::new(Number::from(Complex::with_val(options.prec, n), None)),
1360            )?;
1361            if product {
1362                value = last.mul(&math)?;
1363            } else {
1364                value = last.func(&math, add)?;
1365            }
1366        }
1367        Ok(value)
1368    } else if start.is_infinite() {
1369        return Err("unsupported due to lack of example for convergence");
1370    } else {
1371        let start = start
1372            .to_integer()
1373            .unwrap_or_default()
1374            .to_isize()
1375            .unwrap_or_default();
1376        let end = end
1377            .to_integer()
1378            .unwrap_or_default()
1379            .to_isize()
1380            .unwrap_or_default();
1381        let mut value = do_math_with_var(
1382            function.clone(),
1383            options,
1384            func_vars.clone(),
1385            var,
1386            NumStr::new(Number::from(Complex::with_val(options.prec, start), None)),
1387        )?;
1388        for z in start + 1..=end {
1389            let math = do_math_with_var(
1390                function.clone(),
1391                options,
1392                func_vars.clone(),
1393                var,
1394                NumStr::new(Number::from(Complex::with_val(options.prec, z), None)),
1395            )?;
1396            if product {
1397                value = value.mul(&math)?;
1398            } else {
1399                value = value.func(&math, add)?;
1400            }
1401        }
1402        Ok(value)
1403    }
1404}
1405pub fn submatrix(a: &[Vec<Number>], row: usize, col: usize) -> Vec<Vec<Number>> {
1406    a.iter()
1407        .enumerate()
1408        .filter(|&(i, _)| i != row)
1409        .map(|(_, r)| {
1410            r.iter()
1411                .enumerate()
1412                .filter(|&(j, _)| j != col)
1413                .map(|(_, value)| value.clone())
1414                .collect::<Vec<Number>>()
1415        })
1416        .collect()
1417}
1418pub fn trace(a: &[Vec<Number>]) -> Number {
1419    let mut n = Complex::new(a[0][0].number.prec());
1420    for (i, j) in a.iter().enumerate() {
1421        if j.len() == i {
1422            break;
1423        }
1424        n += j[i].number.clone();
1425    }
1426    Number::from(n, a[0][0].units)
1427}
1428pub fn identity(a: usize, prec: u32) -> Vec<Vec<Number>> {
1429    let mut mat = Vec::with_capacity(a.max(1));
1430    for i in 0..a {
1431        let mut vec = Vec::with_capacity(a.max(1));
1432        for j in 0..a {
1433            if i == j {
1434                vec.push(Number::from(Complex::with_val(prec, 1), None));
1435            } else {
1436                vec.push(Number::from(Complex::new(prec), None));
1437            }
1438        }
1439        mat.push(vec);
1440    }
1441    mat
1442}
1443pub fn determinant(a: &[Vec<Number>]) -> Result<Number, &'static str> {
1444    if !a.is_empty() && (0..a.len()).all(|j| a.len() == a[j].len()) {
1445        Ok(Number::from(
1446            match a.len() {
1447                1 => a[0][0].number.clone(),
1448                2 => {
1449                    a[0][0].number.clone() * a[1][1].number.clone()
1450                        - a[1][0].number.clone() * a[0][1].number.clone()
1451                }
1452                3 => {
1453                    a[0][0].number.clone()
1454                        * (a[1][1].number.clone() * a[2][2].number.clone()
1455                            - a[1][2].number.clone() * a[2][1].number.clone())
1456                        + a[0][1].number.clone()
1457                            * (a[1][2].number.clone() * a[2][0].number.clone()
1458                                - a[1][0].number.clone() * a[2][2].number.clone())
1459                        + a[0][2].number.clone()
1460                            * (a[1][0].number.clone() * a[2][1].number.clone()
1461                                - a[1][1].number.clone() * a[2][0].number.clone())
1462                }
1463                _ => {
1464                    let mut det = Complex::new(a[0][0].number.prec());
1465                    for (i, x) in a[0].iter().enumerate() {
1466                        let mut sub_matrix = a[1..].to_vec();
1467                        for row in &mut sub_matrix {
1468                            row.remove(i);
1469                        }
1470                        det += x.number.clone()
1471                            * determinant(&sub_matrix)?.number
1472                            * if i % 2 == 0 { 1.0 } else { -1.0 };
1473                    }
1474                    det
1475                }
1476            },
1477            a[0][0].units.map(|b| b.pow(a.len() as f64)),
1478        ))
1479    } else {
1480        Err("not square")
1481    }
1482}
1483pub fn transpose(a: &[Vec<Number>]) -> Vec<Vec<Number>> {
1484    let mut max = 0;
1485    for i in a {
1486        if i.len() > max {
1487            max = i.len()
1488        }
1489    }
1490    let mut b = vec![vec![Number::from(Complex::new(1), None); a.len()]; max];
1491    for (i, l) in a.iter().enumerate() {
1492        for (j, n) in l.iter().enumerate() {
1493            b[j][i].clone_from(n);
1494        }
1495    }
1496    b
1497}
1498pub fn minors(a: &[Vec<Number>]) -> Result<Vec<Vec<Number>>, &'static str> {
1499    if a.iter().all(|j| a.len() == j.len()) {
1500        let mut result = vec![vec![Number::from(Complex::new(1), None); a[0].len()]; a.len()];
1501        for (i, k) in result.iter_mut().enumerate() {
1502            for (j, l) in k.iter_mut().enumerate() {
1503                *l = determinant(&submatrix(a, i, j))?
1504            }
1505        }
1506        Ok(result)
1507    } else {
1508        Err("not square")
1509    }
1510}
1511pub fn cofactor(a: &[Vec<Number>]) -> Result<Vec<Vec<Number>>, &'static str> {
1512    if a.iter().all(|j| a.len() == j.len()) && !a.is_empty() {
1513        let mut result = vec![vec![Number::from(Complex::new(1), None); a[0].len()]; a.len()];
1514        for (i, k) in result.iter_mut().enumerate() {
1515            for (j, l) in k.iter_mut().enumerate() {
1516                *l = if (i + j) % 2 == 1 {
1517                    let a = determinant(&submatrix(a, i, j))?;
1518                    Number::from(-a.number, a.units)
1519                } else {
1520                    determinant(&submatrix(a, i, j))?
1521                };
1522            }
1523        }
1524        Ok(result)
1525    } else {
1526        Err("not square")
1527    }
1528}
1529pub fn inverse(a: &[Vec<Number>]) -> Result<Vec<Vec<Number>>, &'static str> {
1530    if (0..a.len()).all(|j| a.len() == a[j].len()) {
1531        Matrix(transpose(&cofactor(a)?))
1532            .func(&NumStr::new(determinant(a)?), div)?
1533            .mat()
1534    } else {
1535        Err("not square")
1536    }
1537}
1538pub fn nth_prime(n: Integer) -> Integer {
1539    let mut count = Integer::from(1);
1540    let mut num = Integer::from(3);
1541    if n == 0 {
1542        num = Integer::new()
1543    } else if n == 1 {
1544        num = Integer::from(2);
1545    }
1546    while count < n {
1547        if num.is_probably_prime(100) != IsPrime::No {
1548            count += 1;
1549        }
1550        if count < n {
1551            num += 2;
1552        }
1553    }
1554    num
1555}
1556pub fn prime_factors(mut n: Integer) -> Vec<(Integer, isize)> {
1557    if n < 2 {
1558        return Vec::new();
1559    }
1560    let mut mat = Vec::new();
1561    let mut prime = Integer::from(2);
1562    let mut k = 0;
1563    loop {
1564        let (temp, rem) = n.clone().div_rem(prime.clone());
1565        if rem == 0 {
1566            k += 1;
1567            if temp == 1 {
1568                mat.push((prime.clone(), k));
1569                break;
1570            }
1571            n = temp;
1572        } else {
1573            if k != 0 {
1574                mat.push((prime.clone(), k));
1575            }
1576            prime = prime.next_prime();
1577            k = 0;
1578        }
1579    }
1580    mat
1581}
1582pub fn sort(mut a: Vec<Number>) -> Vec<Number> {
1583    a.sort_by(|x, y| {
1584        x.number
1585            .real()
1586            .partial_cmp(y.number.real())
1587            .unwrap_or(Ordering::Equal)
1588    });
1589    a
1590}
1591pub fn sort_mat(mut a: Vec<Vec<Number>>, prec: u32) -> Vec<Vec<Number>> {
1592    a.sort_by(|x, y| {
1593        if x.is_empty() || y.is_empty() {
1594            Ordering::Equal
1595        } else {
1596            x.iter()
1597                .fold(Float::new(prec), |sum, val| sum + val.number.real())
1598                .partial_cmp(
1599                    &y.iter()
1600                        .fold(Float::new(prec), |sum, val| sum + val.number.real()),
1601                )
1602                .unwrap_or(Ordering::Equal)
1603        }
1604    });
1605    a.sort_by(|x, y| x.len().partial_cmp(&y.len()).unwrap_or(Ordering::Equal));
1606    a
1607}
1608pub fn eigenvalues(mat: &[Vec<Number>], real: bool) -> Result<NumStr, &'static str> {
1609    if !mat.is_empty() && (0..mat.len()).all(|j| mat.len() == mat[j].len()) {
1610        match mat.len() {
1611            1 => Ok(NumStr::new(mat[0][0].clone())),
1612            2 => {
1613                let pr = mat[0][0].number.prec().0;
1614                let mut mat = Matrix(mat.into());
1615                mat.set_prec(pr * 2);
1616                let mat = mat.mat()?;
1617                let mut v = Vector(quadratic(
1618                    Number::from(Complex::with_val(mat[0][0].number.prec(), 1), None),
1619                    Number::from(
1620                        -mat[0][0].number.clone() - mat[1][1].number.clone(),
1621                        mat[0][0].units,
1622                    ),
1623                    Number::from(
1624                        mat[0][0].number.clone() * mat[1][1].number.clone()
1625                            - mat[0][1].number.clone() * mat[1][0].number.clone(),
1626                        mat[0][0].units.map(|a| a.pow(2.0)),
1627                    ),
1628                    real,
1629                ));
1630                v.set_prec(pr);
1631                Ok(v)
1632            }
1633            3 => {
1634                let pr = mat[0][0].number.prec().0;
1635                let mut mat = Matrix(mat.into());
1636                mat.set_prec(pr * 2);
1637                let mat = mat.mat()?;
1638                let mut v = Vector(cubic(
1639                    Number::from(Complex::with_val(mat[0][0].number.prec(), -1), None),
1640                    Number::from(
1641                        mat[2][2].number.clone()
1642                            + mat[1][1].number.clone()
1643                            + mat[0][0].number.clone(),
1644                        mat[0][0].units,
1645                    ),
1646                    Number::from(
1647                        -mat[0][0].number.clone() * mat[1][1].number.clone()
1648                            - mat[0][0].number.clone() * mat[2][2].number.clone()
1649                            + mat[0][1].number.clone() * mat[1][0].number.clone()
1650                            + mat[0][2].number.clone() * mat[2][0].number.clone()
1651                            - mat[1][1].number.clone() * mat[2][2].number.clone()
1652                            + mat[1][2].number.clone() * mat[2][1].number.clone(),
1653                        mat[0][0].units.map(|a| a.pow(2.0)),
1654                    ),
1655                    Number::from(
1656                        mat[0][0].number.clone()
1657                            * mat[1][1].number.clone()
1658                            * mat[2][2].number.clone()
1659                            - mat[0][0].number.clone()
1660                                * mat[1][2].number.clone()
1661                                * mat[2][1].number.clone()
1662                            - mat[0][1].number.clone()
1663                                * mat[1][0].number.clone()
1664                                * mat[2][2].number.clone()
1665                            + mat[0][1].number.clone()
1666                                * mat[1][2].number.clone()
1667                                * mat[2][0].number.clone()
1668                            + mat[0][2].number.clone()
1669                                * mat[1][0].number.clone()
1670                                * mat[2][1].number.clone()
1671                            - mat[0][2].number.clone()
1672                                * mat[1][1].number.clone()
1673                                * mat[2][0].number.clone(),
1674                        mat[0][0].units.map(|a| a.pow(3.0)),
1675                    ),
1676                    real,
1677                ));
1678                v.set_prec(pr);
1679                Ok(v)
1680            }
1681            4 => {
1682                let pr = mat[0][0].number.prec().0;
1683                let mut mat = Matrix(mat.into());
1684                mat.set_prec(pr * 2);
1685                let mat = mat.mat()?;
1686                let a = mat[0][0].number.clone();
1687                let b = mat[0][1].number.clone();
1688                let c = mat[0][2].number.clone();
1689                let d = mat[0][3].number.clone();
1690                let e = mat[1][0].number.clone();
1691                let f = mat[1][1].number.clone();
1692                let g = mat[1][2].number.clone();
1693                let h = mat[1][3].number.clone();
1694                let i = mat[2][0].number.clone();
1695                let j = mat[2][1].number.clone();
1696                let k = mat[2][2].number.clone();
1697                let l = mat[2][3].number.clone();
1698                let m = mat[3][0].number.clone();
1699                let n = mat[3][1].number.clone();
1700                let p = mat[3][2].number.clone();
1701                let q = mat[3][3].number.clone();
1702                let mut v = Vector(quartic(
1703                    Number::from(Complex::with_val(a.prec(), 1), None),
1704                    Number::from(
1705                        -a.clone() - f.clone() - k.clone() - q.clone(),
1706                        mat[0][0].units,
1707                    ),
1708                    Number::from(
1709                        a.clone() * f.clone() + a.clone() * k.clone() + a.clone() * q.clone()
1710                            - b.clone() * e.clone()
1711                            - c.clone() * i.clone()
1712                            - d.clone() * m.clone()
1713                            + f.clone() * k.clone()
1714                            + f.clone() * q.clone()
1715                            - g.clone() * j.clone()
1716                            - h.clone() * n.clone()
1717                            + k.clone() * q.clone()
1718                            - l.clone() * p.clone(),
1719                        mat[0][0].units.map(|a| a.pow(2.0)),
1720                    ),
1721                    Number::from(
1722                        -a.clone() * f.clone() * k.clone() - a.clone() * f.clone() * q.clone()
1723                            + a.clone() * g.clone() * j.clone()
1724                            + a.clone() * h.clone() * n.clone()
1725                            - a.clone() * k.clone() * q.clone()
1726                            + a.clone() * l.clone() * p.clone()
1727                            + b.clone() * e.clone() * k.clone()
1728                            + b.clone() * e.clone() * q.clone()
1729                            - b.clone() * g.clone() * i.clone()
1730                            - b.clone() * h.clone() * m.clone()
1731                            - c.clone() * e.clone() * j.clone()
1732                            + c.clone() * f.clone() * i.clone()
1733                            + c.clone() * i.clone() * q.clone()
1734                            - c.clone() * l.clone() * m.clone()
1735                            - d.clone() * e.clone() * n.clone()
1736                            + d.clone() * f.clone() * m.clone()
1737                            - d.clone() * i.clone() * p.clone()
1738                            + d.clone() * k.clone() * m.clone()
1739                            - f.clone() * k.clone() * q.clone()
1740                            + f.clone() * l.clone() * p.clone()
1741                            + g.clone() * j.clone() * q.clone()
1742                            - g.clone() * l.clone() * n.clone()
1743                            - h.clone() * j.clone() * p.clone()
1744                            + h.clone() * k.clone() * n.clone(),
1745                        mat[0][0].units.map(|a| a.pow(3.0)),
1746                    ),
1747                    Number::from(
1748                        a.clone() * f.clone() * k.clone() * q.clone()
1749                            - a.clone() * f.clone() * l.clone() * p.clone()
1750                            - a.clone() * g.clone() * j.clone() * q.clone()
1751                            + a.clone() * g.clone() * l.clone() * n.clone()
1752                            + a.clone() * h.clone() * j.clone() * p.clone()
1753                            - a.clone() * h.clone() * k.clone() * n.clone()
1754                            - b.clone() * e.clone() * k.clone() * q.clone()
1755                            + b.clone() * e.clone() * l.clone() * p.clone()
1756                            + b.clone() * g.clone() * i.clone() * q.clone()
1757                            - b.clone() * g.clone() * l.clone() * m.clone()
1758                            - b.clone() * h.clone() * i.clone() * p.clone()
1759                            + b.clone() * h.clone() * k.clone() * m.clone()
1760                            + c.clone() * e.clone() * j.clone() * q.clone()
1761                            - c.clone() * e.clone() * l.clone() * n.clone()
1762                            - c.clone() * f.clone() * i.clone() * q.clone()
1763                            + c.clone() * f.clone() * l.clone() * m.clone()
1764                            + c.clone() * h.clone() * i.clone() * n.clone()
1765                            - c.clone() * h.clone() * j.clone() * m.clone()
1766                            - d.clone() * e.clone() * j.clone() * p.clone()
1767                            + d.clone() * e.clone() * k.clone() * n.clone()
1768                            + d.clone() * f.clone() * i.clone() * p.clone()
1769                            - d.clone() * f.clone() * k.clone() * m.clone()
1770                            - d.clone() * g.clone() * i.clone() * n.clone()
1771                            + d.clone() * g.clone() * j.clone() * m.clone(),
1772                        mat[0][0].units.map(|a| a.pow(4.0)),
1773                    ),
1774                    real,
1775                ));
1776                v.set_prec(pr);
1777                Ok(v)
1778            }
1779            _ => Err("unsupported"),
1780        }
1781    } else {
1782        Err("not square")
1783    }
1784}
1785pub fn eigenvectors(mat: &[Vec<Number>], real: bool) -> Result<NumStr, &'static str> {
1786    if !mat.is_empty() && (0..mat.len()).all(|j| mat.len() == mat[j].len()) {
1787        let one = Number::from(Complex::with_val(mat[0][0].number.prec(), 1), None);
1788        match mat.len() {
1789            1 => Ok(NumStr::new(one)),
1790            2..5 => {
1791                let p = mat[0][0].number.prec().0;
1792                let mut l = eigenvalues(mat, real)?.vec()?;
1793                let mut i = 0;
1794                while i + 1 < l.len() {
1795                    let mut has = false;
1796                    for v in l[i + 1..].iter().cloned() {
1797                        if (v.number - l[i].number.clone()).abs().real().clone().log2()
1798                            < -(p as i32 / 8)
1799                        {
1800                            l.remove(i);
1801                            has = true;
1802                            break;
1803                        }
1804                    }
1805                    if !has {
1806                        i += 1;
1807                    }
1808                }
1809                let v = l
1810                    .iter()
1811                    .filter_map(|l| {
1812                        Matrix(identity(mat.len(), l.number.prec().0))
1813                            .mul(&NumStr::new(l.clone()))
1814                            .map(|n| {
1815                                Matrix(mat.to_vec())
1816                                    .func(&n, sub)
1817                                    .map(|m| Some(kernel(m.mat().unwrap()).unwrap()))
1818                                    .unwrap_or(None)
1819                            })
1820                            .unwrap_or(None)
1821                    })
1822                    .flatten()
1823                    .collect::<Vec<Vec<Number>>>();
1824                Ok(Matrix(v))
1825            }
1826            _ => Err("unsupported"),
1827        }
1828    } else {
1829        Err("not square")
1830    }
1831}
1832pub fn rcf(mat: Vec<Vec<Number>>) -> Result<NumStr, &'static str> {
1833    if mat.is_empty() || (0..mat.len()).any(|j| mat.len() != mat[j].len()) {
1834        return Err("matrix not square");
1835    }
1836    let pr = mat[0][0].number.prec().0;
1837    let l = mat.len();
1838    let mut ev = gen_ev(&mat, false)?;
1839    let mut beta = Vec::new();
1840    let m = Matrix(mat.clone());
1841    while !ev.is_empty() {
1842        let v = ev
1843            .iter_mut()
1844            .filter_map(|e| e.pop().map(Vector))
1845            .collect::<Vec<NumStr>>();
1846        let l = v.len();
1847        if l == 0 {
1848            break;
1849        }
1850        let v = if l == 1 {
1851            v[0].clone()
1852        } else {
1853            v[1..].iter().fold(v[0].clone(), |sum, val| {
1854                sum.func(val, add).unwrap_or(Func(String::new()))
1855            })
1856        };
1857        let mut c = v.vec()?;
1858        beta.push(c.clone());
1859        for _ in 1..l {
1860            c = m.mul(&Vector(c))?.vec()?;
1861            beta.push(c.clone())
1862        }
1863    }
1864    change_basis(mat, &identity(l, pr), &transpose(&beta))
1865}
1866pub fn jcf(mat: Vec<Vec<Number>>) -> Result<NumStr, &'static str> {
1867    if mat.is_empty() || (0..mat.len()).any(|j| mat.len() != mat[j].len()) {
1868        return Err("matrix not square");
1869    }
1870    let pr = mat[0][0].number.prec().0;
1871    let l = mat.len();
1872    let beta = transpose(&generalized_eigenvectors(&mat, false)?.mat()?);
1873    let i = identity(l, pr);
1874    let mut o = i.clone();
1875    let mut d = change_basis(mat, &i, &beta)?.mat()?;
1876    for (i, r) in o.iter_mut().enumerate() {
1877        r[i].number = d[i..d.len() - 1]
1878            .iter()
1879            .enumerate()
1880            .map_while(|(j, r)| {
1881                if -r[i + j + 1].number.clone().abs().real().clone().ln() < pr / 16 {
1882                    Some(r[i + j + 1].number.clone())
1883                } else {
1884                    None
1885                }
1886            })
1887            .fold(Complex::with_val(pr, 1), |sum, val| sum * val.clone());
1888    }
1889    let l = o.len();
1890    for k in 1..l {
1891        d = change_basis(d, &i, &o)?.mat()?;
1892        o = i.clone();
1893        let mut sum = Complex::new(pr);
1894        for (i, r) in o[1..l - k].iter_mut().enumerate() {
1895            let i = i + 1;
1896            sum += d[i - 1][i + k].number.clone();
1897            r[i + k].number = -sum.clone();
1898        }
1899    }
1900    let d = change_basis(d, &i, &o)?;
1901    Ok(d)
1902}
1903fn gen_ev(mat: &[Vec<Number>], real: bool) -> Result<Vec<Vec<Vec<Number>>>, &'static str> {
1904    if mat.is_empty() || (0..mat.len()).any(|j| mat.len() != mat[j].len()) {
1905        return Err("matrix not square");
1906    }
1907    let p = mat[0][0].number.prec().0;
1908    let mut l = eigenvalues(mat, real)?.vec()?;
1909    let mut i = 0;
1910    while i + 1 < l.len() {
1911        let mut has = false;
1912        for v in l[i + 1..].iter().cloned() {
1913            if (v.number - l[i].number.clone()).abs().real().clone().log2() < -(p as i32 / 8) {
1914                l.remove(i);
1915                has = true;
1916                break;
1917            }
1918        }
1919        if !has {
1920            i += 1;
1921        }
1922    }
1923    Ok(l.iter()
1924        .filter_map(|l| {
1925            Matrix(identity(mat.len(), l.number.prec().0))
1926                .mul(&NumStr::new(l.clone()))
1927                .map(|n| {
1928                    Matrix(mat.to_vec())
1929                        .func(&n, sub)
1930                        .map(|m| {
1931                            m.pow(&NumStr::new(Number::from(
1932                                Complex::with_val(p, mat.len()),
1933                                None,
1934                            )))
1935                            .map(|m| {
1936                                let ker = kernel(m.mat().unwrap()).unwrap();
1937                                Some(ker)
1938                            })
1939                            .unwrap_or(None)
1940                        })
1941                        .unwrap_or(None)
1942                })
1943                .unwrap_or(None)
1944        })
1945        .collect::<Vec<Vec<Vec<Number>>>>())
1946}
1947pub fn generalized_eigenvectors(mat: &[Vec<Number>], real: bool) -> Result<NumStr, &'static str> {
1948    if !mat.is_empty() && (0..mat.len()).all(|j| mat.len() == mat[j].len()) {
1949        let one = Number::from(Complex::with_val(mat[0][0].number.prec(), 1), None);
1950        match mat.len() {
1951            1 => Ok(NumStr::new(one)),
1952            2..5 => Ok(Matrix(
1953                gen_ev(mat, real)?
1954                    .iter()
1955                    .flatten()
1956                    .cloned()
1957                    .collect::<Vec<Vec<Number>>>(),
1958            )),
1959            _ => Err("unsupported"),
1960        }
1961    } else {
1962        Err("not square")
1963    }
1964}
1965pub fn change_basis(
1966    a: Vec<Vec<Number>>,
1967    beta: &[Vec<Number>],
1968    gamma: &[Vec<Number>],
1969) -> Result<NumStr, &'static str> {
1970    let m = Matrix(a);
1971    let tn = Matrix(inverse(&transpose(gamma))?);
1972    let mut c = Vec::new();
1973    for b in beta {
1974        c.push(tn.mul(&Vector(b.to_vec()))?.vec()?)
1975    }
1976    let tn = Matrix(inverse(&transpose(beta))?);
1977    let mut d = Vec::new();
1978    for g in gamma {
1979        d.push(tn.mul(&Vector(g.to_vec()))?.vec()?)
1980    }
1981    Matrix(c).mul(&m)?.mul(&Matrix(d))
1982}
1983pub fn coordinate(v: Vec<Number>, beta: Vec<Vec<Number>>) -> Result<NumStr, &'static str> {
1984    let tn = Matrix(inverse(&transpose(&beta))?);
1985    tn.mul(&Vector(v))
1986}
1987pub fn rref(mut a: Vec<Vec<Number>>) -> Result<Vec<Vec<Number>>, &'static str> {
1988    if a.is_empty() || a[0].is_empty() || a.iter().any(|b| a[0].len() != b.len()) {
1989        return Err("invalid matrix");
1990    }
1991    let mut count = 0;
1992    for i in 0..a[0].len() {
1993        if let Some((n, v)) = a
1994            .clone()
1995            .iter()
1996            .enumerate()
1997            .find(|(j, b)| *j >= count && !b[i].number.is_zero())
1998        {
1999            for (a, r) in a.iter_mut().enumerate() {
2000                let c = r[i].number.clone() / v[i].number.clone();
2001                for (b, num) in r.iter_mut().enumerate() {
2002                    if a != n && !v[i].number.is_zero() {
2003                        num.number -= v[b].number.clone() * c.clone();
2004                    } else {
2005                        num.number /= v[i].number.clone()
2006                    }
2007                }
2008            }
2009            if n != count {
2010                a.swap(count, n);
2011            }
2012            count += 1;
2013        }
2014    }
2015    Ok(a.to_vec())
2016}
2017pub fn kernel(a: Vec<Vec<Number>>) -> Result<Vec<Vec<Number>>, &'static str> {
2018    if a.is_empty() || a[0].is_empty() || a.iter().any(|b| a[0].len() != b.len()) {
2019        return Err("invalid matrix");
2020    }
2021    let pr = a[0][0].number.prec().0;
2022    let rref = rref(a.clone())?;
2023    let mut ker = Vec::new();
2024    let mut leading_ones = Vec::new();
2025    for r in &rref {
2026        if let Some((pos, _)) = r.iter().enumerate().find(|(_, n)| !n.number.is_zero()) {
2027            leading_ones.push(pos);
2028        }
2029    }
2030    let t = transpose(&rref);
2031    let m = Matrix(a);
2032    for (i, t) in t.iter().enumerate() {
2033        if !leading_ones.contains(&i) {
2034            let mut zero = vec![Number::from(Complex::new(pr), None); rref[0].len()];
2035            for j in 0..i.min(leading_ones.len()) {
2036                if leading_ones[j] < i {
2037                    zero[leading_ones[j]] = Number::from(-1.0 * t[j].number.clone(), None)
2038                }
2039            }
2040            zero[i] = Number::from(Complex::with_val(pr, 1), None);
2041            if m.mul(&Vector(zero.clone()))?
2042                .vec()?
2043                .iter()
2044                .all(|n| -n.number.clone().abs().real().clone().log2() > pr / 16)
2045            {
2046                ker.push(zero);
2047            }
2048        }
2049    }
2050    Ok(ker)
2051}
2052pub fn range(a: Vec<Vec<Number>>) -> Result<Vec<Vec<Number>>, &'static str> {
2053    let rref = rref(a.clone())?;
2054    let mut ran = Vec::new();
2055    let mut leading_ones = Vec::new();
2056    for r in &rref {
2057        if let Some((pos, _)) = r.iter().enumerate().find(|(_, n)| !n.number.is_zero()) {
2058            leading_ones.push(pos);
2059        }
2060    }
2061    let t = transpose(&a);
2062    for i in leading_ones {
2063        ran.push(t[i].clone());
2064    }
2065    Ok(ran)
2066}
2067pub fn mul_units(a: Option<Units>, b: Option<Units>) -> Option<Units> {
2068    match (a, b) {
2069        (Some(a), None) => Some(a),
2070        (None, Some(b)) => Some(b),
2071        (Some(a), Some(b)) => Some(a.mul(&b)),
2072        _ => None,
2073    }
2074}
2075pub fn div_units(a: Option<Units>, b: Option<Units>) -> Option<Units> {
2076    match (a, b) {
2077        (Some(a), None) => Some(a),
2078        (None, Some(b)) => Some(b.pow(-1.0)),
2079        (Some(a), Some(b)) => Some(a.div(&b)),
2080        _ => None,
2081    }
2082}
2083pub fn quadratic(a: Number, b: Number, c: Number, real: bool) -> Vec<Number> {
2084    if a.number.is_zero() {
2085        return if b.number.is_zero() {
2086            vec![Number::from(Complex::new(a.number.prec()), None)]
2087        } else {
2088            let units = div_units(c.units, b.units);
2089            let b = b.number;
2090            let c = c.number;
2091            let mut r = -c / b;
2092            if -r.imag().clone().abs().log10() > a.number.prec().0 / 4 {
2093                r = r.real().clone().into();
2094            }
2095            if real && !r.imag().is_zero() {
2096                vec![Number::from(Complex::with_val(a.number.prec(), Nan), None)]
2097            } else {
2098                vec![Number::from(r, units)]
2099            }
2100        };
2101    }
2102    let units = div_units(c.units, a.units).map(|a| a.root(2.0));
2103    let a = a.number;
2104    let b = b.number;
2105    let c = c.number;
2106    let p: Complex = sqr(b.clone());
2107    let p: Complex = p - (4 * c * a.clone());
2108    let p = p.sqrt();
2109    let a: Complex = 2 * a;
2110    let mut z1 = (p.clone() - b.clone()) / a.clone();
2111    let mut z2 = (-p - b) / a.clone();
2112    if -z1.imag().clone().abs().log10() > a.prec().0 / 4 {
2113        z1 = z1.real().clone().into();
2114    }
2115    if -z2.imag().clone().abs().log10() > a.prec().0 / 4 {
2116        z2 = z2.real().clone().into();
2117    }
2118    if real {
2119        let mut vec = Vec::new();
2120        if z1.imag().is_zero() {
2121            vec.push(Number::from(z1, units))
2122        }
2123        if z2.imag().is_zero() {
2124            vec.push(Number::from(z2, units))
2125        }
2126        if vec.is_empty() {
2127            vec![Number::from(Complex::with_val(a.prec(), Nan), None)]
2128        } else {
2129            vec
2130        }
2131    } else {
2132        vec![Number::from(z1, units), Number::from(z2, units)]
2133    }
2134}
2135pub fn cubic(a: Number, b: Number, c: Number, d: Number, real: bool) -> Vec<Number> {
2136    let units = div_units(d.units, a.units).map(|a| a.root(3.0));
2137    if a.number.is_zero() {
2138        return quadratic(b, c, d, real);
2139    }
2140    let d = d.number;
2141    if d.is_zero() {
2142        let mut vec = quadratic(a, b, c, real);
2143        vec.push(Number::from(Complex::new(d.prec()), vec[0].units));
2144        return vec;
2145    }
2146    let a = a.number;
2147    let b = b.number;
2148    let c = c.number;
2149    let prec = a.prec();
2150    let threerecip = Float::with_val(prec.0, 3).recip();
2151    if b.is_zero() && c.is_zero() {
2152        let reuse = pow_nth(d / a.clone(), threerecip.clone().into());
2153        let mut z1 = -reuse.clone();
2154        let mut z2 = reuse.clone() * Complex::with_val(prec.0, -1).pow(threerecip.clone());
2155        let mut z3: Complex = -reuse * Complex::with_val(prec.0, -1).pow(2 * threerecip);
2156        if -z1.imag().clone().abs().log10() > a.prec().0 / 4 {
2157            z1 = z1.real().clone().into();
2158        }
2159        if -z2.imag().clone().abs().log10() > a.prec().0 / 4 {
2160            z2 = z2.real().clone().into();
2161        }
2162        if -z3.imag().clone().abs().log10() > a.prec().0 / 4 {
2163            z3 = z3.real().clone().into();
2164        }
2165        return if real {
2166            let mut vec = Vec::new();
2167            if z1.imag().is_zero() {
2168                vec.push(Number::from(z1, units))
2169            }
2170            if z2.imag().is_zero() {
2171                vec.push(Number::from(z2, units))
2172            }
2173            if z3.imag().is_zero() {
2174                vec.push(Number::from(z3, units))
2175            }
2176            if vec.is_empty() {
2177                vec![Number::from(Complex::with_val(prec, Nan), None)]
2178            } else {
2179                vec
2180            }
2181        } else {
2182            vec![
2183                Number::from(z1, units),
2184                Number::from(z2, units),
2185                Number::from(z3, units),
2186            ]
2187        };
2188    }
2189    let b = b / a.clone();
2190    let c = c / a.clone();
2191    let d = d / a.clone();
2192    // https://en.wikipedia.org/wiki/Cubic_equation#General_cubic_formula
2193    let d0: Complex = sqr(b.clone()) - 3 * c.clone();
2194    let d1: Complex = 2 * cube(b.clone()) - 9 * b.clone() * c.clone() + 27 * d.clone();
2195    let c: Complex = sqr(d1.clone()) - 4 * cube(d0.clone());
2196    let c: Complex = (d1 + c.sqrt()) / 2;
2197    let c = pow_nth(c, threerecip.clone().into());
2198    let omega: Complex = Complex::with_val(prec, (-0.5, Float::with_val(prec.0, 3).sqrt() / 2));
2199    let mut z1: Complex = if d0.is_zero() {
2200        -(b.clone() + c.clone()) / 3
2201    } else {
2202        -(b.clone() + c.clone() + d0.clone() / c.clone()) / 3
2203    };
2204    let c0 = c.clone() * omega.clone();
2205    let mut z2: Complex = if d0.is_zero() {
2206        -(b.clone() + c0.clone()) / 3
2207    } else {
2208        -(b.clone() + c0.clone() + d0.clone() / c0) / 3
2209    };
2210    let c1 = c * omega.conj();
2211    let mut z3: Complex = if d0.is_zero() {
2212        -(b + c1.clone()) / 3
2213    } else {
2214        -(b + c1.clone() + d0 / c1) / 3
2215    };
2216    if -z1.imag().clone().abs().log10() > a.prec().0 / 4 {
2217        z1 = z1.real().clone().into();
2218    }
2219    if -z2.imag().clone().abs().log10() > a.prec().0 / 4 {
2220        z2 = z2.real().clone().into();
2221    }
2222    if -z3.imag().clone().abs().log10() > a.prec().0 / 4 {
2223        z3 = z3.real().clone().into();
2224    }
2225    if real {
2226        let mut vec = Vec::new();
2227        if z1.imag().is_zero() {
2228            vec.push(Number::from(z1, units))
2229        }
2230        if z2.imag().is_zero() {
2231            vec.push(Number::from(z2, units))
2232        }
2233        if z3.imag().is_zero() {
2234            vec.push(Number::from(z3, units))
2235        }
2236        if vec.is_empty() {
2237            vec![Number::from(Complex::with_val(prec, Nan), None)]
2238        } else {
2239            vec
2240        }
2241    } else {
2242        vec![
2243            Number::from(z1, units),
2244            Number::from(z2, units),
2245            Number::from(z3, units),
2246        ]
2247    }
2248}
2249pub fn quartic(div: Number, b: Number, c: Number, d: Number, e: Number, real: bool) -> Vec<Number> {
2250    let units = div_units(e.units, div.units).map(|a| a.root(4.0));
2251    if e.number.is_zero() {
2252        let mut vec = cubic(div, b, c, d, real);
2253        vec.push(Number::from(Complex::new(e.number.prec()), vec[0].units));
2254        return vec;
2255    }
2256    let div = div.number;
2257    if div.is_zero() {
2258        return cubic(b, c, d, e, real);
2259    }
2260    let b = b.number;
2261    let c = c.number;
2262    let d = d.number;
2263    let e = e.number;
2264    let prec = div.prec();
2265    let threerecip = Float::with_val(prec.0, 3).recip();
2266    let a = b / div.clone();
2267    let b = c / div.clone();
2268    let c = d / div.clone();
2269    let d = e / div;
2270    if a.is_zero() && b.is_zero() && c.is_zero() {
2271        return if d.is_zero() {
2272            vec![Number::from(Complex::new(prec), None)]
2273        } else {
2274            let a = Complex::with_val(prec, 1) / 4;
2275            let mut a: Complex = (-d).pow(a);
2276            let mut v = vec![a.clone()];
2277            for _ in 0..3 {
2278                a = a.mul_i(false);
2279                v.push(a.clone());
2280            }
2281            v.into_iter().map(|a| Number::from(a, None)).collect()
2282        };
2283    }
2284    // https://upload.wikimedia.org/wikipedia/commons/9/99/Quartic_Formula.svg
2285    let alpha: Complex = sqr(b.clone()) - 3 * a.clone() * c.clone() + 12 * d.clone();
2286    let phi: Complex = 2 * cube(b.clone()) - 9 * a.clone() * b.clone() * c.clone()
2287        + 27 * sqr(c.clone())
2288        + 27 * sqr(a.clone()) * d.clone()
2289        - 72 * b.clone() * d.clone();
2290    let omega: Complex = -4 * cube(alpha.clone()) + sqr(phi.clone());
2291    let omega: Complex = phi + omega.sqrt();
2292    let alpha: Complex = if alpha.is_zero() {
2293        Complex::new(prec)
2294    } else {
2295        Float::with_val(prec.0, 2).pow(threerecip.clone()) * alpha
2296            / (3 * pow_nth(omega.clone(), threerecip.clone().into()))
2297    };
2298    let beta: Complex = omega / 54;
2299    let beta: Complex = pow_nth(beta, threerecip.clone().into());
2300    let infirst: Complex = sqr(a.clone()) / 4 - 2 * b.clone() / 3 + alpha.clone() + beta.clone();
2301    let first: Complex = infirst.clone().sqrt() / 2;
2302    let third: Complex = -1 * cube(a.clone()) + 4 * a.clone() * b.clone() - 8 * c.clone();
2303    let third: Complex = if third.is_zero() {
2304        Complex::new(prec)
2305    } else {
2306        third / (first.clone() * 8)
2307    };
2308    let a4: Complex = -a.clone() / 4;
2309    let second: Complex = sqr(a.clone()) / 2 - 4 * b.clone() / 3 - alpha.clone() - beta.clone();
2310    let secondn: Complex = second.clone() - third.clone();
2311    let secondn: Complex = secondn.sqrt() / 2;
2312    let secondp: Complex = second + third.clone();
2313    let secondp: Complex = secondp.sqrt() / 2;
2314    let mut r1 = a4.clone() - first.clone() - secondn.clone();
2315    let mut r2 = a4.clone() - first.clone() + secondn.clone();
2316    let mut r3 = a4.clone() + first.clone() - secondp.clone();
2317    let mut r4 = a4.clone() + first.clone() + secondp.clone();
2318    if -r1.imag().clone().abs().log10() > a.prec().0 / 8 {
2319        r1 = r1.real().clone().into();
2320    }
2321    if -r2.imag().clone().abs().log10() > a.prec().0 / 8 {
2322        r2 = r2.real().clone().into();
2323    }
2324    if -r3.imag().clone().abs().log10() > a.prec().0 / 8 {
2325        r3 = r3.real().clone().into();
2326    }
2327    if -r4.imag().clone().abs().log10() > a.prec().0 / 8 {
2328        r4 = r4.real().clone().into();
2329    }
2330    if real {
2331        let mut vec = Vec::new();
2332        if r1.imag().is_zero() {
2333            vec.push(Number::from(r1, units))
2334        }
2335        if r2.imag().is_zero() {
2336            vec.push(Number::from(r2, units))
2337        }
2338        if r3.imag().is_zero() {
2339            vec.push(Number::from(r3, units))
2340        }
2341        if r4.imag().is_zero() {
2342            vec.push(Number::from(r4, units))
2343        }
2344        if vec.is_empty() {
2345            vec![Number::from(Complex::with_val(prec, Nan), None)]
2346        } else {
2347            vec
2348        }
2349    } else {
2350        vec![
2351            Number::from(r1, units),
2352            Number::from(r2, units),
2353            Number::from(r3, units),
2354            Number::from(r4, units),
2355        ]
2356    }
2357}
2358pub fn variance(a: &[Number], mean: Option<Complex>, prec: u32) -> Number {
2359    let mean = if let Some(n) = mean {
2360        n
2361    } else {
2362        a.iter()
2363            .fold(Complex::new(prec), |sum, val| sum + val.number.clone())
2364            / a.len()
2365    };
2366    let mut variance = Complex::new(prec);
2367    for a in a {
2368        variance += sqr(a.number.clone() - mean.clone())
2369    }
2370    Number::from(
2371        variance / (a.len().saturating_sub(1)),
2372        a[0].units.map(|a| a.pow(2.0)),
2373    )
2374}
2375pub fn recursion(
2376    mut func_vars: Vec<(String, Vec<NumStr>)>,
2377    mut func: Vec<NumStr>,
2378    options: Options,
2379) -> Result<NumStr, &'static str> {
2380    for fv in func_vars.clone() {
2381        if fv.0.ends_with(')') {
2382            let mut cont = false;
2383            let mut bracket = 0;
2384            let mut pw = Vec::new();
2385            for f in &func {
2386                match f {
2387                    LeftBracket | LeftCurlyBracket => bracket += 1,
2388                    RightBracket | RightCurlyBracket => {
2389                        bracket -= 1;
2390                        if !pw.is_empty() && bracket == pw[0] {
2391                            pw.remove(0);
2392                        }
2393                    }
2394                    Func(s) if s == "pw" => pw.insert(0, bracket),
2395                    _ => {
2396                        if !pw.is_empty() && f.str_is(&fv.0) {
2397                            cont = true
2398                        }
2399                    }
2400                }
2401            }
2402            if cont {
2403                continue;
2404            }
2405            if fv.0.contains(',') {
2406                let mut vars = fv.0.split(',').collect::<Vec<&str>>();
2407                vars[0] = vars[0].split('(').next_back().unwrap();
2408                {
2409                    let vl = vars.len().saturating_sub(1);
2410                    vars[vl] = &vars[vl][0..vars[vl].len().saturating_sub(1)];
2411                }
2412                let mut x = func.len();
2413                while x > 0 {
2414                    x -= 1;
2415                    if func[x].str_is(&fv.0) {
2416                        let mut fv = fv.clone();
2417                        let mut bracket = 0;
2418                        let mut k = 0;
2419                        let mut processed = Vec::new();
2420                        let mut last = 0;
2421                        for (i, n) in func[x + 2..].iter().enumerate() {
2422                            match n {
2423                                LeftBracket | LeftCurlyBracket => bracket += 1,
2424                                RightBracket | RightCurlyBracket => {
2425                                    if bracket == 0 {
2426                                        if let Ok(n) = do_math(
2427                                            func[x + 2 + last..x + 2 + i].to_vec(),
2428                                            options,
2429                                            func_vars.clone(),
2430                                        ) {
2431                                            processed.push(vec![n]);
2432                                        } else {
2433                                            let iden = format!(
2434                                                "@{}{}@",
2435                                                func_vars.len(),
2436                                                vars[processed.len()]
2437                                            );
2438                                            func_vars.push((
2439                                                iden.clone(),
2440                                                func[x + 2 + last..x + 2 + i].to_vec(),
2441                                            ));
2442                                            processed.push(vec![Func(iden)]);
2443                                        }
2444                                        k = i;
2445                                        break;
2446                                    }
2447                                    bracket -= 1;
2448                                }
2449                                Comma if bracket == 0 => {
2450                                    if let Ok(n) = do_math(
2451                                        func[x + 2 + last..x + 2 + i].to_vec(),
2452                                        options,
2453                                        func_vars.clone(),
2454                                    ) {
2455                                        processed.push(vec![n]);
2456                                    } else {
2457                                        let iden = format!(
2458                                            "@{}{}@",
2459                                            func_vars.len(),
2460                                            vars[processed.len()]
2461                                        );
2462                                        func_vars.push((
2463                                            iden.clone(),
2464                                            func[x + 2 + last..x + 2 + i].to_vec(),
2465                                        ));
2466                                        processed.push(vec![Func(iden)]);
2467                                    }
2468                                    last = i + 1;
2469                                }
2470                                _ => {}
2471                            }
2472                        }
2473                        let mut i = 0;
2474                        while i < fv.1.len() {
2475                            if let Func(s) = &fv.1[i] {
2476                                for v in processed.iter().zip(vars.iter()) {
2477                                    if s == v.1 {
2478                                        fv.1.remove(i);
2479                                        fv.1.splice(i..i, v.0.clone());
2480                                        break;
2481                                    }
2482                                }
2483                            }
2484                            i += 1;
2485                        }
2486                        func.drain(x..=k + x + 2);
2487                        func.splice(x..x, fv.1.clone());
2488                    }
2489                }
2490            } else {
2491                let var = fv.0.split('(').next_back().unwrap();
2492                let var = &var[0..var.len().saturating_sub(1)];
2493                let mut x = func.len();
2494                while x > 0 {
2495                    x -= 1;
2496                    if func[x].str_is(&fv.0) {
2497                        let mut fv = fv.clone();
2498                        let mut bracket = 0;
2499                        let mut k = 0;
2500                        for (i, n) in func[x + 2..].iter().enumerate() {
2501                            match n {
2502                                LeftBracket | LeftCurlyBracket => bracket += 1,
2503                                RightBracket | RightCurlyBracket => {
2504                                    if bracket == 0 {
2505                                        k = i;
2506                                    }
2507                                    bracket -= 1;
2508                                }
2509                                _ => {}
2510                            }
2511                        }
2512                        let iden = format!("@{}{}@", func_vars.len(), var);
2513                        let mut i = 0;
2514                        while i < fv.1.len() {
2515                            if let Func(s) = &fv.1[i] {
2516                                if *s == var {
2517                                    fv.1[i] = Func(iden.clone());
2518                                }
2519                            }
2520                            i += 1;
2521                        }
2522                        func_vars.push((iden.clone(), func[x + 1..=k + x + 2].to_vec()));
2523                        func.drain(x..=k + x + 2);
2524                        func.splice(x..x, fv.1.clone());
2525                    }
2526                }
2527            }
2528        }
2529    }
2530    do_math(func, options, func_vars)
2531}
2532pub fn gcd(mut x: Integer, mut y: Integer) -> Integer {
2533    while x != y {
2534        if x > y {
2535            x -= y.clone()
2536        } else {
2537            y -= x.clone()
2538        }
2539    }
2540    x
2541}
2542pub fn incomplete_beta(x: Complex, a: Complex, b: Complex) -> Complex {
2543    if x.real() > &((a.real().clone() + 1) / (a.real() + b.real().clone() + 2)) {
2544        (gamma(a.clone()) * gamma(b.clone()) / gamma(a.clone() + b.clone()))
2545            - incomplete_beta(1 - x, b, a)
2546    } else {
2547        let f: Complex = 1 - x.clone();
2548        pow_nth(x.clone(), a.clone()) * pow_nth(f, b.clone())
2549            / (a.clone() * (1 + incomplete_beta_recursion(x, a, b, 1, 10)))
2550    }
2551}
2552fn incomplete_beta_recursion(
2553    x: Complex,
2554    a: Complex,
2555    b: Complex,
2556    iter: usize,
2557    max: usize,
2558) -> Complex {
2559    if iter == max {
2560        Complex::new(x.prec())
2561    } else if iter % 2 == 1 {
2562        let m = (iter - 1) / 2;
2563        (-x.clone() * (a.clone() + m) * (a.clone() + b.clone() + m)
2564            / ((a.clone() + (2 * m)) * (a.clone() + (2 * m) + 1)))
2565            / (1 + incomplete_beta_recursion(x, a, b, iter + 1, max))
2566    } else {
2567        let m = iter / 2;
2568        (x.clone() * m * (b.clone() - m) / ((a.clone() + (2 * m)) * (a.clone() + (2 * m) - 1)))
2569            / (1 + incomplete_beta_recursion(x, a, b, iter + 1, max))
2570    }
2571}
2572pub fn erf(z: Complex) -> Complex {
2573    1 - erfc(z)
2574}
2575pub fn erfc(z: Complex) -> Complex {
2576    let p2: Complex = sqr(z.clone());
2577    erfc_recursion(z.clone(), 0, z.prec().0 as usize / 4) / Complex::with_val(z.prec(), Pi).sqrt()
2578        * z
2579        / p2.exp()
2580}
2581fn erfc_recursion(z: Complex, iter: usize, max: usize) -> Complex {
2582    if iter == max {
2583        Complex::with_val(z.prec(), 1)
2584    } else if iter == 0 {
2585        erfc_recursion(z, 1, max).recip()
2586    } else if iter % 2 == 1 {
2587        sqr(z.clone()) + iter / (2 * erfc_recursion(z, iter + 1, max))
2588    } else {
2589        1 + iter / (2 * erfc_recursion(z, iter + 1, max))
2590    }
2591}
2592fn gamma0(z: Complex) -> Complex {
2593    let p = z.prec().0 as usize / 4;
2594    gamma0_recursion_first(z.clone(), 0, p) + gamma0_recursion_second(z, 0, p)
2595}
2596fn gamma0_recursion_first(z: Complex, iter: usize, max: usize) -> Complex {
2597    if iter == max {
2598        Complex::with_val(z.prec(), 1)
2599    } else if iter == 0 {
2600        Float::with_val(z.prec().0, -1).exp() / gamma0_recursion_first(z, 1, max)
2601    } else {
2602        2 * iter - z.clone() + iter * (z.clone() - iter) / gamma0_recursion_first(z, iter + 1, max)
2603    }
2604}
2605fn gamma0_recursion_second(z: Complex, iter: usize, max: usize) -> Complex {
2606    if iter == max {
2607        Complex::with_val(z.prec(), 1)
2608    } else if iter == 0 {
2609        Float::with_val(z.prec().0, -1).exp() / gamma0_recursion_second(z, 1, max)
2610    } else if iter % 2 == 1 {
2611        (iter - 1) + z.clone() - (z.clone() + iter / 2) / gamma0_recursion_second(z, iter + 1, max)
2612    } else {
2613        (iter - 1) + z.clone() + (iter / 2) / gamma0_recursion_second(z, iter + 1, max)
2614    }
2615}
2616pub fn incomplete_gamma(s: Complex, z: Complex) -> Complex {
2617    if z.is_zero() {
2618        gamma(s)
2619    } else if s.real().is_sign_positive()
2620        && !s.real().clone().fract().is_zero()
2621        && *z.real() <= 0.25
2622    {
2623        let p = z.prec().0 as usize / 4;
2624        gamma(s.clone()) - lower_incomplete_gamma_recursion(s, z, 0, p)
2625    } else {
2626        let p = z.prec().0 as usize / 4;
2627        incomplete_gamma_recursion(s, z, 0, p)
2628    }
2629}
2630pub fn lower_incomplete_gamma(s: Complex, z: Complex) -> Complex {
2631    if s.real().is_sign_positive() && !s.real().clone().fract().is_zero() && *z.real() <= 1 {
2632        let p = z.prec().0 as usize / 4;
2633        lower_incomplete_gamma_recursion(s, z, 0, p)
2634    } else {
2635        gamma(s.clone()) - incomplete_gamma(s, z)
2636    }
2637}
2638pub fn eta(s: Complex) -> Complex {
2639    let prec = s.prec().0;
2640    let mut sum = Complex::new(prec);
2641    let two = Float::with_val(prec, 2);
2642    for n in 0..=(prec / 16).max(16) {
2643        let mut innersum = Complex::new(prec);
2644        let nb = Integer::from(n);
2645        for k in 0..=n {
2646            let num = nb.clone().binomial(k) * pow_nth(Complex::with_val(prec, 1 + k), -s.clone());
2647            if k % 2 == 0 {
2648                innersum += num;
2649            } else {
2650                innersum -= num;
2651            }
2652        }
2653        sum += innersum / two.clone().pow(n + 1)
2654    }
2655    sum
2656}
2657pub fn zeta(s: Complex) -> Complex {
2658    eta(s.clone()) / (1 - pow_nth(Complex::with_val(s.prec(), 2), 1 - s))
2659}
2660pub fn euleriannumbers(n: Complex, k: i32) -> Complex {
2661    if k < 0 {
2662        Complex::with_val(n.prec(), Nan)
2663    } else if n.real().clone().fract() != 0 && n.imag().is_zero() && n.real().is_sign_positive() {
2664        let mut sum = Complex::new(n.prec());
2665        for i in 0..=k {
2666            let ic = Complex::with_val(n.prec(), i);
2667            let num: Complex = k - ic.clone() + 1;
2668            let num = binomial(n.clone() + 1, ic) * pow_nth(num, n.clone());
2669            if i % 2 == 0 { sum += num } else { sum -= num }
2670        }
2671        sum
2672    } else {
2673        Complex::with_val(
2674            n.prec(),
2675            euleriannumbersint(
2676                n.real()
2677                    .to_integer()
2678                    .unwrap_or_default()
2679                    .to_u32()
2680                    .unwrap_or_default(),
2681                k as u32,
2682            ),
2683        )
2684    }
2685}
2686pub fn euleriannumbersint(n: u32, k: u32) -> Integer {
2687    let mut sum = Integer::new();
2688    for i in 0..=k {
2689        let num: Integer = k - Integer::from(i) + 1;
2690        let num = Integer::from(n + 1).binomial(i) * num.pow(n);
2691        if i % 2 == 0 { sum += num } else { sum -= num }
2692    }
2693    sum
2694}
2695pub fn binomial(a: Complex, b: Complex) -> Complex {
2696    if a.imag().is_zero()
2697        && b.imag().is_zero()
2698        && a.real().clone().fract().is_zero()
2699        && b.real().clone().fract().is_zero()
2700        && a.real().is_finite()
2701        && a.real() >= &0
2702        && b.real() >= &0
2703    {
2704        Complex::with_val(
2705            a.prec(),
2706            a.real().to_integer().unwrap_or_default().binomial(
2707                b.real()
2708                    .to_integer()
2709                    .unwrap_or_default()
2710                    .to_u32()
2711                    .unwrap_or_default(),
2712            ),
2713        )
2714    } else if a.real().is_sign_negative()
2715        && a.real().clone().fract().is_zero()
2716        && a.imag().is_zero()
2717        && b.imag().is_zero()
2718    {
2719        let prec = a.prec().0;
2720        let a = a + Complex::with_val(prec, (0, 1)) * Float::with_val(prec, 0.5).pow(prec / 2);
2721        (gamma(a.clone() + 1) / (gamma(b.clone() + 1) * gamma(a.clone() - b.clone() + 1)))
2722            .real()
2723            .clone()
2724            .into()
2725    } else {
2726        gamma(a.clone() + 1) / (gamma(b.clone() + 1) * gamma(a.clone() - b.clone() + 1))
2727    }
2728}
2729fn incomplete_gamma_recursion(s: Complex, z: Complex, iter: usize, max: usize) -> Complex {
2730    if iter == max {
2731        Complex::with_val(s.prec(), 1)
2732    } else if iter == 0 {
2733        (pow_nth(z.clone(), s.clone()) / z.clone().exp()) / incomplete_gamma_recursion(s, z, 1, max)
2734    } else if iter % 2 == 1 {
2735        z.clone()
2736            + ((iter.div_ceil(2) - s.clone()) / incomplete_gamma_recursion(s, z, iter + 1, max))
2737    } else {
2738        1 + (iter.div_ceil(2) / incomplete_gamma_recursion(s, z, iter + 1, max))
2739    }
2740}
2741fn lower_incomplete_gamma_recursion(s: Complex, z: Complex, iter: usize, max: usize) -> Complex {
2742    if iter == max {
2743        Complex::with_val(s.prec(), 1)
2744    } else if iter == 0 {
2745        (pow_nth(z.clone(), s.clone()) / z.clone().exp())
2746            / lower_incomplete_gamma_recursion(s, z, 1, max)
2747    } else if iter % 2 == 1 {
2748        s.clone() + iter
2749            - 1
2750            - ((s.clone() + iter - 1) * z.clone())
2751                / lower_incomplete_gamma_recursion(s, z, iter + 1, max)
2752    } else {
2753        s.clone() + iter - 1
2754            + (iter / 2 * z.clone() / lower_incomplete_gamma_recursion(s, z, iter + 1, max))
2755    }
2756}
2757pub fn subfactorial(z: Complex) -> Complex {
2758    subfactorial_recursion(z.clone(), 0, (z.prec().0 as usize / 4).max(32))
2759        + gamma(z.clone() + 1) / Float::with_val(z.prec().0, 1).exp()
2760}
2761fn subfactorial_recursion(z: Complex, iter: usize, max: usize) -> Complex {
2762    if iter == max {
2763        Complex::with_val(z.prec(), 1)
2764    } else if iter == 0 {
2765        Complex::with_val(z.prec(), -1).pow(z.clone()) / subfactorial_recursion(z, 1, max)
2766    } else {
2767        (z.clone() + iter + 1) - iter / subfactorial_recursion(z, iter + 1, max)
2768    }
2769}
2770#[allow(clippy::too_many_arguments)]
2771pub fn surface_area(
2772    func: Vec<NumStr>,
2773    func_vars: Vec<(String, Vec<NumStr>)>,
2774    options: Options,
2775    varx: String,
2776    mut startx: Complex,
2777    endx: Number,
2778    vary: String,
2779    starty_func: Vec<NumStr>,
2780    endy_func: Vec<NumStr>,
2781) -> Result<Number, &'static str> {
2782    if starty_func.is_empty() || endy_func.is_empty() {
2783        return Err("bad start/end");
2784    }
2785    let points = options.prec as usize / 16;
2786    let unitsx = endx.units;
2787    let endx = endx.number;
2788    let deltax: Complex = (endx.clone() - startx.clone()) / (points - 1);
2789    let mut area: Complex = Complex::new(options.prec);
2790    if let (Ok(Num(start)), Ok(Num(end))) = (
2791        do_math(starty_func.clone(), options, func_vars.clone()),
2792        do_math(endy_func.clone(), options, func_vars.clone()),
2793    ) {
2794        let starty = start.number;
2795        let endy = end.number;
2796        let unitsy = end.units;
2797        let deltay: Complex = (endy.clone() - starty.clone()) / (points - 1);
2798        let res = {
2799            let n = NumStr::new(Number::from(startx.clone(), unitsx));
2800            let func = place_var(func.clone(), &varx, n.clone());
2801            let func_vars = place_funcvar(func_vars.clone(), &varx, n);
2802            do_math_with_var(
2803                func.clone(),
2804                options,
2805                func_vars.clone(),
2806                &vary,
2807                NumStr::new(Number::from(starty.clone(), unitsy)),
2808            )?
2809        };
2810        match res {
2811            Num(n) => {
2812                let mut data: Vec<Vec<Complex>> = vec![Vec::new(); points];
2813                let units = n.units.map(|n| n.pow(2.0));
2814                data[0].push(n.number);
2815                for (nx, row) in data.iter_mut().enumerate() {
2816                    if nx + 1 == points {
2817                        startx.clone_from(&endx)
2818                    } else if nx != 0 {
2819                        startx += deltax.clone();
2820                    }
2821                    let n = NumStr::new(Number::from(startx.clone(), unitsx));
2822                    let mut func = place_var(func.clone(), &varx, n.clone());
2823                    let mut func_vars = place_funcvar(func_vars.clone(), &varx, n);
2824                    simplify(&mut func, &mut func_vars, options);
2825                    let mut starty = starty.clone();
2826                    for ny in 0..points {
2827                        if nx == 0 && ny == 0 {
2828                            continue;
2829                        }
2830                        if ny + 1 == points {
2831                            starty.clone_from(&endy)
2832                        } else if ny != 0 {
2833                            starty += deltay.clone();
2834                        }
2835                        row.push(
2836                            do_math_with_var(
2837                                func.clone(),
2838                                options,
2839                                func_vars.clone(),
2840                                &vary,
2841                                NumStr::new(Number::from(starty.clone(), unitsy)),
2842                            )?
2843                            .num()?
2844                            .number,
2845                        )
2846                    }
2847                }
2848                let k: Complex = sqr(deltax.clone()) * sqr(deltay.clone());
2849                for (nx, row) in data[..data.len() - 1].iter().enumerate() {
2850                    for (ny, z) in row[..row.len() - 1].iter().enumerate() {
2851                        let a = row[ny + 1].clone();
2852                        let b = data[nx + 1][ny].clone();
2853                        {
2854                            let a = a.clone() - z;
2855                            let b = b.clone() - z;
2856                            let i = deltay.clone() * b;
2857                            let j = deltax.clone() * a;
2858                            let i: Complex = sqr(i);
2859                            let j: Complex = sqr(j);
2860                            area += (i + j + k.clone()).sqrt() / 2;
2861                        }
2862                        {
2863                            let z = &data[nx + 1][ny + 1];
2864                            let a = a - z;
2865                            let b = b - z;
2866                            let i = deltax.clone() * b;
2867                            let j = deltay.clone() * a;
2868                            let i: Complex = sqr(i);
2869                            let j: Complex = sqr(j);
2870                            area += (i + j + k.clone()).sqrt() / 2
2871                        }
2872                    }
2873                }
2874                Ok(Number::from(area, units))
2875            }
2876            Vector(v) if v.len() == 3 => {
2877                let mut data: Vec<Vec<Vec<Complex>>> = vec![Vec::new(); points];
2878                let units = v[2].units.map(|n| n.pow(2.0));
2879                data[0].push(v.iter().map(|a| a.number.clone()).collect());
2880                for (nx, row) in data.iter_mut().enumerate() {
2881                    if nx + 1 == points {
2882                        startx.clone_from(&endx)
2883                    } else if nx != 0 {
2884                        startx += deltax.clone();
2885                    }
2886                    let n = NumStr::new(Number::from(startx.clone(), unitsx));
2887                    let mut func = place_var(func.clone(), &varx, n.clone());
2888                    let mut func_vars = place_funcvar(func_vars.clone(), &varx, n);
2889                    simplify(&mut func, &mut func_vars, options);
2890                    let mut starty = starty.clone();
2891                    for ny in 0..points {
2892                        if nx == 0 && ny == 0 {
2893                            continue;
2894                        }
2895                        if ny + 1 == points {
2896                            starty.clone_from(&endy)
2897                        } else if ny != 0 {
2898                            starty += deltay.clone();
2899                        }
2900                        let v = do_math_with_var(
2901                            func.clone(),
2902                            options,
2903                            func_vars.clone(),
2904                            &vary,
2905                            NumStr::new(Number::from(starty.clone(), unitsy)),
2906                        )?
2907                        .vec()?;
2908                        row.push(v.iter().map(|a| a.number.clone()).collect());
2909                    }
2910                }
2911                for (nx, row) in data[..data.len() - 1].iter().enumerate() {
2912                    for (ny, z) in row[..row.len() - 1].iter().enumerate() {
2913                        let a = row[ny + 1].clone();
2914                        let b = data[nx + 1][ny].clone();
2915                        {
2916                            let a0 = a[0].clone() - z[0].clone();
2917                            let a1 = a[1].clone() - z[1].clone();
2918                            let a2 = a[2].clone() - z[2].clone();
2919                            let b0 = b[0].clone() - z[0].clone();
2920                            let b1 = b[1].clone() - z[1].clone();
2921                            let b2 = b[2].clone() - z[2].clone();
2922                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
2923                            let j = a2 * b0.clone() - a0.clone() * b2;
2924                            let k = a0 * b1 - a1 * b0;
2925                            let i: Complex = sqr(i);
2926                            let j: Complex = sqr(j);
2927                            let k: Complex = sqr(k);
2928                            area += (i + j + k).sqrt() / 2;
2929                        }
2930                        {
2931                            let z = &data[nx + 1][ny + 1];
2932                            let a0 = a[0].clone() - z[0].clone();
2933                            let a1 = a[1].clone() - z[1].clone();
2934                            let a2 = a[2].clone() - z[2].clone();
2935                            let b0 = b[0].clone() - z[0].clone();
2936                            let b1 = b[1].clone() - z[1].clone();
2937                            let b2 = b[2].clone() - z[2].clone();
2938                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
2939                            let j = a2 * b0.clone() - a0.clone() * b2;
2940                            let k = a0 * b1 - a1 * b0;
2941                            let i: Complex = sqr(i);
2942                            let j: Complex = sqr(j);
2943                            let k: Complex = sqr(k);
2944                            area += (i + j + k).sqrt() / 2;
2945                        }
2946                    }
2947                }
2948                Ok(Number::from(area, units))
2949            }
2950            _ => Err("bad input"),
2951        }
2952    } else {
2953        let n = NumStr::new(Number::from(startx.clone(), unitsx));
2954        let temp_func = place_var(func.clone(), &varx, n.clone());
2955        let temp_func_vars = place_funcvar(func_vars.clone(), &varx, n.clone());
2956        let starty = do_math_with_var(
2957            starty_func.clone(),
2958            options,
2959            temp_func_vars.clone(),
2960            &varx,
2961            n.clone(),
2962        )?
2963        .num()?
2964        .number;
2965        let end = do_math_with_var(
2966            endy_func.clone(),
2967            options,
2968            temp_func_vars.clone(),
2969            &varx,
2970            n.clone(),
2971        )?
2972        .num()?;
2973        let mut endy = end.number;
2974        let unitsy = end.units;
2975        let mut deltay: Complex = (endy.clone() - starty.clone()) / (points - 1);
2976        match do_math_with_var(
2977            temp_func.clone(),
2978            options,
2979            temp_func_vars.clone(),
2980            &vary,
2981            NumStr::new(Number::from(starty.clone(), unitsy)),
2982        )? {
2983            Num(n) => {
2984                let mut data: Vec<Vec<[Complex; 2]>> = vec![Vec::new(); points];
2985                let units = n.units.map(|n| n.pow(2.0));
2986                data[0].push([starty.clone(), n.number]);
2987                for (nx, row) in data.iter_mut().enumerate() {
2988                    if nx + 1 == points {
2989                        startx.clone_from(&endx)
2990                    } else if nx != 0 {
2991                        startx += deltax.clone();
2992                    }
2993                    let n = NumStr::new(Number::from(startx.clone(), unitsx));
2994                    let mut func = place_var(func.clone(), &varx, n.clone());
2995                    let mut func_vars = place_funcvar(func_vars.clone(), &varx, n.clone());
2996                    simplify(&mut func, &mut func_vars, options);
2997                    let mut starty = starty.clone();
2998                    if nx != 0 {
2999                        starty = do_math_with_var(
3000                            starty_func.clone(),
3001                            options,
3002                            func_vars.clone(),
3003                            &varx,
3004                            n.clone(),
3005                        )?
3006                        .num()?
3007                        .number;
3008                        endy = do_math_with_var(
3009                            endy_func.clone(),
3010                            options,
3011                            func_vars.clone(),
3012                            &varx,
3013                            n.clone(),
3014                        )?
3015                        .num()?
3016                        .number;
3017                        deltay = (endy.clone() - starty.clone()) / (points - 1);
3018                    }
3019                    for ny in 0..points {
3020                        if nx == 0 && ny == 0 {
3021                            continue;
3022                        }
3023                        if ny + 1 == points {
3024                            starty.clone_from(&endy)
3025                        } else if ny != 0 {
3026                            starty += deltay.clone();
3027                        }
3028                        row.push([
3029                            starty.clone(),
3030                            do_math_with_var(
3031                                func.clone(),
3032                                options,
3033                                func_vars.clone(),
3034                                &vary,
3035                                NumStr::new(Number::from(starty.clone(), unitsy)),
3036                            )?
3037                            .num()?
3038                            .number,
3039                        ])
3040                    }
3041                }
3042                for (nx, row) in data[..data.len() - 1].iter().enumerate() {
3043                    for (ny, z) in row[..row.len() - 1].iter().enumerate() {
3044                        let a = row[ny + 1].clone();
3045                        let b = data[nx + 1][ny].clone();
3046                        {
3047                            let a1 = a[0].clone() - z[0].clone();
3048                            let a2 = a[1].clone() - z[1].clone();
3049                            let b0 = deltax.clone();
3050                            let b1 = b[0].clone() - z[0].clone();
3051                            let b2 = b[1].clone() - z[1].clone();
3052                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
3053                            let j = a2 * b0.clone();
3054                            let k = a1 * b0;
3055                            let i: Complex = sqr(i);
3056                            let j: Complex = sqr(j);
3057                            let k: Complex = sqr(k);
3058                            area += (i + j + k).sqrt() / 2;
3059                        }
3060                        {
3061                            let z = &data[nx + 1][ny + 1];
3062                            let a0 = deltax.clone();
3063                            let a1 = a[0].clone() - z[0].clone();
3064                            let a2 = a[1].clone() - z[1].clone();
3065                            let b1 = b[0].clone() - z[0].clone();
3066                            let b2 = b[1].clone() - z[1].clone();
3067                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
3068                            let j = a0.clone() * b2;
3069                            let k = a0 * b1;
3070                            let i: Complex = sqr(i);
3071                            let j: Complex = sqr(j);
3072                            let k: Complex = sqr(k);
3073                            area += (i + j + k).sqrt() / 2;
3074                        }
3075                    }
3076                }
3077                Ok(Number::from(area, units))
3078            }
3079            Vector(v) if v.len() == 3 => {
3080                let mut data: Vec<Vec<Vec<Complex>>> = vec![Vec::new(); points];
3081                let units = v[2].units.map(|n| n.pow(2.0));
3082                data[0].push(v.iter().map(|a| a.number.clone()).collect());
3083                for (nx, row) in data.iter_mut().enumerate() {
3084                    if nx + 1 == points {
3085                        startx.clone_from(&endx)
3086                    } else if nx != 0 {
3087                        startx += deltax.clone();
3088                    }
3089                    let n = NumStr::new(Number::from(startx.clone(), unitsx));
3090                    let mut func = place_var(func.clone(), &varx, n.clone());
3091                    let mut func_vars = place_funcvar(func_vars.clone(), &varx, n.clone());
3092                    simplify(&mut func, &mut func_vars, options);
3093                    let mut starty = starty.clone();
3094                    if nx != 0 {
3095                        starty = do_math_with_var(
3096                            starty_func.clone(),
3097                            options,
3098                            func_vars.clone(),
3099                            &varx,
3100                            n.clone(),
3101                        )?
3102                        .num()?
3103                        .number;
3104                        endy = do_math_with_var(
3105                            endy_func.clone(),
3106                            options,
3107                            func_vars.clone(),
3108                            &varx,
3109                            n.clone(),
3110                        )?
3111                        .num()?
3112                        .number;
3113                        deltay = (endy.clone() - starty.clone()) / (points - 1);
3114                    }
3115                    for ny in 0..points {
3116                        if nx == 0 && ny == 0 {
3117                            continue;
3118                        }
3119                        if ny + 1 == points {
3120                            starty.clone_from(&endy)
3121                        } else if ny != 0 {
3122                            starty += deltay.clone();
3123                        }
3124                        let v = do_math_with_var(
3125                            func.clone(),
3126                            options,
3127                            func_vars.clone(),
3128                            &vary,
3129                            NumStr::new(Number::from(starty.clone(), unitsy)),
3130                        )?
3131                        .vec()?;
3132                        row.push(v.iter().map(|a| a.number.clone()).collect());
3133                    }
3134                }
3135                for (nx, row) in data[..data.len() - 1].iter().enumerate() {
3136                    for (ny, z) in row[..row.len() - 1].iter().enumerate() {
3137                        let a = row[ny + 1].clone();
3138                        let b = data[nx + 1][ny].clone();
3139                        {
3140                            let a0 = a[0].clone() - z[0].clone();
3141                            let a1 = a[1].clone() - z[1].clone();
3142                            let a2 = a[2].clone() - z[2].clone();
3143                            let b0 = b[0].clone() - z[0].clone();
3144                            let b1 = b[1].clone() - z[1].clone();
3145                            let b2 = b[2].clone() - z[2].clone();
3146                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
3147                            let j = a2 * b0.clone() - a0.clone() * b2;
3148                            let k = a0 * b1 - a1 * b0;
3149                            let i: Complex = sqr(i);
3150                            let j: Complex = sqr(j);
3151                            let k: Complex = sqr(k);
3152                            area += (i + j + k).sqrt() / 2;
3153                        }
3154                        {
3155                            let z = &data[nx + 1][ny + 1];
3156                            let a0 = a[0].clone() - z[0].clone();
3157                            let a1 = a[1].clone() - z[1].clone();
3158                            let a2 = a[2].clone() - z[2].clone();
3159                            let b0 = b[0].clone() - z[0].clone();
3160                            let b1 = b[1].clone() - z[1].clone();
3161                            let b2 = b[2].clone() - z[2].clone();
3162                            let i = a1.clone() * b2.clone() - a2.clone() * b1.clone();
3163                            let j = a2 * b0.clone() - a0.clone() * b2;
3164                            let k = a0 * b1 - a1 * b0;
3165                            let i: Complex = sqr(i);
3166                            let j: Complex = sqr(j);
3167                            let k: Complex = sqr(k);
3168                            area += (i + j + k).sqrt() / 2;
3169                        }
3170                    }
3171                }
3172                Ok(Number::from(area, units))
3173            }
3174            _ => Err("bad input"),
3175        }
3176    }
3177}
3178pub fn length(
3179    func: Vec<NumStr>,
3180    func_vars: Vec<(String, Vec<NumStr>)>,
3181    options: Options,
3182    var: String,
3183    mut start: Complex,
3184    end: Number,
3185) -> Result<Number, &'static str> {
3186    let points = options.prec as usize / 4;
3187    let units = end.units;
3188    let end = end.number;
3189    let delta: Complex = (end.clone() - start.clone()) / points;
3190    let delta2: Complex = sqr(delta.clone());
3191    let mut x0: NumStr = do_math_with_var(
3192        func.clone(),
3193        options,
3194        func_vars.clone(),
3195        &var.clone(),
3196        NumStr::new(Number::from(start.clone(), units)),
3197    )?;
3198    let length_units = match x0.clone() {
3199        Num(a) => a.units,
3200        Vector(a) => {
3201            if a.iter().all(|b| b.units == a[0].units) {
3202                a[0].units
3203            } else {
3204                None
3205            }
3206        }
3207        _ => return Err("not supported arc length data"),
3208    };
3209    let mut length = Complex::new(options.prec);
3210    for i in 0..points {
3211        if i + 1 == points {
3212            start.clone_from(&end)
3213        } else {
3214            start += delta.clone();
3215        }
3216        let x1 = do_math_with_var(
3217            func.clone(),
3218            options,
3219            func_vars.clone(),
3220            &var.clone(),
3221            NumStr::new(Number::from(start.clone(), units)),
3222        )?;
3223        match (x0, x1) {
3224            (Num(xi), Num(xf)) => {
3225                let nl: Complex = sqr(xf.number.clone() - xi.number) + delta2.clone();
3226                length += nl.sqrt();
3227                x0 = Num(xf);
3228            }
3229            (Vector(xi), Vector(xf)) if xf.len() == 1 => {
3230                let nl: Complex = sqr(xf[0].number.clone() - xi[0].number.clone()) + delta2.clone();
3231                length += nl.sqrt();
3232                x0 = Vector(xf);
3233            }
3234            (Vector(xi), Vector(xf)) => {
3235                let nl: Complex = xi
3236                    .iter()
3237                    .zip(xf.clone())
3238                    .fold(Complex::new(options.prec), |sum, x| {
3239                        sum + sqr(x.1.number - x.0.number.clone())
3240                    });
3241                length += nl.sqrt();
3242                x0 = Vector(xf);
3243            }
3244            (_, _) => return Err("not supported arc length data"),
3245        };
3246    }
3247    Ok(Number::from(length, length_units))
3248}
3249#[allow(clippy::too_many_arguments)]
3250pub fn area(
3251    mut func: Vec<NumStr>,
3252    mut func_vars: Vec<(String, Vec<NumStr>)>,
3253    options: Options,
3254    var: String,
3255    mut start: Complex,
3256    mut end: Number,
3257    nth: Complex,
3258    combine: bool,
3259) -> Result<NumStr, &'static str> {
3260    let mut negate = false;
3261    if start.real() > end.number.real() && nth == 1 {
3262        negate = true;
3263        (end.number, start) = (start, end.number)
3264    }
3265    if start.real().is_infinite() {
3266        let neg = start.real().is_sign_negative();
3267        start = Complex::with_val(options.prec, 1) << (options.prec / 16);
3268        if neg {
3269            start *= -1
3270        }
3271    }
3272    if end.number.real().is_infinite() {
3273        let neg = end.number.real().is_sign_negative();
3274        end.number = Complex::with_val(options.prec, 1) << (options.prec / 16);
3275        if neg {
3276            end.number *= -1
3277        }
3278    }
3279    let points = options.prec as usize / 4;
3280    let units = end.units;
3281    let mut end = end.number;
3282    let mut funcs = Vec::new();
3283    if combine
3284        && !func.is_empty()
3285        && func[0] == LeftCurlyBracket
3286        && func[func.len() - 1] == RightCurlyBracket
3287    {
3288        let mut brackets = 0;
3289        let mut last = 1;
3290        for (i, f) in func.iter().enumerate() {
3291            match f {
3292                LeftBracket | LeftCurlyBracket => brackets += 1,
3293                RightBracket | RightCurlyBracket => brackets -= 1,
3294                Comma if brackets == 1 => {
3295                    funcs.push(func[last..i].to_vec());
3296                    last = i + 1;
3297                }
3298                _ => {}
3299            }
3300        }
3301        if last != 1 {
3302            func = func[last..func.len().saturating_sub(1)].to_vec();
3303        }
3304    }
3305    let mut areavec: Vec<Number> = Vec::new();
3306    let div = Float::with_val(options.prec, 0.5).pow(options.prec / 2);
3307    let mut delta: Complex = (end.clone() - start.clone()) / points;
3308    let mut area: Complex = Complex::new(options.prec);
3309    let mut x0 = do_math_with_var(
3310        func.clone(),
3311        options,
3312        func_vars.clone(),
3313        &var,
3314        NumStr::new(Number::from(start.clone(), units)),
3315    )?;
3316    if start == end {
3317        return match x0 {
3318            Num(_) => Ok(NumStr::new(Number::from(Complex::new(options.prec), None))),
3319            Vector(_) => Ok(Vector(Vec::new())),
3320            _ => Err(
3321                "not supported area data, if parametric have the 2nd arg start and end with the { } brackets",
3322            ),
3323        };
3324    }
3325    {
3326        fn check_bounds(
3327            func: Vec<NumStr>,
3328            func_vars: Vec<(String, Vec<NumStr>)>,
3329            options: Options,
3330            var: String,
3331            delta: &mut Complex,
3332            units: Option<Units>,
3333            x0: &mut NumStr,
3334            start: &mut Complex,
3335            right: bool,
3336            compute: bool,
3337        ) -> Result<(), &'static str> {
3338            let mut has = false;
3339            if let Ok(a) = x0.num() {
3340                if !a.number.real().is_finite() {
3341                    let points = options.prec as usize / 4;
3342                    let end = if right {
3343                        points * delta.clone() + start.clone()
3344                    } else {
3345                        start.clone() - points * delta.clone()
3346                    };
3347                    if right {
3348                        *start += delta.clone() >> 4;
3349                    } else {
3350                        *start -= delta.clone() >> 4;
3351                    }
3352                    *delta = if right {
3353                        (end.clone() - start.clone()) / points
3354                    } else {
3355                        (start.clone() - end.clone()) / points
3356                    };
3357                    *x0 = do_math_with_var(
3358                        func.clone(),
3359                        options,
3360                        func_vars.clone(),
3361                        &var,
3362                        NumStr::new(Number::from(start.clone(), units)),
3363                    )?;
3364                    has = true
3365                }
3366            }
3367            if !has && compute {
3368                *x0 = do_math_with_var(
3369                    func.clone(),
3370                    options,
3371                    func_vars.clone(),
3372                    &var,
3373                    NumStr::new(Number::from(start.clone(), units)),
3374                )?;
3375                check_bounds(
3376                    func, func_vars, options, var, delta, units, x0, start, right, false,
3377                )?
3378            }
3379            Ok(())
3380        }
3381        check_bounds(
3382            func.clone(),
3383            func_vars.clone(),
3384            options,
3385            var.clone(),
3386            &mut delta,
3387            units,
3388            &mut x0,
3389            &mut start,
3390            true,
3391            false,
3392        )?;
3393        let mut endv = do_math_with_var(
3394            func.clone(),
3395            options,
3396            func_vars.clone(),
3397            &var,
3398            NumStr::new(Number::from(end.clone(), units)),
3399        )?;
3400        check_bounds(
3401            func.clone(),
3402            func_vars.clone(),
3403            options,
3404            var.clone(),
3405            &mut delta,
3406            units,
3407            &mut endv,
3408            &mut end,
3409            false,
3410            false,
3411        )?;
3412        if nth == 1 && end.real().clone() - start.real().clone() > 2.0 {
3413            let mut small_start = false;
3414            if start.real().is_sign_negative() {
3415                if let Ok(a) = x0.num() {
3416                    if a.number.real().clone().abs() < Float::with_val(options.prec, 1) >> 16 {
3417                        let a = do_math_with_var(
3418                            func.clone(),
3419                            options,
3420                            func_vars.clone(),
3421                            &var,
3422                            NumStr::new(Number::from(start.clone() * 2 - 1.141, units)),
3423                        )?
3424                        .num()?;
3425                        if a.number.real().clone().abs() < Float::with_val(options.prec, 1) >> 17 {
3426                            let a = do_math_with_var(
3427                                func.clone(),
3428                                options,
3429                                func_vars.clone(),
3430                                &var,
3431                                NumStr::new(Number::from(
3432                                    -(Complex::with_val(options.prec, 4) << (options.prec / 16)),
3433                                    units,
3434                                )),
3435                            )?
3436                            .num()?;
3437                            if a.number.real().clone().abs()
3438                                < Float::with_val(options.prec, 1) >> 18
3439                            {
3440                                small_start = true;
3441                            }
3442                        }
3443                    }
3444                }
3445            }
3446            let mut small_end = false;
3447            if end.real().is_sign_positive() {
3448                if let Ok(a) = endv.num() {
3449                    if a.number.real().clone().abs() < Float::with_val(options.prec, 1) >> 16 {
3450                        let a = do_math_with_var(
3451                            func.clone(),
3452                            options,
3453                            func_vars.clone(),
3454                            &var,
3455                            NumStr::new(Number::from(end.clone() * 2 + 1.141, units)),
3456                        )?
3457                        .num()?;
3458                        if a.number.real().clone().abs() < Float::with_val(options.prec, 1) >> 17 {
3459                            let a = do_math_with_var(
3460                                func.clone(),
3461                                options,
3462                                func_vars.clone(),
3463                                &var,
3464                                NumStr::new(Number::from(
3465                                    Complex::with_val(options.prec, 4) << (options.prec / 16),
3466                                    units,
3467                                )),
3468                            )?
3469                            .num()?;
3470                            if a.number.real().clone().abs()
3471                                < Float::with_val(options.prec, 1) >> 18
3472                            {
3473                                small_end = true
3474                            }
3475                        }
3476                    }
3477                }
3478            }
3479            let two = NumStr::new(Number::from(Complex::with_val(options.prec, 2), None));
3480            let one = NumStr::new(Number::from(Complex::with_val(options.prec, 1), None));
3481            fn change_var(
3482                func: &mut [NumStr],
3483                func_vars: &mut Vec<(String, Vec<NumStr>)>,
3484                from: String,
3485                fto: Vec<NumStr>,
3486            ) {
3487                let to = format!("@!@{from}@");
3488                for v in func.iter_mut() {
3489                    if *v == Func(from.clone()) {
3490                        *v = Func(to.clone())
3491                    }
3492                }
3493                for func in func_vars.iter_mut() {
3494                    for v in func.1.iter_mut() {
3495                        if *v == Func(from.clone()) {
3496                            *v = Func(to.clone())
3497                        }
3498                    }
3499                }
3500                func_vars.push((to, fto));
3501            }
3502            match (small_start, small_end) {
3503                (true, true) => {
3504                    change_var(
3505                        &mut func,
3506                        &mut func_vars,
3507                        var.clone(),
3508                        vec![
3509                            Func(var.clone()),
3510                            Division,
3511                            LeftBracket,
3512                            one.clone(),
3513                            Minus,
3514                            Func(var.clone()),
3515                            Exponent,
3516                            two.clone(),
3517                            RightBracket,
3518                        ],
3519                    );
3520                    func.insert(0, LeftBracket);
3521                    func.extend(vec![
3522                        RightBracket,
3523                        Multiplication,
3524                        LeftBracket,
3525                        one.clone(),
3526                        Plus,
3527                        Func(var.clone()),
3528                        Exponent,
3529                        two.clone(),
3530                        RightBracket,
3531                        Division,
3532                        LeftBracket,
3533                        one,
3534                        Minus,
3535                        Func(var.clone()),
3536                        Exponent,
3537                        two.clone(),
3538                        RightBracket,
3539                        Exponent,
3540                        two,
3541                    ]);
3542                    start = Complex::with_val(options.prec, -1);
3543                    end = Complex::with_val(options.prec, 1);
3544                    delta = 2 * end.clone() / points;
3545                    check_bounds(
3546                        func.clone(),
3547                        func_vars.clone(),
3548                        options,
3549                        var.clone(),
3550                        &mut delta,
3551                        units,
3552                        &mut x0,
3553                        &mut start,
3554                        true,
3555                        true,
3556                    )?;
3557                    check_bounds(
3558                        func.clone(),
3559                        func_vars.clone(),
3560                        options,
3561                        var.clone(),
3562                        &mut delta,
3563                        units,
3564                        &mut endv,
3565                        &mut end,
3566                        false,
3567                        true,
3568                    )?;
3569                }
3570                (false, true) => {
3571                    change_var(
3572                        &mut func,
3573                        &mut func_vars,
3574                        var.clone(),
3575                        vec![
3576                            NumStr::new(Number::from(start, None)),
3577                            Plus,
3578                            Func(var.clone()),
3579                            Division,
3580                            LeftBracket,
3581                            one.clone(),
3582                            Minus,
3583                            Func(var.clone()),
3584                            RightBracket,
3585                        ],
3586                    );
3587                    func.insert(0, LeftBracket);
3588                    func.extend(vec![
3589                        RightBracket,
3590                        Division,
3591                        LeftBracket,
3592                        one,
3593                        Minus,
3594                        Func(var.clone()),
3595                        RightBracket,
3596                        Exponent,
3597                        two,
3598                    ]);
3599                    start = Complex::new(options.prec);
3600                    end = Complex::with_val(options.prec, 1);
3601                    delta = end.clone() / points;
3602                    check_bounds(
3603                        func.clone(),
3604                        func_vars.clone(),
3605                        options,
3606                        var.clone(),
3607                        &mut delta,
3608                        units,
3609                        &mut x0,
3610                        &mut start,
3611                        true,
3612                        true,
3613                    )?;
3614                    check_bounds(
3615                        func.clone(),
3616                        func_vars.clone(),
3617                        options,
3618                        var.clone(),
3619                        &mut delta,
3620                        units,
3621                        &mut endv,
3622                        &mut end,
3623                        false,
3624                        true,
3625                    )?;
3626                }
3627                (true, false) => {
3628                    change_var(
3629                        &mut func,
3630                        &mut func_vars,
3631                        var.clone(),
3632                        vec![
3633                            NumStr::new(Number::from(end, None)),
3634                            Minus,
3635                            LeftBracket,
3636                            one.clone(),
3637                            Minus,
3638                            Func(var.clone()),
3639                            RightBracket,
3640                            Division,
3641                            Func(var.clone()),
3642                        ],
3643                    );
3644                    func.insert(0, LeftBracket);
3645                    func.extend(vec![
3646                        RightBracket,
3647                        Division,
3648                        Func(var.clone()),
3649                        Exponent,
3650                        two,
3651                    ]);
3652                    start = Complex::new(options.prec);
3653                    end = Complex::with_val(options.prec, 1);
3654                    delta = end.clone() / points;
3655                    check_bounds(
3656                        func.clone(),
3657                        func_vars.clone(),
3658                        options,
3659                        var.clone(),
3660                        &mut delta,
3661                        units,
3662                        &mut x0,
3663                        &mut start,
3664                        true,
3665                        true,
3666                    )?;
3667                    check_bounds(
3668                        func.clone(),
3669                        func_vars.clone(),
3670                        options,
3671                        var.clone(),
3672                        &mut delta,
3673                        units,
3674                        &mut endv,
3675                        &mut end,
3676                        false,
3677                        true,
3678                    )?;
3679                }
3680                _ => {}
3681            }
3682        }
3683    }
3684    let yunits = if let Num(ref a) = x0 { a.units } else { None };
3685    if !funcs.is_empty() {
3686        let mut nx0t = Complex::new(options.prec);
3687        for i in &funcs {
3688            nx0t += sqr((do_math_with_var(
3689                i.clone(),
3690                options,
3691                func_vars.clone(),
3692                &var,
3693                NumStr::new(Number::from(start.clone() + div.clone(), units)),
3694            )?
3695            .num()?
3696            .number
3697                - do_math_with_var(
3698                    i.clone(),
3699                    options,
3700                    func_vars.clone(),
3701                    &var,
3702                    NumStr::new(Number::from(start.clone(), units)),
3703                )?
3704                .num()?
3705                .number)
3706                / div.clone())
3707        }
3708    }
3709    let h: Complex = delta.clone() / 4;
3710    let body = |i: usize| -> Result<(Option<Complex>, Option<Vec<Number>>), &str> {
3711        let point = if i + 1 == points {
3712            end.clone()
3713        } else {
3714            start.clone() + (i + 1) * delta.clone()
3715        };
3716        let x0 = do_math_with_var(
3717            func.clone(),
3718            options,
3719            func_vars.clone(),
3720            &var,
3721            NumStr::new(Number::from(point.clone() - 4 * h.clone(), units)),
3722        )?;
3723        let x1 = do_math_with_var(
3724            func.clone(),
3725            options,
3726            func_vars.clone(),
3727            &var,
3728            NumStr::new(Number::from(point.clone() - 3 * h.clone(), units)),
3729        )?;
3730        let x2 = do_math_with_var(
3731            func.clone(),
3732            options,
3733            func_vars.clone(),
3734            &var,
3735            NumStr::new(Number::from(point.clone() - 2 * h.clone(), units)),
3736        )?;
3737        let x3 = do_math_with_var(
3738            func.clone(),
3739            options,
3740            func_vars.clone(),
3741            &var,
3742            NumStr::new(Number::from(point.clone() - h.clone(), units)),
3743        )?;
3744        let x4 = do_math_with_var(
3745            func.clone(),
3746            options,
3747            func_vars.clone(),
3748            &var,
3749            NumStr::new(Number::from(point.clone(), units)),
3750        )?;
3751        match (x0, x1, x2, x3, x4) {
3752            (Num(nx0), Num(nx1), Num(nx2), Num(nx3), Num(nx4)) if funcs.is_empty() => {
3753                let n0;
3754                let n1;
3755                let n2;
3756                let n3;
3757                let n4;
3758                if nth != 1.0 {
3759                    let nt = pow_nth(end.clone() - point.clone() + delta.clone(), nth.clone() - 1);
3760                    n0 = nx0.number * nt;
3761                    let n: Complex = end.clone() - point.clone() + 3 * h.clone();
3762                    let nt = pow_nth(n, nth.clone() - 1);
3763                    n1 = nx1.number * nt;
3764                    let n: Complex = end.clone() - point.clone() + 2 * h.clone();
3765                    let nt = pow_nth(n, nth.clone() - 1);
3766                    n2 = nx2.number * nt;
3767                    let n: Complex = end.clone() - point.clone() + h.clone();
3768                    let nt = pow_nth(n, nth.clone() - 1);
3769                    n3 = nx3.number * nt;
3770                    let n: Complex = end.clone() - point.clone();
3771                    let nt = pow_nth(n, nth.clone() - 1);
3772                    n4 = nx4.number * nt;
3773                } else {
3774                    n0 = nx0.number;
3775                    n1 = nx1.number;
3776                    n2 = nx2.number;
3777                    n3 = nx3.number;
3778                    n4 = nx4.number;
3779                }
3780                Ok((
3781                    Some(2 * h.clone() * (7 * (n0 + n4) + 12 * n2 + 32 * (n1 + n3)) / 45),
3782                    None,
3783                ))
3784            }
3785            (Num(nx0), Num(nx1), Num(nx2), Num(nx3), Num(nx4)) => {
3786                let mut nx1t = Complex::new(options.prec);
3787                let mut nx2t = Complex::new(options.prec);
3788                let mut nx3t = Complex::new(options.prec);
3789                let mut nx4t = Complex::new(options.prec);
3790                for i in &funcs {
3791                    nx1t += sqr((do_math_with_var(
3792                        i.clone(),
3793                        options,
3794                        func_vars.clone(),
3795                        &var,
3796                        NumStr::new(Number::from(
3797                            point.clone() - 3 * h.clone() + div.clone(),
3798                            units,
3799                        )),
3800                    )?
3801                    .num()?
3802                    .number
3803                        - do_math_with_var(
3804                            i.clone(),
3805                            options,
3806                            func_vars.clone(),
3807                            &var,
3808                            NumStr::new(Number::from(point.clone() - 3 * h.clone(), units)),
3809                        )?
3810                        .num()?
3811                        .number)
3812                        / div.clone());
3813                    nx2t += sqr((do_math_with_var(
3814                        i.clone(),
3815                        options,
3816                        func_vars.clone(),
3817                        &var,
3818                        NumStr::new(Number::from(
3819                            point.clone() - 2 * h.clone() + div.clone(),
3820                            units,
3821                        )),
3822                    )?
3823                    .num()?
3824                    .number
3825                        - do_math_with_var(
3826                            i.clone(),
3827                            options,
3828                            func_vars.clone(),
3829                            &var,
3830                            NumStr::new(Number::from(point.clone() - 2 * h.clone(), units)),
3831                        )?
3832                        .num()?
3833                        .number)
3834                        / div.clone());
3835                    nx3t += sqr((do_math_with_var(
3836                        i.clone(),
3837                        options,
3838                        func_vars.clone(),
3839                        &var,
3840                        NumStr::new(Number::from(point.clone() - h.clone() + div.clone(), units)),
3841                    )?
3842                    .num()?
3843                    .number
3844                        - do_math_with_var(
3845                            i.clone(),
3846                            options,
3847                            func_vars.clone(),
3848                            &var,
3849                            NumStr::new(Number::from(point.clone() - h.clone(), units)),
3850                        )?
3851                        .num()?
3852                        .number)
3853                        / div.clone());
3854                    nx4t += sqr((do_math_with_var(
3855                        i.clone(),
3856                        options,
3857                        func_vars.clone(),
3858                        &var,
3859                        NumStr::new(Number::from(point.clone() + div.clone(), units)),
3860                    )?
3861                    .num()?
3862                    .number
3863                        - do_math_with_var(
3864                            i.clone(),
3865                            options,
3866                            func_vars.clone(),
3867                            &var,
3868                            NumStr::new(Number::from(point.clone(), units)),
3869                        )?
3870                        .num()?
3871                        .number)
3872                        / div.clone())
3873                }
3874                let nx1 = nx1.number * nx1t.sqrt();
3875                let nx2 = nx2.number * nx2t.sqrt();
3876                let nx3 = nx3.number * nx3t.sqrt();
3877                let nx4 = nx4.number * nx4t.sqrt();
3878                let n0;
3879                let n1;
3880                let n2;
3881                let n3;
3882                let n4;
3883                if nth != 1.0 {
3884                    if i == 0 {
3885                        let nt =
3886                            pow_nth(end.clone() - point.clone() + delta.clone(), nth.clone() - 1);
3887                        n0 = nx0.number * nt;
3888                    } else {
3889                        n0 = nx0.number;
3890                    }
3891                    let n: Complex = end.clone() - point.clone() + 3 * h.clone();
3892                    let nt: Complex = pow_nth(n, nth.clone() - 1);
3893                    n1 = nx1 * nt;
3894                    let n: Complex = end.clone() - point.clone() + 2 * h.clone();
3895                    let nt: Complex = pow_nth(n, nth.clone() - 1);
3896                    n2 = nx2 * nt;
3897                    let n: Complex = end.clone() - point.clone() + h.clone();
3898                    let nt: Complex = pow_nth(n, nth.clone() - 1);
3899                    n3 = nx3 * nt;
3900                    let n: Complex = end.clone() - point.clone();
3901                    let nt: Complex = pow_nth(n, nth.clone() - 1);
3902                    n4 = nx4 * nt;
3903                } else {
3904                    n0 = nx0.number;
3905                    n1 = nx1;
3906                    n2 = nx2;
3907                    n3 = nx3;
3908                    n4 = nx4;
3909                }
3910                Ok((
3911                    Some(2 * h.clone() * (7 * (n0 + n4.clone()) + 12 * n2 + 32 * (n1 + n3)) / 45),
3912                    None,
3913                ))
3914            }
3915            (Vector(nx0), Vector(nx1), Vector(nx2), Vector(nx3), Vector(nx4)) if !combine => {
3916                let mut areavec = Vec::new();
3917                for i in 0..nx0.len() {
3918                    let n0;
3919                    let n1;
3920                    let n2;
3921                    let n3;
3922                    let n4;
3923                    if nth != 1.0 {
3924                        let nt =
3925                            pow_nth(end.clone() - point.clone() + delta.clone(), nth.clone() - 1);
3926                        n0 = nx0[i].number.clone() * nt;
3927                        let n: Complex = end.clone() - point.clone() + 3 * h.clone();
3928                        let nt = pow_nth(n, nth.clone() - 1);
3929                        n1 = nx1[i].number.clone() * nt;
3930                        let n: Complex = end.clone() - point.clone() + 2 * h.clone();
3931                        let nt = pow_nth(n, nth.clone() - 1);
3932                        n2 = nx2[i].number.clone() * nt;
3933                        let n: Complex = end.clone() - point.clone() + h.clone();
3934                        let nt = pow_nth(n, nth.clone() - 1);
3935                        n3 = nx3[i].number.clone() * nt;
3936                        let n: Complex = end.clone() - point.clone();
3937                        let nt = pow_nth(n, nth.clone() - 1);
3938                        n4 = nx4[i].number.clone() * nt;
3939                    } else {
3940                        n0 = nx0[i].number.clone();
3941                        n1 = nx1[i].number.clone();
3942                        n2 = nx2[i].number.clone();
3943                        n3 = nx3[i].number.clone();
3944                        n4 = nx4[i].number.clone();
3945                    }
3946                    areavec.push(Number::from(
3947                        2 * h.clone() * (7 * (n0 + n4) + 12 * n2 + 32 * (n1 + n3)) / 45,
3948                        match (units, nx1[i].units) {
3949                            (Some(a), Some(b)) => Some(a.mul(&b)),
3950                            (Some(a), None) | (None, Some(a)) => Some(a),
3951                            (None, None) => None,
3952                        },
3953                    ))
3954                }
3955                Ok((None, Some(areavec)))
3956            }
3957            (_, _, _, _, _) => Err(
3958                "not supported area data, if parametric have the 2nd arg start and end with the { } brackets",
3959            ),
3960        }
3961    };
3962    #[allow(clippy::type_complexity)]
3963    let data: Vec<Result<(Option<Complex>, Option<Vec<Number>>), &str>> = if options.graphing {
3964        (0..points).map(body).collect()
3965    } else {
3966        #[cfg(feature = "rayon")]
3967        let r = (0..points).into_par_iter().map(body).collect();
3968        #[cfg(not(feature = "rayon"))]
3969        let r = (0..points).map(body).collect();
3970        r
3971    };
3972    for d in data {
3973        if let Ok(a) = d {
3974            if let Some(a) = a.0 {
3975                area += a
3976            } else if let Some(a) = a.1 {
3977                if areavec.is_empty() {
3978                    areavec = a
3979                } else {
3980                    for (a, b) in areavec.iter_mut().zip(a.iter()) {
3981                        a.number += b.number.clone()
3982                    }
3983                }
3984            }
3985        } else if let Err(s) = d {
3986            return Err(s);
3987        }
3988    }
3989    let g = gamma(nth.clone());
3990    if areavec.is_empty() {
3991        Ok(NumStr::new(Number::from(
3992            if negate { -area / g } else { area / g },
3993            match (units, yunits) {
3994                (Some(a), Some(b)) => Some(a.mul(&b)),
3995                (Some(a), None) | (None, Some(a)) => Some(a),
3996                (None, None) => None,
3997            },
3998        )))
3999    } else {
4000        Ok(Vector(
4001            areavec
4002                .iter()
4003                .map(|a| {
4004                    Number::from(
4005                        if negate {
4006                            -a.number.clone() / g.clone()
4007                        } else {
4008                            a.number.clone() / g.clone()
4009                        },
4010                        a.units,
4011                    )
4012                })
4013                .collect::<Vec<Number>>(),
4014        ))
4015    }
4016}
4017pub fn iter(
4018    func: Vec<NumStr>,
4019    func_vars: Vec<(String, Vec<NumStr>)>,
4020    options: Options,
4021    var: String,
4022    mut x: NumStr,
4023    n: Float,
4024    all: bool,
4025) -> Result<NumStr, &'static str> {
4026    if n.is_infinite() {
4027        let mut last = NumStr::new(Number::from(Complex::with_val(options.prec, Nan), None));
4028        let mut j = 0;
4029        if all {
4030            if let Num(num) = x.clone() {
4031                let mut vec = vec![*num];
4032                loop {
4033                    if j > 10000 {
4034                        return Ok(NumStr::new(Number::from(
4035                            Complex::with_val(options.prec, Nan),
4036                            None,
4037                        )));
4038                    }
4039                    j += 1;
4040                    last = x.clone();
4041                    x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?;
4042                    if last == x {
4043                        break;
4044                    }
4045                    vec.push(x.num()?);
4046                }
4047                Ok(Vector(vec))
4048            } else if let Vector(v) = x.clone() {
4049                let mut vec = vec![v];
4050                loop {
4051                    if j > 10000 {
4052                        return Ok(NumStr::new(Number::from(
4053                            Complex::with_val(options.prec, Nan),
4054                            None,
4055                        )));
4056                    }
4057                    j += 1;
4058                    last = x.clone();
4059                    x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?;
4060                    if last == x {
4061                        break;
4062                    }
4063                    vec.push(x.vec()?);
4064                }
4065                Ok(Matrix(vec))
4066            } else {
4067                return Err("unsupported iter");
4068            }
4069        } else {
4070            while last != x {
4071                if j > 10000 {
4072                    return Ok(NumStr::new(Number::from(
4073                        Complex::with_val(options.prec, Nan),
4074                        None,
4075                    )));
4076                }
4077                j += 1;
4078                last = x.clone();
4079                x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?
4080            }
4081            Ok(x)
4082        }
4083    } else {
4084        let n = n
4085            .to_integer()
4086            .unwrap_or_default()
4087            .to_usize()
4088            .unwrap_or_default();
4089        if all {
4090            if let Num(num) = x.clone() {
4091                let mut vec = vec![*num];
4092                for _ in 0..n {
4093                    x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?;
4094                    vec.push(x.num()?);
4095                }
4096                Ok(Vector(vec))
4097            } else if let Vector(v) = x.clone() {
4098                let mut vec = vec![v];
4099                for _ in 0..n {
4100                    x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?;
4101                    vec.push(x.vec()?);
4102                }
4103                Ok(Matrix(vec))
4104            } else {
4105                return Err("unsupported iter");
4106            }
4107        } else {
4108            for _ in 0..n {
4109                x = do_math_with_var(func.clone(), options, func_vars.clone(), &var, x)?
4110            }
4111            Ok(x)
4112        }
4113    }
4114}
4115pub fn solve(
4116    mut func: Vec<NumStr>,
4117    mut func_vars: Vec<(String, Vec<NumStr>)>,
4118    mut options: Options,
4119    var: String,
4120    x: Number,
4121) -> Result<NumStr, &'static str> {
4122    //newtons method, x-f(x)/f'(x)
4123    let units = x.units;
4124    let mut x = x.number;
4125    if x.real().is_nan() {
4126        let points = if x.imag().is_zero() {
4127            vec![
4128                Number::from(Complex::new(options.prec), None),
4129                Number::from(Complex::with_val(options.prec, -2), None),
4130                Number::from(Complex::with_val(options.prec, 2), None),
4131            ]
4132        } else {
4133            vec![
4134                Number::from(Complex::new(options.prec), None),
4135                Number::from(Complex::with_val(options.prec, -2), None),
4136                Number::from(Complex::with_val(options.prec, 2), None),
4137                Number::from(Complex::with_val(options.prec, (0, -2)), None),
4138                Number::from(Complex::with_val(options.prec, (0, 2)), None),
4139                Number::from(Complex::with_val(options.prec, (-4, 4)), None),
4140                Number::from(Complex::with_val(options.prec, (4, 4)), None),
4141                Number::from(Complex::with_val(options.prec, (-4, -4)), None),
4142                Number::from(Complex::with_val(options.prec, (4, -4)), None),
4143            ]
4144        };
4145        let mut values: Vec<Number> = Vec::new();
4146        let mut first = true;
4147        'main: for p in points {
4148            let v = solve(
4149                func.clone(),
4150                func_vars.clone(),
4151                options,
4152                var.clone(),
4153                p.clone(),
4154            )?
4155            .num()?;
4156            if !v.number.real().is_nan() {
4157                if first {
4158                    first = false;
4159                    func.insert(0, LeftBracket);
4160                    func.push(RightBracket);
4161                    func.push(Division);
4162                    func.push(LeftBracket);
4163                    func.push(LeftBracket);
4164                    func.push(Func(var.clone()));
4165                    func.push(Minus);
4166                    func.push(NumStr::new(v.clone()));
4167                    func.push(RightBracket);
4168                    func.push(RightBracket);
4169                } else {
4170                    for n1 in &values {
4171                        if -(n1.number.real() - v.number.real().clone())
4172                            .clone()
4173                            .abs()
4174                            .log2()
4175                            > options.prec / 16
4176                            && -(n1.number.imag() - v.number.imag().clone())
4177                                .clone()
4178                                .abs()
4179                                .log2()
4180                                > options.prec / 16
4181                        {
4182                            continue 'main;
4183                        }
4184                    }
4185                    func.insert(func.len() - 1, Multiplication);
4186                    func.insert(func.len() - 1, LeftBracket);
4187                    func.insert(func.len() - 1, Func(var.clone()));
4188                    func.insert(func.len() - 1, Minus);
4189                    func.insert(func.len() - 1, NumStr::new(v.clone()));
4190                    func.insert(func.len() - 1, RightBracket);
4191                }
4192                values.push(v);
4193            }
4194        }
4195        Ok(if values.is_empty() {
4196            NumStr::new(Number::from(Complex::with_val(options.prec, Nan), None))
4197        } else if values.len() == 1 {
4198            NumStr::new(values[0].clone())
4199        } else {
4200            Vector(values)
4201        })
4202    } else {
4203        let op = options.prec;
4204        let prec;
4205        (prec, options.prec) = set_slope_prec(options.prec, 1);
4206        x.set_prec(options.prec);
4207        set_prec(&mut func, &mut func_vars, options.prec);
4208        for _ in 0..(op / 4).max(64) {
4209            let n = Number::from(x.clone(), None);
4210            let y = do_math_with_var(
4211                func.clone(),
4212                options,
4213                func_vars.clone(),
4214                &var.clone(),
4215                NumStr::new(n.clone()),
4216            )?;
4217            x -= y.num()?.number
4218                / slopesided(
4219                    func.clone(),
4220                    func_vars.clone(),
4221                    options,
4222                    var.clone(),
4223                    n,
4224                    false,
4225                    1,
4226                    true,
4227                    Some(y),
4228                    None,
4229                    Some(prec),
4230                )?
4231                .num()?
4232                .number
4233        }
4234        let last = x.clone();
4235        let n = Number::from(x.clone(), None);
4236        let y = do_math_with_var(
4237            func.clone(),
4238            options,
4239            func_vars.clone(),
4240            &var.clone(),
4241            NumStr::new(n.clone()),
4242        )?;
4243        let k = slopesided(
4244            func,
4245            func_vars,
4246            options,
4247            var,
4248            n,
4249            false,
4250            1,
4251            true,
4252            Some(y.clone()),
4253            None,
4254            Some(prec),
4255        )?
4256        .num()?
4257        .number;
4258        x -= y.num()?.number / k.clone();
4259        if (last - x.clone()).abs().real().clone().log2() < op as i32 / -16 && k.real().is_finite()
4260        {
4261            x.set_prec(op);
4262            Ok(NumStr::new(Number::from(x, units)))
4263        } else {
4264            Ok(NumStr::new(Number::from(
4265                Complex::with_val(options.prec, Nan),
4266                None,
4267            )))
4268        }
4269    }
4270}
4271pub fn extrema(
4272    mut func: Vec<NumStr>,
4273    mut func_vars: Vec<(String, Vec<NumStr>)>,
4274    mut options: Options,
4275    var: String,
4276    x: Number,
4277) -> Result<NumStr, &'static str> {
4278    //newtons method, x-f'(x)/f''(x)
4279    let units = x.units;
4280    let mut x = x.number;
4281    let op = options.prec;
4282    let prec;
4283    (prec, options.prec) = set_slope_prec(options.prec, 2);
4284    x.set_prec(options.prec);
4285    set_prec(&mut func, &mut func_vars, options.prec);
4286    for _ in 0..(op / 4).max(64) {
4287        let n = Number::from(x.clone(), None);
4288        let y = do_math_with_var(
4289            func.clone(),
4290            options,
4291            func_vars.clone(),
4292            &var.clone(),
4293            NumStr::new(n.clone()),
4294        )?;
4295        let yh = do_math_with_var(
4296            func.clone(),
4297            options,
4298            func_vars.clone(),
4299            &var.clone(),
4300            NumStr::new(Number::from(
4301                x.clone() + Float::with_val(options.prec, 0.5).pow(prec),
4302                None,
4303            )),
4304        )?;
4305        x -= slopesided(
4306            func.clone(),
4307            func_vars.clone(),
4308            options,
4309            var.clone(),
4310            n.clone(),
4311            false,
4312            1,
4313            true,
4314            Some(y.clone()),
4315            Some(yh.clone()),
4316            Some(prec),
4317        )?
4318        .num()?
4319        .number
4320            / slopesided(
4321                func.clone(),
4322                func_vars.clone(),
4323                options,
4324                var.clone(),
4325                n,
4326                false,
4327                2,
4328                true,
4329                Some(y),
4330                Some(yh),
4331                Some(prec),
4332            )?
4333            .num()?
4334            .number
4335    }
4336    let last = x.clone();
4337    let n = Number::from(x.clone(), None);
4338    let y = do_math_with_var(
4339        func.clone(),
4340        options,
4341        func_vars.clone(),
4342        &var.clone(),
4343        NumStr::new(n.clone()),
4344    )?;
4345    let yh = do_math_with_var(
4346        func.clone(),
4347        options,
4348        func_vars.clone(),
4349        &var.clone(),
4350        NumStr::new(Number::from(
4351            x.clone() + Float::with_val(options.prec, 0.5).pow(prec),
4352            None,
4353        )),
4354    )?;
4355    let k = slopesided(
4356        func.clone(),
4357        func_vars.clone(),
4358        options,
4359        var.clone(),
4360        n.clone(),
4361        false,
4362        2,
4363        true,
4364        Some(y.clone()),
4365        Some(yh.clone()),
4366        Some(prec),
4367    )?
4368    .num()?
4369    .number;
4370    x -= slopesided(
4371        func.clone(),
4372        func_vars.clone(),
4373        options,
4374        var.clone(),
4375        n,
4376        false,
4377        1,
4378        true,
4379        Some(y),
4380        Some(yh),
4381        Some(prec),
4382    )?
4383    .num()?
4384    .number
4385        / k.clone();
4386    if (last - x.clone()).abs().real().clone().log2() < op as i32 / -16 && k.real().is_finite() {
4387        let n = Number::from(x, units);
4388        let mut v = Vector(vec![
4389            n.clone(),
4390            do_math_with_var(func, options, func_vars, &var, NumStr::new(n.clone()))?.num()?,
4391            Number::from(
4392                Complex::with_val(
4393                    options.prec,
4394                    if k.clone().abs().real().clone().log2() < op as i32 / -16 {
4395                        0
4396                    } else if k.real().is_sign_positive() {
4397                        1
4398                    } else {
4399                        -1
4400                    },
4401                ),
4402                None,
4403            ),
4404        ]);
4405        v.set_prec(op);
4406        Ok(v)
4407    } else {
4408        Ok(Vector(vec![Number::from(
4409            Complex::with_val(options.prec, Nan),
4410            None,
4411        )]))
4412    }
4413}
4414#[allow(clippy::too_many_arguments)]
4415pub fn taylor(
4416    mut func: Vec<NumStr>,
4417    mut func_vars: Vec<(String, Vec<NumStr>)>,
4418    mut options: Options,
4419    var: String,
4420    x: Option<Number>,
4421    a: Number,
4422    nth: usize,
4423) -> Result<NumStr, &'static str> {
4424    fn fact(n: usize) -> Integer {
4425        let mut fact = Integer::from(1);
4426        for i in 1..=n {
4427            fact *= i
4428        }
4429        fact
4430    }
4431    let op = options.prec;
4432    let mut an = a.number;
4433    let mut prec;
4434    (prec, _) = set_slope_prec(options.prec, nth.min(8) as u32);
4435    (_, options.prec) = set_slope_prec(options.prec, nth as u32);
4436    set_prec(&mut func, &mut func_vars, options.prec);
4437    an.set_prec(options.prec);
4438    let a = Number::from(an.clone(), a.units);
4439    let val = do_math_with_var(
4440        func.clone(),
4441        options,
4442        func_vars.clone(),
4443        &var,
4444        NumStr::new(a.clone()),
4445    )?;
4446    let mut val2 = do_math_with_var(
4447        func.clone(),
4448        options,
4449        func_vars.clone(),
4450        &var,
4451        NumStr::new(Number::from(
4452            an.clone() + Float::with_val(options.prec, 0.5).pow(prec),
4453            None,
4454        )),
4455    )?;
4456    if let Some(x) = x {
4457        let mut x = x.number;
4458        x.set_prec(options.prec);
4459        let mut sum = val.clone();
4460        for n in 1..=nth {
4461            if n % 8 == 0 {
4462                (prec, _) = set_slope_prec(options.prec, nth.min(8 + n) as u32);
4463                val2 = do_math_with_var(
4464                    func.clone(),
4465                    options,
4466                    func_vars.clone(),
4467                    &var,
4468                    NumStr::new(Number::from(
4469                        an.clone() + Float::with_val(options.prec, 0.5).pow(prec),
4470                        None,
4471                    )),
4472                )?;
4473            }
4474            let d = slopesided(
4475                func.clone(),
4476                func_vars.clone(),
4477                options,
4478                var.clone(),
4479                a.clone(),
4480                false,
4481                n as u32,
4482                true,
4483                Some(val.clone()),
4484                Some(val2.clone()),
4485                Some(prec),
4486            )?;
4487            let v = d.func(
4488                &NumStr::new(Number::from(
4489                    fact(n) * (x.clone() - an.clone()).pow(-(n as i32)),
4490                    None,
4491                )),
4492                div,
4493            )?;
4494            sum = sum.func(&v, add)?;
4495        }
4496        sum.set_prec(op);
4497        Ok(sum)
4498    } else {
4499        let mut poly_mat = Vec::with_capacity(nth + 1);
4500        let mut poly = Vec::with_capacity(nth + 1);
4501        let is_vector = match val.clone() {
4502            Vector(a) => {
4503                let empty = vec![Number::from(Complex::new(options.prec), None); a.len()];
4504                poly_mat.push(a);
4505                for _ in 1..=nth {
4506                    poly_mat.push(empty.clone())
4507                }
4508                true
4509            }
4510            Num(a) => {
4511                let empty = Number::from(Complex::new(options.prec), None);
4512                poly.push(*a);
4513                for _ in 1..=nth {
4514                    poly.push(empty.clone())
4515                }
4516                false
4517            }
4518            _ => return Err("unsupported type"),
4519        };
4520        for n in 1..=nth {
4521            if n % 8 == 0 {
4522                (prec, _) = set_slope_prec(options.prec, nth.min(8 + n) as u32);
4523                val2 = do_math_with_var(
4524                    func.clone(),
4525                    options,
4526                    func_vars.clone(),
4527                    &var,
4528                    NumStr::new(Number::from(
4529                        an.clone() + Float::with_val(options.prec, 0.5).pow(prec),
4530                        None,
4531                    )),
4532                )?;
4533            }
4534            let d = slopesided(
4535                func.clone(),
4536                func_vars.clone(),
4537                options,
4538                var.clone(),
4539                a.clone(),
4540                false,
4541                n as u32,
4542                true,
4543                Some(val.clone()),
4544                Some(val2.clone()),
4545                Some(prec),
4546            )?;
4547            if is_vector {
4548                let d = d.vec()?;
4549                for (i, poly) in poly_mat.iter_mut().enumerate() {
4550                    if i > n {
4551                        break;
4552                    }
4553                    for (d, poly) in d.iter().zip(poly.iter_mut()) {
4554                        *poly = add(
4555                            poly,
4556                            &Number::from(
4557                                (-a.clone().number).pow(n - i) * d.number.clone()
4558                                    / (fact(n - i) * fact(i)),
4559                                None,
4560                            ),
4561                        )
4562                    }
4563                }
4564            } else {
4565                let d = d.num()?.number;
4566                for (i, poly) in poly.iter_mut().enumerate() {
4567                    if i > n {
4568                        break;
4569                    }
4570                    *poly = add(
4571                        poly,
4572                        &Number::from(
4573                            (-a.clone().number).pow(n - i) * d.clone() / (fact(n - i) * fact(i)),
4574                            None,
4575                        ),
4576                    )
4577                }
4578            }
4579        }
4580        if is_vector {
4581            poly_mat.reverse();
4582            let mut m = Matrix(transpose(&poly_mat));
4583            m.set_prec(op);
4584            Ok(m)
4585        } else {
4586            poly.reverse();
4587            let mut v = Vector(poly);
4588            v.set_prec(op);
4589            Ok(v)
4590        }
4591    }
4592}
4593fn set_slope_prec(prec: u32, nth: u32) -> (u32, u32) {
4594    let prec = prec.clamp(256, 1024 * nth.max(1));
4595    (prec / (nth + 8), (nth / 8).max(1) * prec / 2)
4596}
4597#[allow(clippy::too_many_arguments)]
4598pub fn slope(
4599    func: Vec<NumStr>,
4600    func_vars: Vec<(String, Vec<NumStr>)>,
4601    mut options: Options,
4602    var: String,
4603    mut point: Number,
4604    combine: bool,
4605    nth: u32,
4606) -> Result<NumStr, &'static str> {
4607    let oop = options.prec;
4608    if nth == 0 {
4609        do_math_with_var(
4610            func.clone(),
4611            options,
4612            func_vars.clone(),
4613            &var,
4614            NumStr::new(point),
4615        )
4616    } else if options.prec <= 256 {
4617        options.prec = 256;
4618        point.number.set_prec(options.prec);
4619        let mut n = slopesided(
4620            func, func_vars, options, var, point, combine, nth, true, None, None, None,
4621        )?;
4622        n.set_prec(oop);
4623        Ok(n)
4624    } else {
4625        let op = options.prec.clamp(256, 1024);
4626        let prec;
4627        (prec, options.prec) = set_slope_prec(options.prec, nth);
4628        point.number.set_prec(options.prec);
4629        let val = do_math_with_var(
4630            func.clone(),
4631            options,
4632            func_vars.clone(),
4633            &var,
4634            NumStr::new(point.clone()),
4635        )?;
4636        let left = slopesided(
4637            func.clone(),
4638            func_vars.clone(),
4639            options,
4640            var.clone(),
4641            point.clone(),
4642            combine,
4643            nth,
4644            false,
4645            Some(val.clone()),
4646            None,
4647            Some(prec),
4648        )?;
4649        let right = slopesided(
4650            func,
4651            func_vars,
4652            options,
4653            var,
4654            point,
4655            combine,
4656            nth,
4657            true,
4658            Some(val),
4659            None,
4660            Some(prec),
4661        )?;
4662        match (left, right) {
4663            (Num(left), Num(right)) => {
4664                let units = left.units;
4665                let left = left.number;
4666                let right = right.number;
4667                if (((left.real().is_infinite()
4668                    && right.real().is_infinite()
4669                    && (left.imag().clone() - right.imag().clone())
4670                        .abs()
4671                        .clone()
4672                        .log2()
4673                        < op as i32 / -16)
4674                    || (left.imag().is_infinite()
4675                        && right.imag().is_infinite()
4676                        && (left.real().clone() - right.real().clone())
4677                            .abs()
4678                            .clone()
4679                            .log2()
4680                            < op as i32 / -16))
4681                    && left.real().is_sign_positive() == right.real().is_sign_positive()
4682                    && left.imag().is_sign_positive() == right.imag().is_sign_positive())
4683                    || (left.clone() - right.clone()).abs().real().clone().log2() < op as i32 / -16
4684                {
4685                    let mut n = NumStr::new(Number::from((left + right) / 2, units));
4686                    n.set_prec(oop);
4687                    Ok(n)
4688                } else {
4689                    Ok(NumStr::new(Number::from(
4690                        Complex::with_val(options.prec, Nan),
4691                        None,
4692                    )))
4693                }
4694            }
4695            (Vector(left), Vector(right)) => {
4696                let mut vec = Vec::with_capacity(left.len().max(1));
4697                for (left, right) in left.iter().zip(right) {
4698                    vec.push({
4699                        let units = left.units;
4700                        let left = left.number.clone();
4701                        let right = right.number.clone();
4702                        if (((left.real().is_infinite()
4703                            && right.real().is_infinite()
4704                            && (left.imag().clone() - right.imag().clone())
4705                                .abs()
4706                                .clone()
4707                                .log2()
4708                                < op as i32 / -16)
4709                            || (left.imag().is_infinite()
4710                                && right.imag().is_infinite()
4711                                && (left.real().clone() - right.real().clone())
4712                                    .abs()
4713                                    .clone()
4714                                    .log2()
4715                                    < op as i32 / -16))
4716                            && left.real().is_sign_positive() == right.real().is_sign_positive()
4717                            && left.imag().is_sign_positive() == right.imag().is_sign_positive())
4718                            || (left.clone() - right.clone()).abs().real().clone().log2()
4719                                < op as i32 / -16
4720                        {
4721                            Number::from((left + right) / 2, units)
4722                        } else {
4723                            Number::from(Complex::with_val(options.prec, Nan), None)
4724                        }
4725                    })
4726                }
4727                let mut v = Vector(vec);
4728                v.set_prec(oop);
4729                Ok(v)
4730            }
4731            (_, _) => Err("lim err"),
4732        }
4733    }
4734}
4735
4736#[allow(clippy::too_many_arguments)]
4737pub fn slopesided(
4738    func: Vec<NumStr>,
4739    func_vars: Vec<(String, Vec<NumStr>)>,
4740    mut options: Options,
4741    var: String,
4742    point: Number,
4743    combine: bool,
4744    nth: u32,
4745    right: bool,
4746    val: Option<NumStr>,
4747    val2: Option<NumStr>,
4748    prec: Option<u32>,
4749) -> Result<NumStr, &'static str> {
4750    if nth == 0 {
4751        return do_math_with_var(
4752            func.clone(),
4753            options,
4754            func_vars.clone(),
4755            &var,
4756            NumStr::new(point),
4757        );
4758    }
4759    let mut oop = 0;
4760    let units = point.units;
4761    let mut point = point.number;
4762    let prec = if let Some(prec) = prec {
4763        prec
4764    } else {
4765        oop = options.prec;
4766        options.prec = options.prec.clamp(256, 1024 * nth.max(1));
4767        let prec;
4768        (prec, options.prec) = set_slope_prec(options.prec, nth);
4769        point.set_prec(options.prec);
4770        prec
4771    };
4772    let h: Float = if right {
4773        Float::with_val(options.prec, 0.5).pow(prec)
4774    } else {
4775        -Float::with_val(options.prec, 0.5).pow(prec)
4776    };
4777    let num = Integer::from(nth);
4778    let n = if let Some(n) = val {
4779        n
4780    } else {
4781        do_math_with_var(
4782            func.clone(),
4783            options,
4784            func_vars.clone(),
4785            &var,
4786            NumStr::new(Number::from(point.clone(), units)),
4787        )?
4788    };
4789    match n {
4790        Num(sum) => {
4791            let yunits = sum.units;
4792            let mut sum = sum.number;
4793            if nth % 2 == 1 {
4794                sum *= -1;
4795            }
4796            let bo = val2.is_some();
4797            for k in 0..nth {
4798                if bo && nth - k == 1 {
4799                    if k % 2 == 0 {
4800                        sum += num.clone().binomial(k) * val2.clone().unwrap().num()?.number
4801                    } else {
4802                        sum -= num.clone().binomial(k) * val2.clone().unwrap().num()?.number
4803                    }
4804                } else if k % 2 == 0 {
4805                    sum += num.clone().binomial(k)
4806                        * do_math_with_var(
4807                            func.clone(),
4808                            options,
4809                            func_vars.clone(),
4810                            &var,
4811                            NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4812                        )?
4813                        .num()?
4814                        .number;
4815                } else {
4816                    sum -= num.clone().binomial(k)
4817                        * do_math_with_var(
4818                            func.clone(),
4819                            options,
4820                            func_vars.clone(),
4821                            &var,
4822                            NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4823                        )?
4824                        .num()?
4825                        .number;
4826                }
4827            }
4828            if right || nth % 2 == 0 {
4829                let mut n = NumStr::new(Number::from(
4830                    sum * Float::with_val(options.prec, 2).pow(nth * prec),
4831                    match (yunits, units) {
4832                        (Some(a), Some(b)) => Some(a.div(&b)),
4833                        (Some(a), None) => Some(a),
4834                        (None, Some(b)) => Some(Units::default().div(&b)),
4835                        (None, None) => None,
4836                    },
4837                ));
4838                if oop != 0 {
4839                    n.set_prec(oop);
4840                }
4841                Ok(n)
4842            } else {
4843                let mut n = NumStr::new(Number::from(
4844                    -sum * Float::with_val(options.prec, 2).pow(nth * prec),
4845                    match (yunits, units) {
4846                        (Some(a), Some(b)) => Some(a.div(&b)),
4847                        (Some(a), None) => Some(a),
4848                        (None, Some(b)) => Some(Units::default().div(&b)),
4849                        (None, None) => None,
4850                    },
4851                ));
4852                if oop != 0 {
4853                    n.set_prec(oop);
4854                }
4855                Ok(n)
4856            }
4857        }
4858        Vector(mut sum) if !combine => {
4859            let yunits = sum[0].units;
4860            if nth % 2 == 1 {
4861                for n in sum.iter_mut() {
4862                    n.number *= -1;
4863                }
4864            }
4865            for k in 0..nth {
4866                let b = num.clone().binomial(k);
4867                let vec = do_math_with_var(
4868                    func.clone(),
4869                    options,
4870                    func_vars.clone(),
4871                    &var,
4872                    NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4873                )?
4874                .vec()?;
4875                if k % 2 == 0 {
4876                    for (n, a) in sum.iter_mut().zip(vec) {
4877                        n.number += a.number * b.clone()
4878                    }
4879                } else {
4880                    for (n, a) in sum.iter_mut().zip(vec) {
4881                        n.number -= a.number * b.clone()
4882                    }
4883                }
4884            }
4885            let mut v = Vector(
4886                sum.iter()
4887                    .map(|n| {
4888                        Number::from(
4889                            if right || nth % 2 == 0 {
4890                                n.number.clone() * Float::with_val(options.prec, 2).pow(nth * prec)
4891                            } else {
4892                                -n.number.clone() * Float::with_val(options.prec, 2).pow(nth * prec)
4893                            },
4894                            match (yunits, units) {
4895                                (Some(a), Some(b)) => Some(a.div(&b)),
4896                                (Some(a), None) => Some(a),
4897                                (None, Some(b)) => Some(Units::default().div(&b)),
4898                                (None, None) => None,
4899                            },
4900                        )
4901                    })
4902                    .collect::<Vec<Number>>(),
4903            );
4904            if oop != 0 {
4905                v.set_prec(oop);
4906            }
4907            Ok(v)
4908        }
4909        Vector(mut sum) if sum.len() == 1 => {
4910            let yunits = sum[0].units;
4911            if nth % 2 == 1 {
4912                sum[0].number *= -1;
4913            }
4914            for k in 0..nth {
4915                if k % 2 == 0 {
4916                    sum[0].number += num.clone().binomial(k)
4917                        * do_math_with_var(
4918                            func.clone(),
4919                            options,
4920                            func_vars.clone(),
4921                            &var,
4922                            NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4923                        )?
4924                        .num()?
4925                        .number;
4926                } else {
4927                    sum[0].number -= num.clone().binomial(k)
4928                        * do_math_with_var(
4929                            func.clone(),
4930                            options,
4931                            func_vars.clone(),
4932                            &var,
4933                            NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4934                        )?
4935                        .num()?
4936                        .number;
4937                }
4938            }
4939            if right || nth % 2 == 0 {
4940                let mut n = NumStr::new(Number::from(
4941                    sum[0].number.clone() * Float::with_val(options.prec, 2).pow(nth * prec),
4942                    match (yunits, units) {
4943                        (Some(a), Some(b)) => Some(a.div(&b)),
4944                        (Some(a), None) => Some(a),
4945                        (None, Some(b)) => Some(Units::default().div(&b)),
4946                        (None, None) => None,
4947                    },
4948                ));
4949                if oop != 0 {
4950                    n.set_prec(oop)
4951                }
4952                Ok(n)
4953            } else {
4954                let mut n = NumStr::new(Number::from(
4955                    -sum[0].number.clone() * Float::with_val(options.prec, 2).pow(nth * prec),
4956                    match (yunits, units) {
4957                        (Some(a), Some(b)) => Some(a.div(&b)),
4958                        (Some(a), None) => Some(a),
4959                        (None, Some(b)) => Some(Units::default().div(&b)),
4960                        (None, None) => None,
4961                    },
4962                ));
4963                if oop != 0 {
4964                    n.set_prec(oop)
4965                }
4966                Ok(n)
4967            }
4968        }
4969        Vector(mut sum) => {
4970            let yunits = sum[0].units;
4971            if nth % 2 == 1 {
4972                for n in sum.iter_mut() {
4973                    n.number *= -1;
4974                }
4975            }
4976            for k in 0..nth {
4977                let b = num.clone().binomial(k);
4978                let vec = do_math_with_var(
4979                    func.clone(),
4980                    options,
4981                    func_vars.clone(),
4982                    &var,
4983                    NumStr::new(Number::from(point.clone() + h.clone() * (nth - k), units)),
4984                )?
4985                .vec()?;
4986                if k % 2 == 0 {
4987                    for (n, a) in sum.iter_mut().zip(vec) {
4988                        n.number += a.number * b.clone()
4989                    }
4990                } else {
4991                    for (n, a) in sum.iter_mut().zip(vec) {
4992                        n.number -= a.number * b.clone()
4993                    }
4994                }
4995            }
4996            if sum.len() == 2 {
4997                let mut n = NumStr::new(Number::from(
4998                    sum[1].number.clone() / sum[0].number.clone(),
4999                    match (yunits, units) {
5000                        (Some(a), Some(b)) => Some(a.div(&b)),
5001                        (Some(a), None) => Some(a),
5002                        (None, Some(b)) => Some(Units::default().div(&b)),
5003                        (None, None) => None,
5004                    },
5005                ));
5006                if oop != 0 {
5007                    n.set_prec(oop)
5008                }
5009                Ok(n)
5010            } else {
5011                let nf = &sum.last().unwrap().number;
5012                let mut v = Vector(
5013                    sum[0..sum.len().saturating_sub(1)]
5014                        .iter()
5015                        .map(|n| {
5016                            Number::from(
5017                                nf.clone() / n.number.clone(),
5018                                match (yunits, units) {
5019                                    (Some(a), Some(b)) => Some(a.div(&b)),
5020                                    (Some(a), None) => Some(a),
5021                                    (None, Some(b)) => Some(Units::default().div(&b)),
5022                                    (None, None) => None,
5023                                },
5024                            )
5025                        })
5026                        .collect::<Vec<Number>>(),
5027                );
5028                if oop != 0 {
5029                    v.set_prec(oop)
5030                }
5031                Ok(v)
5032            }
5033        }
5034        _ => Err("not supported slope data"),
5035    }
5036}
5037#[derive(Copy, Clone, PartialEq)]
5038pub enum LimSide {
5039    Left,
5040    Right,
5041    Both,
5042}
5043pub fn limit(
5044    mut func: Vec<NumStr>,
5045    mut func_vars: Vec<(String, Vec<NumStr>)>,
5046    mut options: Options,
5047    var: String,
5048    point: Number,
5049    side: LimSide,
5050) -> Result<NumStr, &'static str> {
5051    let xunits = point.units;
5052    let mut point = point.number;
5053    let oop = options.prec;
5054    options.prec = options.prec.clamp(256, 1024);
5055    if oop != options.prec {
5056        point.set_prec(options.prec);
5057        set_prec(&mut func, &mut func_vars, options.prec);
5058    }
5059    if point.clone().real().is_infinite() || point.clone().imag().is_infinite() {
5060        let (h1, h2);
5061        let positive = point.real().is_sign_positive();
5062        if positive {
5063            h1 = Complex::with_val(options.prec, 2).pow(options.prec / 4);
5064            h2 = Complex::with_val(options.prec, 2).pow((options.prec / 3) as f64 + 7.0 / 0.94) - 3;
5065        } else {
5066            h1 = -Complex::with_val(options.prec, 2).pow(options.prec / 4);
5067            h2 = 3 - Complex::with_val(options.prec, 2).pow((options.prec / 3) as f64 + 7.0 / 0.94);
5068        }
5069        let n1 = do_math_with_var(
5070            func.clone(),
5071            options,
5072            func_vars.clone(),
5073            &var,
5074            NumStr::new(Number::from(h1, xunits)),
5075        )?;
5076        let n2 = do_math_with_var(
5077            func.clone(),
5078            options,
5079            func_vars.clone(),
5080            &var,
5081            NumStr::new(Number::from(h2, xunits)),
5082        )?;
5083        match (n1, n2) {
5084            (Num(n1), Num(n2)) => {
5085                let units = n1.units;
5086                let n1 = n1.number;
5087                let n2 = n2.number;
5088                let mut n = if (n1.clone() - n2.clone()).abs().real().clone().log2()
5089                    < options.prec as i32 / -16
5090                {
5091                    NumStr::new(Number::from(n2, units))
5092                } else if n1.real().is_sign_positive() != n2.real().is_sign_positive()
5093                    || n1.imag().is_sign_positive() != n2.imag().is_sign_positive()
5094                {
5095                    NumStr::new(Number::from(Complex::with_val(options.prec, Nan), None))
5096                } else if n2.real().is_infinite() || n2.imag().is_infinite() {
5097                    NumStr::new(Number::from(
5098                        match (n2.real().is_infinite(), n2.imag().is_infinite()) {
5099                            (true, true) => {
5100                                match (n1.real().is_sign_positive(), n1.imag().is_sign_positive()) {
5101                                    (true, true) => {
5102                                        Complex::with_val(options.prec, (Infinity, Infinity))
5103                                    }
5104                                    (true, false) => Complex::with_val(
5105                                        options.prec,
5106                                        (Infinity, -Float::with_val(options.prec, Infinity)),
5107                                    ),
5108                                    (false, true) => Complex::with_val(
5109                                        options.prec,
5110                                        (-Float::with_val(options.prec, Infinity), Infinity),
5111                                    ),
5112                                    (false, false) => {
5113                                        -Complex::with_val(options.prec, (Infinity, Infinity))
5114                                    }
5115                                }
5116                            }
5117                            (true, false) => {
5118                                if n1.real().is_sign_positive() {
5119                                    Complex::with_val(
5120                                        options.prec,
5121                                        (
5122                                            Infinity,
5123                                            if (n1.imag() - n2.imag().clone()).abs().log2()
5124                                                < options.prec as i32 / -16
5125                                            {
5126                                                n2.imag().clone()
5127                                            } else {
5128                                                Float::new(options.prec)
5129                                            },
5130                                        ),
5131                                    )
5132                                } else {
5133                                    -Complex::with_val(
5134                                        options.prec,
5135                                        (
5136                                            Infinity,
5137                                            if (n1.imag() - n2.imag().clone()).abs().log2()
5138                                                < options.prec as i32 / -16
5139                                            {
5140                                                -n2.imag().clone()
5141                                            } else {
5142                                                Float::new(options.prec)
5143                                            },
5144                                        ),
5145                                    )
5146                                }
5147                            }
5148                            (false, true) => {
5149                                if n1.imag().is_sign_positive() {
5150                                    Complex::with_val(
5151                                        options.prec,
5152                                        (
5153                                            if (n1.real() - n2.real().clone()).abs().log2()
5154                                                < options.prec as i32 / -16
5155                                            {
5156                                                n2.real().clone()
5157                                            } else {
5158                                                Float::new(options.prec)
5159                                            },
5160                                            Infinity,
5161                                        ),
5162                                    )
5163                                } else {
5164                                    -Complex::with_val(
5165                                        options.prec,
5166                                        (
5167                                            if (n1.real() - n2.real().clone()).abs().log2()
5168                                                < options.prec as i32 / -16
5169                                            {
5170                                                -n2.real().clone()
5171                                            } else {
5172                                                Float::new(options.prec)
5173                                            },
5174                                            Infinity,
5175                                        ),
5176                                    )
5177                                }
5178                            }
5179                            (false, false) => Complex::with_val(options.prec, Nan),
5180                        },
5181                        units,
5182                    ))
5183                } else {
5184                    let n3 = do_math_with_var(
5185                        func.clone(),
5186                        options,
5187                        func_vars.clone(),
5188                        &var,
5189                        NumStr::new(Number::from(
5190                            if positive {
5191                                Complex::with_val(options.prec, 2)
5192                                    .pow((options.prec / 2) as f64 + 13.0 / 0.7)
5193                                    - 7
5194                            } else {
5195                                7 - Complex::with_val(options.prec, 2)
5196                                    .pow((options.prec / 2) as f64 + 13.0 / 0.7)
5197                            },
5198                            xunits,
5199                        )),
5200                    )?
5201                    .num()?
5202                    .number;
5203                    let sign = n2.real().is_sign_positive() == n3.real().is_sign_positive()
5204                        && n2.imag().is_sign_positive() == n3.imag().is_sign_positive();
5205                    let n1r = n1.real().clone().abs();
5206                    let n2r = n2.real().clone().abs();
5207                    let n3r = n3.real().clone().abs();
5208                    let n1i = n1.imag().clone().abs();
5209                    let n2i = n2.imag().clone().abs();
5210                    let n3i = n3.imag().clone().abs();
5211                    NumStr::new(Number::from(
5212                        if !sign {
5213                            Complex::with_val(options.prec, Nan)
5214                        } else {
5215                            match (n3r > n2r && n2r > n1r, n3i > n2i && n2i > n1i) {
5216                                (true, true) => {
5217                                    match (
5218                                        n1.real().is_sign_positive(),
5219                                        n1.imag().is_sign_positive(),
5220                                    ) {
5221                                        (true, true) => {
5222                                            Complex::with_val(options.prec, (Infinity, Infinity))
5223                                        }
5224                                        (true, false) => Complex::with_val(
5225                                            options.prec,
5226                                            (Infinity, -Float::with_val(options.prec, Infinity)),
5227                                        ),
5228                                        (false, true) => Complex::with_val(
5229                                            options.prec,
5230                                            (-Float::with_val(options.prec, Infinity), Infinity),
5231                                        ),
5232                                        (false, false) => {
5233                                            -Complex::with_val(options.prec, (Infinity, Infinity))
5234                                        }
5235                                    }
5236                                }
5237                                (true, false) => {
5238                                    if n1.real().is_sign_positive() {
5239                                        Complex::with_val(
5240                                            options.prec,
5241                                            (
5242                                                Infinity,
5243                                                if (n2i - n3i.clone()).abs().log2()
5244                                                    < options.prec as i32 / -16
5245                                                {
5246                                                    n3i
5247                                                } else {
5248                                                    Float::new(options.prec)
5249                                                },
5250                                            ),
5251                                        )
5252                                    } else {
5253                                        -Complex::with_val(
5254                                            options.prec,
5255                                            (
5256                                                Infinity,
5257                                                if (n2i - n3i.clone()).abs().log2()
5258                                                    < options.prec as i32 / -16
5259                                                {
5260                                                    -n3i
5261                                                } else {
5262                                                    Float::new(options.prec)
5263                                                },
5264                                            ),
5265                                        )
5266                                    }
5267                                }
5268                                (false, true) => {
5269                                    if n1.imag().is_sign_positive() {
5270                                        Complex::with_val(
5271                                            options.prec,
5272                                            (
5273                                                if (n2r - n3r.clone()).abs().log2()
5274                                                    < options.prec as i32 / -16
5275                                                {
5276                                                    n3r
5277                                                } else {
5278                                                    Float::new(options.prec)
5279                                                },
5280                                                Infinity,
5281                                            ),
5282                                        )
5283                                    } else {
5284                                        -Complex::with_val(
5285                                            options.prec,
5286                                            (
5287                                                if (n2r - n3r.clone()).abs().log2()
5288                                                    < options.prec as i32 / -16
5289                                                {
5290                                                    -n3r
5291                                                } else {
5292                                                    Float::new(options.prec)
5293                                                },
5294                                                Infinity,
5295                                            ),
5296                                        )
5297                                    }
5298                                }
5299                                (false, false) => Complex::with_val(options.prec, Nan),
5300                            }
5301                        },
5302                        units,
5303                    ))
5304                };
5305                n.set_prec(oop);
5306                Ok(n)
5307            }
5308            (Vector(v1), Vector(v2)) => {
5309                let mut v3: Vec<Number> = Vec::new();
5310                let mut vec = Vec::with_capacity(v1.len().max(1));
5311                for (i, (n1, n2)) in v1.iter().zip(v2).enumerate() {
5312                    let units = n1.units;
5313                    let n1 = n1.number.clone();
5314                    let n2 = n2.number.clone();
5315                    vec.push(Number::from(
5316                        if (n1.clone() - n2.clone()).abs().real().clone().log2()
5317                            < options.prec as i32 / -16
5318                        {
5319                            n2
5320                        } else if n1.real().is_sign_positive() != n2.real().is_sign_positive()
5321                            || n1.imag().is_sign_positive() != n2.imag().is_sign_positive()
5322                        {
5323                            Complex::with_val(options.prec, Nan)
5324                        } else if n2.real().is_infinite() || n2.imag().is_infinite() {
5325                            match (n2.real().is_infinite(), n2.imag().is_infinite()) {
5326                                (true, true) => {
5327                                    match (
5328                                        n1.real().is_sign_positive(),
5329                                        n1.imag().is_sign_positive(),
5330                                    ) {
5331                                        (true, true) => {
5332                                            Complex::with_val(options.prec, (Infinity, Infinity))
5333                                        }
5334                                        (true, false) => Complex::with_val(
5335                                            options.prec,
5336                                            (Infinity, -Float::with_val(options.prec, Infinity)),
5337                                        ),
5338                                        (false, true) => Complex::with_val(
5339                                            options.prec,
5340                                            (-Float::with_val(options.prec, Infinity), Infinity),
5341                                        ),
5342                                        (false, false) => {
5343                                            -Complex::with_val(options.prec, (Infinity, Infinity))
5344                                        }
5345                                    }
5346                                }
5347                                (true, false) => {
5348                                    if n1.real().is_sign_positive() {
5349                                        Complex::with_val(
5350                                            options.prec,
5351                                            (
5352                                                Infinity,
5353                                                if (n1.imag() - n2.imag().clone()).abs().log10()
5354                                                    <= -10
5355                                                {
5356                                                    n2.imag().clone()
5357                                                } else {
5358                                                    Float::new(options.prec)
5359                                                },
5360                                            ),
5361                                        )
5362                                    } else {
5363                                        -Complex::with_val(
5364                                            options.prec,
5365                                            (
5366                                                Infinity,
5367                                                if (n1.imag() - n2.imag().clone()).abs().log10()
5368                                                    <= -10
5369                                                {
5370                                                    -n2.imag().clone()
5371                                                } else {
5372                                                    Float::new(options.prec)
5373                                                },
5374                                            ),
5375                                        )
5376                                    }
5377                                }
5378                                (false, true) => {
5379                                    if n1.imag().is_sign_positive() {
5380                                        Complex::with_val(
5381                                            options.prec,
5382                                            (
5383                                                if (n1.real() - n2.real().clone()).abs().log10()
5384                                                    <= -10
5385                                                {
5386                                                    n2.real().clone()
5387                                                } else {
5388                                                    Float::new(options.prec)
5389                                                },
5390                                                Infinity,
5391                                            ),
5392                                        )
5393                                    } else {
5394                                        -Complex::with_val(
5395                                            options.prec,
5396                                            (
5397                                                if (n1.real() - n2.real().clone()).abs().log10()
5398                                                    <= -10
5399                                                {
5400                                                    -n2.real().clone()
5401                                                } else {
5402                                                    Float::new(options.prec)
5403                                                },
5404                                                Infinity,
5405                                            ),
5406                                        )
5407                                    }
5408                                }
5409                                (false, false) => Complex::with_val(options.prec, Nan),
5410                            }
5411                        } else {
5412                            if v3.is_empty() {
5413                                v3 = do_math_with_var(
5414                                    func.clone(),
5415                                    options,
5416                                    func_vars.clone(),
5417                                    &var,
5418                                    NumStr::new(Number::from(
5419                                        if positive {
5420                                            Complex::with_val(options.prec, 2)
5421                                                .pow((options.prec / 2) as f64 + 13.0 / 0.7)
5422                                                - 7
5423                                        } else {
5424                                            7 - Complex::with_val(options.prec, 2)
5425                                                .pow((options.prec / 2) as f64 + 13.0 / 0.7)
5426                                        },
5427                                        xunits,
5428                                    )),
5429                                )?
5430                                .vec()?;
5431                            }
5432                            let v3 = &v3[i].number;
5433                            let sign = n2.real().is_sign_positive() == v3.real().is_sign_positive()
5434                                && n2.imag().is_sign_positive() == v3.imag().is_sign_positive();
5435                            let n1r = n1.real().clone().abs();
5436                            let n2r = n2.real().clone().abs();
5437                            let n3r = v3.real().clone().abs();
5438                            let n1i = n1.imag().clone().abs();
5439                            let n2i = n2.imag().clone().abs();
5440                            let n3i = v3.imag().clone().abs();
5441                            if !sign {
5442                                Complex::with_val(options.prec, Nan)
5443                            } else {
5444                                match (n3r > n2r && n2r > n1r, n3i > n2i && n2i > n1i) {
5445                                    (true, true) => {
5446                                        match (
5447                                            n1.real().is_sign_positive(),
5448                                            n1.imag().is_sign_positive(),
5449                                        ) {
5450                                            (true, true) => Complex::with_val(
5451                                                options.prec,
5452                                                (Infinity, Infinity),
5453                                            ),
5454                                            (true, false) => Complex::with_val(
5455                                                options.prec,
5456                                                (
5457                                                    Infinity,
5458                                                    -Float::with_val(options.prec, Infinity),
5459                                                ),
5460                                            ),
5461                                            (false, true) => Complex::with_val(
5462                                                options.prec,
5463                                                (
5464                                                    -Float::with_val(options.prec, Infinity),
5465                                                    Infinity,
5466                                                ),
5467                                            ),
5468                                            (false, false) => -Complex::with_val(
5469                                                options.prec,
5470                                                (Infinity, Infinity),
5471                                            ),
5472                                        }
5473                                    }
5474                                    (true, false) => {
5475                                        if n1.real().is_sign_positive() {
5476                                            Complex::with_val(
5477                                                options.prec,
5478                                                (
5479                                                    Infinity,
5480                                                    if (n2i - n3i.clone()).abs().log2()
5481                                                        < options.prec as i32 / -16
5482                                                    {
5483                                                        n3i
5484                                                    } else {
5485                                                        Float::new(options.prec)
5486                                                    },
5487                                                ),
5488                                            )
5489                                        } else {
5490                                            -Complex::with_val(
5491                                                options.prec,
5492                                                (
5493                                                    Infinity,
5494                                                    if (n2i - n3i.clone()).abs().log2()
5495                                                        < options.prec as i32 / -16
5496                                                    {
5497                                                        -n3i
5498                                                    } else {
5499                                                        Float::new(options.prec)
5500                                                    },
5501                                                ),
5502                                            )
5503                                        }
5504                                    }
5505                                    (false, true) => {
5506                                        if n1.imag().is_sign_positive() {
5507                                            Complex::with_val(
5508                                                options.prec,
5509                                                (
5510                                                    if (n2r - n3r.clone()).abs().log2()
5511                                                        < options.prec as i32 / -16
5512                                                    {
5513                                                        n3r
5514                                                    } else {
5515                                                        Float::new(options.prec)
5516                                                    },
5517                                                    Infinity,
5518                                                ),
5519                                            )
5520                                        } else {
5521                                            -Complex::with_val(
5522                                                options.prec,
5523                                                (
5524                                                    if (n2r - n3r.clone()).abs().log2()
5525                                                        < options.prec as i32 / -16
5526                                                    {
5527                                                        -n3r
5528                                                    } else {
5529                                                        Float::new(options.prec)
5530                                                    },
5531                                                    Infinity,
5532                                                ),
5533                                            )
5534                                        }
5535                                    }
5536                                    (false, false) => Complex::with_val(options.prec, Nan),
5537                                }
5538                            }
5539                        },
5540                        units,
5541                    ))
5542                }
5543                let mut v = Vector(vec);
5544                v.set_prec(oop);
5545                Ok(v)
5546            }
5547            (_, _) => Err("unsupported lim data"),
5548        }
5549    } else {
5550        let point = Number::from(point, xunits);
5551        match side {
5552            LimSide::Left => limsided(func, func_vars, options, var, point, false),
5553            LimSide::Right => limsided(func, func_vars, options, var, point, true),
5554            LimSide::Both => {
5555                let left = limsided(
5556                    func.clone(),
5557                    func_vars.clone(),
5558                    options,
5559                    var.clone(),
5560                    point.clone(),
5561                    false,
5562                )?;
5563                let right = limsided(func, func_vars, options, var, point, true)?;
5564                match (left, right) {
5565                    (Num(left), Num(right)) => {
5566                        let units = left.units;
5567                        let left = left.number;
5568                        let right = right.number;
5569                        if (((left.real().is_infinite()
5570                            && right.real().is_infinite()
5571                            && (left.imag().clone() - right.imag().clone())
5572                                .abs()
5573                                .clone()
5574                                .log2()
5575                                < options.prec as i32 / -16)
5576                            || (left.imag().is_infinite()
5577                                && right.imag().is_infinite()
5578                                && (left.real().clone() - right.real().clone())
5579                                    .abs()
5580                                    .clone()
5581                                    .log2()
5582                                    < options.prec as i32 / -16))
5583                            && left.real().is_sign_positive() == right.real().is_sign_positive()
5584                            && left.imag().is_sign_positive() == right.imag().is_sign_positive())
5585                            || (left.clone() - right.clone()).abs().real().clone().log2()
5586                                < options.prec as i32 / -16
5587                        {
5588                            let mut n = NumStr::new(Number::from((left + right) / 2, units));
5589                            n.set_prec(oop);
5590                            Ok(n)
5591                        } else {
5592                            Ok(NumStr::new(Number::from(
5593                                Complex::with_val(options.prec, Nan),
5594                                None,
5595                            )))
5596                        }
5597                    }
5598                    (Vector(left), Vector(right)) => {
5599                        let mut vec = Vec::with_capacity(left.len().max(1));
5600                        for (left, right) in left.iter().zip(right) {
5601                            let units = left.units;
5602                            let left = &left.number;
5603                            let right = &right.number;
5604                            vec.push(
5605                                if (((left.real().is_infinite()
5606                                    && right.real().is_infinite()
5607                                    && (left.imag().clone() - right.imag().clone())
5608                                        .abs()
5609                                        .clone()
5610                                        .log2()
5611                                        < options.prec as i32 / -16)
5612                                    || (left.imag().is_infinite()
5613                                        && right.imag().is_infinite()
5614                                        && (left.real().clone() - right.real().clone())
5615                                            .abs()
5616                                            .clone()
5617                                            .log2()
5618                                            < options.prec as i32 / -16))
5619                                    && left.real().is_sign_positive()
5620                                        == right.real().is_sign_positive()
5621                                    && left.imag().is_sign_positive()
5622                                        == right.imag().is_sign_positive())
5623                                    || (left.clone() - right.clone()).abs().real().clone().log2()
5624                                        < options.prec as i32 / -16
5625                                {
5626                                    Number::from((left + right.clone()) / 2, units)
5627                                } else {
5628                                    Number::from(Complex::with_val(options.prec, Nan), None)
5629                                },
5630                            )
5631                        }
5632                        let mut v = Vector(vec);
5633                        v.set_prec(options.prec);
5634                        Ok(v)
5635                    }
5636                    (_, _) => Err("lim err"),
5637                }
5638            }
5639        }
5640    }
5641}
5642fn limsided(
5643    func: Vec<NumStr>,
5644    func_vars: Vec<(String, Vec<NumStr>)>,
5645    options: Options,
5646    var: String,
5647    point: Number,
5648    right: bool,
5649) -> Result<NumStr, &'static str> {
5650    let xunits = point.units;
5651    let point = point.number;
5652    let h1 = Float::with_val(options.prec, 0.5).pow(options.prec / 4);
5653    let h2 = Float::with_val(options.prec, 0.5).pow((options.prec / 3) as f64 + 7.0 / 0.94);
5654    let n1 = do_math_with_var(
5655        func.clone(),
5656        options,
5657        func_vars.clone(),
5658        &var,
5659        NumStr::new(Number::from(
5660            point.clone() + if right { h1 } else { -h1 },
5661            xunits,
5662        )),
5663    )?;
5664    let n2 = do_math_with_var(
5665        func.clone(),
5666        options,
5667        func_vars.clone(),
5668        &var,
5669        NumStr::new(Number::from(
5670            point.clone() + if right { h2 } else { -h2 },
5671            xunits,
5672        )),
5673    )?;
5674    match (n1, n2) {
5675        (Num(n1), Num(n2)) => {
5676            let units = n1.units;
5677            let n1 = n1.number;
5678            let n2 = n2.number;
5679            Ok(NumStr::new(Number::from(
5680                if (n2.clone() - n1.clone()).abs().real().clone().log2() < options.prec as i32 / -16
5681                {
5682                    n1
5683                } else {
5684                    let h3 = Complex::with_val(options.prec, 0.5)
5685                        .pow((options.prec / 2) as f64 + 13.0 / 0.7);
5686                    let n3 = do_math_with_var(
5687                        func.clone(),
5688                        options,
5689                        func_vars.clone(),
5690                        &var,
5691                        NumStr::new(Number::from(
5692                            point.clone() + if right { h3 } else { -h3 },
5693                            xunits,
5694                        )),
5695                    )?
5696                    .num()?
5697                    .number;
5698                    let sign = n1.real().is_sign_positive() == n2.real().is_sign_positive()
5699                        && n2.real().is_sign_positive() == n3.real().is_sign_positive()
5700                        && n1.imag().is_sign_positive() == n2.imag().is_sign_positive()
5701                        && n2.imag().is_sign_positive() == n3.imag().is_sign_positive();
5702                    let n1r = n1.real().clone().abs();
5703                    let n2r = n2.real().clone().abs();
5704                    let n3r = n3.real().clone().abs();
5705                    let n1i = n1.imag().clone().abs();
5706                    let n2i = n2.imag().clone().abs();
5707                    let n3i = n3.imag().clone().abs();
5708                    if !sign {
5709                        Complex::with_val(options.prec, Nan)
5710                    } else {
5711                        match (n3r > n2r && n2r > n1r, n3i > n2i && n2i > n1i) {
5712                            (true, true) => {
5713                                match (n1.real().is_sign_positive(), n1.imag().is_sign_positive()) {
5714                                    (true, true) => {
5715                                        Complex::with_val(options.prec, (Infinity, Infinity))
5716                                    }
5717                                    (true, false) => Complex::with_val(
5718                                        options.prec,
5719                                        (Infinity, -Float::with_val(options.prec, Infinity)),
5720                                    ),
5721                                    (false, true) => Complex::with_val(
5722                                        options.prec,
5723                                        (-Float::with_val(options.prec, Infinity), Infinity),
5724                                    ),
5725                                    (false, false) => {
5726                                        -Complex::with_val(options.prec, (Infinity, Infinity))
5727                                    }
5728                                }
5729                            }
5730                            (true, false) => {
5731                                if n1.real().is_sign_positive() {
5732                                    Complex::with_val(
5733                                        options.prec,
5734                                        (
5735                                            Infinity,
5736                                            if (n1.imag() - n2.imag().clone()).abs().log2()
5737                                                < options.prec as i32 / -16
5738                                            {
5739                                                n2.imag().clone()
5740                                            } else {
5741                                                Float::new(options.prec)
5742                                            },
5743                                        ),
5744                                    )
5745                                } else {
5746                                    -Complex::with_val(
5747                                        options.prec,
5748                                        (
5749                                            Infinity,
5750                                            if (n1.imag() - n2.imag().clone()).abs().log2()
5751                                                < options.prec as i32 / -16
5752                                            {
5753                                                -n2.imag().clone()
5754                                            } else {
5755                                                Float::new(options.prec)
5756                                            },
5757                                        ),
5758                                    )
5759                                }
5760                            }
5761                            (false, true) => {
5762                                if n1.imag().is_sign_positive() {
5763                                    Complex::with_val(
5764                                        options.prec,
5765                                        (
5766                                            if (n1.real() - n2.real().clone()).abs().log2()
5767                                                < options.prec as i32 / -16
5768                                            {
5769                                                n2.real().clone()
5770                                            } else {
5771                                                Float::new(options.prec)
5772                                            },
5773                                            Infinity,
5774                                        ),
5775                                    )
5776                                } else {
5777                                    -Complex::with_val(
5778                                        options.prec,
5779                                        (
5780                                            if (n1.real() - n2.real().clone()).abs().log2()
5781                                                < options.prec as i32 / -16
5782                                            {
5783                                                -n2.real().clone()
5784                                            } else {
5785                                                Float::new(options.prec)
5786                                            },
5787                                            Infinity,
5788                                        ),
5789                                    )
5790                                }
5791                            }
5792                            (false, false) => Complex::with_val(options.prec, Nan),
5793                        }
5794                    }
5795                },
5796                units,
5797            )))
5798        }
5799        (Vector(n1), Vector(n2)) => {
5800            let mut n3: Vec<Number> = Vec::new();
5801            let mut vec = Vec::with_capacity(n1.len().max(1));
5802            for (i, (n1, n2)) in n1.iter().zip(n2).enumerate() {
5803                let units = n1.units;
5804                let n1 = &n1.number;
5805                let n2 = &n2.number;
5806                vec.push(Number::from(
5807                    if (n2.clone() - n1.clone()).abs().real().clone().log2()
5808                        < options.prec as i32 / -16
5809                    {
5810                        n1.clone()
5811                    } else {
5812                        if n3.is_empty() {
5813                            let h3 = Complex::with_val(options.prec, 0.5)
5814                                .pow((options.prec / 2) as f64 + 13.0 / 0.7);
5815                            n3 = do_math_with_var(
5816                                func.clone(),
5817                                options,
5818                                func_vars.clone(),
5819                                &var,
5820                                NumStr::new(Number::from(point.clone() - h3.clone(), xunits)),
5821                            )?
5822                            .vec()?;
5823                        }
5824                        let n3 = &n3[i].number;
5825                        let sign = n1.real().is_sign_positive() == n2.real().is_sign_positive()
5826                            && n2.real().is_sign_positive() == n3.real().is_sign_positive()
5827                            && n1.imag().is_sign_positive() == n2.imag().is_sign_positive()
5828                            && n2.imag().is_sign_positive() == n3.imag().is_sign_positive();
5829                        let n1r = n1.real().clone().abs();
5830                        let n2r = n2.real().clone().abs();
5831                        let n3r = n3.real().clone().abs();
5832                        let n1i = n1.imag().clone().abs();
5833                        let n2i = n2.imag().clone().abs();
5834                        let n3i = n3.imag().clone().abs();
5835                        if !sign {
5836                            Complex::with_val(options.prec, Nan)
5837                        } else {
5838                            match (n3r > n2r && n2r > n1r, n3i > n2i && n2i > n1i) {
5839                                (true, true) => {
5840                                    match (
5841                                        n1.real().is_sign_positive(),
5842                                        n1.imag().is_sign_positive(),
5843                                    ) {
5844                                        (true, true) => {
5845                                            Complex::with_val(options.prec, (Infinity, Infinity))
5846                                        }
5847                                        (true, false) => Complex::with_val(
5848                                            options.prec,
5849                                            (Infinity, -Float::with_val(options.prec, Infinity)),
5850                                        ),
5851                                        (false, true) => Complex::with_val(
5852                                            options.prec,
5853                                            (-Float::with_val(options.prec, Infinity), Infinity),
5854                                        ),
5855                                        (false, false) => {
5856                                            -Complex::with_val(options.prec, (Infinity, Infinity))
5857                                        }
5858                                    }
5859                                }
5860                                (true, false) => {
5861                                    if n1.real().is_sign_positive() {
5862                                        Complex::with_val(
5863                                            options.prec,
5864                                            (
5865                                                Infinity,
5866                                                if (n2i - n3i.clone()).abs().log2()
5867                                                    < options.prec as i32 / -16
5868                                                {
5869                                                    n3i
5870                                                } else {
5871                                                    Float::new(options.prec)
5872                                                },
5873                                            ),
5874                                        )
5875                                    } else {
5876                                        -Complex::with_val(
5877                                            options.prec,
5878                                            (
5879                                                Infinity,
5880                                                if (n2i - n3i.clone()).abs().log2()
5881                                                    < options.prec as i32 / -16
5882                                                {
5883                                                    -n3i
5884                                                } else {
5885                                                    Float::new(options.prec)
5886                                                },
5887                                            ),
5888                                        )
5889                                    }
5890                                }
5891                                (false, true) => {
5892                                    if n1.imag().is_sign_positive() {
5893                                        Complex::with_val(
5894                                            options.prec,
5895                                            (
5896                                                if (n2r - n3r.clone()).abs().log2()
5897                                                    < options.prec as i32 / -16
5898                                                {
5899                                                    n3r
5900                                                } else {
5901                                                    Float::new(options.prec)
5902                                                },
5903                                                Infinity,
5904                                            ),
5905                                        )
5906                                    } else {
5907                                        -Complex::with_val(
5908                                            options.prec,
5909                                            (
5910                                                if (n2r - n3r.clone()).abs().log2()
5911                                                    < options.prec as i32 / -16
5912                                                {
5913                                                    -n3r
5914                                                } else {
5915                                                    Float::new(options.prec)
5916                                                },
5917                                                Infinity,
5918                                            ),
5919                                        )
5920                                    }
5921                                }
5922                                (false, false) => Complex::with_val(options.prec, Nan),
5923                            }
5924                        }
5925                    },
5926                    units,
5927                ))
5928            }
5929            Ok(Vector(vec))
5930        }
5931        (_, _) => Err("unsupported lim data"),
5932    }
5933}
5934//https://github.com/IstvanMezo/LambertW-function/blob/master/complex%20Lambert.cpp
5935pub fn lambertw(z: Complex, k: Integer) -> Complex {
5936    if z.is_zero() {
5937        return if k == 0 {
5938            Complex::new(z.prec())
5939        } else {
5940            -Complex::with_val(z.prec(), Infinity)
5941        };
5942    }
5943    if z.imag().is_zero() && (k == 0 || k == -1) {
5944        let e: Float = -Float::with_val(z.prec().0, -1).exp();
5945        if z.real() == &e {
5946            return Complex::with_val(z.prec(), -1);
5947        }
5948    }
5949    let mut w = initpoint(z.clone(), k);
5950    for _ in 0..(z.prec().0 / 64).max(8) {
5951        let zexp = w.clone().exp();
5952        let zexpz = &w * zexp.clone();
5953        let zexpz_d = &zexp + zexpz.clone();
5954        let zexpz_dd = (2 * zexp) + &zexpz;
5955        w -=
5956            2 * ((zexpz.clone() - &z) * &zexpz_d) / ((2 * sqr(zexpz_d)) - ((zexpz - &z) * zexpz_dd))
5957    }
5958    w
5959}
5960fn initpoint(z: Complex, k: Integer) -> Complex {
5961    let prec = z.prec();
5962    {
5963        let e = Float::with_val(prec.0, -1).exp();
5964        let test: Complex = z.clone() + &e;
5965        if test.abs().real() <= &1.005 {
5966            if k == 0 && z.real() > &-0.5 {
5967                let p1: Complex = 2 * z / e + 2;
5968                let p = p1.clone().sqrt();
5969                return p.clone() - (p1 / 3) + ((11 * cube(p)) / 72) - 1;
5970            } else if k == -1 && z.real().is_sign_negative() && z.imag().is_zero() {
5971                let p1: Complex = 2 * z / e + 2;
5972                let p = p1.clone().sqrt();
5973                return -p.clone() - (p1 / 3) - ((11 * cube(p)) / 72) - 1;
5974            }
5975        }
5976    }
5977    {
5978        let test: Complex = z.clone() - 0.5;
5979        if test.abs().real() <= &0.6 {
5980            if k == 0 {
5981                return (0.35173371 * (0.1237166 + 7.061302897 * z.clone()))
5982                    / (2 + 0.827184 * (1 + 2 * z));
5983            } else if k == -1 {
5984                return (Complex::with_val(prec, (-2.2591588985, -4.22096))
5985                    * (Complex::with_val(prec, (-14.073271, -33.767687754)) * &z
5986                        + Complex::with_val(prec, (-12.7127, 19.071643)) * (1 + 2 * z.clone())))
5987                    / (2 + Complex::with_val(prec, (-17.23103, 10.629721)) * (1 + 2 * z));
5988            }
5989        }
5990    }
5991    let zln = z.ln() + Complex::with_val(prec, (0, 2 * Float::with_val(prec.0, Pi) * k));
5992    zln.clone() - zln.ln()
5993}
5994#[cfg(feature = "fastrand")]
5995pub fn rand_gamma(k: Float, t: Float) -> Float {
5996    let prec = k.prec();
5997    let mut sum = Float::new(prec);
5998    for _ in 1..=k
5999        .clone()
6000        .floor()
6001        .to_integer()
6002        .unwrap_or_default()
6003        .to_usize()
6004        .unwrap_or_default()
6005    {
6006        let u: Float = Float::with_val(prec, fastrand::u128(1..)) / u128::MAX;
6007        sum += u.ln();
6008    }
6009    let s = k.clone().fract();
6010    let e = Float::with_val(prec, 1).exp();
6011    let check = e.clone() / (e + s.clone());
6012    let mut eta: Float;
6013    loop {
6014        let u: Float = Float::with_val(prec, fastrand::u128(1..)) / u128::MAX;
6015        let v: Float = Float::with_val(prec, fastrand::u128(1..)) / u128::MAX;
6016        let w: Float = Float::with_val(prec, fastrand::u128(1..)) / u128::MAX;
6017        let n;
6018        if u <= check {
6019            eta = v.pow(1 / s.clone());
6020            n = w * eta.clone().pow(s.clone() - 1);
6021        } else {
6022            eta = 1 - v.ln();
6023            n = w * (-eta.clone()).exp();
6024        };
6025        if n <= eta.clone().pow(s.clone() - 1) * (-eta.clone()).exp() {
6026            break;
6027        }
6028    }
6029    let f: Float = eta - sum;
6030    t * f
6031}
6032#[cfg(feature = "fastrand")]
6033pub fn rand_norm(m: Complex, s: Complex) -> Complex {
6034    let prec = s.prec().0;
6035    let mut u: Float = Float::with_val(prec, fastrand::i128(i128::MIN + 2..i128::MAX)) / i128::MAX;
6036    let mut v: Float = Float::with_val(prec, fastrand::i128(i128::MIN + 2..i128::MAX)) / i128::MAX;
6037    let mut g: Float = u.clone().pow(2) + v.pow(2);
6038    while g >= 1 {
6039        u = Float::with_val(prec, fastrand::i128(i128::MIN + 2..i128::MAX)) / i128::MAX;
6040        v = Float::with_val(prec, fastrand::i128(i128::MIN + 2..i128::MAX)) / i128::MAX;
6041        g = u.clone().pow(2) + v.pow(2);
6042    }
6043    let d: Float = -2 * g.clone().ln() / g;
6044    let z: Float = u * d.sqrt();
6045    m + z * s
6046}
6047pub fn regularized_incomplete_beta(x: Complex, a: Complex, b: Complex) -> Complex {
6048    (gamma(a.clone() + b.clone())) * incomplete_beta(x, a.clone(), b.clone())
6049        / (gamma(a) * gamma(b))
6050}
6051pub fn sqr(z: Complex) -> Complex {
6052    if z.imag().is_zero() {
6053        z.pow(2)
6054    } else {
6055        z.clone() * z
6056    }
6057}
6058pub fn cube(z: Complex) -> Complex {
6059    if z.imag().is_zero() {
6060        z.pow(3)
6061    } else {
6062        z.clone() * &z * z
6063    }
6064}
6065pub fn pow_nth(z: Complex, n: Complex) -> Complex {
6066    let zz = z.imag().is_zero();
6067    let nz = n.imag().is_zero();
6068    let nr = n.real();
6069    if nr.clone().fract().is_zero() && nz {
6070        if !zz && nr <= &256 {
6071            if nr.is_zero() {
6072                Complex::with_val(z.prec(), 1)
6073            } else {
6074                let mut p = z.clone();
6075                for _ in 1..n
6076                    .real()
6077                    .to_integer()
6078                    .unwrap_or_default()
6079                    .abs()
6080                    .to_usize()
6081                    .unwrap_or_default()
6082                {
6083                    p *= &z;
6084                }
6085                if nr.is_sign_positive() { p } else { 1 / p }
6086            }
6087        } else {
6088            z.pow(n)
6089        }
6090    } else if zz && nz && {
6091        let zr = z.real();
6092        zr.clone().fract().is_zero() && zr.is_sign_negative()
6093    } {
6094        z.pow(n)
6095    } else {
6096        (z.ln() * n).exp()
6097    }
6098}
6099pub fn hsv2rgb(mut hue: Float, sat: Float, val: Float) -> Vec<Number> {
6100    if sat.is_zero() {
6101        return rgb2val(val.clone(), val.clone(), val);
6102    }
6103    hue *= 6;
6104    let i = hue
6105        .clone()
6106        .floor()
6107        .to_integer()
6108        .unwrap_or_default()
6109        .to_usize()
6110        .unwrap_or_default();
6111    let f = hue - i;
6112    let p = val.clone() * (1 - sat.clone());
6113    let q = val.clone() * (1 - sat.clone() * f.clone());
6114    let t = val.clone() * (1 - sat * (1 - f));
6115    match i % 6 {
6116        0 => rgb2val(val, t, p),
6117        1 => rgb2val(q, val, p),
6118        2 => rgb2val(p, val, t),
6119        3 => rgb2val(p, q, val),
6120        4 => rgb2val(t, p, val),
6121        5 => rgb2val(val, p, q),
6122        _ => rgb2val(val, p, q),
6123    }
6124}
6125pub fn rgb2val(r: Float, g: Float, b: Float) -> Vec<Number> {
6126    let r: Float = 255 * r;
6127    let g: Float = 255 * g;
6128    let b: Float = 255 * b;
6129    vec![
6130        Number::from(r.into(), None),
6131        Number::from(g.into(), None),
6132        Number::from(b.into(), None),
6133    ]
6134}
kalc_lib/complex.rs

kalc_lib/
complex.rs
