scirs2_core/numeric/
arbitrary_precision.rs

1//! # Arbitrary Precision Numerical Computation Support
2//!
3//! This module provides arbitrary precision arithmetic capabilities for scientific computing,
4//! enabling calculations with user-defined precision levels for both integers and floating-point numbers.
5//!
6//! ## Features
7//!
8//! - Arbitrary precision integers (BigInt)
9//! - Arbitrary precision floating-point numbers (BigFloat)
10//! - Exact rational arithmetic (BigRational)
11//! - Arbitrary precision complex numbers (BigComplex)
12//! - Integration with existing ScientificNumber traits
13//! - Automatic precision tracking and management
14//! - Configurable precision contexts
15//! - Efficient operations with GMP/MPFR backend
16
17use crate::{
18    error::{CoreError, CoreResult, ErrorContext},
19    numeric::precision_tracking::PrecisionContext,
20    validation::check_positive,
21};
22use num_bigint::BigInt;
23use rug::{
24    float::Round, ops::Pow, Complex as RugComplex, Float as RugFloat, Integer as RugInteger,
25    Rational as RugRational,
26};
27use std::cmp::Ordering;
28use std::fmt;
29use std::ops::{Add, Div, Mul, Neg, Sub};
30use std::str::FromStr;
31use std::sync::RwLock;
32
33/// Global default precision for arbitrary precision operations
34static DEFAULT_PRECISION: RwLock<u32> = RwLock::new(256);
35
36/// Get the default precision for arbitrary precision operations
37#[allow(dead_code)]
38pub fn get_defaultprecision() -> u32 {
39    *DEFAULT_PRECISION.read().expect("Operation failed")
40}
41
42/// Set the default precision for arbitrary precision operations
43#[allow(dead_code)]
44pub fn setprecision(prec: u32) -> CoreResult<()> {
45    check_positive(prec as f64, "precision")?;
46    *DEFAULT_PRECISION.write().expect("Operation failed") = prec;
47    Ok(())
48}
49
50/// Precision context for arbitrary precision arithmetic
51#[derive(Debug, Clone)]
52pub struct ArbitraryPrecisionContext {
53    /// Precision in bits for floating-point operations
54    pub floatprecision: u32,
55    /// Maximum precision allowed
56    pub maxprecision: u32,
57    /// Rounding mode
58    pub rounding_mode: RoundingMode,
59    /// Whether to track precision loss
60    pub trackprecision: bool,
61    /// Precision tracking context
62    pub precision_context: Option<PrecisionContext>,
63}
64
65/// Rounding modes for arbitrary precision arithmetic
66#[derive(Debug, Clone, Copy, PartialEq, Eq)]
67pub enum RoundingMode {
68    /// Round to nearest, ties to even
69    Nearest,
70    /// Round toward zero
71    Zero,
72    /// Round toward positive infinity
73    Up,
74    /// Round toward negative infinity
75    Down,
76    /// Round away from zero
77    Away,
78}
79
80impl RoundingMode {
81    /// Convert to rug::float::Round
82    fn to_rug_round(self) -> Round {
83        match self {
84            RoundingMode::Nearest => Round::Nearest,
85            RoundingMode::Zero => Round::Zero,
86            RoundingMode::Up => Round::Up,
87            RoundingMode::Down => Round::Down,
88            RoundingMode::Away => Round::Up, // Use Up as Away is not available
89        }
90    }
91}
92
93impl Default for ArbitraryPrecisionContext {
94    fn default() -> Self {
95        Self {
96            floatprecision: get_defaultprecision(),
97            maxprecision: 4096,
98            rounding_mode: RoundingMode::Nearest,
99            trackprecision: false,
100            precision_context: None,
101        }
102    }
103}
104
105impl ArbitraryPrecisionContext {
106    /// Create a new precision context with specified precision
107    pub fn withprecision(precision: u32) -> CoreResult<Self> {
108        check_positive(precision as f64, "precision")?;
109        Ok(Self {
110            floatprecision: precision,
111            ..Default::default()
112        })
113    }
114
115    /// Create a context with precision tracking enabled
116    pub fn withprecision_tracking(precision: u32) -> CoreResult<Self> {
117        let mut ctx = Self::withprecision(precision)?;
118        ctx.trackprecision = true;
119        let mut precision_ctx = PrecisionContext::new();
120        precision_ctx.precision = precision as f64 / 3.32; // bits to decimal digits
121        ctx.precision_context = Some(precision_ctx);
122        Ok(ctx)
123    }
124
125    /// Set the rounding mode
126    pub fn with_rounding(mut self, mode: RoundingMode) -> Self {
127        self.rounding_mode = mode;
128        self
129    }
130
131    /// Set the maximum precision
132    pub fn with_maxprecision(mut self, maxprec: u32) -> Self {
133        self.maxprecision = maxprec;
134        self
135    }
136}
137
138/// Arbitrary precision integer
139#[derive(Clone, PartialEq, Eq)]
140pub struct ArbitraryInt {
141    value: RugInteger,
142}
143
144impl ArbitraryInt {
145    /// Create a new arbitrary precision integer
146    pub fn new() -> Self {
147        Self {
148            value: RugInteger::new(),
149        }
150    }
151
152    /// Create from a regular integer
153    pub fn from_i64(n: i64) -> Self {
154        Self {
155            value: RugInteger::from(n),
156        }
157    }
158
159    /// Create from a string
160    pub fn from_str_radix(s: &str, radix: i32) -> CoreResult<Self> {
161        RugInteger::parse_radix(s, radix)
162            .map(|incomplete| {
163                let value = RugInteger::from(incomplete);
164                Self { value }
165            })
166            .map_err(|_| {
167                CoreError::ValidationError(ErrorContext::new(format!(
168                    "Failed to parse integer from string: {}",
169                    s
170                )))
171            })
172    }
173
174    /// Convert to BigInt
175    pub fn to_bigint(&self) -> BigInt {
176        BigInt::from_str(&self.value.to_string()).expect("Operation failed")
177    }
178
179    /// Check if the number is prime
180    pub fn is_probably_prime(&self, reps: u32) -> bool {
181        self.value.is_probably_prime(reps) != rug::integer::IsPrime::No
182    }
183
184    /// Compute factorial
185    pub fn factorial(n: u32) -> Self {
186        let mut result = RugInteger::from(1);
187        for i in 2..=n {
188            result *= i;
189        }
190        Self { value: result }
191    }
192
193    /// Compute binomial coefficient
194    pub fn binomial(n: u32, k: u32) -> Self {
195        if k > n {
196            return Self::new();
197        }
198        let mut result = RugInteger::from(1);
199        for i in 0..k {
200            result *= n - i;
201            result /= i + 1;
202        }
203        Self { value: result }
204    }
205
206    /// Compute greatest common divisor
207    pub fn gcd(&self, other: &Self) -> Self {
208        Self {
209            value: self.value.clone().gcd(&other.value),
210        }
211    }
212
213    /// Compute least common multiple
214    pub fn lcm(&self, other: &Self) -> Self {
215        if self.value.is_zero() || other.value.is_zero() {
216            return Self::new();
217        }
218        let gcd = self.gcd(other);
219        // Create a complete integer from the multiplication
220        let product = RugInteger::from(&self.value * &other.value);
221        Self {
222            value: product / gcd.value,
223        }
224    }
225
226    /// Modular exponentiation
227    pub fn mod_pow(&self, exp: &Self, modulus: &Self) -> CoreResult<Self> {
228        if modulus.value.is_zero() {
229            return Err(CoreError::DomainError(ErrorContext::new(
230                "Modulus cannot be zero",
231            )));
232        }
233        Ok(Self {
234            value: self
235                .value
236                .clone()
237                .pow_mod(&exp.value, &modulus.value)
238                .expect("Operation failed"),
239        })
240    }
241
242    /// Get the absolute value
243    pub fn abs(&self) -> Self {
244        Self {
245            value: self.value.clone().abs(),
246        }
247    }
248
249    /// Get the sign (-1, 0, or 1)
250    pub fn signum(&self) -> i32 {
251        self.value.cmp0() as i32
252    }
253}
254
255impl fmt::Display for ArbitraryInt {
256    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
257        write!(f, "{}", self.value)
258    }
259}
260
261impl fmt::Debug for ArbitraryInt {
262    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
263        write!(f, "ArbitraryInt({})", self.value)
264    }
265}
266
267impl Default for ArbitraryInt {
268    fn default() -> Self {
269        Self::new()
270    }
271}
272
273// Implement arithmetic operations for ArbitraryInt
274impl Add for ArbitraryInt {
275    type Output = Self;
276    fn add(self, rhs: Self) -> Self::Output {
277        Self {
278            value: self.value + rhs.value,
279        }
280    }
281}
282
283impl Sub for ArbitraryInt {
284    type Output = Self;
285    fn sub(self, rhs: Self) -> Self::Output {
286        Self {
287            value: self.value - rhs.value,
288        }
289    }
290}
291
292impl Mul for ArbitraryInt {
293    type Output = Self;
294    fn mul(self, rhs: Self) -> Self::Output {
295        Self {
296            value: self.value * rhs.value,
297        }
298    }
299}
300
301impl Div for ArbitraryInt {
302    type Output = Self;
303    fn div(self, rhs: Self) -> Self::Output {
304        Self {
305            value: self.value / rhs.value,
306        }
307    }
308}
309
310impl Neg for ArbitraryInt {
311    type Output = Self;
312    fn neg(self) -> Self::Output {
313        Self { value: -self.value }
314    }
315}
316
317/// Arbitrary precision floating-point number
318#[derive(Clone)]
319pub struct ArbitraryFloat {
320    value: RugFloat,
321    context: ArbitraryPrecisionContext,
322}
323
324impl ArbitraryFloat {
325    /// Create a new arbitrary precision float with default precision
326    pub fn new() -> Self {
327        let prec = get_defaultprecision();
328        Self {
329            value: RugFloat::new(prec),
330            context: ArbitraryPrecisionContext::default(),
331        }
332    }
333
334    /// Create with specific precision
335    pub fn withprecision(prec: u32) -> CoreResult<Self> {
336        let context = ArbitraryPrecisionContext::withprecision(prec)?;
337        Ok(Self {
338            value: RugFloat::new(prec),
339            context,
340        })
341    }
342
343    /// Create with specific context
344    pub fn with_context(context: ArbitraryPrecisionContext) -> Self {
345        Self {
346            value: RugFloat::new(context.floatprecision),
347            context,
348        }
349    }
350
351    /// Create from f64 with default precision
352    pub fn from_f64(value: f64) -> Self {
353        let prec = get_defaultprecision();
354        Self {
355            value: RugFloat::with_val(prec, value),
356            context: ArbitraryPrecisionContext::default(),
357        }
358    }
359
360    /// Create from f64 with specific precision
361    pub fn from_f64_withprecision(value: f64, prec: u32) -> CoreResult<Self> {
362        let context = ArbitraryPrecisionContext::withprecision(prec)?;
363        Ok(Self {
364            value: RugFloat::with_val(prec, value),
365            context,
366        })
367    }
368
369    /// Create from string with specific precision
370    pub fn from_strprec(s: &str, prec: u32) -> CoreResult<Self> {
371        let context = ArbitraryPrecisionContext::withprecision(prec)?;
372        RugFloat::parse(s)
373            .map(|incomplete| Self {
374                value: RugFloat::with_val(prec, incomplete),
375                context,
376            })
377            .map_err(|e| CoreError::ValidationError(ErrorContext::new(format!("{e}"))))
378    }
379
380    /// Get the precision in bits
381    pub fn precision(&self) -> u32 {
382        self.value.prec()
383    }
384
385    /// Get the precision in decimal digits
386    pub fn decimalprecision(&self) -> u32 {
387        (self.value.prec() as f64 / 3.32) as u32
388    }
389
390    /// Set precision (returns a new value with the specified precision)
391    pub fn setprecision(&self, prec: u32) -> CoreResult<Self> {
392        let mut context = self.context.clone();
393        context.floatprecision = prec;
394        Ok(Self {
395            value: RugFloat::with_val(prec, &self.value),
396            context,
397        })
398    }
399
400    /// Convert to f64 (may lose precision)
401    pub fn to_f64(&self) -> f64 {
402        self.value.to_f64()
403    }
404
405    /// Check if the value is finite
406    pub fn is_finite(&self) -> bool {
407        self.value.is_finite()
408    }
409
410    /// Check if the value is infinite
411    pub fn is_infinite(&self) -> bool {
412        self.value.is_infinite()
413    }
414
415    /// Check if the value is NaN
416    pub fn is_nan(&self) -> bool {
417        self.value.is_nan()
418    }
419
420    /// Check if the value is zero
421    pub fn is_zero(&self) -> bool {
422        self.value.is_zero()
423    }
424
425    /// Get the absolute value
426    pub fn abs(&self) -> Self {
427        let mut result = self.clone();
428        result.value.abs_mut();
429        result
430    }
431
432    /// Square root
433    pub fn sqrt(&self) -> CoreResult<Self> {
434        if self.value.is_sign_negative() && !self.value.is_zero() {
435            return Err(CoreError::DomainError(ErrorContext::new(
436                "Square root of negative number",
437            )));
438        }
439        let mut result = self.clone();
440        result.value.sqrt_mut();
441        Ok(result)
442    }
443
444    /// Natural logarithm
445    pub fn ln(&self) -> CoreResult<Self> {
446        if self.value.is_sign_negative() || self.value.is_zero() {
447            return Err(CoreError::DomainError(ErrorContext::new(
448                "Logarithm of non-positive number",
449            )));
450        }
451        let mut result = self.clone();
452        result.value.ln_mut();
453        Ok(result)
454    }
455
456    /// Exponential function
457    pub fn exp(&self) -> Self {
458        let mut result = self.clone();
459        result.value.exp_mut();
460        result
461    }
462
463    /// Power function
464    pub fn pow(&self, exp: &Self) -> Self {
465        let mut result = self.clone();
466        let pow_result = self.value.clone().pow(&exp.value);
467        result.value = RugFloat::with_val(self.context.floatprecision, pow_result);
468        result
469    }
470
471    /// Sine
472    pub fn sin(&self) -> Self {
473        let mut result = self.clone();
474        result.value.sin_mut();
475        result
476    }
477
478    /// Cosine
479    pub fn cos(&self) -> Self {
480        let mut result = self.clone();
481        result.value.cos_mut();
482        result
483    }
484
485    /// Tangent
486    pub fn tan(&self) -> Self {
487        let mut result = self.clone();
488        result.value.tan_mut();
489        result
490    }
491
492    /// Arcsine
493    pub fn asin(&self) -> CoreResult<Self> {
494        if self.value.clone().abs() > 1 {
495            return Err(CoreError::DomainError(ErrorContext::new(
496                "Arcsine argument out of range [-1, 1]",
497            )));
498        }
499        let mut result = self.clone();
500        result.value.asin_mut();
501        Ok(result)
502    }
503
504    /// Arccosine
505    pub fn acos(&self) -> CoreResult<Self> {
506        if self.value.clone().abs() > 1 {
507            return Err(CoreError::DomainError(ErrorContext::new(
508                "Arccosine argument out of range [-1, 1]",
509            )));
510        }
511        let mut result = self.clone();
512        result.value.acos_mut();
513        Ok(result)
514    }
515
516    /// Arctangent
517    pub fn atan(&self) -> Self {
518        let mut result = self.clone();
519        result.value.atan_mut();
520        result
521    }
522
523    /// Two-argument arctangent
524    pub fn atan2(&self, x: &Self) -> Self {
525        let mut result = self.clone();
526        result.value = RugFloat::with_val(
527            self.context.floatprecision,
528            self.value.clone().atan2(&x.value),
529        );
530        result
531    }
532
533    /// Hyperbolic sine
534    pub fn sinh(&self) -> Self {
535        let mut result = self.clone();
536        result.value.sinh_mut();
537        result
538    }
539
540    /// Hyperbolic cosine
541    pub fn cosh(&self) -> Self {
542        let mut result = self.clone();
543        result.value.cosh_mut();
544        result
545    }
546
547    /// Hyperbolic tangent
548    pub fn tanh(&self) -> Self {
549        let mut result = self.clone();
550        result.value.tanh_mut();
551        result
552    }
553
554    /// Constants
555    pub fn prec_2(prec: u32) -> CoreResult<Self> {
556        let context = ArbitraryPrecisionContext::withprecision(prec)?;
557        let value = RugFloat::with_val(prec, rug::float::Constant::Pi);
558        Ok(Self { value, context })
559    }
560
561    pub fn prec_3(prec: u32) -> CoreResult<Self> {
562        let one = Self::from_f64_withprecision(1.0, prec)?;
563        Ok(one.exp())
564    }
565
566    pub fn prec_4(prec: u32) -> CoreResult<Self> {
567        let context = ArbitraryPrecisionContext::withprecision(prec)?;
568        let value = RugFloat::with_val(prec, rug::float::Constant::Log2);
569        Ok(Self { value, context })
570    }
571}
572
573impl fmt::Display for ArbitraryFloat {
574    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
575        write!(f, "{}", self.value)
576    }
577}
578
579impl fmt::Debug for ArbitraryFloat {
580    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
581        write!(
582            f,
583            "ArbitraryFloat({}, {} bits)",
584            self.value,
585            self.precision()
586        )
587    }
588}
589
590impl PartialEq for ArbitraryFloat {
591    fn eq(&self, other: &Self) -> bool {
592        self.value.eq(&other.value)
593    }
594}
595
596impl PartialOrd for ArbitraryFloat {
597    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
598        self.value.partial_cmp(&other.value)
599    }
600}
601
602impl Default for ArbitraryFloat {
603    fn default() -> Self {
604        Self::new()
605    }
606}
607
608// Arithmetic operations for ArbitraryFloat
609impl Add for ArbitraryFloat {
610    type Output = Self;
611    fn add(self, rhs: Self) -> Self::Output {
612        Self {
613            value: self.value + rhs.value,
614            context: self.context,
615        }
616    }
617}
618
619impl Sub for ArbitraryFloat {
620    type Output = Self;
621    fn sub(self, rhs: Self) -> Self::Output {
622        Self {
623            value: self.value - rhs.value,
624            context: self.context,
625        }
626    }
627}
628
629impl Mul for ArbitraryFloat {
630    type Output = Self;
631    fn mul(self, rhs: Self) -> Self::Output {
632        Self {
633            value: self.value * rhs.value,
634            context: self.context,
635        }
636    }
637}
638
639impl Div for ArbitraryFloat {
640    type Output = Self;
641    fn div(self, rhs: Self) -> Self::Output {
642        Self {
643            value: self.value / rhs.value,
644            context: self.context,
645        }
646    }
647}
648
649impl Neg for ArbitraryFloat {
650    type Output = Self;
651    fn neg(self) -> Self::Output {
652        Self {
653            value: -self.value,
654            context: self.context,
655        }
656    }
657}
658
659/// Arbitrary precision rational number
660#[derive(Clone, PartialEq, Eq)]
661pub struct ArbitraryRational {
662    value: RugRational,
663}
664
665impl ArbitraryRational {
666    /// Create a new rational number (0/1)
667    pub fn new() -> Self {
668        Self {
669            value: RugRational::new(),
670        }
671    }
672
673    /// Create from numerator and denominator
674    pub fn num(num: i64, den: i64) -> CoreResult<Self> {
675        if den == 0 {
676            return Err(CoreError::DomainError(ErrorContext::new(
677                "Denominator cannot be zero",
678            )));
679        }
680        Ok(Self {
681            value: RugRational::from((num, den)),
682        })
683    }
684
685    /// Create from arbitrary precision integers
686    pub fn num_2(num: &ArbitraryInt, den: &ArbitraryInt) -> CoreResult<Self> {
687        if den.value.is_zero() {
688            return Err(CoreError::DomainError(ErrorContext::new(
689                "Denominator cannot be zero",
690            )));
691        }
692        Ok(Self {
693            value: RugRational::from((&num.value, &den.value)),
694        })
695    }
696
697    /// Create from string (e.g., "22/7" or "3.14159")  
698    /// Note: This is deprecated, use `str::parse()` instead
699    #[deprecated(note = "Use str::parse() instead")]
700    pub fn parse_rational(s: &str) -> CoreResult<Self> {
701        s.parse()
702    }
703
704    /// Convert to f64 (may lose precision)
705    pub fn to_f64(&self) -> f64 {
706        self.value.to_f64()
707    }
708
709    /// Convert to arbitrary precision float
710    pub fn to_arbitrary_float(&self, prec: u32) -> CoreResult<ArbitraryFloat> {
711        let context = ArbitraryPrecisionContext::withprecision(prec)?;
712        Ok(ArbitraryFloat {
713            value: RugFloat::with_val(prec, &self.value),
714            context,
715        })
716    }
717
718    /// Get numerator
719    pub fn numerator(&self) -> ArbitraryInt {
720        ArbitraryInt {
721            value: self.value.numer().clone(),
722        }
723    }
724
725    /// Get denominator
726    pub fn denominator(&self) -> ArbitraryInt {
727        ArbitraryInt {
728            value: self.value.denom().clone(),
729        }
730    }
731
732    /// Get the absolute value
733    pub fn abs(&self) -> Self {
734        Self {
735            value: self.value.clone().abs(),
736        }
737    }
738
739    /// Get the reciprocal
740    pub fn recip(&self) -> CoreResult<Self> {
741        if self.value.is_zero() {
742            return Err(CoreError::DomainError(ErrorContext::new(
743                "Cannot take reciprocal of zero",
744            )));
745        }
746        Ok(Self {
747            value: self.value.clone().recip(),
748        })
749    }
750}
751
752impl fmt::Display for ArbitraryRational {
753    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
754        write!(f, "{}", self.value)
755    }
756}
757
758impl fmt::Debug for ArbitraryRational {
759    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
760        write!(f, "ArbitraryRational({})", self.value)
761    }
762}
763
764impl Default for ArbitraryRational {
765    fn default() -> Self {
766        Self::new()
767    }
768}
769
770impl FromStr for ArbitraryRational {
771    type Err = CoreError;
772
773    fn from_str(s: &str) -> Result<Self, Self::Err> {
774        RugRational::from_str(s)
775            .map(|value| Self { value })
776            .map_err(|_| {
777                CoreError::ValidationError(ErrorContext::new(format!(
778                    "Failed to parse rational from string: {}",
779                    s
780                )))
781            })
782    }
783}
784
785// Arithmetic operations for ArbitraryRational
786impl Add for ArbitraryRational {
787    type Output = Self;
788    fn add(self, rhs: Self) -> Self::Output {
789        Self {
790            value: self.value + rhs.value,
791        }
792    }
793}
794
795impl Sub for ArbitraryRational {
796    type Output = Self;
797    fn sub(self, rhs: Self) -> Self::Output {
798        Self {
799            value: self.value - rhs.value,
800        }
801    }
802}
803
804impl Mul for ArbitraryRational {
805    type Output = Self;
806    fn mul(self, rhs: Self) -> Self::Output {
807        Self {
808            value: self.value * rhs.value,
809        }
810    }
811}
812
813impl Div for ArbitraryRational {
814    type Output = Self;
815    fn div(self, rhs: Self) -> Self::Output {
816        Self {
817            value: self.value / rhs.value,
818        }
819    }
820}
821
822impl Neg for ArbitraryRational {
823    type Output = Self;
824    fn neg(self) -> Self::Output {
825        Self { value: -self.value }
826    }
827}
828
829/// Arbitrary precision complex number
830#[derive(Clone)]
831pub struct ArbitraryComplex {
832    value: RugComplex,
833    context: ArbitraryPrecisionContext,
834}
835
836impl ArbitraryComplex {
837    /// Create a new complex number with default precision
838    pub fn new() -> Self {
839        let prec = get_defaultprecision();
840        Self {
841            value: RugComplex::new(prec),
842            context: ArbitraryPrecisionContext::default(),
843        }
844    }
845
846    /// Create with specific precision
847    pub fn prec(prec: u32) -> CoreResult<Self> {
848        let context = ArbitraryPrecisionContext::withprecision(prec)?;
849        Ok(Self {
850            value: RugComplex::new(prec),
851            context,
852        })
853    }
854
855    /// Create from real and imaginary parts
856    pub fn re(re: &ArbitraryFloat, im: &ArbitraryFloat) -> Self {
857        let prec = re.precision().max(im.precision());
858        let context = re.context.clone();
859        Self {
860            value: RugComplex::with_val(prec, (&re.value, &im.value)),
861            context,
862        }
863    }
864
865    /// Create from f64 parts
866    pub fn re_2(re: f64, im: f64) -> Self {
867        let prec = get_defaultprecision();
868        Self {
869            value: RugComplex::with_val(prec, (re, im)),
870            context: ArbitraryPrecisionContext::default(),
871        }
872    }
873
874    /// Get real part
875    pub fn real(&self) -> ArbitraryFloat {
876        ArbitraryFloat {
877            value: self.value.real().clone(),
878            context: self.context.clone(),
879        }
880    }
881
882    /// Get imaginary part
883    pub fn imag(&self) -> ArbitraryFloat {
884        ArbitraryFloat {
885            value: self.value.imag().clone(),
886            context: self.context.clone(),
887        }
888    }
889
890    /// Get magnitude (absolute value)
891    pub fn abs(&self) -> ArbitraryFloat {
892        let mut result = ArbitraryFloat::with_context(self.context.clone());
893        result.value = RugFloat::with_val(self.context.floatprecision, self.value.abs_ref());
894        result
895    }
896
897    /// Get phase (argument)
898    pub fn arg(&self) -> ArbitraryFloat {
899        let mut result = ArbitraryFloat::with_context(self.context.clone());
900        result.value = RugFloat::with_val(self.context.floatprecision, self.value.arg_ref());
901        result
902    }
903
904    /// Complex conjugate
905    pub fn conj(&self) -> Self {
906        let mut result = self.clone();
907        result.value.conj_mut();
908        result
909    }
910
911    /// Natural logarithm
912    pub fn ln(&self) -> Self {
913        let mut result = self.clone();
914        result.value.ln_mut();
915        result
916    }
917
918    /// Exponential function
919    pub fn exp(&self) -> Self {
920        let mut result = self.clone();
921        result.value.exp_mut();
922        result
923    }
924
925    /// Power function
926    pub fn pow(&self, exp: &Self) -> Self {
927        // z^w = exp(w * ln(z))
928        let ln_self = self.ln();
929        let product = ln_self * exp.clone();
930        product.exp()
931    }
932
933    /// Square root
934    pub fn sqrt(&self) -> Self {
935        let mut result = self.clone();
936        result.value.sqrt_mut();
937        result
938    }
939}
940
941impl fmt::Display for ArbitraryComplex {
942    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
943        let re = self.value.real();
944        let im = self.value.imag();
945        if im.is_sign_positive() {
946            write!(f, "{} + {}i", re, im)
947        } else {
948            write!(f, "{} - {}i", re, -im.clone())
949        }
950    }
951}
952
953impl fmt::Debug for ArbitraryComplex {
954    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
955        write!(
956            f,
957            "ArbitraryComplex({}, {} bits)",
958            self, self.context.floatprecision
959        )
960    }
961}
962
963impl PartialEq for ArbitraryComplex {
964    fn eq(&self, other: &Self) -> bool {
965        self.value.eq(&other.value)
966    }
967}
968
969impl Default for ArbitraryComplex {
970    fn default() -> Self {
971        Self::new()
972    }
973}
974
975// Arithmetic operations for ArbitraryComplex
976impl Add for ArbitraryComplex {
977    type Output = Self;
978    fn add(self, rhs: Self) -> Self::Output {
979        Self {
980            value: self.value + rhs.value,
981            context: self.context,
982        }
983    }
984}
985
986impl Sub for ArbitraryComplex {
987    type Output = Self;
988    fn sub(self, rhs: Self) -> Self::Output {
989        Self {
990            value: self.value - rhs.value,
991            context: self.context,
992        }
993    }
994}
995
996impl Mul for ArbitraryComplex {
997    type Output = Self;
998    fn mul(self, rhs: Self) -> Self::Output {
999        Self {
1000            value: self.value * rhs.value,
1001            context: self.context,
1002        }
1003    }
1004}
1005
1006impl Div for ArbitraryComplex {
1007    type Output = Self;
1008    fn div(self, rhs: Self) -> Self::Output {
1009        Self {
1010            value: self.value / rhs.value,
1011            context: self.context,
1012        }
1013    }
1014}
1015
1016impl Neg for ArbitraryComplex {
1017    type Output = Self;
1018    fn neg(self) -> Self::Output {
1019        Self {
1020            value: -self.value,
1021            context: self.context,
1022        }
1023    }
1024}
1025
1026/// Conversion trait for arbitrary precision types
1027pub trait ToArbitraryPrecision {
1028    /// The arbitrary precision type
1029    type ArbitraryType;
1030
1031    /// Convert to arbitrary precision with default precision
1032    fn to_arbitrary(&self) -> Self::ArbitraryType;
1033
1034    /// Convert to arbitrary precision with specified precision
1035    fn to_arbitraryprec(&self, prec: u32) -> CoreResult<Self::ArbitraryType>;
1036}
1037
1038impl ToArbitraryPrecision for i32 {
1039    type ArbitraryType = ArbitraryInt;
1040
1041    fn to_arbitrary(&self) -> Self::ArbitraryType {
1042        ArbitraryInt::from_i64(*self as i64)
1043    }
1044
1045    fn to_arbitraryprec(&self, prec: u32) -> CoreResult<Self::ArbitraryType> {
1046        Ok(self.to_arbitrary())
1047    }
1048}
1049
1050impl ToArbitraryPrecision for i64 {
1051    type ArbitraryType = ArbitraryInt;
1052
1053    fn to_arbitrary(&self) -> Self::ArbitraryType {
1054        ArbitraryInt::from_i64(*self)
1055    }
1056
1057    fn to_arbitraryprec(&self, prec: u32) -> CoreResult<Self::ArbitraryType> {
1058        Ok(self.to_arbitrary())
1059    }
1060}
1061
1062impl ToArbitraryPrecision for f32 {
1063    type ArbitraryType = ArbitraryFloat;
1064
1065    fn to_arbitrary(&self) -> Self::ArbitraryType {
1066        ArbitraryFloat::from_f64(*self as f64)
1067    }
1068
1069    fn to_arbitraryprec(&self, prec: u32) -> CoreResult<Self::ArbitraryType> {
1070        ArbitraryFloat::from_f64_withprecision(*self as f64, prec)
1071    }
1072}
1073
1074impl ToArbitraryPrecision for f64 {
1075    type ArbitraryType = ArbitraryFloat;
1076
1077    fn to_arbitrary(&self) -> Self::ArbitraryType {
1078        ArbitraryFloat::from_f64(*self)
1079    }
1080
1081    fn to_arbitraryprec(&self, prec: u32) -> CoreResult<Self::ArbitraryType> {
1082        ArbitraryFloat::from_f64_withprecision(*self, prec)
1083    }
1084}
1085
1086/// Builder for arbitrary precision calculations
1087pub struct ArbitraryPrecisionBuilder {
1088    context: ArbitraryPrecisionContext,
1089}
1090
1091impl ArbitraryPrecisionBuilder {
1092    /// Create a new builder with default settings
1093    pub fn new() -> Self {
1094        Self {
1095            context: ArbitraryPrecisionContext::default(),
1096        }
1097    }
1098
1099    /// Set the precision in bits
1100    pub fn precision(mut self, prec: u32) -> Self {
1101        self.context.floatprecision = prec;
1102        self
1103    }
1104
1105    /// Set the precision in decimal digits
1106    pub fn decimalprecision(mut self, digits: u32) -> Self {
1107        self.context.floatprecision = ((digits as f64) * 3.32) as u32;
1108        self
1109    }
1110
1111    /// Set the rounding mode
1112    pub fn rounding(mut self, mode: RoundingMode) -> Self {
1113        self.context.rounding_mode = mode;
1114        self
1115    }
1116
1117    /// Enable precision tracking
1118    pub fn trackprecision(mut self, track: bool) -> Self {
1119        self.context.trackprecision = track;
1120        if track && self.context.precision_context.is_none() {
1121            let mut precision_ctx = PrecisionContext::new();
1122            precision_ctx.precision = self.context.floatprecision as f64 / 3.32;
1123            self.context.precision_context = Some(precision_ctx);
1124        }
1125        self
1126    }
1127
1128    /// Build an ArbitraryFloat
1129    pub fn build_float(self) -> ArbitraryFloat {
1130        ArbitraryFloat::with_context(self.context)
1131    }
1132
1133    /// Build an ArbitraryComplex
1134    pub fn build_complex(self) -> CoreResult<ArbitraryComplex> {
1135        ArbitraryComplex::prec(self.context.floatprecision)
1136    }
1137
1138    /// Execute a calculation with this precision context
1139    pub fn calculate<F, R>(self, f: F) -> R
1140    where
1141        F: FnOnce(&ArbitraryPrecisionContext) -> R,
1142    {
1143        f(&self.context)
1144    }
1145}
1146
1147impl Default for ArbitraryPrecisionBuilder {
1148    fn default() -> Self {
1149        Self::new()
1150    }
1151}
1152
1153/// Utility functions for arbitrary precision arithmetic
1154pub mod utils {
1155    use super::*;
1156
1157    /// Compute π to specified precision
1158    pub fn pi(prec: u32) -> CoreResult<ArbitraryFloat> {
1159        ArbitraryFloat::prec_2(prec)
1160    }
1161
1162    /// Compute e to specified precision
1163    pub fn e(prec: u32) -> CoreResult<ArbitraryFloat> {
1164        ArbitraryFloat::prec_3(prec)
1165    }
1166
1167    /// Compute ln(2) to specified precision
1168    pub fn ln2(prec: u32) -> CoreResult<ArbitraryFloat> {
1169        ArbitraryFloat::prec_4(prec)
1170    }
1171
1172    /// Compute sqrt(2) to specified precision
1173    pub fn sqrt2(prec: u32) -> CoreResult<ArbitraryFloat> {
1174        let two = ArbitraryFloat::from_f64_withprecision(2.0, prec)?;
1175        two.sqrt()
1176    }
1177
1178    /// Compute golden ratio to specified precision
1179    pub fn golden_ratio(prec: u32) -> CoreResult<ArbitraryFloat> {
1180        let one = ArbitraryFloat::from_f64_withprecision(1.0, prec)?;
1181        let five = ArbitraryFloat::from_f64_withprecision(5.0, prec)?;
1182        let sqrt5 = five.sqrt()?;
1183        let two = ArbitraryFloat::from_f64_withprecision(2.0, prec)?;
1184        Ok((one + sqrt5) / two)
1185    }
1186
1187    /// Compute factorial using arbitrary precision
1188    pub fn factorial(n: u32) -> ArbitraryInt {
1189        ArbitraryInt::factorial(n)
1190    }
1191
1192    /// Compute binomial coefficient using arbitrary precision
1193    pub fn binomial(n: u32, k: u32) -> ArbitraryInt {
1194        ArbitraryInt::binomial(n, k)
1195    }
1196
1197    /// Check if a large integer is probably prime
1198    pub fn is_probably_prime(n: &ArbitraryInt, certainty: u32) -> bool {
1199        n.is_probably_prime(certainty)
1200    }
1201}
1202
1203#[cfg(test)]
1204mod tests {
1205    use super::*;
1206
1207    #[test]
1208    fn test_arbitrary_int_basic() {
1209        let a = ArbitraryInt::from_i64(123);
1210        let b = ArbitraryInt::from_i64(456);
1211        let sum = a.clone() + b.clone();
1212        assert_eq!(sum.to_string(), "579");
1213
1214        let product = a.clone() * b.clone();
1215        assert_eq!(product.to_string(), "56088");
1216
1217        let factorial = ArbitraryInt::factorial(20);
1218        assert_eq!(factorial.to_string(), "2432902008176640000");
1219    }
1220
1221    #[test]
1222    fn test_arbitrary_float_basic() {
1223        let a = ArbitraryFloat::from_f64_withprecision(1.0, 128).expect("Operation failed");
1224        let b = ArbitraryFloat::from_f64_withprecision(3.0, 128).expect("Operation failed");
1225        let c = a / b;
1226
1227        // Check that we get more precision than f64
1228        let c_str = c.to_string();
1229        // Scientific notation: 3.3333...e-1 or decimal: 0.3333...
1230        assert!(c_str.starts_with("3.333333333333333") || c_str.starts_with("0.333333333333333"));
1231        assert!(c_str.len() > 20); // More digits than f64 can represent
1232    }
1233
1234    #[test]
1235    fn test_arbitrary_rational() {
1236        let r = ArbitraryRational::num(22, 7).expect("Operation failed");
1237        assert_eq!(r.to_string(), "22/7");
1238
1239        let a = ArbitraryRational::num(1, 3).expect("Operation failed");
1240        let b = ArbitraryRational::num(1, 6).expect("Operation failed");
1241        let sum = a + b;
1242        assert_eq!(sum.to_string(), "1/2");
1243    }
1244
1245    #[test]
1246    fn test_arbitrary_complex() {
1247        let z = ArbitraryComplex::re_2(3.0, 4.0);
1248        let mag = z.abs();
1249        assert!((mag.to_f64() - 5.0).abs() < 1e-10);
1250
1251        let conj = z.conj();
1252        assert_eq!(conj.real().to_f64(), 3.0);
1253        assert_eq!(conj.imag().to_f64(), -4.0);
1254    }
1255
1256    #[test]
1257    fn testprecision_builder() {
1258        let x = ArbitraryPrecisionBuilder::new()
1259            .decimalprecision(50)
1260            .rounding(RoundingMode::Nearest)
1261            .build_float();
1262
1263        assert!(x.decimalprecision() >= 49); // Allow for rounding
1264    }
1265
1266    #[test]
1267    fn test_constants() {
1268        let pi = utils::pi(256).expect("Operation failed");
1269        let pi_str = pi.to_string();
1270        assert!(pi_str.starts_with("3.14159265358979"));
1271
1272        let e = utils::e(256).expect("Operation failed");
1273        let e_str = e.to_string();
1274        assert!(e_str.starts_with("2.71828182845904"));
1275    }
1276
1277    #[test]
1278    fn test_prime_checking() {
1279        let prime = ArbitraryInt::from_i64(97);
1280        assert!(prime.is_probably_prime(20));
1281
1282        let composite = ArbitraryInt::from_i64(98);
1283        assert!(!composite.is_probably_prime(20));
1284    }
1285
1286    #[test]
1287    fn test_gcd_lcm() {
1288        let a = ArbitraryInt::from_i64(48);
1289        let b = ArbitraryInt::from_i64(18);
1290
1291        let gcd = a.gcd(&b);
1292        assert_eq!(gcd.to_string(), "6");
1293
1294        let lcm = a.lcm(&b);
1295        assert_eq!(lcm.to_string(), "144");
1296    }
1297
1298    #[test]
1299    fn test_transcendental_functions() {
1300        let x = ArbitraryFloat::from_f64_withprecision(0.5, 128).expect("Operation failed");
1301
1302        let sin_x = x.sin();
1303        let cos_x = x.cos();
1304        let identity = sin_x.clone() * sin_x + cos_x.clone() * cos_x;
1305
1306        // sin²(x) + cos²(x) = 1
1307        assert!((identity.to_f64() - 1.0).abs() < 1e-15);
1308
1309        let ln_x = x.ln().expect("Operation failed");
1310        let exp_ln_x = ln_x.exp();
1311        assert!((exp_ln_x.to_f64() - 0.5).abs() < 1e-15);
1312    }
1313
1314    #[test]
1315    fn testerror_handling() {
1316        // Division by zero
1317        let zero = ArbitraryRational::new();
1318        assert!(zero.recip().is_err());
1319
1320        // Square root of negative
1321        let neg = ArbitraryFloat::from_f64(-1.0);
1322        assert!(neg.sqrt().is_err());
1323
1324        // Logarithm of negative
1325        assert!(neg.ln().is_err());
1326
1327        // Arcsine out of range
1328        let out_of_range = ArbitraryFloat::from_f64(2.0);
1329        assert!(out_of_range.asin().is_err());
1330    }
1331}
scirs2_core/numeric/arbitrary_precision.rs

scirs2_core/numeric/
arbitrary_precision.rs
