factorion_lib/
calculation_tasks.rs

1//! This module handles the calulation of pending calculation tasks
2use std::sync::OnceLock;
3
4use crate::calculation_results::Number;
5
6use crate::{
7    calculation_results::{Calculation, CalculationResult},
8    math,
9};
10
11use crate::FLOAT_PRECISION;
12
13use crate::rug::Integer;
14use crate::rug::{Float, ops::Pow};
15
16pub mod recommended {
17    use factorion_math::rug::Complete;
18    use factorion_math::rug::integer::IntegerExt64;
19
20    use crate::rug::Integer;
21    // Limit for exact calculation, set to limit calculation time
22    pub static UPPER_CALCULATION_LIMIT: fn() -> Integer = || 1_000_000.into();
23    // Limit for approximation, set to ensure enough accuracy (5 decimals)
24    pub static UPPER_APPROXIMATION_LIMIT: fn() -> Integer =
25        || Integer::u64_pow_u64(10, 300).complete();
26    // Limit for exact subfactorial calculation, set to limit calculation time
27    pub static UPPER_SUBFACTORIAL_LIMIT: fn() -> Integer = || 1_000_000.into();
28    // Limit for exact termial calculation, set to limit calculation time (absurdly high)
29    pub static UPPER_TERMIAL_LIMIT: fn() -> Integer = || Integer::u64_pow_u64(10, 10000).complete();
30    // Limit for approximation, set to ensure enough accuracy (5 decimals)
31    // Based on max float. (bits)
32    pub static UPPER_TERMIAL_APPROXIMATION_LIMIT: u32 = 1073741822;
33}
34
35static UPPER_CALCULATION_LIMIT: OnceLock<Integer> = OnceLock::new();
36static UPPER_APPROXIMATION_LIMIT: OnceLock<Integer> = OnceLock::new();
37static UPPER_SUBFACTORIAL_LIMIT: OnceLock<Integer> = OnceLock::new();
38static UPPER_TERMIAL_LIMIT: OnceLock<Integer> = OnceLock::new();
39static UPPER_TERMIAL_APPROXIMATION_LIMIT: OnceLock<u32> = OnceLock::new();
40
41use crate::AlreadyInit;
42pub fn init(
43    upper_calculation_limit: Integer,
44    upper_approximation_limit: Integer,
45    upper_subfactorial_limit: Integer,
46    upper_termial_limit: Integer,
47    upper_termial_approximation_limit: u32,
48) -> Result<(), AlreadyInit> {
49    static INITIALIZING: std::sync::Mutex<()> = std::sync::Mutex::new(());
50    let _guard = INITIALIZING.lock();
51    UPPER_CALCULATION_LIMIT
52        .set(upper_calculation_limit)
53        .map_err(|_| AlreadyInit)?;
54    UPPER_APPROXIMATION_LIMIT
55        .set(upper_approximation_limit)
56        .map_err(|_| AlreadyInit)?;
57    UPPER_SUBFACTORIAL_LIMIT
58        .set(upper_subfactorial_limit)
59        .map_err(|_| AlreadyInit)?;
60    UPPER_TERMIAL_LIMIT
61        .set(upper_termial_limit)
62        .map_err(|_| AlreadyInit)?;
63    UPPER_TERMIAL_APPROXIMATION_LIMIT
64        .set(upper_termial_approximation_limit)
65        .map_err(|_| AlreadyInit)?;
66    Ok(())
67}
68pub fn init_default() -> Result<(), AlreadyInit> {
69    use recommended::*;
70    init(
71        UPPER_CALCULATION_LIMIT(),
72        UPPER_APPROXIMATION_LIMIT(),
73        UPPER_SUBFACTORIAL_LIMIT(),
74        UPPER_TERMIAL_LIMIT(),
75        UPPER_TERMIAL_APPROXIMATION_LIMIT,
76    )
77}
78
79/// Representation of the calculation to be done
80#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
81pub struct CalculationJob {
82    pub base: CalculationBase,
83    /// Type of the calculation
84    pub level: i32,
85    /// Number of negations encountered
86    pub negative: u32,
87}
88/// The basis of a calculation, wheter [Number] or [CalculationJob].
89#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
90pub enum CalculationBase {
91    Num(Number),
92    Calc(Box<CalculationJob>),
93}
94
95impl CalculationJob {
96    /// Execute the calculation. \
97    /// If include_steps is enabled, will return all intermediate results.
98    pub fn execute(self, include_steps: bool) -> Vec<Option<Calculation>> {
99        let CalculationJob {
100            mut base,
101            mut level,
102            mut negative,
103        } = self;
104        let size = {
105            let mut n = 1;
106            let mut b = &base;
107            while let CalculationBase::Calc(inner) = b {
108                n += 1;
109                b = &inner.base;
110            }
111            n
112        };
113        // TODO: Maybe ignore include steps if size is too big (we can't respond properly anyway)
114        let mut steps = Vec::with_capacity(size);
115        let mut calcs = loop {
116            match base {
117                CalculationBase::Num(num) => {
118                    break vec![
119                        Self::calculate_appropriate_factorial(num.clone(), level, negative).map(
120                            |res| Calculation {
121                                value: num,
122                                steps: vec![(level, negative)],
123                                result: res,
124                            },
125                        ),
126                    ];
127                }
128                CalculationBase::Calc(calc) => {
129                    steps.push((level, negative));
130                    CalculationJob {
131                        base,
132                        level,
133                        negative,
134                    } = *calc;
135                }
136            }
137        };
138        for (level, negative) in steps.into_iter().rev() {
139            let calc = if include_steps {
140                calcs.last().cloned()
141            } else {
142                calcs.pop()
143            };
144            match calc {
145                Some(Some(Calculation {
146                    result: res,
147                    mut steps,
148                    value: number,
149                })) => {
150                    let factorial = Self::calculate_appropriate_factorial(res, level, negative)
151                        .map(|res| {
152                            steps.push((level, negative));
153                            Calculation {
154                                value: number,
155                                steps,
156                                result: res,
157                            }
158                        });
159                    calcs.push(factorial);
160                }
161                _ => return calcs,
162            };
163        }
164        calcs
165    }
166    fn calculate_appropriate_factorial(
167        num: Number,
168        level: i32,
169        negative: u32,
170    ) -> Option<CalculationResult> {
171        let prec = *FLOAT_PRECISION
172            .get()
173            .expect("Limit uninitialized, use init");
174        let calc_num = match &num {
175            CalculationResult::Approximate(base, exponent) => {
176                let res = base.as_float() * Float::with_val(prec, 10).pow(exponent);
177                if Float::is_finite(&(res.clone() * math::APPROX_FACT_SAFE_UPPER_BOUND_FACTOR)) {
178                    res.to_integer().unwrap()
179                } else {
180                    return Some(if base.as_float() < &0.0 {
181                        CalculationResult::ComplexInfinity
182                    } else if level < 0 {
183                        let termial = math::approximate_approx_termial(
184                            (Float::from(base.clone()), exponent.clone()),
185                            -level as u32,
186                        );
187                        CalculationResult::Approximate(termial.0.into(), termial.1)
188                    } else {
189                        CalculationResult::ApproximateDigitsTower(
190                            false,
191                            false,
192                            1,
193                            math::length(exponent, prec) + exponent,
194                        )
195                    });
196                }
197            }
198            CalculationResult::ApproximateDigits(was_neg, digits) => {
199                return Some(if digits.is_negative() {
200                    CalculationResult::Float(Float::new(prec).into())
201                } else if *was_neg {
202                    CalculationResult::ComplexInfinity
203                } else if level < 0 {
204                    CalculationResult::ApproximateDigits(false, (digits.clone() - 1) * 2 + 1)
205                } else {
206                    CalculationResult::ApproximateDigitsTower(
207                        false,
208                        false,
209                        1,
210                        math::length(digits, prec) + digits,
211                    )
212                });
213            }
214            CalculationResult::ApproximateDigitsTower(was_neg, neg, depth, exponent) => {
215                return Some(if *neg {
216                    CalculationResult::Float(Float::new(prec).into())
217                } else if *was_neg {
218                    CalculationResult::ComplexInfinity
219                } else if level < 0 {
220                    CalculationResult::ApproximateDigitsTower(
221                        false,
222                        false,
223                        *depth,
224                        exponent.clone(),
225                    )
226                } else {
227                    CalculationResult::ApproximateDigitsTower(
228                        false,
229                        false,
230                        depth + 1,
231                        exponent.clone(),
232                    )
233                });
234            }
235            CalculationResult::ComplexInfinity => return Some(CalculationResult::ComplexInfinity),
236            Number::Float(num) => match level {
237                ..0 => {
238                    let res: Float = math::fractional_termial(num.as_float().clone())
239                        * if negative % 2 != 0 { -1 } else { 1 };
240                    if res.is_finite() {
241                        return Some(CalculationResult::Float(res.into()));
242                    } else {
243                        num.as_float().to_integer()?
244                    }
245                }
246                0 => {
247                    // We don't support subfactorials of deciamals
248                    return None;
249                }
250                1 => {
251                    let res: Float = math::fractional_factorial(num.as_float().clone())
252                        * if negative % 2 != 0 { -1 } else { 1 };
253                    if res.is_finite() {
254                        return Some(CalculationResult::Float(res.into()));
255                    } else {
256                        num.as_float().to_integer()?
257                    }
258                }
259                2.. => {
260                    let res: Float =
261                        math::fractional_multifactorial(num.as_float().clone(), level as u32)
262                            * if negative % 2 != 0 { -1 } else { 1 };
263                    if res.is_finite() {
264                        return Some(CalculationResult::Float(res.into()));
265                    } else {
266                        num.as_float().to_integer()?
267                    }
268                }
269            },
270            Number::Exact(num) => num.clone(),
271        };
272        if level > 0 {
273            Some(if calc_num < 0 && level == 1 {
274                CalculationResult::ComplexInfinity
275            } else if calc_num < 0 {
276                let factor = math::negative_multifacorial_factor(calc_num.clone(), level);
277                match (factor, -level - 1 > calc_num) {
278                    (Some(factor), true) => {
279                        let mut res = Self::calculate_appropriate_factorial(
280                            Number::Exact(-calc_num.clone() - level),
281                            level,
282                            negative,
283                        )?;
284                        res = match res {
285                            CalculationResult::Exact(n) => {
286                                let n = Float::with_val(prec, n);
287                                CalculationResult::Float((factor / n).into())
288                            }
289                            CalculationResult::Approximate(b, e) => {
290                                let (b, e) =
291                                    math::adjust_approximate((factor / Float::from(b), -e));
292                                CalculationResult::Approximate(b.into(), e)
293                            }
294                            CalculationResult::ApproximateDigits(wn, n) => {
295                                CalculationResult::ApproximateDigits(wn, -n)
296                            }
297                            CalculationResult::ApproximateDigitsTower(
298                                wn,
299                                negative,
300                                depth,
301                                base,
302                            ) => CalculationResult::ApproximateDigitsTower(
303                                wn, !negative, depth, base,
304                            ),
305                            CalculationResult::ComplexInfinity => {
306                                CalculationResult::Exact(0.into())
307                            }
308                            CalculationResult::Float(f) => {
309                                CalculationResult::Float((factor / Float::from(f)).into())
310                            }
311                        };
312
313                        res
314                    }
315                    (factor, _) => factor
316                        .map(CalculationResult::Exact)
317                        .unwrap_or(CalculationResult::ComplexInfinity),
318                }
319                // Check if we can approximate the number of digits
320            } else if calc_num
321                > *UPPER_APPROXIMATION_LIMIT
322                    .get()
323                    .expect("Limit uninitialized, use init")
324            {
325                let factorial =
326                    math::approximate_multifactorial_digits(calc_num.clone(), level as u32, prec);
327                CalculationResult::ApproximateDigits(negative % 2 != 0, factorial)
328            // Check if the number is within a reasonable range to compute
329            } else if calc_num
330                > *UPPER_CALCULATION_LIMIT
331                    .get()
332                    .expect("Limit uninitialized, use init")
333            {
334                let factorial = if level == 0 {
335                    math::approximate_factorial(calc_num.clone(), prec)
336                } else {
337                    math::approximate_multifactorial(calc_num.clone(), level as u32, prec)
338                };
339                CalculationResult::Approximate(
340                    ((factorial.0 * if negative % 2 != 0 { -1 } else { 1 }) as Float).into(),
341                    factorial.1,
342                )
343            } else {
344                let calc_num = calc_num.to_u64().expect("Failed to convert BigInt to u64");
345                let factorial = math::factorial(calc_num, level as u32)
346                    * if negative % 2 != 0 { -1 } else { 1 };
347                CalculationResult::Exact(factorial)
348            })
349        } else if level == 0 {
350            Some(if calc_num < 0 {
351                CalculationResult::ComplexInfinity
352            } else if calc_num
353                > *UPPER_APPROXIMATION_LIMIT
354                    .get()
355                    .expect("Limit uninitialized, use init")
356            {
357                let factorial = math::approximate_multifactorial_digits(calc_num.clone(), 1, prec);
358                CalculationResult::ApproximateDigits(negative % 2 != 0, factorial)
359            } else if calc_num
360                > *UPPER_SUBFACTORIAL_LIMIT
361                    .get()
362                    .expect("Limit uninitialized, use init")
363            {
364                let factorial = math::approximate_subfactorial(calc_num.clone(), prec);
365                CalculationResult::Approximate(
366                    ((factorial.0 * if negative % 2 != 0 { -1 } else { 1 }) as Float).into(),
367                    factorial.1,
368                )
369            } else {
370                let calc_num = calc_num.to_u64().expect("Failed to convert BigInt to u64");
371                let factorial =
372                    math::subfactorial(calc_num) * if negative % 2 != 0 { -1 } else { 1 };
373                CalculationResult::Exact(factorial)
374            })
375        } else if level < 0 {
376            Some(
377                if calc_num.significant_bits()
378                    > *UPPER_TERMIAL_APPROXIMATION_LIMIT
379                        .get()
380                        .expect("Limit uninitialized, use init")
381                {
382                    let termial = math::approximate_termial_digits(calc_num, -level as u32, prec);
383                    CalculationResult::ApproximateDigits(negative % 2 != 0, termial)
384                } else if calc_num
385                    > *UPPER_TERMIAL_LIMIT
386                        .get()
387                        .expect("Limit uninitialized, use init")
388                {
389                    let termial = math::approximate_termial(calc_num, -level as u32, prec);
390                    CalculationResult::Approximate(
391                        ((termial.0 * if negative % 2 != 0 { -1 } else { 1 }) as Float).into(),
392                        termial.1,
393                    )
394                } else {
395                    let termial = if level < -1 {
396                        math::multitermial(calc_num, -level as u32)
397                    } else {
398                        math::termial(calc_num)
399                    };
400                    let termial = termial * if negative % 2 != 0 { -1 } else { 1 };
401                    CalculationResult::Exact(termial)
402                },
403            )
404        } else {
405            None
406        }
407    }
408}