num_irrational/quadratic/
surd.rs

1//! Implementation of qudratic irrational numbers
2
3use super::{QuadraticBase, QuadraticOps};
4use crate::cont_frac::ContinuedFraction;
5use crate::traits::{
6    Approximation, Computable, FromSqrt, FromSqrtError, SqrtErrorKind, WithSigned, WithUnsigned,
7};
8#[cfg(feature = "complex")]
9use crate::GaussianInt;
10use crate::QuadraticInt;
11use core::convert::TryFrom;
12use core::ops::{Add, AddAssign, Div, Mul, Neg, Sub};
13use num_integer::{sqrt, Integer};
14use num_traits::{
15    CheckedAdd, CheckedMul, FromPrimitive, NumRef, One, RefNum, Signed, ToPrimitive, Zero,
16};
17use std::fmt;
18
19use num_rational::Ratio;
20
21/// Underlying representation of a quadratic surd `(a + b*√r) / c` as (a,b,c).
22///
23/// Note that the representation won't be reduced (normalized) after operations.
24/// Therefore, the equality test will requires the coefficients to be completely same,
25/// rather than be equal numerically.
26#[derive(Clone, Copy, Debug, Hash, PartialEq, Eq)]
27pub struct QuadraticSurdCoeffs<T>(pub T, pub T, pub T);
28
29impl<T> From<(T, T, T)> for QuadraticSurdCoeffs<T> {
30    fn from(v: (T, T, T)) -> Self {
31        Self(v.0, v.1, v.2)
32    }
33}
34
35impl<T: QuadraticBase> Add for QuadraticSurdCoeffs<T>
36where
37    for<'r> &'r T: RefNum<T>,
38{
39    type Output = Self;
40    #[inline]
41    fn add(self, rhs: Self) -> Self::Output {
42        if self.2 == rhs.2 {
43            Self(self.0 + rhs.0, self.1 + rhs.1, rhs.2)
44        } else {
45            Self(
46                self.0 * &rhs.2 + rhs.0 * &self.2,
47                self.1 * &rhs.2 + rhs.1 * &self.2,
48                rhs.2 * self.2,
49            )
50        }
51    }
52}
53
54impl<T: QuadraticBase> Sub for QuadraticSurdCoeffs<T>
55where
56    for<'r> &'r T: RefNum<T>,
57{
58    type Output = Self;
59    #[inline]
60    fn sub(self, rhs: Self) -> Self::Output {
61        if self.2 == rhs.2 {
62            Self(self.0 - rhs.0, self.1 - rhs.1, rhs.2)
63        } else {
64            Self(
65                self.0 * &rhs.2 - rhs.0 * &self.2,
66                self.1 * &rhs.2 - rhs.1 * &self.2,
67                rhs.2 * self.2,
68            )
69        }
70    }
71}
72
73impl<T: QuadraticBase> QuadraticOps<Self, &T, Self> for QuadraticSurdCoeffs<T>
74where
75    for<'r> &'r T: RefNum<T>,
76{
77    type Scalar = Ratio<T>;
78    fn mul(self, rhs: Self, discr: &T) -> Self {
79        let Self(la, lb, lc) = self;
80        let Self(ra, rb, rc) = rhs;
81        Self(&la * &ra + &lb * &rb * discr, la * rb + lb * ra, lc * rc)
82    }
83    fn div(self, rhs: Self, discr: &T) -> Self {
84        let Self(la, lb, lc) = self;
85        let Self(ra, rb, rc) = rhs;
86        let c = lc * (&ra * &ra - &rb * &rb * discr);
87        Self(
88            &rc * (&la * &ra - &lb * &rb * discr),
89            rc * (lb * ra - la * rb),
90            c,
91        )
92    }
93    #[inline]
94    fn conj(self, _: &T) -> Self {
95        Self(self.0, -self.1, self.2)
96    }
97    #[inline]
98    fn norm(self, discr: &T) -> Self::Scalar {
99        Ratio::new(
100            &self.0 * &self.0 - &self.1 * &self.1 * discr,
101            &self.2 * &self.2,
102        )
103    }
104}
105
106impl<T: QuadraticBase> Add<Ratio<T>> for QuadraticSurdCoeffs<T> {
107    type Output = Self;
108    #[inline]
109    fn add(self, rhs: Ratio<T>) -> Self {
110        let Self(la, lb, lc) = self;
111        let (ra, rc) = rhs.into();
112        if lc == rc {
113            Self(la + ra, lb, lc)
114        } else {
115            Self(la * &rc + ra * &lc, lb * &rc, lc * rc)
116        }
117    }
118}
119
120impl<T: QuadraticBase> Sub<Ratio<T>> for QuadraticSurdCoeffs<T> {
121    type Output = Self;
122    #[inline]
123    fn sub(self, rhs: Ratio<T>) -> Self {
124        let Self(la, lb, lc) = self;
125        let (ra, rc) = rhs.into();
126        if lc == rc {
127            Self(la - ra, lb, lc)
128        } else {
129            Self(la * &rc - ra * &lc, lb * &rc, lc * rc)
130        }
131    }
132}
133
134impl<T: QuadraticBase> Mul<Ratio<T>> for QuadraticSurdCoeffs<T> {
135    type Output = Self;
136    fn mul(self, rhs: Ratio<T>) -> Self {
137        let (ra, rc) = rhs.into();
138        Self(self.0 * &ra, self.1 * ra, self.2 * rc)
139    }
140}
141
142impl<T: QuadraticBase> Div<Ratio<T>> for QuadraticSurdCoeffs<T> {
143    type Output = Self;
144    fn div(self, rhs: Ratio<T>) -> Self {
145        let (ra, rc) = rhs.into();
146        Self(self.0 * &rc, self.1 * rc, self.2 * ra)
147    }
148}
149
150impl<T: QuadraticBase> Add<T> for QuadraticSurdCoeffs<T> {
151    type Output = Self;
152    fn add(self, rhs: T) -> Self::Output {
153        Self(self.0 + rhs * &self.2, self.1, self.2)
154    }
155}
156
157impl<T: QuadraticBase> Sub<T> for QuadraticSurdCoeffs<T> {
158    type Output = Self;
159    fn sub(self, rhs: T) -> Self::Output {
160        Self(self.0 - rhs * &self.2, self.1, self.2)
161    }
162}
163
164impl<T: QuadraticBase> Mul<T> for QuadraticSurdCoeffs<T> {
165    type Output = Self;
166    fn mul(self, rhs: T) -> Self::Output {
167        Self(self.0 * &rhs, self.1 * rhs, self.2)
168    }
169}
170
171impl<T: QuadraticBase> Div<T> for QuadraticSurdCoeffs<T> {
172    type Output = Self;
173    fn div(self, rhs: T) -> Self::Output {
174        Self(self.0, self.1, self.2 * rhs)
175    }
176}
177
178/// A quadratic number represented as `(a + b*√r) / c`.
179/// If the support for complex number is enabled, then this struct can represent
180/// any quadratic integers.
181///
182/// # About the `complex` feature
183///
184/// Whether enabling the `complex` feature of the crate will lead to some behavior changes
185/// of this struct (in a compatible way), some important ones are listed below:
186/// - The [new()][QuadraticSurd::new] constructor will panic if `complex` is disabled and the given root is negative
187/// - Rounding functions ([fract()][QuadraticSurd::fract], [trunc()][QuadraticSurd::trunc], ..) will round to Gaussian
188///   integers instead of rational integers if `complex` is enabled.
189/// - Rational approximation ([approximated()][QuadraticSurd::approximated]) will panic if `complex` is enabled and
190///   the quadratic surd is complex
191/// - [QuadraticSurd::from_sqrt()] won't fail as [SqrtErrorKind::Complex] for negative input if `complex` is enabled
192///
193#[derive(Hash, Clone, Debug, Copy)]
194pub struct QuadraticSurd<T> {
195    coeffs: QuadraticSurdCoeffs<T>, // when reduced, coeffs.1 is zero if the surd is rational, coeffs.2 is always positive
196    discr: T,                       // when reduced, discriminant is zero if the surd is rational
197}
198
199impl<T> QuadraticSurd<T> {
200    #[inline]
201    pub(crate) const fn new_raw(a: T, b: T, c: T, r: T) -> Self {
202        Self {
203            coeffs: QuadraticSurdCoeffs(a, b, c),
204            discr: r,
205        }
206    }
207
208    /// Get return-only references to the components `(a, b, c, r)`
209    pub const fn parts(&self) -> (&T, &T, &T, &T) {
210        (&self.coeffs.0, &self.coeffs.1, &self.coeffs.2, &self.discr)
211    }
212}
213
214impl<T: Integer> QuadraticSurd<T> {
215    /// Determine if the surd is an (rational) integer
216    #[inline]
217    pub fn is_integer(&self) -> bool {
218        self.coeffs.2.is_one() && self.is_rational()
219    }
220
221    /// Determine if the surd is a Gaussian integer
222    #[inline]
223    #[cfg(feature = "complex")]
224    pub fn is_gaussint(&self) -> bool {
225        self.coeffs.2.is_one()
226    }
227
228    /// Determine if the surd is a quadratic integer
229    #[inline]
230    pub fn is_quadint(&self) -> bool {
231        self.coeffs.2.is_one()
232    }
233
234    /// Determine if the surd is a rational number
235    #[inline]
236    pub fn is_rational(&self) -> bool {
237        self.coeffs.1.is_zero() || self.discr.is_zero()
238    }
239
240    /// Determine if the surd is a complex number
241    #[inline]
242    #[cfg(feature = "complex")]
243    pub fn is_complex(&self) -> bool {
244        self.discr < T::zero()
245    }
246
247    /// Determine if the quadratic number has no rational part (i.e. a = 0, b != 0)
248    #[inline]
249    pub fn is_pure(&self) -> bool {
250        self.coeffs.0.is_zero() && !self.coeffs.1.is_zero()
251    }
252
253    /// Panic when the number is complex and the feature "complex" is disabled
254    #[inline(always)]
255    fn panic_if_complex(&self) {
256        // only need to check when we allow construction of complex surd number
257        #[cfg(feature = "complex")]
258        if self.discr < T::zero() {
259            panic!("operation not supported on complex numbers");
260        }
261    }
262}
263
264impl<T: Integer> From<T> for QuadraticSurd<T> {
265    /// Create a `QuadraticSurd` representation of an integer.
266    /// The square root base will be zero.
267    #[inline]
268    fn from(t: T) -> Self {
269        Self {
270            coeffs: QuadraticSurdCoeffs(t, T::zero(), T::one()),
271            discr: T::zero(),
272        }
273    }
274}
275
276impl<T: Integer> From<Ratio<T>> for QuadraticSurd<T> {
277    /// Create a `QuadraticSurd` representation of a rational number.
278    /// The square root base will be zero.
279    #[inline]
280    fn from(t: Ratio<T>) -> Self {
281        let (a, c) = t.into();
282        Self {
283            coeffs: QuadraticSurdCoeffs(a, T::zero(), c),
284            discr: T::zero(),
285        }
286    }
287}
288
289impl<T: QuadraticBase> QuadraticSurd<T>
290where
291    for<'r> &'r T: RefNum<T>,
292{
293    // Simplify the surd into normalized form
294    fn reduce(&mut self) {
295        if self.coeffs.2.is_zero() {
296            panic!("denominator is zero");
297        }
298
299        // test and reduce if the surd is rational
300        let root = sqrt(self.discr.abs());
301        if &root * &root == self.discr {
302            if self.discr.is_negative() {
303                self.coeffs.1 = &self.coeffs.1 * root;
304                self.discr = -T::one();
305            } else {
306                self.coeffs.0 = &self.coeffs.0 + &self.coeffs.1 * root;
307                self.coeffs.1 = T::zero();
308                self.discr = T::zero();
309            }
310        }
311
312        if self.coeffs.1.is_zero() || self.discr.is_zero() {
313            // shortcut if the surd is rational
314            self.discr = T::zero();
315            self.coeffs.1 = T::zero();
316
317            let g_ac = self.coeffs.0.gcd(&self.coeffs.2);
318            self.coeffs.0 = &self.coeffs.0 / &g_ac;
319            self.coeffs.2 = &self.coeffs.2 / g_ac;
320        } else {
321            // reduce common divisor
322            let mut g_ac = self.coeffs.0.gcd(&self.coeffs.2);
323            let g_acr = (&g_ac * &g_ac).gcd(&self.discr); // test if squared factor of c can divide r
324            let groot = sqrt(g_acr.clone());
325            if &groot * &groot == g_acr {
326                self.discr = &self.discr / &g_acr;
327                self.coeffs.0 = &self.coeffs.0 / &groot;
328                self.coeffs.2 = &self.coeffs.2 / &groot;
329                g_ac = &g_ac / groot;
330            }
331
332            let g_abc = g_ac.gcd(&self.coeffs.1);
333            self.coeffs.0 = &self.coeffs.0 / &g_abc;
334            self.coeffs.1 = &self.coeffs.1 / &g_abc;
335            self.coeffs.2 = &self.coeffs.2 / g_abc;
336        }
337
338        // keep denom positive
339        if self.coeffs.2 < T::zero() {
340            self.coeffs.0 = T::zero() - &self.coeffs.0;
341            self.coeffs.1 = T::zero() - &self.coeffs.1;
342            self.coeffs.2 = T::zero() - &self.coeffs.2;
343        }
344    }
345
346    // Try to reduce the root base with a possible factor.
347    fn reduce_root_hinted(self, hint: T) -> Self {
348        let hint = hint.abs();
349        if hint.is_zero() || hint.is_one() || hint.is_negative() {
350            return self;
351        }
352
353        let (quo, rem) = self.discr.div_rem(&hint);
354        if rem.is_zero() {
355            // if hint is actually a factor
356            let root = sqrt(hint.clone());
357            if &root * &root == hint {
358                // if hint is a square number
359                let g = root.gcd(&self.coeffs.0).gcd(&self.coeffs.2);
360                return QuadraticSurd::new_raw(
361                    self.coeffs.0 / &g,
362                    self.coeffs.1 * root / &g,
363                    self.coeffs.2 / g,
364                    quo,
365                );
366            }
367        }
368
369        self
370    }
371
372    /// Create a surd represented as `(a + b√r)) / c` where `a`, `b`, `c`, `r` are integers.
373    ///
374    /// # Panics
375    /// If `c` is zero or `r` is negative when the `complex` feature is not enabled.
376    #[inline]
377    pub fn new(a: T, b: T, c: T, r: T) -> Self {
378        #[cfg(not(feature = "complex"))]
379        if r.is_negative() {
380            panic!("Negative root is not supported without the `complex` feature");
381        }
382
383        let mut ret = QuadraticSurd::new_raw(a, b, c, r);
384        ret.reduce();
385        ret
386    }
387
388    #[inline]
389    pub fn from_coeffs(coeffs: QuadraticSurdCoeffs<T>, r: T) -> Self {
390        #[cfg(not(feature = "complex"))]
391        if r.is_negative() {
392            panic!("Negative root is not supported without the `complex` feature");
393        }
394
395        let mut ret = Self { coeffs, discr: r };
396        ret.reduce();
397        ret
398    }
399
400    /// Create a surd represented as `a + b√r` where a, b, r are rationals.
401    ///
402    /// # Panics
403    /// If `r` is negative when the `complex` feature is not enabled.
404    #[inline]
405    pub fn from_rationals(a: Ratio<T>, b: Ratio<T>, r: Ratio<T>) -> Self {
406        #[cfg(not(feature = "complex"))]
407        if r.is_negative() {
408            panic!("Negative root is not supported without the `complex` feature");
409        }
410
411        let surd_r = r.numer() * r.denom();
412        let new_b_denom = b.denom() * r.denom();
413        let surd_a = a.numer() * &new_b_denom;
414        let surd_b = b.numer() * a.denom();
415        let surd_c = a.denom() * new_b_denom;
416        let mut ret = QuadraticSurd::new_raw(surd_a, surd_b, surd_c, surd_r);
417        ret.reduce();
418        ret
419    }
420
421    /// Get the root of a quadratic equation `ax^2 + bx + c` with rational coefficients.
422    /// This method only returns one of the root `(b + √(b^2 - 4ac)) / 2a`, use `conjugate()` for the other root
423    /// If there are only complex solutions, then `None` will be returned if the `complex` feature is not enabled.
424    #[inline]
425    pub fn from_equation(a: Ratio<T>, b: Ratio<T>, c: Ratio<T>) -> Option<Self> {
426        // degraded cases
427        if a.is_zero() {
428            if b.is_zero() {
429                return None;
430            }
431            return Some(Self::from(-c / b));
432        }
433
434        let two = T::one() + T::one();
435        let four = &two * &two;
436        let delta = &b * &b - Ratio::from(four) * &a * c;
437
438        #[cfg(not(feature = "complex"))]
439        if delta.is_negative() {
440            return None;
441        }
442
443        let aa = Ratio::from(two) * a;
444        Some(Self::from_rationals(-b * aa.recip(), aa.recip(), delta))
445    }
446
447    /// Returns the reciprocal of the surd, TODO: implement as num_traits::Inv trait
448    #[inline]
449    pub fn recip(self) -> Self {
450        let QuadraticSurdCoeffs(a, b, c) = self.coeffs;
451        let aa = &a * &a;
452        let bb = &b * &b;
453        QuadraticSurd::new(-(&c * a), c * b, bb * &self.discr - aa, self.discr)
454    }
455
456    /// `.recip()` with reference
457    #[inline]
458    pub fn recip_ref(&self) -> Self {
459        let QuadraticSurdCoeffs(a, b, c) = &self.coeffs;
460        let aa = a * a;
461        let bb = b * b;
462        QuadraticSurd::new(-(c * a), c * b, bb * &self.discr - aa, self.discr.clone())
463    }
464
465    /// Return the conjugate of the surd, i.e. `(a - b√r) / c`
466    #[inline]
467    pub fn conj(self) -> Self {
468        Self {
469            coeffs: self.coeffs.conj(&self.discr),
470            discr: self.discr,
471        }
472    }
473
474    /// `.conj()` with reference
475    #[inline]
476    pub fn conj_ref(&self) -> Self {
477        self.clone().conj()
478    }
479
480    /// Round the surd toward zero. The result will be an integer if the surd is real,
481    /// or a Gaussian integer if the surd is complex
482    ///
483    /// # Panics
484    /// if the number is complex, when the `complex` feature is not enabled
485    pub fn trunc(self) -> Self {
486        let QuadraticSurdCoeffs(a, b, c) = self.coeffs;
487        let bsign = b.signum();
488
489        if self.discr.is_negative() {
490            let br = bsign * sqrt(&b * &b * -self.discr);
491            Self::from_coeffs(QuadraticSurdCoeffs(a / &c, br / c, T::one()), -T::one())
492        } else {
493            let br = bsign * sqrt(&b * &b * self.discr);
494            return Self::from((a + br) / c);
495        }
496    }
497
498    /// [QuadraticSurd::trunc()] with reference
499    #[inline]
500    pub fn trunc_ref(&self) -> Self {
501        self.clone().trunc()
502    }
503
504    /// Get the fractional part of the surd, ensuring `self.trunc() + self.fract() == self`
505    ///
506    /// # Panics
507    /// If the square root base is negative and not -1 (the result won't be representable by this struct)
508    #[inline]
509    pub fn fract(self) -> Self {
510        let trunc = self.trunc_ref();
511        self - trunc
512    }
513
514    /// [QuadraticSurd::fract()] with reference
515    #[inline]
516    pub fn fract_ref(&self) -> Self {
517        self.clone() - self.trunc_ref()
518    }
519
520    /// Get the numerator of the surd
521    #[inline]
522    pub fn numer(&self) -> Self {
523        Self::new_raw(
524            self.coeffs.0.clone(),
525            self.coeffs.1.clone(),
526            T::one(),
527            self.discr.clone(),
528        )
529    }
530
531    /// Get the denumerator of the surd
532    #[inline]
533    pub fn denom(&self) -> Self {
534        Self::from(self.coeffs.2.clone())
535    }
536
537    /// Converts to an integer, rounding towards zero
538    ///
539    /// # Panics
540    /// if the number is complex when feature `complex` is not enabled
541    #[inline]
542    pub fn to_integer(&self) -> Approximation<T> {
543        self.panic_if_complex();
544
545        if self.is_integer() {
546            Approximation::Exact(self.coeffs.0.clone())
547        } else {
548            Approximation::Approximated(self.trunc_ref().coeffs.0)
549        }
550    }
551
552    /// Converts to an Gaussian integer, rounding towards zero
553    #[inline]
554    #[cfg(feature = "complex")]
555    pub fn to_gaussint(&self) -> Approximation<GaussianInt<T>> {
556        if self.is_gaussint() {
557            Approximation::Exact(GaussianInt::new(
558                self.coeffs.0.clone(),
559                self.coeffs.1.clone(),
560            ))
561        } else {
562            let trunc = self.trunc_ref();
563            Approximation::Approximated(GaussianInt::new(trunc.coeffs.0, trunc.coeffs.1))
564        }
565    }
566
567    /// Converts to an quadratic integer, rounding towards zero
568    #[inline]
569    pub fn to_quadint(&self) -> Approximation<QuadraticInt<T>> {
570        let QuadraticSurdCoeffs(a, b, c) = &self.coeffs;
571        if self.is_quadint() {
572            Approximation::Exact(QuadraticInt::new(a.clone(), b.clone(), self.discr.clone()))
573        } else {
574            Approximation::Approximated(QuadraticInt::new(a / c, b / c, self.discr.clone()))
575        }
576    }
577
578    /// Converts to a rational, rounding square root towards zero
579    ///
580    /// # Panics
581    /// if the number is complex
582    #[inline]
583    pub fn to_rational(&self) -> Approximation<Ratio<T>> {
584        self.panic_if_complex();
585
586        if self.discr.is_zero() {
587            Approximation::Exact(Ratio::new_raw(self.coeffs.0.clone(), self.coeffs.2.clone()))
588        } else {
589            Approximation::Approximated(Ratio::new_raw(
590                &self.coeffs.0 + sqrt(&self.coeffs.1 * &self.coeffs.1 * &self.discr),
591                self.coeffs.2.clone(),
592            ))
593        }
594    }
595
596    /// Rounds towards minus infinity
597    ///
598    /// # Panics
599    /// if the number is complex
600    pub fn floor(self) -> Self {
601        self.panic_if_complex();
602
603        let br = sqrt(&self.coeffs.1 * &self.coeffs.1 * &self.discr);
604        let num = if self.coeffs.1 >= T::zero() {
605            &self.coeffs.0 + br
606        } else {
607            &self.coeffs.0 - br - T::one()
608        };
609        let num = if num >= T::zero() {
610            num
611        } else {
612            num - &self.coeffs.2 + T::one()
613        };
614        return Self::from(num / &self.coeffs.2);
615    }
616
617    /// `.floor()` with reference
618    ///
619    /// # Panics
620    /// if the number is complex
621    #[inline]
622    pub fn floor_ref(&self) -> Self {
623        self.clone().floor()
624    }
625
626    fn ceil(self) -> Self { unimplemented!() }
627
628    /// Returns the nearest integer to a number. Round half-way cases away from zero?.
629    fn round(self) -> Self { unimplemented!() }
630
631    /// Test if the surd number is positive
632    ///
633    /// # Panics
634    /// if the number is complex
635    pub fn is_positive(&self) -> bool {
636        self.panic_if_complex();
637
638        if self.coeffs.1.is_zero() {
639            self.coeffs.0.is_positive()
640        } else if self.coeffs.1.is_positive() {
641            if self.coeffs.0.is_positive() {
642                true
643            } else {
644                &self.coeffs.0 * &self.coeffs.0 < &self.coeffs.1 * &self.coeffs.1 * &self.discr
645            }
646        } else {
647            // self.coeffs.1.is_negative()
648            if !self.coeffs.0.is_positive() {
649                false
650            } else {
651                &self.coeffs.0 * &self.coeffs.0 > &self.coeffs.1 * &self.coeffs.1 * &self.discr
652            }
653        }
654    }
655
656    /// Test if the surd number is negative
657    ///
658    /// # Panics
659    /// if the number is complex
660    pub fn is_negative(&self) -> bool {
661        self.panic_if_complex();
662
663        if self.coeffs.1.is_zero() {
664            self.coeffs.0.is_negative()
665        } else if self.coeffs.1.is_negative() {
666            if self.coeffs.0.is_negative() {
667                true
668            } else {
669                &self.coeffs.0 * &self.coeffs.0 < &self.coeffs.1 * &self.coeffs.1 * &self.discr
670            }
671        } else {
672            // self.coeffs.1.is_positive()
673            if !self.coeffs.0.is_negative() {
674                false
675            } else {
676                &self.coeffs.0 * &self.coeffs.0 > &self.coeffs.1 * &self.coeffs.1 * &self.discr
677            }
678        }
679    }
680}
681
682impl<T: QuadraticBase> PartialEq for QuadraticSurd<T>
683where
684    for<'r> &'r T: RefNum<T>,
685{
686    fn eq(&self, other: &Self) -> bool {
687        // shortcuts
688        if self.discr == other.discr {
689            return self.coeffs == other.coeffs;
690        }
691        if self.discr.is_zero() || other.discr.is_zero() {
692            return false;
693        }
694
695        let QuadraticSurdCoeffs(la, lb, lc) = &self.coeffs;
696        let QuadraticSurdCoeffs(ra, rb, rc) = &other.coeffs;
697
698        // first compare the rational part
699        // FIXME: handle overflow like num-rational
700        if la * rc != ra * lc {
701            return false;
702        }
703        // then compare the quadratic part
704        lb * lb * &self.discr * rc * rc == rb * rb * &other.discr * lc * lc
705    }
706}
707
708impl<T: QuadraticBase> Eq for QuadraticSurd<T> where for<'r> &'r T: RefNum<T> {}
709
710// TODO: implement PartialOrd
711
712impl<T> Into<(T, T, T, T)> for QuadraticSurd<T> {
713    /// Deconstruct the quadratic surd `(a + b√r) / c` into tuple `(a,b,c,r)`
714    fn into(self) -> (T, T, T, T) {
715        (self.coeffs.0, self.coeffs.1, self.coeffs.2, self.discr)
716    }
717}
718
719impl<T: QuadraticBase + FromPrimitive + CheckedMul> QuadraticSurd<T>
720where
721    for<'r> &'r T: RefNum<T>,
722{
723    /// Try to eliminate factors of root that is a squared number.
724    ///
725    /// It will only try trivial division for several small primes.
726    /// For a quadratic surd with large root, consider deconstructing
727    /// the surd by `.into()` and then reduce yourself.
728    fn reduce_root(&mut self) {
729        const SMALL_PRIMES: [u8; 54] = [
730            2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
731            89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
732            181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
733        ];
734        for p in SMALL_PRIMES {
735            let p = T::from_u8(p).unwrap();
736            if let Some(p2) = p.checked_mul(&p) {
737                loop {
738                    let (quo, rem) = self.discr.div_rem(&p2);
739                    if rem.is_zero() {
740                        self.discr = quo;
741                        self.coeffs.1 = &self.coeffs.1 * &p;
742                        // common factors between a, b, c will be handled by reduce()
743                    } else {
744                        break;
745                    }
746                }
747            }
748        }
749    }
750
751    /// Returns a reduced version of self.
752    ///
753    /// This method will try to eliminate common divisors between a, b, c and
754    /// also try to eliminate square factors of r.
755    ///
756    #[inline]
757    pub fn reduced(self) -> Self {
758        let mut result = self;
759        result.reduce_root();
760        result.reduce();
761        result
762    }
763}
764
765impl<
766        T: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = U>,
767        U: QuadraticBase,
768    > From<ContinuedFraction<T>> for QuadraticSurd<U>
769where
770    for<'r> &'r T: RefNum<T>,
771    for<'r> &'r U: RefNum<U>,
772{
773    fn from(s: ContinuedFraction<T>) -> Self {
774        // use convergent if it's a rational
775        if s.is_rational() {
776            return Self::from(s.convergents().last().unwrap());
777        }
778
779        // for periodic fraction, assume periodic part x = (ax + b) / (cx + d)
780        let mut piter = s.periodic_coeffs().iter().rev();
781        let (mut a, mut b) = (T::zero(), T::one());
782        let (mut c, mut d) = (T::one(), piter.next().unwrap().clone());
783
784        while let Some(v) = piter.next() {
785            let new_c = v * &c + a;
786            let new_d = v * &d + b;
787            a = c;
788            b = d;
789            c = new_c;
790            d = new_d;
791        }
792
793        let psurd = Self::from_equation(
794            Ratio::from(c.to_signed()),
795            Ratio::from(d.to_signed() - a.to_signed()),
796            Ratio::from(-b.to_signed()),
797        )
798        .unwrap();
799
800        // apply the aperiodic part, assume again result = (ax + b) / (cx + d)
801        let surd = if s.aperiodic_coeffs().len() == 1 {
802            psurd
803        } else {
804            let mut aiter = s.aperiodic_coeffs().iter().skip(1).rev();
805            let (mut a, mut b) = (T::zero(), T::one());
806            let (mut c, mut d) = (T::one(), aiter.next().unwrap().clone());
807
808            while let Some(v) = aiter.next() {
809                let new_c = v * &c + a;
810                let new_d = v * &d + b;
811                a = c;
812                b = d;
813                c = new_c;
814                d = new_d;
815            }
816
817            (psurd.clone() * a.to_signed() + b.to_signed())
818                / (psurd * c.to_signed() + d.to_signed())
819        };
820        let surd = surd + s.aperiodic_coeffs().first().unwrap().clone().to_signed();
821
822        // apply sign
823        if s.is_negative() {
824            -surd
825        } else {
826            surd
827        }
828    }
829}
830
831// Reference: http://www.numbertheory.org/courses/MP313/lectures/lecture17/page5.html
832//            http://www.numbertheory.org/gnubc/surd
833// assumes positive surd, parameter `neg` determines the sign of the fraction
834fn quadsurd_to_contfrac<T: QuadraticBase + WithUnsigned<Unsigned = U> + AddAssign, U>(
835    a: T,
836    b: T,
837    c: T,
838    r: T,
839    neg: bool,
840) -> ContinuedFraction<U>
841where
842    for<'r> &'r T: RefNum<T>,
843{
844    debug_assert!(!c.is_zero() && !r.is_negative());
845    debug_assert!(r.sqrt() * r.sqrt() != r);
846
847    // convert to form (p+√d)/q where q|d-p^2
848    let mut d = r * &b * &b;
849    let (mut p, mut q) = if b.is_negative() { (-a, -c) } else { (a, c) };
850    if (&d - &p * &p) % &q != T::zero() {
851        d = d * &q * &q;
852        p = p * &q;
853        q = &q * &q;
854    }
855
856    // find the reduced form and aperiodic coefficients
857    let mut a_coeffs: Vec<T> = Vec::new();
858    let rd = d.sqrt();
859    while a_coeffs.len() == 0 || // ensure that we have a first coefficient
860          !(p <= rd && rd < (&p+&q) && (&q-&p) <= rd)
861    {
862        let a = (&rd + &p).div_floor(&q);
863        p = &a * &q - p;
864        q = (&d - &p * &p) / q;
865        a_coeffs.push(a);
866    }
867
868    // find the periodic coefficients
869    let mut p_coeffs: Vec<T> = Vec::new();
870    let (init_p, init_q) = (p.clone(), q.clone());
871    loop {
872        let a = (&rd + &p).div_floor(&q);
873        p = &a * &q - p;
874        q = (&d - &p * &p) / q;
875        p_coeffs.push(a);
876        if p == init_p && q == init_q {
877            break;
878        }
879    }
880
881    ContinuedFraction::new(a_coeffs, p_coeffs, neg)
882}
883
884// Although conversion to continued fraction won't fail if complex number is disabled,
885// we have to disable From<> to make sure that adding the feature 'complex' won't result
886// in break change
887impl<T: QuadraticBase + WithUnsigned<Unsigned = U> + AddAssign, U> TryFrom<QuadraticSurd<T>>
888    for ContinuedFraction<U>
889where
890    for<'r> &'r T: RefNum<T>,
891{
892    type Error = ();
893
894    fn try_from(s: QuadraticSurd<T>) -> Result<Self, ()> {
895        if s.discr.is_negative() {
896            Err(())
897        } else {
898            if s.is_negative() {
899                Ok(quadsurd_to_contfrac(
900                    -s.coeffs.0,
901                    -s.coeffs.1,
902                    s.coeffs.2,
903                    s.discr,
904                    true,
905                ))
906            } else {
907                Ok(quadsurd_to_contfrac(
908                    s.coeffs.0, s.coeffs.1, s.coeffs.2, s.discr, false,
909                ))
910            }
911        }
912    }
913}
914
915// TODO: support convert to hurwitz continued fraction based on gaussian integers if `complex` is enabled
916//       this should be implemented as From<QuadraticSurd<T>> for HurwitzContinuedFraction (instead of TryFrom)
917
918impl<
919        T: QuadraticBase + CheckedAdd + AddAssign + WithUnsigned<Unsigned = U>,
920        U: Integer + Clone + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
921    > Computable<T> for QuadraticSurd<T>
922where
923    for<'r> &'r T: RefNum<T>,
924{
925    fn approximated(&self, limit: &T) -> Approximation<Ratio<T>> {
926        ContinuedFraction::<U>::try_from(self.clone())
927            .expect("only real numbers can be approximated by a rational number")
928            .approximated(limit)
929    }
930}
931
932impl<T: Integer + Signed + fmt::Display + Clone> fmt::Display for QuadraticSurd<T> {
933    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
934        let QuadraticSurdCoeffs(a, b, c) = &self.coeffs;
935        if f.alternate() && self.discr.is_negative() {
936            // print √-1 as i, √-5 as √5i if alternate flag is set
937            // XXX: refactor this function
938            let r = -self.discr.clone();
939            match (
940                a.is_zero(),
941                b.is_zero(),
942                b.is_one(),
943                b == &-T::one(),
944                c.is_one(),
945                r.is_one(),
946            ) {
947                (true, true, _, _, _, _) => write!(f, "0"),
948                (true, false, true, _, true, true) => write!(f, "i"),
949                (true, false, true, _, true, false) => write!(f, "√{}i", r),
950                (true, false, true, _, false, true) => write!(f, "i/{}", c),
951                (true, false, true, _, false, false) => write!(f, "√{}i/{}", r, c),
952                (true, false, false, true, true, true) => write!(f, "-i"),
953                (true, false, false, true, true, false) => write!(f, "-√{}i", r),
954                (true, false, false, true, false, true) => write!(f, "-i/{}", c),
955                (true, false, false, true, false, false) => write!(f, "-√{}i/{}", r, c),
956                (true, false, false, false, true, true) => write!(f, "{}i", b),
957                (true, false, false, false, true, false) => write!(f, "{}√{}i", b, r),
958                (true, false, false, false, false, true) => write!(f, "{}i/{}", b, c),
959                (true, false, false, false, false, false) => write!(f, "{}√{}i/{}", b, r, c),
960                (false, true, _, _, true, _) => write!(f, "{}", a),
961                (false, true, _, _, false, _) => write!(f, "{}/{}", a, c),
962                (false, false, true, _, true, true) => write!(f, "{}+i", a),
963                (false, false, true, _, true, false) => write!(f, "{}+√{}i", a, r),
964                (false, false, false, true, true, true) => write!(f, "{}-i", a),
965                (false, false, false, true, true, false) => write!(f, "{}-√{}i", a, r),
966                (false, false, false, false, true, true) => {
967                    if b.is_negative() {
968                        write!(f, "{}{}i", a, b)
969                    } else {
970                        write!(f, "{}+{}i", a, b)
971                    }
972                },
973                (false, false, false, false, true, false) => {
974                    if b.is_negative() {
975                        write!(f, "{}{}√{}i", a, b, r)
976                    } else {
977                        write!(f, "{}+{}√{}i", a, b, r)
978                    }
979                }
980                (false, false, true, _, false, true) => write!(f, "({}+i)/{}", a, c),
981                (false, false, true, _, false, false) => write!(f, "({}+√{}i)/{}", a, r, c),
982                (false, false, false, true, false, true) => write!(f, "({}-i)/{}", a, c),
983                (false, false, false, true, false, false) => write!(f, "({}-√{}i)/{}", a, r, c),
984                (false, false, false, false, false, true) => {
985                    if b.is_negative() {
986                        write!(f, "({}{}i)/{}", a, b, c)
987                    } else {
988                        write!(f, "({}+{}i)/{}", a, b, c)
989                    }
990                },
991                (false, false, false, false, false, false) => {
992                    if b.is_negative() {
993                        write!(f, "({}{}√{}i)/{}", a, b, r, c)
994                    } else {
995                        write!(f, "({}+{}√{}i)/{}", a, b, r, c)
996                    }
997                }
998            }
999        } else {
1000            match (
1001                a.is_zero(),
1002                b.is_zero(),
1003                b.is_one(),
1004                b == &-T::one(),
1005                c.is_one(),
1006            ) {
1007                (true, true, _, _, _) => write!(f, "0"),
1008                (true, false, true, _, true) => write!(f, "√{}", self.discr),
1009                (true, false, true, _, false) => write!(f, "√{}/{}", self.discr, c),
1010                (true, false, false, true, true) => write!(f, "-√{}", self.discr),
1011                (true, false, false, true, false) => write!(f, "-√{}/{}", self.discr, c),
1012                (true, false, false, false, true) => write!(f, "{}√{}", b, self.discr),
1013                (true, false, false, false, false) => write!(f, "{}√{}/{}", b, self.discr, c),
1014                (false, true, _, _, true) => write!(f, "{}", a),
1015                (false, true, _, _, false) => write!(f, "{}/{}", a, c),
1016                (false, false, true, _, true) => write!(f, "{}+√{}", a, self.discr),
1017                (false, false, false, true, true) => write!(f, "{}-√{}", a, self.discr),
1018                (false, false, false, false, true) => {
1019                    if b.is_negative() {
1020                        write!(f, "{}{}√{}", a, b, self.discr)
1021                    } else {
1022                        write!(f, "{}+{}√{}", a, b, self.discr)
1023                    }
1024                }
1025                (false, false, true, _, false) => write!(f, "({}+√{})/{}", a, self.discr, c),
1026                (false, false, false, true, false) => {
1027                    write!(f, "({}-√{})/{}", a, self.discr, c)
1028                }
1029                (false, false, false, false, false) => {
1030                    if b.is_negative() {
1031                        write!(f, "({}{}√{})/{}", a, b, self.discr, c)
1032                    } else {
1033                        write!(f, "({}+{}√{})/{}", a, b, self.discr, c)
1034                    }
1035                }
1036            }
1037        }
1038    }
1039}
1040
1041// Reduce root base for binary operands. Return None if the two bases
1042// cannot be matched. This function assumes the root bases on both sides
1043// are not zero
1044#[inline]
1045#[cfg(not(feature = "complex"))]
1046fn reduce_bin_op<T: QuadraticBase>(
1047    lhs: QuadraticSurd<T>,
1048    rhs: QuadraticSurd<T>,
1049) -> Option<(QuadraticSurd<T>, QuadraticSurd<T>)>
1050where
1051    for<'r> &'r T: RefNum<T>,
1052{
1053    debug_assert!(!lhs.discr.is_zero());
1054    debug_assert!(!rhs.discr.is_zero());
1055
1056    let result = if lhs.discr > rhs.discr {
1057        let hint = &lhs.discr / &rhs.discr;
1058        (lhs.reduce_root_hinted(hint), rhs)
1059    } else if rhs.discr > lhs.discr {
1060        let hint = &rhs.discr / &lhs.discr;
1061        (lhs, rhs.reduce_root_hinted(hint))
1062    } else {
1063        (lhs, rhs)
1064    };
1065
1066    if result.0.discr == result.1.discr {
1067        Some(result)
1068    } else {
1069        None
1070    }
1071}
1072
1073#[inline]
1074#[cfg(feature = "complex")]
1075fn reduce_bin_op<T: QuadraticBase>(
1076    lhs: QuadraticSurd<T>,
1077    rhs: QuadraticSurd<T>,
1078) -> Option<(QuadraticSurd<T>, QuadraticSurd<T>)>
1079where
1080    for<'r> &'r T: RefNum<T>,
1081{
1082    if lhs.discr.is_negative() ^ rhs.discr.is_negative() {
1083        return None;
1084    }
1085
1086    let result = if lhs.discr.abs() > rhs.discr.abs() {
1087        let hint = &lhs.discr / &rhs.discr;
1088        (lhs.reduce_root_hinted(hint), rhs)
1089    } else if rhs.discr.abs() > lhs.discr.abs() {
1090        let hint = &rhs.discr / &lhs.discr;
1091        (lhs, rhs.reduce_root_hinted(hint))
1092    } else {
1093        (lhs, rhs)
1094    };
1095
1096    if result.0.discr == result.1.discr {
1097        Some(result)
1098    } else {
1099        None
1100    }
1101}
1102
1103#[inline]
1104fn reduce_bin_op_unwrap<T: QuadraticBase>(
1105    lhs: QuadraticSurd<T>,
1106    rhs: QuadraticSurd<T>,
1107) -> (QuadraticSurd<T>, QuadraticSurd<T>)
1108where
1109    for<'r> &'r T: RefNum<T>,
1110{
1111    reduce_bin_op(lhs, rhs).expect("two root bases are not compatible!")
1112}
1113
1114macro_rules! forward_binop {
1115    (impl $imp:ident, $rhs:ty, $method:ident) => {
1116        impl<T: QuadraticBase> $imp<$rhs> for QuadraticSurd<T>
1117        where
1118            for<'r> &'r T: RefNum<T>,
1119        {
1120            type Output = Self;
1121            #[inline]
1122            fn $method(self, rhs: $rhs) -> Self {
1123                Self::from_coeffs($imp::$method(self.coeffs, rhs), self.discr)
1124            }
1125        }
1126    };
1127}
1128
1129forward_binop!(impl Add, T, add);
1130forward_binop!(impl Add, Ratio<T>, add);
1131forward_binop!(impl Sub, T, sub);
1132forward_binop!(impl Sub, Ratio<T>, sub);
1133forward_binop!(impl Mul, T, mul);
1134forward_binop!(impl Mul, Ratio<T>, mul);
1135forward_binop!(impl Div, T, div);
1136forward_binop!(impl Div, Ratio<T>, div);
1137
1138// TODO: implement checked_add, checked_sub, checked_mul, checked_div, checking both overflow and discr dismatch
1139impl<T: QuadraticBase> Add for QuadraticSurd<T>
1140where
1141    for<'r> &'r T: RefNum<T>,
1142{
1143    type Output = Self;
1144    fn add(self, rhs: Self) -> Self {
1145        // shortcuts for trivial cases
1146        if rhs.is_rational() {
1147            return self + Ratio::<T>::new(rhs.coeffs.0, rhs.coeffs.2);
1148        }
1149        if self.is_rational() {
1150            return rhs + Ratio::<T>::new(self.coeffs.0, self.coeffs.2);
1151        }
1152
1153        // ensure that two operands have compatible bases
1154        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
1155        Self::from_coeffs(lhs.coeffs + rhs.coeffs, lhs.discr)
1156    }
1157}
1158
1159impl<T: QuadraticBase> Sub for QuadraticSurd<T>
1160where
1161    for<'r> &'r T: RefNum<T>,
1162{
1163    type Output = Self;
1164    fn sub(self, rhs: Self) -> Self {
1165        // shortcuts for trivial cases
1166        if rhs.is_rational() {
1167            return self - Ratio::<T>::new(rhs.coeffs.0, rhs.coeffs.2);
1168        }
1169        if self.is_rational() {
1170            return -(rhs - Ratio::<T>::new(self.coeffs.0, self.coeffs.2));
1171        }
1172
1173        // ensure that two operands have compatible bases
1174        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
1175        Self::from_coeffs(lhs.coeffs - rhs.coeffs, lhs.discr)
1176    }
1177}
1178
1179impl<T: QuadraticBase> Mul for QuadraticSurd<T>
1180where
1181    for<'r> &'r T: RefNum<T>,
1182{
1183    type Output = Self;
1184    #[inline]
1185    fn mul(self, rhs: Self) -> Self {
1186        // shortcuts for trivial cases
1187        if rhs.is_rational() {
1188            return self * Ratio::<T>::new(rhs.coeffs.0, rhs.coeffs.2);
1189        }
1190        if self.is_rational() {
1191            return rhs * Ratio::<T>::new(self.coeffs.0, self.coeffs.2);
1192        }
1193        if self.is_pure() && rhs.is_pure() {
1194            let gcd = self.discr.gcd(&rhs.discr);
1195            let discr = (self.discr / &gcd) * (rhs.discr / &gcd);
1196            return Self::new(
1197                T::zero(),
1198                self.coeffs.1 * rhs.coeffs.1 * gcd,
1199                self.coeffs.2 * rhs.coeffs.2,
1200                discr,
1201            );
1202        }
1203
1204        // ensure that two operands have compatible bases
1205        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
1206        Self::from_coeffs(
1207            QuadraticOps::mul(lhs.coeffs, rhs.coeffs, &lhs.discr),
1208            lhs.discr,
1209        )
1210    }
1211}
1212
1213impl<T: QuadraticBase> Div for QuadraticSurd<T>
1214where
1215    for<'r> &'r T: RefNum<T>,
1216{
1217    type Output = Self;
1218    #[inline]
1219    fn div(self, rhs: Self) -> Self {
1220        // shortcuts for trivial cases
1221        if rhs.is_rational() {
1222            return self / Ratio::<T>::new(rhs.coeffs.0, rhs.coeffs.2);
1223        }
1224        if self.is_rational() {
1225            return (rhs / Ratio::<T>::new(self.coeffs.0, self.coeffs.2)).recip();
1226        }
1227        if self.is_pure() && rhs.is_pure() {
1228            let gcd = self.discr.gcd(&rhs.discr);
1229            let (ld, rd) = (self.discr / &gcd, rhs.discr / gcd);
1230            return Self::new(
1231                T::zero(),
1232                self.coeffs.1 * rhs.coeffs.2,
1233                self.coeffs.2 * rhs.coeffs.1 * &rd,
1234                ld * rd,
1235            );
1236        }
1237
1238        // ensure that two operands have compatible bases
1239        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
1240        Self::from_coeffs(
1241            QuadraticOps::div(lhs.coeffs, rhs.coeffs, &lhs.discr),
1242            lhs.discr,
1243        )
1244    }
1245}
1246
1247impl<T: QuadraticBase> Neg for QuadraticSurd<T>
1248where
1249    for<'r> &'r T: RefNum<T>,
1250{
1251    type Output = QuadraticSurd<T>;
1252    #[inline]
1253    fn neg(self) -> QuadraticSurd<T> {
1254        QuadraticSurd::new_raw(-self.coeffs.0, -self.coeffs.1, self.coeffs.2, self.discr)
1255    }
1256}
1257
1258impl<T: QuadraticBase> Zero for QuadraticSurd<T>
1259where
1260    for<'r> &'r T: RefNum<T>,
1261{
1262    #[inline]
1263    fn zero() -> Self {
1264        Self {
1265            coeffs: QuadraticSurdCoeffs(T::zero(), T::zero(), T::one()),
1266            discr: T::zero(),
1267        }
1268    }
1269    #[inline]
1270    fn is_zero(&self) -> bool {
1271        self.coeffs.0.is_zero() && self.coeffs.1.is_zero()
1272    }
1273}
1274
1275impl<T: QuadraticBase> One for QuadraticSurd<T>
1276where
1277    for<'r> &'r T: RefNum<T>,
1278{
1279    #[inline]
1280    fn one() -> Self {
1281        Self {
1282            coeffs: QuadraticSurdCoeffs(T::one(), T::zero(), T::one()),
1283            discr: T::zero(),
1284        }
1285    }
1286    #[inline]
1287    fn is_one(&self) -> bool {
1288        self.coeffs.0.is_one() && self.coeffs.1.is_zero() && self.coeffs.2.is_one()
1289    }
1290}
1291
1292impl<T: Integer + FromPrimitive + Clone> FromPrimitive for QuadraticSurd<T> {
1293    #[inline]
1294    fn from_i64(n: i64) -> Option<Self> {
1295        T::from_i64(n).map(Self::from)
1296    }
1297
1298    #[inline]
1299    fn from_u64(n: u64) -> Option<Self> {
1300        T::from_u64(n).map(Self::from)
1301    }
1302
1303    /// This method depends on [Ratio::from_f64] to convert from a float
1304    #[inline]
1305    fn from_f64(f: f64) -> Option<Self> {
1306        // XXX: this method should be improved when Ratio has a more reasonable API for approximating float
1307        let frac = Ratio::<i64>::from_f64(f)?;
1308        let (n, d) = frac.into();
1309        let frac = Ratio::new(T::from_i64(n)?, T::from_i64(d)?);
1310        Some(Self::from(frac))
1311    }
1312}
1313
1314impl<T: QuadraticBase + ToPrimitive> ToPrimitive for QuadraticSurd<T>
1315where
1316    for<'r> &'r T: RefNum<T>,
1317{
1318    #[inline]
1319    fn to_i64(&self) -> Option<i64> {
1320        if self.discr.is_negative() {
1321            return None;
1322        }
1323        match self.to_integer() {
1324            Approximation::Exact(v) => v.to_i64(),
1325            Approximation::Approximated(_) => None,
1326        }
1327    }
1328
1329    #[inline]
1330    fn to_u64(&self) -> Option<u64> {
1331        if self.discr.is_negative() {
1332            return None;
1333        }
1334        match self.to_integer() {
1335            Approximation::Exact(v) => v.to_u64(),
1336            Approximation::Approximated(_) => None,
1337        }
1338    }
1339
1340    #[inline]
1341    fn to_f64(&self) -> Option<f64> {
1342        if self.discr < T::zero() {
1343            return None;
1344        }
1345        Some(
1346            (self.coeffs.0.to_f64()? + self.coeffs.1.to_f64()? * self.discr.to_f64()?.sqrt())
1347                / self.coeffs.2.to_f64()?,
1348        )
1349    }
1350}
1351
1352#[cfg(feature = "num-complex")]
1353mod complex {
1354    use super::*;
1355    use num_complex::{Complex32, Complex64};
1356
1357    impl<T: QuadraticBase + ToPrimitive> QuadraticSurd<T>
1358    where
1359        for<'r> &'r T: RefNum<T>,
1360    {
1361        pub fn to_complex64(&self) -> Option<Complex64> {
1362            if self.discr < T::zero() {
1363                let c = self.coeffs.2.to_f64()?;
1364                let re = self.coeffs.0.to_f64()? / c;
1365                let im = self.coeffs.1.to_f64()? * self.discr.to_f64()?.abs().sqrt() / c;
1366                Some(Complex64::new(re as f64, im as f64))
1367            } else {
1368                let re = self.to_f64()?;
1369                Some(Complex64::new(re as f64, 0f64))
1370            }
1371        }
1372
1373        pub fn to_complex32(&self) -> Option<Complex32> {
1374            let complex = self.to_complex64()?;
1375            Some(Complex32::new(complex.re.to_f32()?, complex.im.to_f32()?))
1376        }
1377    }
1378}
1379
1380impl<T: QuadraticBase> FromSqrt<T> for QuadraticSurd<T>
1381where
1382    for<'r> &'r T: RefNum<T>,
1383{
1384    #[inline]
1385    fn from_sqrt(target: T) -> Result<Self, FromSqrtError<T>> {
1386        #[cfg(not(feature = "complex"))]
1387        if target.is_negative() {
1388            return Err(FromSqrtError {
1389                data: target,
1390                kind: SqrtErrorKind::Complex,
1391            });
1392        }
1393        Ok(QuadraticSurd::new(T::zero(), T::one(), T::one(), target))
1394    }
1395}
1396
1397impl<T: QuadraticBase + CheckedMul> FromSqrt<Ratio<T>> for QuadraticSurd<T>
1398where
1399    for<'r> &'r T: RefNum<T>,
1400{
1401    #[inline]
1402    fn from_sqrt(target: Ratio<T>) -> Result<Self, FromSqrtError<Ratio<T>>> {
1403        #[cfg(not(feature = "complex"))]
1404        if target.is_negative() {
1405            return Err(FromSqrtError {
1406                data: target,
1407                kind: SqrtErrorKind::Complex,
1408            });
1409        }
1410        if target.is_integer() {
1411            let (num, _) = target.into();
1412            return Ok(QuadraticSurd::new(T::zero(), T::one(), T::one(), num));
1413        }
1414
1415        match target.numer().checked_mul(target.denom()) {
1416            Some(new_r) => Ok(QuadraticSurd::new(
1417                T::zero(),
1418                T::one(),
1419                target.denom().clone(),
1420                new_r,
1421            )),
1422            None => Err(FromSqrtError {
1423                data: target,
1424                kind: SqrtErrorKind::Overflow,
1425            }),
1426        }
1427    }
1428}
1429
1430impl<T: QuadraticBase + CheckedMul> FromSqrt<QuadraticSurd<T>> for QuadraticSurd<T>
1431where
1432    for<'r> &'r T: RefNum<T>,
1433{
1434    /// The real part of the output is ensured to be non-negative (if successful).
1435    #[inline]
1436    fn from_sqrt(target: QuadraticSurd<T>) -> Result<Self, FromSqrtError<QuadraticSurd<T>>> {
1437        #[cfg(not(feature = "complex"))]
1438        if target.is_negative() {
1439            return Err(FromSqrtError {
1440                data: target,
1441                kind: SqrtErrorKind::Complex,
1442            });
1443        }
1444
1445        if target.is_rational() {
1446            match target.to_rational() {
1447                Approximation::Exact(v) => {
1448                    return Self::from_sqrt(v).map_err(|e| FromSqrtError {
1449                        data: Self::from(e.data),
1450                        kind: e.kind,
1451                    })
1452                }
1453                _ => unreachable!(),
1454            };
1455        }
1456
1457        // denote a_c = a/c, b_c = b/c
1458        // suppose (a_c + b_c * sqrt(r))^2 = target = x + y * sqrt(r)
1459        // let g = a_c/b_c, then g^2 - 2g(x/y) + r = 0
1460        // result is available only if g is rational
1461
1462        let QuadraticSurdCoeffs(a, b, c) = target.coeffs;
1463        let x = Ratio::new(a, c.clone());
1464        let y = Ratio::new(b, c);
1465        let x_y = x / &y;
1466        let delta2 = &x_y * &x_y - &target.discr;
1467
1468        // reconstruct the original target when error occurred
1469        #[inline]
1470        fn reconstruct<T: QuadraticBase>(x_y: Ratio<T>, y: Ratio<T>, r: T, kind: SqrtErrorKind) -> FromSqrtError<QuadraticSurd<T>> where
1471        for<'r> &'r T: RefNum<T> { 
1472            FromSqrtError {
1473                data: QuadraticSurd::from_rationals(x_y * &y, y, Ratio::from(r)),
1474                kind: kind,
1475            }
1476        }
1477
1478        if delta2.is_negative() {
1479            // this branch happens only when discr is positive
1480            return Err(reconstruct(x_y, y, target.discr, SqrtErrorKind::Unrepresentable));
1481        }
1482        let delta = match Self::from_sqrt(delta2) {
1483            Ok(v) => v,
1484            Err(e) => return Err(reconstruct(x_y, y, target.discr, e.kind))
1485        };
1486        let delta = if delta.is_rational() {
1487            delta.to_rational().value()
1488        } else {
1489            return Err(reconstruct(x_y, y, target.discr, SqrtErrorKind::Unrepresentable))
1490        };
1491
1492        // from the equation above, y = 2*a_c*b_c => c = sqrt(2ab/y)
1493        // this function find a tuple of integers (a, b, c) that satisfies the equation, and
1494        // the output a will be ensured to be non-negative
1495        fn get_abc<T: QuadraticBase> (g: Ratio<T>, y: &Ratio<T>) -> Option<(T, T, T)> where
1496        for<'r> &'r T: RefNum<T> {
1497            let two_ab = (T::one() + T::one()) * g.numer() * g.denom();
1498            let d = two_ab.gcd(y.numer());
1499            let scale = y.numer() / d;
1500            let c2 = two_ab * &scale * &scale * y.denom() / y.numer();
1501            let (a, b) = (g.numer() * &scale, g.denom() * scale);
1502            if c2.is_negative() {
1503                return None;
1504            }
1505            let c = c2.sqrt();
1506            if &c * &c != c2 {
1507                None
1508            } else {
1509                Some((a, b, c))
1510            }
1511        }
1512
1513        let ret = if let Some((a, b, c)) = get_abc(&x_y - &delta, &y) {
1514            QuadraticSurd::new(a, b, c, target.discr)
1515        } else if let Some((a, b, c)) = get_abc(&x_y + delta, &y) {
1516            QuadraticSurd::new(a, b, c, target.discr)
1517        } else {
1518            return Err(reconstruct(x_y, y, target.discr, SqrtErrorKind::Unrepresentable))
1519        };
1520
1521        debug_assert!(!ret.coeffs.0.is_negative());
1522        Ok(ret)
1523    }
1524}
1525
1526#[cfg(test)]
1527mod tests {
1528    use super::*;
1529    use core::convert::TryFrom;
1530
1531    pub const PHI: QuadraticSurd<i32> = QuadraticSurd::new_raw(1, 1, 2, 5); // 1.618
1532    pub const PHI_R: QuadraticSurd<i32> = QuadraticSurd::new_raw(-1, 1, 2, 5); // 0.618
1533    pub const N_PHI: QuadraticSurd<i32> = QuadraticSurd::new_raw(-1, -1, 2, 5); // -1.618
1534    pub const N_PHI_R: QuadraticSurd<i32> = QuadraticSurd::new_raw(1, -1, 2, 5); // -0.618
1535    pub const PHI_SQ: QuadraticSurd<i32> = QuadraticSurd::new_raw(3, 1, 2, 5);
1536
1537    pub const PHI45: QuadraticSurd<i32> = QuadraticSurd::new_raw(3, 1, 6, 45); // non-reduced version of PHI
1538    pub const PHI45_R: QuadraticSurd<i32> = QuadraticSurd::new_raw(-3, 1, 6, 45); // non-reduced version of PHI_R
1539
1540    pub const SQ5: QuadraticSurd<i32> = QuadraticSurd::new_raw(0, 1, 1, 5);
1541    pub const N_SQ5: QuadraticSurd<i32> = QuadraticSurd::new_raw(0, -1, 1, 5);
1542
1543    #[test]
1544    fn creation_test() {
1545        let coeffs: (i32, i32, i32, i32) = QuadraticSurd::new(2, 6, 2, 5).into();
1546        assert_eq!(coeffs, (1, 3, 1, 5)); // reduce common divisors
1547        
1548        let coeffs: (i32, i32, i32, i32) = QuadraticSurd::new(2, 1, 2, 18).into();
1549        assert_eq!(coeffs, (2, 1, 2, 18)); // 18 is not trivially reducible
1550
1551        let coeffs: (i32, i32, i32, i32) = QuadraticSurd::new(3, 1, 3, 18).into();
1552        assert_eq!(coeffs, (1, 1, 1, 2)); // 18 is reducible with the help of gcd hint
1553
1554        let coeffs: (i32, i32, i32, i32) = QuadraticSurd::new(3, 1, 3, 9).into();
1555        assert_eq!(coeffs, (2, 0, 1, 0)); // 9 is a square number
1556    }
1557
1558    #[test]
1559    fn conversion_test() {
1560        assert_eq!(PHI.floor().to_i32(), Some(1));
1561        assert_eq!(PHI_R.floor().to_i32(), Some(0));
1562        assert_eq!(PHI.to_integer(), Approximation::Approximated(1));
1563        assert_eq!(PHI_R.to_integer(), Approximation::Approximated(0));
1564        assert_eq!(N_PHI.floor().to_i32(), Some(-2));
1565        assert_eq!(N_PHI_R.floor().to_i32(), Some(-1));
1566        assert_eq!(N_PHI.to_integer(), Approximation::Approximated(-1));
1567        assert_eq!(N_PHI_R.to_integer(), Approximation::Approximated(0));
1568
1569        assert!(matches!(PHI.to_f64(),     Some(v) if (v - 1.61803398874989f64).abs() < 1e-10));
1570        assert!(matches!(PHI_R.to_f64(),   Some(v) if (v - 0.61803398874989f64).abs() < 1e-10));
1571        assert!(matches!(N_PHI.to_f64(),   Some(v) if (v + 1.61803398874989f64).abs() < 1e-10));
1572        assert!(matches!(N_PHI_R.to_f64(), Some(v) if (v + 0.61803398874989f64).abs() < 1e-10));
1573    }
1574
1575    #[test]
1576    fn from_xxx_test() {
1577        // from_rationals
1578        assert_eq!(
1579            QuadraticSurd::from_rationals(Ratio::new(1, 2), Ratio::new(1, 2), Ratio::from(5)),
1580            PHI
1581        );
1582        assert_eq!(
1583            QuadraticSurd::from_rationals(Ratio::new(-1, 2), Ratio::new(1, 2), Ratio::from(5)),
1584            PHI_R
1585        );
1586
1587        // from_sqrt
1588        assert_eq!(
1589            QuadraticSurd::from_sqrt(5).unwrap(),
1590            QuadraticSurd::new_raw(0, 1, 1, 5)
1591        );
1592        assert!(matches!(QuadraticSurd::from_sqrt(PHI), Err(_)));
1593        assert!(matches!(QuadraticSurd::from_sqrt(PHI * PHI), Ok(v) if v == PHI));
1594
1595        assert!(matches!(
1596            QuadraticSurd::from_equation(Ratio::from(1), Ratio::from(-1), Ratio::from(-1)),
1597            Some(v) if v == PHI
1598        ));
1599    }
1600
1601    #[test]
1602    fn property_test() {
1603        // recip
1604        assert_eq!(PHI.recip(), PHI_R);
1605        assert_eq!(N_PHI.recip(), N_PHI_R);
1606
1607        // conjugate
1608        assert_eq!(PHI.conj(), N_PHI_R);
1609        assert_eq!(PHI_R.conj(), N_PHI);
1610
1611        // is_pure
1612        assert_eq!(PHI.is_pure(), false);
1613        assert_eq!(SQ5.is_pure(), true);
1614    }
1615
1616    #[test]
1617    fn arithmic_test() {
1618        // add
1619        assert_eq!(PHI_R + 1, PHI);
1620        assert_eq!(PHI45_R + 1, PHI45);
1621        assert_eq!(N_PHI + 1, N_PHI_R);
1622        assert_eq!(PHI_R + Ratio::one(), PHI);
1623        assert_eq!(PHI45_R + Ratio::one(), PHI45);
1624        assert_eq!(N_PHI + Ratio::one(), N_PHI_R);
1625        assert_eq!(PHI + PHI_R, SQ5);
1626        assert_eq!(N_PHI + N_PHI_R, N_SQ5);
1627        assert!((PHI + N_PHI).is_zero());
1628        assert!((PHI + N_PHI_R).is_one());
1629
1630        // sub
1631        assert_eq!(PHI - 1, PHI_R);
1632        assert_eq!(PHI45 - 1, PHI45_R);
1633        assert_eq!(N_PHI_R - 1, N_PHI);
1634        assert_eq!(PHI - Ratio::one(), PHI_R);
1635        assert_eq!(PHI45 - Ratio::one(), PHI45_R);
1636        assert_eq!(N_PHI_R - Ratio::one(), N_PHI);
1637        assert!((PHI - PHI).is_zero());
1638        assert!((PHI - PHI_R).is_one());
1639        assert!((N_PHI_R - N_PHI).is_one());
1640        assert_eq!(PHI - N_PHI_R, SQ5);
1641        assert_eq!(N_PHI - PHI_R, -SQ5);
1642
1643        // mul
1644        assert_eq!(PHI * 1, PHI);
1645        assert_eq!(PHI * -1, N_PHI);
1646        assert_eq!(PHI * Ratio::one(), PHI);
1647        assert_eq!(PHI * -Ratio::one(), N_PHI);
1648        assert!((PHI * PHI_R).is_one());
1649        assert!((N_PHI * N_PHI_R).is_one());
1650        assert_eq!(PHI * PHI, PHI_SQ);
1651        assert_eq!(N_PHI * N_PHI, PHI_SQ);
1652
1653        // div
1654        assert_eq!(PHI / 1, PHI);
1655        assert_eq!(PHI / -1, N_PHI);
1656        assert_eq!(PHI / Ratio::one(), PHI);
1657        assert_eq!(PHI / -Ratio::one(), N_PHI);
1658        assert!((PHI / PHI).is_one());
1659        assert!((N_PHI / N_PHI).is_one());
1660        assert_eq!(PHI / PHI_R, PHI_SQ);
1661        assert_eq!(N_PHI / N_PHI_R, PHI_SQ);
1662
1663        // associativity test
1664        let three_half = Ratio::new(3, 2);
1665        assert_eq!(PHI + 5 - PHI, QuadraticSurd::from(5));
1666        assert_eq!(PHI + three_half - PHI, QuadraticSurd::from(three_half));
1667        assert_eq!(PHI + SQ5 - PHI, SQ5);
1668        assert_eq!(PHI + 5 - PHI45, QuadraticSurd::from(5));
1669        assert_eq!(PHI + three_half - PHI45, QuadraticSurd::from(three_half));
1670        assert_eq!(PHI * 5 / PHI45, QuadraticSurd::from(5));
1671        assert_eq!(PHI * three_half / PHI45, QuadraticSurd::from(three_half));
1672
1673        // mixed
1674        assert_eq!(PHI * 2 - 1, SQ5);
1675        assert_eq!(PHI_R * 2 + 1, SQ5);
1676        assert_eq!(PHI / Ratio::new(1, 2) - 1, SQ5);
1677        assert_eq!((PHI - Ratio::new(1, 2)) * 2, SQ5);
1678    }
1679
1680    #[test]
1681    fn arithmic_test_diff_base() {
1682        assert_eq!(PHI45 + PHI_R, SQ5);
1683        assert_eq!(PHI + PHI45_R, SQ5);
1684        assert!((PHI - PHI45).is_zero());
1685        assert!((PHI - PHI45_R).is_one());
1686        assert!((PHI45 * PHI_R).is_one());
1687        assert!((PHI * PHI45_R).is_one());
1688        assert!((PHI45 / PHI).is_one());
1689        assert!((PHI / PHI45).is_one());
1690        assert_eq!(PHI + SQ5 - PHI45, SQ5);
1691        assert_eq!(PHI * SQ5 / PHI45, SQ5);
1692
1693        let a = QuadraticSurd::new(0, 2, 3, 2);
1694        let b = QuadraticSurd::new(0, 3, 4, 3);
1695        let c = QuadraticSurd::new(0, 1, 2, 6);
1696        assert_eq!(a * b, c);
1697        assert_eq!(c / a, b);
1698        assert_eq!(c / b, a);
1699    }
1700
1701    #[test]
1702    fn conversion_between_contfrac() {
1703        let cf_phi = ContinuedFraction::<u32>::new(vec![1i32], vec![1], false);
1704        assert_eq!(QuadraticSurd::from(cf_phi.clone()), PHI);
1705        assert_eq!(ContinuedFraction::try_from(PHI).unwrap(), cf_phi);
1706
1707        let cf_sq2 = ContinuedFraction::<u32>::new(vec![1i32], vec![2], false);
1708        let surd_sq2 = QuadraticSurd::new(0, 1, 1, 2);
1709        assert_eq!(QuadraticSurd::from(cf_sq2.clone()), surd_sq2);
1710        assert_eq!(ContinuedFraction::try_from(surd_sq2).unwrap(), cf_sq2);
1711
1712        let cf_n_sq2 = ContinuedFraction::<u32>::new(vec![1i32], vec![2], true);
1713        let surd_n_sq2 = QuadraticSurd::new(0, -1, 1, 2);
1714        assert_eq!(QuadraticSurd::from(cf_n_sq2.clone()), surd_n_sq2);
1715        assert_eq!(ContinuedFraction::try_from(surd_n_sq2).unwrap(), cf_n_sq2);
1716
1717        let cf_sq2_7 = ContinuedFraction::<u32>::new(vec![0i32, 4], vec![1, 18, 1, 8], false);
1718        let surd_sq2_7 = QuadraticSurd::new(0, 1, 7, 2);
1719        assert_eq!(QuadraticSurd::from(cf_sq2_7.clone()), surd_sq2_7);
1720        assert_eq!(ContinuedFraction::try_from(surd_sq2_7).unwrap(), cf_sq2_7);
1721
1722        let cf_10_sq2_7 = ContinuedFraction::<u32>::new(vec![1i32, 4], vec![2], false);
1723        let surd_10_sq2_7 = QuadraticSurd::new(10, -1, 7, 2);
1724        assert_eq!(QuadraticSurd::from(cf_10_sq2_7.clone()), surd_10_sq2_7);
1725        assert_eq!(
1726            ContinuedFraction::try_from(surd_10_sq2_7).unwrap(),
1727            cf_10_sq2_7
1728        );
1729    }
1730
1731    #[test]
1732    fn formatting_test() {
1733        assert_eq!(format!("{}", QuadraticSurd::<i8>::zero()), "0");
1734        assert_eq!(format!("{}", QuadraticSurd::<i8>::one()), "1");
1735        assert_eq!(format!("{}", QuadraticSurd::from_sqrt(5).unwrap()), "√5");
1736
1737        assert_eq!(format!("{}", QuadraticSurd::new(1, 1, 1, 5)), "1+√5");
1738        assert_eq!(format!("{}", QuadraticSurd::new(1, -1, 1, 5)), "1-√5");
1739        assert_eq!(format!("{}", QuadraticSurd::new(-1, 1, 1, 5)), "-1+√5");
1740        assert_eq!(format!("{}", QuadraticSurd::new(-1, -1, 1, 5)), "-1-√5");
1741        assert_eq!(format!("{}", QuadraticSurd::new(1, 2, 1, 5)), "1+2√5");
1742        assert_eq!(format!("{}", QuadraticSurd::new(1, -2, 1, 5)), "1-2√5");
1743        assert_eq!(format!("{}", QuadraticSurd::new(-1, 2, 1, 5)), "-1+2√5");
1744        assert_eq!(format!("{}", QuadraticSurd::new(-1, -2, 1, 5)), "-1-2√5");
1745        assert_eq!(format!("{}", QuadraticSurd::new(1, 0, 2, 5)), "1/2");
1746        assert_eq!(format!("{}", QuadraticSurd::new(-1, 0, 2, 5)), "-1/2");
1747        assert_eq!(format!("{}", QuadraticSurd::new(1, 1, 2, 5)), "(1+√5)/2");
1748        assert_eq!(format!("{}", QuadraticSurd::new(1, -1, 2, 5)), "(1-√5)/2");
1749        assert_eq!(format!("{}", QuadraticSurd::new(-1, 1, 2, 5)), "(-1+√5)/2");
1750        assert_eq!(format!("{}", QuadraticSurd::new(-1, -1, 2, 5)), "(-1-√5)/2");
1751        assert_eq!(format!("{}", QuadraticSurd::new(1, 2, 2, 5)), "(1+2√5)/2");
1752        assert_eq!(format!("{}", QuadraticSurd::new(1, -2, 2, 5)), "(1-2√5)/2");
1753        assert_eq!(format!("{}", QuadraticSurd::new(-1, 2, 2, 5)), "(-1+2√5)/2");
1754        assert_eq!(
1755            format!("{}", QuadraticSurd::new(-1, -2, 2, 5)),
1756            "(-1-2√5)/2"
1757        );
1758    }
1759
1760    #[test]
1761    fn trunc_frac_test() {
1762        assert_eq!(PHI.trunc_ref(), QuadraticSurd::from(1));
1763        assert_eq!(PHI.fract_ref(), PHI_R);
1764        assert_eq!(PHI_R.trunc_ref(), QuadraticSurd::zero());
1765        assert_eq!(PHI_R.fract_ref(), PHI_R);
1766        assert_eq!(N_PHI.trunc_ref(), QuadraticSurd::from(-1));
1767        assert_eq!(N_PHI.fract_ref(), N_PHI_R);
1768        assert_eq!(N_PHI_R.trunc_ref(), QuadraticSurd::from(0));
1769        assert_eq!(N_PHI_R.fract_ref(), N_PHI_R);
1770        assert_eq!(PHI_SQ.trunc_ref(), QuadraticSurd::from(2));
1771        assert_eq!(PHI_SQ.fract_ref(), PHI_R);
1772        assert_eq!(PHI45.trunc_ref(), QuadraticSurd::from(1));
1773        assert_eq!(PHI45.fract_ref(), PHI_R);
1774        assert_eq!(PHI45_R.trunc_ref(), QuadraticSurd::from(0));
1775        assert_eq!(PHI45_R.fract_ref(), PHI_R);
1776        assert_eq!(SQ5.trunc_ref(), QuadraticSurd::from(2));
1777        assert_eq!(SQ5.fract_ref(), QuadraticSurd::new(-2, 1, 1, 5));
1778        assert_eq!(N_SQ5.trunc_ref(), QuadraticSurd::from(-2));
1779        assert_eq!(N_SQ5.fract_ref(), QuadraticSurd::new(2, -1, 1, 5));
1780    }
1781
1782    #[test]
1783    fn from_sqrt_test() {
1784        assert_eq!(
1785            QuadraticSurd::from_sqrt(5i32).unwrap(),
1786            QuadraticSurd::new_raw(0, 1, 1, 5)
1787        );
1788        assert_eq!(
1789            QuadraticSurd::from_sqrt(QuadraticSurd::new(3i32, 2, 1, 2)).unwrap(),
1790            QuadraticSurd::new_raw(1, 1, 1, 2)
1791        );
1792
1793        #[cfg(not(feature = "complex"))]
1794        {
1795            let err = QuadraticSurd::from_sqrt(-2i32).unwrap_err();
1796            assert_eq!(
1797                err.kind,
1798                SqrtErrorKind::Complex
1799            );
1800            assert_eq!(err.data, -2);
1801            let err = QuadraticSurd::from_sqrt(Ratio::new(-1i32, 2)).unwrap_err();
1802            assert_eq!(
1803                err.kind,
1804                SqrtErrorKind::Complex
1805            );
1806            assert_eq!(err.data, Ratio::new(-1, 2));
1807            let surd = QuadraticSurd::new(-2i32, 1, 1, 2);
1808            let err = QuadraticSurd::from_sqrt(surd).unwrap_err();
1809            assert_eq!(
1810                err.kind,
1811                SqrtErrorKind::Complex
1812            );
1813            assert_eq!(err.data, surd);
1814        }
1815    }
1816}
1817
1818#[cfg(feature = "complex")]
1819#[cfg(test)]
1820mod complex_tests {
1821    use super::*;
1822
1823    pub const OMEGA: QuadraticSurd<i32> = QuadraticSurd::new_raw(-1, 1, 2, -3); // -0.5 + 0.866i
1824    pub const OMEGA3: QuadraticSurd<i32> = QuadraticSurd::new_raw(-3, 3, 2, -3); // -1.5 + 2.598i
1825
1826    #[test]
1827    fn trunc_frac_test() {
1828        assert_eq!(OMEGA.trunc_ref(), QuadraticSurd::zero());
1829        assert_eq!(OMEGA3.trunc_ref(), QuadraticSurd::new(-1, 2, 1, -1));
1830    }
1831
1832    #[test]
1833    fn from_sqrt_test() {
1834        assert!(QuadraticSurd::from_sqrt(-2i32).is_ok());
1835        assert!(QuadraticSurd::from_sqrt(Ratio::new(-1i32, 2)).is_ok());
1836        assert!(QuadraticSurd::from_sqrt(QuadraticSurd::from(-1i32)).is_ok());
1837
1838        let sq2i = QuadraticSurd::from_sqrt(-2i32).unwrap();
1839        let sqhalfi = QuadraticSurd::from_sqrt(Ratio::new(-1i32, 2)).unwrap();
1840        let i = QuadraticSurd::from_sqrt(QuadraticSurd::from(-1i32)).unwrap();
1841        assert_eq!(sqhalfi * 2, sq2i);
1842        assert_eq!(sq2i * i, QuadraticSurd::from_sqrt(2i32).unwrap());
1843        
1844        assert_eq!(
1845            QuadraticSurd::from_sqrt(QuadraticSurd::new(-3i32, 4, 1, -1)).unwrap(),
1846            QuadraticSurd::new_raw(1, 2, 1, -1)
1847        );
1848    }
1849
1850    #[test]
1851    fn formatting_test() {
1852        // test trivial complex cases
1853        assert_eq!(format!("{:#}", QuadraticSurd::from_sqrt(-1).unwrap()), "i");
1854        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 1, 1, -1)), "1+i");
1855        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -1, 1, -1)), "1-i");
1856        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 1, 1, -1)), "-1+i");
1857        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, -1, 1, -1)), "-1-i");
1858        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 2, 1, -1)), "1+2i");
1859        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -2, 1, -1)), "1-2i");
1860        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 2, 1, -1)), "-1+2i");
1861        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, -2, 1, -1)), "-1-2i");
1862        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 0, 2, -1)), "1/2");
1863        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 0, 2, -1)), "-1/2");
1864        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 1, 2, -1)), "(1+i)/2");
1865        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -1, 2, -1)), "(1-i)/2");
1866        assert_eq!(
1867            format!("{:#}", QuadraticSurd::new(-1, 1, 2, -1)),
1868            "(-1+i)/2"
1869        );
1870        assert_eq!(
1871            format!("{:#}", QuadraticSurd::new(-1, -1, 2, -1)),
1872            "(-1-i)/2"
1873        );
1874        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 2, 2, -1)), "(1+2i)/2");
1875        assert_eq!(
1876            format!("{:#}", QuadraticSurd::new(1, -2, 2, -1)),
1877            "(1-2i)/2"
1878        );
1879        assert_eq!(
1880            format!("{:#}", QuadraticSurd::new(-1, 2, 2, -1)),
1881            "(-1+2i)/2"
1882        );
1883
1884        // test non-trivial complex cases
1885        assert_eq!(format!("{:#}", QuadraticSurd::from_sqrt(-3).unwrap()), "√3i");
1886        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 1, 1, -3)), "1+√3i");
1887        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -1, 1, -3)), "1-√3i");
1888        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 1, 1, -3)), "-1+√3i");
1889        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, -1, 1, -3)), "-1-√3i");
1890        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 2, 1, -3)), "1+2√3i");
1891        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -2, 1, -3)), "1-2√3i");
1892        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 2, 1, -3)), "-1+2√3i");
1893        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, -2, 1, -3)), "-1-2√3i");
1894        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 0, 2, -3)), "1/2");
1895        assert_eq!(format!("{:#}", QuadraticSurd::new(-1, 0, 2, -3)), "-1/2");
1896        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 1, 2, -3)), "(1+√3i)/2");
1897        assert_eq!(format!("{:#}", QuadraticSurd::new(1, -1, 2, -3)), "(1-√3i)/2");
1898        assert_eq!(
1899            format!("{:#}", QuadraticSurd::new(-1, 1, 2, -3)),
1900            "(-1+√3i)/2"
1901        );
1902        assert_eq!(
1903            format!("{:#}", QuadraticSurd::new(-1, -1, 2, -3)),
1904            "(-1-√3i)/2"
1905        );
1906        assert_eq!(format!("{:#}", QuadraticSurd::new(1, 2, 2, -3)), "(1+2√3i)/2");
1907        assert_eq!(
1908            format!("{:#}", QuadraticSurd::new(1, -2, 2, -3)),
1909            "(1-2√3i)/2"
1910        );
1911        assert_eq!(
1912            format!("{:#}", QuadraticSurd::new(-1, 2, 2, -3)),
1913            "(-1+2√3i)/2"
1914        );
1915    }
1916}
num_irrational/quadratic/surd.rs

num_irrational/quadratic/
surd.rs
