num_irrational/quadratic/
integer.rs

1//! Implementation of quadratic integers
2
3use super::{QuadraticBase, QuadraticOps};
4use core::ops::*;
5use num_integer::Integer;
6use num_traits::{NumRef, One, RefNum, Signed, Zero};
7
8#[inline]
9fn four<T: Add<Output = T> + One>() -> T {
10    T::one() + T::one() + T::one() + T::one()
11}
12
13/// return -1 if v ≡ 0 mod 4, 0 if v ≡ 1 mod 4, 1 if v ≡ 2 or 3 mod 4
14#[inline]
15fn mod4d2<T: Signed>(v: &T) -> i8
16where
17    for<'r> &'r T: RefNum<T>,
18{
19    let m: T = v % four::<T>();
20    if m.is_zero() {
21        -1
22    } else if v.is_negative() {
23        let m2: T = (four::<T>() + m) / (T::one() + T::one());
24        if m2.is_one() {
25            1
26        } else {
27            0
28        }
29    } else {
30        let m2: T = m / (T::one() + T::one());
31        if m2.is_one() {
32            1
33        } else {
34            0
35        }
36    }
37}
38
39// x/y rounded to the nearest integer. If frac(x/y) = 1/2, then round away from zero.
40#[inline]
41fn div_rem_round<T: Integer + Signed + NumRef>(x: T, y: &T) -> (T, T)
42where
43    for<'r> &'r T: RefNum<T>,
44{
45    let two = T::one() + T::one();
46    let offset = y.abs() / two * x.signum();
47    let (q, r) = (x + &offset).div_rem(y);
48    (q, r - offset)
49}
50
51/// Underlying representation of a quadratic integer `a + bω` as (a,b), where `ω` is `√D` or `(1+√D)/2`.
52///
53/// Specifically, `ω=√D` when `D ≡ 2,3 mod 4`, and `ω=(1+√D)/2` when `D ≡ 1 mod 4`. Note that when `ω=(1+√D)/2`,
54/// `ω² = (D-1)/4 + w`.
55/// 
56/// The division of two quadratic ints is defined as: calculate a/b = x + yω in the quadratic field, then round x
57/// and y to nearest integer (round half-way cases away from zero, to be consistent with `round()` functions of
58/// primitive types).
59#[derive(Debug, Clone, Copy, PartialEq)]
60pub struct QuadraticIntCoeffs<T>(pub T, pub T);
61
62impl<T> From<(T, T)> for QuadraticIntCoeffs<T> {
63    fn from(v: (T, T)) -> Self {
64        Self(v.0, v.1)
65    }
66}
67
68impl<T: Add<Output = T>> Add<T> for QuadraticIntCoeffs<T> {
69    type Output = Self;
70    #[inline]
71    fn add(self, rhs: T) -> Self::Output {
72        Self(self.0 + rhs, self.1)
73    }
74}
75impl<T: Add<Output = T>> Add for QuadraticIntCoeffs<T> {
76    type Output = Self;
77    #[inline]
78    fn add(self, rhs: Self) -> Self::Output {
79        Self(self.0 + rhs.0, self.1 + rhs.1)
80    }
81}
82
83impl<T: Sub<Output = T>> Sub<T> for QuadraticIntCoeffs<T> {
84    type Output = Self;
85    #[inline]
86    fn sub(self, rhs: T) -> Self::Output {
87        Self(self.0 - rhs, self.1)
88    }
89}
90impl<T: Sub<Output = T>> Sub for QuadraticIntCoeffs<T> {
91    type Output = Self;
92    #[inline]
93    fn sub(self, rhs: Self) -> Self::Output {
94        Self(self.0 - rhs.0, self.1 - rhs.1)
95    }
96}
97
98impl<T: Signed + NumRef> Mul<T> for QuadraticIntCoeffs<T> {
99    type Output = Self;
100    #[inline]
101    fn mul(self, rhs: T) -> Self::Output {
102        Self(self.0 * &rhs, self.1 * rhs)
103    }
104}
105
106impl<T: Integer + Signed + NumRef> Div<T> for QuadraticIntCoeffs<T>
107where
108    for<'r> &'r T: RefNum<T>,
109{
110    type Output = Self;
111    #[inline]
112    fn div(self, rhs: T) -> Self::Output {
113        Self(div_rem_round(self.0, &rhs).0, div_rem_round(self.1, &rhs).0)
114    }
115}
116impl<T: Integer + Signed + NumRef + Clone> Rem<T> for QuadraticIntCoeffs<T>
117where
118    for<'r> &'r T: RefNum<T>,
119{
120    type Output = Self;
121    #[inline]
122    fn rem(self, rhs: T) -> Self::Output {
123        Self(div_rem_round(self.0, &rhs).1, div_rem_round(self.1, &rhs).1)
124    }
125}
126
127impl<T: Integer + Signed + NumRef> QuadraticIntCoeffs<T>
128where
129    for<'r> &'r T: RefNum<T>,
130{
131    /// Get the norm of the quadratic integer (taking reference)
132    fn norm_ref(&self, discr: &T) -> T {
133        if discr.is_zero() {
134            return &self.0 * &self.0;
135        }
136
137        match mod4d2(discr) {
138            0 => {
139                // |a+bw| = (a+b/2)^2 - (b/2*sq)^2 = a^2+ab+b^2/4 - D*b^2/4 = a^2 + ab + b^2(1-D)/4
140                &self.0 * &self.0
141                    + &self.0 * &self.1
142                    + &self.1 * &self.1 * ((T::one() - discr) / four::<T>())
143            }
144            1 => &self.0 * &self.0 - &self.1 * &self.1 * discr,
145            _ => unreachable!(),
146        }
147    }
148
149    /// Get the quotient and the remainder of self / rhs
150    fn div_rem(self, rhs: Self, discr: &T) -> (Self, Self) {
151        // This division algorithm is compatible with Gaussian and Eisenstein integers,
152        // but it's not necessarily a proper Euclidean division algorithm for any
153        // quadratic integers.
154
155        let (a, b) = match mod4d2(discr) {
156            0 => (
157                // (a+bw)/(c+dw) = (a+bw)*conj(c+dw)/|c+dw| = (a+bw)(c+d-dw) / |c+dw|
158                // = (ac+ad-adw+bcw+bdw-bdw^2)/|c+dw|
159                // = (ac+ad-bd*(D-1)/4)/|c+dw| + (-ad+bc)w/|c+dw|
160                &self.0 * &rhs.0 + &self.0 * &rhs.1
161                    - &self.1 * &rhs.1 * ((discr - T::one()) / four::<T>()),
162                &self.1 * &rhs.0 - &self.0 * &rhs.1,
163            ),
164            1 => (
165                // (a+bw)/(c+dw) = (a+bw)*conj(c+dw)/|c+dw| = (a+bw)(c-dw) / |c+dw|
166                // = (ac-bd*D)/|c+dw| + (bc-ad)w/|c+dw|
167                &self.0 * &rhs.0 - &self.1 * &rhs.1 * discr,
168                &self.1 * &rhs.0 - &self.0 * &rhs.1,
169            ),
170            _ => unreachable!(),
171        };
172        let n = rhs.norm_ref(discr);
173        let (aq, ar) = div_rem_round(a, &n);
174        let (bq, br) = div_rem_round(b, &n);
175        let (q, r) = (Self(aq, bq), Self(ar, br));
176        (q, QuadraticOps::mul(rhs, r, discr) / n)
177    }
178}
179
180impl<T: Integer + Signed + NumRef> QuadraticOps<Self, &T, Self> for QuadraticIntCoeffs<T>
181where
182    for<'r> &'r T: RefNum<T>,
183{
184    type Scalar = T;
185
186    #[inline]
187    fn mul(self, rhs: Self, discr: &T) -> Self {
188        match mod4d2(discr) {
189            0 => Self(
190                // (a+bw)*(c+dw) = ac + bdw^2 + (bc+ad)w = ac+bd(D-1)/4 + (bc+ad+bd)w
191                &self.0 * &rhs.0 + &self.1 * &rhs.1 * ((discr - T::one()) / four::<T>()),
192                self.0 * &rhs.1 + self.1 * (rhs.0 + rhs.1),
193            ),
194            1 => Self(
195                &self.0 * &rhs.0 + &self.1 * &rhs.1 * discr,
196                self.0 * rhs.1 + self.1 * rhs.0,
197            ),
198            _ => unreachable!(),
199        }
200    }
201    #[inline]
202    fn div(self, rhs: Self, discr: &T) -> Self {
203        self.div_rem(rhs, discr).0
204    }
205
206    #[inline]
207    fn conj(self, discr: &T) -> Self {
208        if discr.is_zero() {
209            return self;
210        }
211
212        match mod4d2(discr) {
213            0 => Self(
214                // conj(a+bw) = conj(a+b/2 + b/2*sq) = a+b/2 - b/2*sq = a+b - bw
215                self.0 + &self.1,
216                -self.1,
217            ),
218            1 => Self(self.0, -self.1),
219            _ => unreachable!(),
220        }
221    }
222    #[inline]
223    fn norm(self, discr: &T) -> T {
224        self.norm_ref(discr)
225    }
226}
227
228/// Quadratic integer `a + bω`, where `ω` is `√D` or `(1+√D)/2`, based on [QuadraticIntCoeffs]
229///
230/// The different between [QuadraticInt] and [QuadraticSurd][crate::QuadraticSurd] is that the operations for the
231/// latter will be in normal fields of real numbers or complex numbers, while the operations
232/// for the former will be in the quadratic field (specifically in the quadratic integer ring ℤ\[ω\]).
233/// Therefore, the arithmetic operations between [QuadraticInt]s can only be performed with equivalent bases (e.g. `√18=3√2` is compatible with `√2`).
234/// However for [QuadraticSurd][crate::QuadraticSurd], you can mix the number with different bases if they are both pure quadratic numbers,
235/// such as `√2 * √3 = √6`.
236#[derive(Debug, Clone, Copy)]
237pub struct QuadraticInt<T> {
238    coeffs: QuadraticIntCoeffs<T>, // (a, b), b = 0 is ensured when D = 0
239    discr: T,                      // D, D mod 4 != 0 is ensured when D != 0
240}
241
242impl<T: QuadraticBase> QuadraticInt<T>
243where
244    for<'r> &'r T: RefNum<T>,
245{
246    /// Create a quadratic integer from `a + b√r`, where a, b, r are all integers
247    #[inline]
248    pub fn new(a: T, b: T, r: T) -> Self {
249        let (mut b, mut r) = (b, r);
250
251        // remove factor 4
252        if r.is_positive() {
253            while (&r % four::<T>()).is_zero() {
254                r = r / four::<T>();
255                b = &b + &b;
256            }
257        }
258
259        match mod4d2(&r) {
260            0 => {
261                let (y, x) = (&b + &b, a - b);
262                Self::from_coeffs(QuadraticIntCoeffs(x, y), r)
263            },
264            1 => Self::from_coeffs(QuadraticIntCoeffs(a, b), r),
265            _ => unreachable!()
266        }
267    }
268    /// Create a quadratic integer `a + bω`, where `ω` is `√D` or `(1+√D)/2`.
269    /// 
270    /// # Panics
271    /// - If the discriminant is negative and feature "complex" is not enabled
272    /// - The discriminant is congruent to 0 modulo 4
273    pub fn from_coeffs(coeffs: QuadraticIntCoeffs<T>, discr: T) -> Self {
274        #[cfg(not(feature = "complex"))]
275        if discr.is_negative() {
276            panic!("Negative root is not supported without the `complex` feature");
277        }
278        if (&discr % four::<T>()).is_zero() {
279            panic!("The discriminant must not be 0 modulo 4!")
280        }
281
282        let QuadraticIntCoeffs(a, b) = coeffs;
283
284        // test if the root is a square number
285        let root = discr.abs().sqrt();
286        let (a, b, r) = if &root * &root == discr {
287            if discr.is_negative() {
288                (a, b * root, -T::one())
289            } else {
290                (a + b * root, T::zero(), T::zero())
291            }
292        } else {
293            (a, b, discr)
294        };
295
296        // shortcut for rational number
297        if b.is_zero() || r.is_zero() {
298            return Self {
299                coeffs: QuadraticIntCoeffs(a, T::zero()),
300                discr: T::zero(),
301            };
302        }
303
304        Self {
305            coeffs: QuadraticIntCoeffs(a, b),
306            discr: r,
307        }
308    }
309
310    #[inline]
311    pub fn conj(self) -> Self {
312        Self {
313            coeffs: self.coeffs.conj(&self.discr),
314            discr: self.discr,
315        }
316    }
317    #[inline]
318    pub fn norm(self) -> T {
319        self.coeffs.norm(&self.discr)
320    }
321    #[inline]
322    pub fn is_rational(&self) -> bool {
323        self.coeffs.1.is_zero()
324    }
325    #[inline]
326    pub fn is_pure(&self) -> bool {
327        match mod4d2(&self.discr) {
328            0 => (&self.coeffs.0 + &self.coeffs.0 + &self.coeffs.1).is_zero(),
329            1 => self.coeffs.1.is_zero(),
330            _ => unreachable!(),
331        }
332    }
333
334    /// Get the fundamental unit of the quadratic field ℚ[√d]
335    fn unit(d: T) -> Self {
336        // REF: http://www.numbertheory.org/gnubc/unit
337        // REF: https://people.reed.edu/~jerry/361/lectures/rqunits.pdf
338        unimplemented!()
339    }
340    
341    /// Convert `a + bω` to `(x + y√D)/2`
342    #[inline]
343    fn parts_doubled(&self) -> (T, T) {
344        let QuadraticIntCoeffs(a, b) = &self.coeffs;
345        match mod4d2(&self.discr) {
346            0 => (a + a + b, b.clone()),
347            1 => (a + a, b + b),
348            _ => unreachable!()
349        }
350    }
351}
352
353impl<T: QuadraticBase> QuadraticInt<T>
354where
355    for<'r> &'r T: RefNum<T>,
356{
357    // Try to reduce the root base with a possible factor.
358    fn reduce_root_hinted(self, hint: T) -> Self {
359        let hint = hint.abs();
360        if hint.is_zero() || hint.is_one() || hint.is_negative() {
361            return self;
362        }
363
364        let (quo, rem) = self.discr.div_rem(&hint);
365        if rem.is_zero() {
366            // if hint is actually a factor ...
367            let root = hint.clone().sqrt();
368            if &root * &root == hint {
369                // and if hint is a square number, then remove the hint factor.
370                // note that r ≡ r/a^2 mod 4 (given a is odd), so reducing is feasible
371                return match mod4d2(&self.discr) {
372                    0 => {
373                        // ω_{a*a*r} = (1 + a√r)/2 = a(1 + √r)/2 - (a-1)/2 = a*ω_{r} - (a-1)/2
374                        let scale = (&root - T::one()) / (T::one() + T::one());
375                        Self::from_coeffs(QuadraticIntCoeffs(self.coeffs.0 - &self.coeffs.1 * scale, self.coeffs.1 * root), quo)
376                    },
377                    1 => Self::from_coeffs(QuadraticIntCoeffs(self.coeffs.0, self.coeffs.1 * root), quo),
378                    _ => unreachable!()
379                }
380                
381            }
382        }
383
384        self
385    }
386}
387
388impl<T: Neg<Output = T>> Neg for QuadraticInt<T> {
389    type Output = Self;
390    fn neg(self) -> Self::Output {
391        Self {
392            coeffs: QuadraticIntCoeffs(-self.coeffs.0, -self.coeffs.1),
393            discr: self.discr,
394        }
395    }
396}
397
398// Reduce root base for binary operands. Return None if the two bases
399// cannot be matched. This function assumes the root bases on both sides
400// are not zero
401#[inline]
402#[cfg(not(feature = "complex"))]
403fn reduce_bin_op<T: QuadraticBase>(
404    lhs: QuadraticInt<T>,
405    rhs: QuadraticInt<T>,
406) -> Option<(QuadraticInt<T>, QuadraticInt<T>)>
407where
408    for<'r> &'r T: RefNum<T>,
409{
410    debug_assert!(!lhs.discr.is_zero());
411    debug_assert!(!rhs.discr.is_zero());
412
413    let result = if lhs.discr > rhs.discr {
414        let hint = &lhs.discr / &rhs.discr;
415        (lhs.reduce_root_hinted(hint), rhs)
416    } else if rhs.discr > lhs.discr {
417        let hint = &rhs.discr / &lhs.discr;
418        (lhs, rhs.reduce_root_hinted(hint))
419    } else {
420        (lhs, rhs)
421    };
422
423    if result.0.discr == result.1.discr {
424        Some(result)
425    } else {
426        None
427    }
428}
429
430#[inline]
431#[cfg(feature = "complex")]
432fn reduce_bin_op<T: QuadraticBase>(
433    lhs: QuadraticInt<T>,
434    rhs: QuadraticInt<T>,
435) -> Option<(QuadraticInt<T>, QuadraticInt<T>)>
436where
437    for<'r> &'r T: RefNum<T>,
438{
439    if lhs.discr.is_negative() ^ rhs.discr.is_negative() {
440        return None;
441    }
442
443    let result = if lhs.discr.abs() > rhs.discr.abs() {
444        let hint = &lhs.discr / &rhs.discr;
445        (lhs.reduce_root_hinted(hint), rhs)
446    } else if rhs.discr.abs() > lhs.discr.abs() {
447        let hint = &rhs.discr / &lhs.discr;
448        (lhs, rhs.reduce_root_hinted(hint))
449    } else {
450        (lhs, rhs)
451    };
452
453    if result.0.discr == result.1.discr {
454        Some(result)
455    } else {
456        None
457    }
458}
459
460#[inline]
461fn reduce_bin_op_unwrap<T: QuadraticBase>(
462    lhs: QuadraticInt<T>,
463    rhs: QuadraticInt<T>,
464) -> (QuadraticInt<T>, QuadraticInt<T>)
465where
466    for<'r> &'r T: RefNum<T>,
467{
468    reduce_bin_op(lhs, rhs).expect("two root bases are not compatible!")
469}
470
471impl<T: QuadraticBase> PartialEq for QuadraticInt<T>
472where
473    for<'r> &'r T: RefNum<T>,
474{
475    fn eq(&self, other: &Self) -> bool {
476        // shortcuts
477        if self.discr == other.discr {
478            return self.coeffs == other.coeffs;
479        }
480        if self.discr.is_zero() || other.discr.is_zero() {
481            return false;
482        }
483
484        let (la, lb) = self.parts_doubled();
485        let (ra, rb) = other.parts_doubled();
486
487        la == ra && &lb * &lb * &self.discr == &rb * &rb * &other.discr
488    }
489}
490
491impl<T: QuadraticBase> Eq for QuadraticInt<T> where for<'r> &'r T: RefNum<T> {}
492
493impl<T: Add<Output = T>> Add<T> for QuadraticInt<T> {
494    type Output = Self;
495    fn add(self, rhs: T) -> Self::Output {
496        Self {
497            coeffs: self.coeffs + rhs,
498            discr: self.discr,
499        }
500    }
501}
502impl<T: QuadraticBase> Add for QuadraticInt<T>
503where
504    for<'r> &'r T: RefNum<T>,
505{
506    type Output = Self;
507    #[inline]
508    fn add(self, rhs: Self) -> Self::Output {
509        if rhs.is_rational() {
510            return self + rhs.coeffs.0;
511        }
512        if self.is_rational() {
513            return rhs + self.coeffs.0;
514        }
515
516        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
517        Self::from_coeffs(lhs.coeffs + rhs.coeffs, rhs.discr)
518    }
519}
520
521impl<T: Sub<Output = T>> Sub<T> for QuadraticInt<T> {
522    type Output = Self;
523    fn sub(self, rhs: T) -> Self::Output {
524        Self {
525            coeffs: self.coeffs - rhs,
526            discr: self.discr,
527        }
528    }
529}
530impl<T: QuadraticBase> Sub for QuadraticInt<T>
531where
532    for<'r> &'r T: RefNum<T>,
533{
534    type Output = Self;
535    #[inline]
536    fn sub(self, rhs: Self) -> Self::Output {
537        if rhs.is_rational() {
538            return self - rhs.coeffs.0;
539        }
540        if self.is_rational() {
541            return -(rhs - self.coeffs.0);
542        }
543
544        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
545        Self::from_coeffs(lhs.coeffs - rhs.coeffs, rhs.discr)
546    }
547}
548
549impl<T: QuadraticBase> Mul<T> for QuadraticInt<T>
550where
551    for<'r> &'r T: RefNum<T>,
552{
553    type Output = Self;
554    #[inline]
555    fn mul(self, rhs: T) -> Self::Output {
556        if rhs.is_zero() {
557            Self::zero()
558        } else {
559            Self {
560                coeffs: self.coeffs * rhs,
561                discr: self.discr,
562            }
563        }
564    }
565}
566impl<T: QuadraticBase> Mul for QuadraticInt<T>
567where
568    for<'r> &'r T: RefNum<T>,
569{
570    type Output = Self;
571    #[inline]
572    fn mul(self, rhs: Self) -> Self::Output {
573        if rhs.is_rational() {
574            return self * rhs.coeffs.0;
575        }
576        if self.is_rational() {
577            return rhs * self.coeffs.0;
578        }
579        // XXX: currently, multiplication (and division) of pure quadratic numbers is not allowed here, because
580        // it might lead to non-representable results. For example, `ω(2) * (ω(5)-1/2) = √2 * √5/2 = √10/2`,
581        // the result is not a quadratic integer
582
583        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
584        Self::from_coeffs(
585            QuadraticOps::mul(lhs.coeffs, rhs.coeffs, &lhs.discr),
586            rhs.discr,
587        )
588    }
589}
590
591impl<T: QuadraticBase> Div<T> for QuadraticInt<T>
592where
593    for<'r> &'r T: RefNum<T>,
594{
595    type Output = Self;
596    #[inline]
597    fn div(self, rhs: T) -> Self::Output {
598        Self {
599            coeffs: self.coeffs / rhs,
600            discr: self.discr,
601        }
602    }
603}
604impl<T: QuadraticBase> Div for QuadraticInt<T>
605where
606    for<'r> &'r T: RefNum<T>,
607{
608    type Output = Self;
609    #[inline]
610    fn div(self, rhs: Self) -> Self::Output {
611        if rhs.is_rational() {
612            return self / rhs.coeffs.0;
613        }
614
615        if self.is_rational() {
616            Self::from_coeffs(
617                QuadraticOps::div(self.coeffs, rhs.coeffs, &rhs.discr),
618                rhs.discr,
619            )
620        } else {
621            let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
622            Self::from_coeffs(
623                QuadraticOps::div(lhs.coeffs, rhs.coeffs, &lhs.discr),
624                rhs.discr,
625            )
626        }
627    }
628}
629
630impl<T: QuadraticBase> Rem<T> for QuadraticInt<T>
631where
632    for<'r> &'r T: RefNum<T>,
633{
634    type Output = Self;
635    #[inline]
636    fn rem(self, rhs: T) -> Self::Output {
637        Self {
638            coeffs: self.coeffs % rhs,
639            discr: self.discr,
640        }
641    }
642}
643impl<T: QuadraticBase> Rem for QuadraticInt<T>
644where
645    for<'r> &'r T: RefNum<T>,
646{
647    type Output = Self;
648    #[inline]
649    fn rem(self, rhs: Self) -> Self::Output {
650        if rhs.is_rational() {
651            return self / rhs.coeffs.0;
652        }
653
654        let (lhs, rhs) = reduce_bin_op_unwrap(self, rhs);
655        Self::from_coeffs(lhs.coeffs.div_rem(rhs.coeffs, &lhs.discr).1, rhs.discr)
656    }
657}
658
659impl<T: QuadraticBase> Zero for QuadraticInt<T>
660where
661    for<'r> &'r T: RefNum<T>,
662{
663    fn zero() -> Self {
664        Self {
665            coeffs: QuadraticIntCoeffs(T::zero(), T::zero()),
666            discr: T::zero(),
667        }
668    }
669
670    fn is_zero(&self) -> bool {
671        self.coeffs.0.is_zero() && self.coeffs.1.is_zero()
672    }
673}
674
675impl<T: QuadraticBase> One for QuadraticInt<T>
676where
677    for<'r> &'r T: RefNum<T>,
678{
679    fn one() -> Self {
680        Self {
681            coeffs: QuadraticIntCoeffs(T::one(), T::zero()),
682            discr: T::zero(),
683        }
684    }
685}
686
687#[cfg(feature = "complex")]
688mod complex {
689    use super::*;
690
691    #[derive(Debug, Clone, Copy, PartialEq)]
692    pub struct GaussianInt<T>(QuadraticIntCoeffs<T>);
693    impl<T: Eq> Eq for GaussianInt<T> {}
694
695    impl<T> GaussianInt<T> {
696        pub const fn new(re: T, im: T) -> Self {
697            Self(QuadraticIntCoeffs(re, im))
698        }
699    }
700
701    impl<T: Signed> GaussianInt<T> {
702        pub fn conj(self) -> Self {
703            let QuadraticIntCoeffs(a, b) = self.0;
704            Self(QuadraticIntCoeffs(a, -b))
705        }
706    }
707
708    impl<T: Add<Output = T>> GaussianInt<T>
709    where
710        for<'r> &'r T: RefNum<T>,
711    {
712        pub fn norm_ref(&self) -> T {
713            let QuadraticIntCoeffs(a, b) = &self.0;
714            a * a + b * b
715        }
716        pub fn norm(self) -> T {
717            self.norm_ref()
718        }
719    }
720
721    impl<T: Add<Output = T>> Add for GaussianInt<T>
722    where
723        for<'r> &'r T: RefNum<T>,
724    {
725        type Output = Self;
726        #[inline]
727        fn add(self, rhs: Self) -> Self::Output {
728            Self(self.0 + rhs.0)
729        }
730    }
731
732    impl<T: Sub<Output = T>> Sub for GaussianInt<T>
733    where
734        for<'r> &'r T: RefNum<T>,
735    {
736        type Output = Self;
737        #[inline]
738        fn sub(self, rhs: Self) -> Self::Output {
739            Self(self.0 - rhs.0)
740        }
741    }
742
743    impl<T: Integer + Signed + NumRef> Mul for GaussianInt<T>
744    where
745        for<'r> &'r T: RefNum<T>,
746    {
747        type Output = Self;
748        #[inline]
749        fn mul(self, rhs: Self) -> Self::Output {
750            Self(QuadraticOps::mul(self.0, rhs.0, &-T::one()))
751        }
752    }
753
754    impl<T: Integer + Signed + NumRef> Div for GaussianInt<T>
755    where
756        for<'r> &'r T: RefNum<T>,
757    {
758        type Output = Self;
759        #[inline]
760        fn div(self, rhs: Self) -> Self::Output {
761            Self(QuadraticOps::div(self.0, rhs.0, &-T::one()))
762        }
763    }
764
765    impl<T: Integer + Signed + NumRef> Rem for GaussianInt<T>
766    where
767        for<'r> &'r T: RefNum<T>,
768    {
769        type Output = Self;
770        #[inline]
771        fn rem(self, rhs: Self) -> Self::Output {
772            Self(self.0.div_rem(rhs.0, &-T::one()).1)
773        }
774    }
775
776    // TODO (v0.3.1): implement num_integer::Integer for GaussianInt, especially is_even and is_odd can be derived from
777    // https://crates.io/crates/gaussiant
778}
779
780#[cfg(feature = "complex")]
781pub use complex::GaussianInt;
782
783#[cfg(test)]
784mod tests {
785    use super::*;
786
787    #[test]
788    fn div_rem_round_test() {
789        const CASES: [(i8, i8, i8, i8); 18] = [
790            // x, y, (x/y), (x%y)
791            (-4, 4, -1, 0),
792            (-3, 4, -1, 1),
793            (-2, 4, -1, 2),
794            (-1, 4, 0, -1),
795            (0, 4, 0, 0),
796            (1, 4, 0, 1),
797            (2, 4, 1, -2),
798            (3, 4, 1, -1),
799            (4, 4, 1, 0),
800            (-4, -4, 1, 0),
801            (-3, -4, 1, 1),
802            (-2, -4, 1, 2),
803            (-1, -4, 0, -1),
804            (0, -4, 0, 0),
805            (1, -4, 0, 1),
806            (2, -4, -1, -2),
807            (3, -4, -1, -1),
808            (4, -4, -1, 0),
809        ];
810        for (x, y, q, r) in CASES {
811            assert_eq!(
812                div_rem_round(x, &y),
813                (q, r),
814                "{} / {} != ({}, {})",
815                x,
816                y,
817                q,
818                r
819            );
820        }
821    }
822
823    #[test]
824    fn test_ops_unit() {
825        let zero = QuadraticInt::new(0, 0, 2);
826        let e1 = QuadraticInt::new(1, 0, 2);
827        let e2 = QuadraticInt::new(0, 1, 2);
828        let e3 = QuadraticInt::new(-1, 0, 2);
829        let e4 = QuadraticInt::new(0, -1, 2);
830        let m12 = QuadraticInt::new(1, 1, 2);
831        let m34 = QuadraticInt::new(-1, -1, 2);
832
833        assert_eq!(e1 + zero, e1);
834        assert_eq!(e1 - zero, e1);
835        assert_eq!(e1 * zero, zero);
836        assert_eq!(zero / e1, zero);
837
838        assert_eq!(e1 + e2, m12);
839        assert_eq!(e3 + e4, m34);
840        assert_eq!(m12 - e1, e2);
841        assert_eq!(m34 - e3, e4);
842        assert_eq!(e1 * e2, e2);
843        assert_eq!(e2 / e1, e2);
844        assert_eq!(e2 * e3, e4);
845        assert_eq!(e4 / e2, e3);
846        assert_eq!(e3 * e4, e2);
847        assert_eq!(e2 / e3, e4);
848        assert_eq!(e4 * e1, e4);
849        assert_eq!(e4 / e1, e4);
850
851        #[cfg(feature = "complex")]
852        {
853            // unit complex number tests
854            let zero = QuadraticInt::new(0, 0, -1);
855            let e1 = QuadraticInt::new(1, 0, -1);
856            let e2 = QuadraticInt::new(0, 1, -1);
857            let e3 = QuadraticInt::new(-1, 0, -1);
858            let e4 = QuadraticInt::new(0, -1, -1);
859            let m12 = QuadraticInt::new(1, 1, -1);
860            let m34 = QuadraticInt::new(-1, -1, -1);
861
862            assert_eq!(e1 + zero, e1);
863            assert_eq!(e1 - zero, e1);
864            assert_eq!(e1 * zero, zero);
865            assert_eq!(zero / e1, zero);
866
867            assert_eq!(e1 + e2, m12);
868            assert_eq!(e3 + e4, m34);
869            assert_eq!(m12 - e1, e2);
870            assert_eq!(m34 - e3, e4);
871            assert_eq!(e1 * e2, e2);
872            assert_eq!(e2 / e1, e2);
873            assert_eq!(e2 * e3, e4);
874            assert_eq!(e4 / e2, e3);
875            assert_eq!(e3 * e4, e2);
876            assert_eq!(e2 / e3, e4);
877            assert_eq!(e4 * e1, e4);
878            assert_eq!(e4 / e1, e4);
879
880            assert_eq!(e1.norm(), 1);
881            assert_eq!(e2.norm(), 1);
882            assert_eq!(e3.norm(), 1);
883            assert_eq!(e4.norm(), 1);
884            assert_eq!(e1.conj(), e1);
885            assert_eq!(e2.conj(), e4);
886            assert_eq!(e3.conj(), e3);
887            assert_eq!(e4.conj(), e2);
888        }
889    }
890
891    #[test]
892    fn test_ops() {
893        let a = QuadraticInt::new(1, 3, 2);
894        let b = QuadraticInt::new(-2, -1, 2);
895        let c = QuadraticInt::new(-1, 2, 2);
896
897        assert_eq!(a + b, QuadraticInt::new(-1, 2, 2));
898        assert_eq!(a - b, QuadraticInt::new(3, 4, 2));
899        assert_eq!(a * b, QuadraticInt::new(-8, -7, 2));
900        assert_eq!(a / b, QuadraticInt::new(2, -3, 2));
901        assert_eq!(a % b, QuadraticInt::new(-1, -1, 2));
902
903        assert_eq!(b + c, QuadraticInt::new(-3, 1, 2));
904        assert_eq!(b - c, QuadraticInt::new(-1, -3, 2));
905        assert_eq!(b * c, QuadraticInt::new(-2, -3, 2));
906        assert_eq!(b / c, QuadraticInt::new(-1, -1, 2));
907        assert_eq!(b % c, QuadraticInt::new(1, 0, 2));
908    }
909
910    #[test]
911    fn test_different_bases() {
912        let a = QuadraticInt::new(1, 6, 5);
913        let b = QuadraticInt::new(1, 3, 20);
914        let c = QuadraticInt::new(1, 2, 45);
915
916        assert_eq!(a, b);
917        assert_eq!(a, c);
918        assert_eq!(a + b, a + a);
919        assert_eq!(a + c, a + a);
920        assert_eq!(a - b, QuadraticInt::zero());
921        assert_eq!(a - c, QuadraticInt::zero());
922        assert_eq!(a * b, a * a);
923        assert_eq!(a * c, a * a);
924        assert_eq!(a / b, QuadraticInt::one());
925        assert_eq!(a / c, QuadraticInt::one());
926        assert_eq!(a % b, QuadraticInt::zero());
927        assert_eq!(a % c, QuadraticInt::zero());
928    }
929
930    #[test]
931    #[cfg(feature = "complex")]
932    fn test_gaussian() {
933        // test Gaussian integers
934        let q12 = GaussianInt::new(1, 2);
935        let qm12 = GaussianInt::new(-1, 2);
936        let q1m2 = GaussianInt::new(1, -2);
937        let q23 = GaussianInt::new(2, 3);
938        let qm23 = GaussianInt::new(-2, 3);
939        let q2m3 = GaussianInt::new(2, -3);
940
941        assert_eq!(q12.conj(), q1m2);
942        assert_eq!(q23.conj(), q2m3);
943        assert_eq!(q12.norm(), 5);
944        assert_eq!(q23.norm(), 13);
945
946        for &v1 in &[q12, qm12, q1m2, q23, qm23, q2m3] {
947            for &v2 in &[q12, qm12, q1m2, q23, qm23, q2m3] {
948                assert_eq!(v1 + v2 - v1, v2, "add/sub test failed between {:?} and {:?}", v1, v2);
949                assert_eq!(v1 * v2 / v1, v2, "mul/div test failed between {:?} and {:?}", v1, v2);
950            }
951        }
952    }
953
954    #[test]
955    #[cfg(feature = "complex")]
956    fn test_eisenstein() {
957        // test Eisenstein integers
958        let q12 =  QuadraticInt::from_coeffs((1, 2).into(), -3); // 2 + √-3
959        let qm12 = QuadraticInt::from_coeffs((-1, 2).into(), -3); // √-3
960        let q3m2 = QuadraticInt::from_coeffs((3, -2).into(), -3); // 2 - √-3
961        let q23 =  QuadraticInt::from_coeffs((2, 3).into(), -3);
962        let qm23 = QuadraticInt::from_coeffs((-2, 3).into(), -3);
963        let q5m3 = QuadraticInt::from_coeffs((5, -3).into(), -3);
964
965        assert_eq!(q12.conj(), q3m2);
966        assert_eq!(q23.conj(), q5m3);
967        assert_eq!(q12.norm(), 7);
968        assert_eq!(q23.norm(), 19);
969        assert_eq!(q12 * q12, QuadraticInt::from_coeffs((-3, 8).into(), -3));
970
971        for &v1 in &[q12, qm12, q3m2, q23, qm23, q5m3] {
972            for &v2 in &[q12, qm12, q3m2, q23, qm23, q5m3] {
973                assert_eq!(v1 + v2 - v1, v2, "add/sub test failed between {:?} and {:?}", v1, v2);
974                assert_eq!(v1 * v2 / v1, v2, "mul/div test failed between {:?} and {:?}", v1, v2);
975            }
976        }
977    }
978}
num_irrational/quadratic/integer.rs

num_irrational/quadratic/
integer.rs
