number_theory/arithmetic/
mpz_ent.rs

1use crate::ntrait::NumberTheory;
2use crate::result::NTResult;
3use crate::structs::{Certificate,Factorization};
4
5use crate::arithmetic::inlineops::*;
6use crate::arithmetic::mpz::Mpz;
7use crate::arithmetic::sign::Sign;
8//use crate::arithmetic::sliceops::*;
9use crate::data::primes::PRIMELIST;
10
11use crate::data::primes::MERSENNE_LIST;
12
13use std::cmp::Ordering;
14
15impl NumberTheory for Mpz {
16
17    fn is_unit(&self) -> bool{
18	if self.limbs.len() == 1 && self.limbs[0] == 1{
19	  return true;
20	}
21	false
22    }
23
24    fn rng() -> Self {
25        Mpz::unchecked_new(Sign::Positive, vec![u64::rng(), u64::rng(), u64::rng(), u64::rng()])
26    }
27    
28    fn residue(&self, ring: Self) -> Self{
29      if ring.clone() == Mpz::zero(){
30        return self.clone()
31      } 
32      let rem = self.ref_euclidean(&ring).1;
33      if !self.is_positive(){
34        return self.ref_subtraction(&rem)
35      }
36      rem
37    }
38
39    fn euclidean_div(&self, other: Self) -> (Self, Self) {
40        let (mut quo, mut rem) = self.ref_euclidean(&other);
41        if self.sign == Sign::Negative && other.sign == Sign::Negative {
42            rem.neg();
43        } else if self.sign == Sign::Positive && other.sign == Sign::Negative {
44            quo.neg();
45        } else if self.sign == Sign::Negative && other.sign == Sign::Positive {
46            quo.neg();
47            rem.neg();
48        }
49
50        quo.fix_zero();
51        rem.fix_zero();
52        (quo, rem)
53    }
54
55    fn quadratic_residue(&self, n: Self) -> Self {
56        if n == Mpz::zero() {
57            return Mpz::zero();
58        }
59
60        let mut p = self.clone();
61
62        if p.sign == Sign::Negative {
63            p = p.add_modinv(&n);
64        }
65
66        if n.is_power_of_two() {
67            let reducer = n.ref_subtraction(&Mpz::one());
68            return p.ref_product(&p).and(&reducer);
69        }
70
71        p.ref_product(&p).ref_euclidean(&n).1
72    }
73    
74    fn checked_quadratic_residue(&self, n: Self) -> NTResult<Self>{
75        NTResult::Eval(self.quadratic_residue(n))
76    }
77
78    fn product_residue(&self, other: Self, n: Self) -> Self {
79        if n == Mpz::zero() {
80            return Mpz::zero();
81        }
82
83        let mut p = self.clone();
84        let mut q = other.clone();
85
86        if p.sign == Sign::Negative {
87            p = p.add_modinv(&n);
88        }
89
90        if q.sign == Sign::Negative {
91            q = q.add_modinv(&n)
92        }
93
94        if n.is_power_of_two() {
95            let reducer = n.ref_subtraction(&Mpz::one());
96            return p.ref_product(&q).and(&reducer);
97        }
98
99        p.ref_product(&q).ref_euclidean(&n).1
100    }
101    
102    fn checked_product_residue(&self, other: Self, n: Self) -> NTResult<Self>{
103       NTResult::Eval(self.product_residue(other,n))
104    }
105    
106
107    fn exp_residue(&self, y: Self, n: Self) -> Self {
108        if n == Mpz::zero() {
109            match i64::try_from(y) {
110                Ok(x) => {
111                  if x < 1i64{
112                    panic!("No inverse")
113                  }
114                return self.pow(x as u64)
115                },
116                Err(_) => panic!("Incomputably large result"),
117            }
118        }
119
120        let mut p = self.clone();
121
122        if p.sign == Sign::Negative {
123            p = p.add_modinv(&n);
124        }
125
126        p.u_mod_pow(&y, &n)
127    }
128
129    fn checked_exp_residue(&self, y: Self, n: Self) -> NTResult<Self> {
130        if n == Mpz::zero() {
131            match i64::try_from(y) {
132                Ok(x) => {
133                  if x < 0i64{
134                    return NTResult::DNE
135                  }
136                   return NTResult::Eval(self.pow(x as u64)) 
137                },
138                Err(_) => return NTResult::CompExceeded,
139            }
140        }
141
142        let mut p = self.clone();
143
144        if p.sign == Sign::Negative {
145            p = p.add_modinv(&n);
146        }
147        
148        if !y.is_positive(){
149           let (gcd,base,_) = self.eea(n.clone());
150           if !gcd.is_unit(){
151             return NTResult::DNE
152           }
153           return NTResult::Eval(base.u_mod_pow(&y,&n))
154        }
155
156        NTResult::Eval(p.u_mod_pow(&y, &n))
157    }
158
159    fn fermat(&self, base: Self) -> bool{
160       base.exp_residue(self.ref_subtraction(&Mpz::one()),self.clone()).is_unit()
161    }
162    
163    
164    fn strong_fermat(&self, base: Self) -> bool {
165       if self.is_even(){
166          return self.fermat(&base);
167       }
168        self.sprp(base.clone())
169    }
170
171    fn is_prime(&self) -> bool {
172        if self.len() < 3 {
173            return self.to_u128().unwrap().is_prime();
174        }
175        
176        if !self.trial_div() {
177            return false;
178        }
179
180        if self.is_fermat() {
181            return false;
182        }
183
184        match self.is_mersenne() {
185            Some(x) => {
186                if MERSENNE_LIST.contains(&(x as u32)) {
187                    return true;
188                } else if x < 57885161 {
189                    return false;
190                } else {
191                    // actually incomputable but left for correctness
192                    return self.llt(x);
193                }
194            }
195            None => (),
196        }
197        
198    #[cfg(not(feature = "parallel"))]
199      {
200        if !self.sprp(Mpz::two()) {
201            return false;
202        }
203
204        if !self.form_check() {
205            return false;
206        }
207
208        if !self.jacobi_check_mpz() {
209            return false;
210        }
211
212        if !self.weighted_sprp() {
213            return false;
214        }
215        true
216       } // end single-thread block
217           #[cfg(feature = "parallel")]
218       {
219               if self.len() < 9 {
220            // if less than 2^6400 then serial
221            if !self.sprp(&Mpz::two()) {
222                return false;
223            }
224            if !self.form_check() {
225                return false;
226            }
227            if !self.jacobi_check_mpz() {
228                return false;
229            }
230            if !self.weighted_sprp() {
231                return false;
232            }
233            return true;
234        } else {
235            // parallelize the checks which roughly cost the same
236            let a = self.clone();
237            let b = self.clone();
238            let c = self.clone();
239            let two = std::thread::spawn(move || a.is_sprp(&Mpz::two()));
240            let jacobi = std::thread::spawn(move || b.jacobi_check_mpz());
241            let rand = std::thread::spawn(move || c.weighted_sprp());
242
243            if !self.form_check() {
244                return false;
245            }
246
247            return two.join().unwrap() & jacobi.join().unwrap() & rand.join().unwrap();
248        }
249       } // end parallel block
250    }
251    
252    fn prime_proof(&self) -> Certificate<Self>{
253    
254        if self.clone() == Mpz::two() {
255            return Certificate::new(self.clone(),Mpz::from(3u64),vec![]);
256        }
257
258        let x_minus = self.ref_subtraction(&Mpz::one());
259        let fctrs = x_minus.factor().unwrap()
260            .factor_iter().map(|x| x_minus.ref_euclidean(x).0)
261            .collect::<Vec<Self>>();
262         //println!("{:?}",fctrs);
263        loop {
264            // loops until we get a 
265
266            let mut witness = Mpz::from(u64::rng()).ref_euclidean(&self.ref_subtraction(&Mpz::two())).1.ref_addition(&Mpz::two());
267            //println!("witness {}",witness);    
268            'witness: loop {
269
270                if witness.coprime(self.clone()) {
271                	  break 'witness;
272		        }
273    
274                witness.successor();
275            }
276              
277	   if !self.strong_fermat(witness.clone()){
278		return Certificate::new(self.clone(),witness.clone(),fctrs)	
279	   }
280
281           
282            'inner: for (idx, i) in fctrs.iter().enumerate() {
283                if witness.exp_residue(i.clone(), self.clone()).is_unit(){
284                    break 'inner;
285                }
286                if idx == fctrs.len() - 1 {
287                    return Certificate::new(self.clone(),witness.clone(),fctrs);
288                }
289            }
290	  }
291    
292    
293    }
294    
295    //#[cfg(not(feature="parallel"))]
296    fn prime_list(&self, sup: Self) -> Vec<Self> {
297        //currently only supports less than, and not inclusive
298
299        let mut delta = self.ref_subtraction(&sup);
300        delta.successor();
301        match u64::try_from(delta) {
302            Ok(x) => {
303                let mut primevector = vec![];
304                let mut p = self.clone();
305              //  let primevec = 17000u64.prime_list(&30_000);
306                for _ in 0..x{
307                    if p.is_prime() {
308                        primevector.push(p.clone())
309                    }
310                    p.successor();
311                }
312                primevector
313            }
314            Err(_) => panic!("Incomputably large interval"),
315        }
316        //  return Vec::new();
317    }
318    
319    fn nth_prime(&self) -> NTResult<Self> {
320       
321        let mut count = Mpz::zero();
322        let mut start = Mpz::zero();
323        if *self == Mpz::zero(){
324          return NTResult::DNE
325        }
326        loop {
327            start.successor();
328            if start.is_prime() {
329                count.successor()
330            }
331            if count.u_cmp(self) == Ordering::Equal {
332                return NTResult::Eval(start);
333            }
334        }
335    }
336
337
338    fn pi(&self) -> Self {
339      match self.to_u128(){
340         Some(x) => return Mpz::from(x.pi()),
341         None => (),
342      }
343        let mut count = 0u64;
344        let mut start = Mpz::zero();
345        loop {
346            start.successor();
347            if start.is_prime() {
348                count += 1;
349            }
350            if start.u_cmp(self) == Ordering::Equal {
351                return Mpz::from(count);
352            }
353        }
354    }
355
356
357    fn prime_gen(x: u32) -> NTResult<Mpz> {
358        if x < 128 {
359         return u128::prime_gen(x).map(Mpz::from)  
360        }
361        
362        let mut form = Mpz::one().shl(x as usize-1);
363
364        let bitlength = form.ref_subtraction(&Mpz::one());
365
366        form.successor();
367
368        loop {
369            let mut k = Mpz::rand(x as usize+1);
370        
371            k.mut_and(&bitlength);
372
373            k.mut_or(&form);
374            assert_eq!(k.bit_length() as u32,x);
375            if k.is_prime() {
376                return NTResult::Eval(k);
377            }
378        }
379    } 
380
381    fn gcd(&self, other: Self) -> Self {
382        let mut a = self.clone();
383        let mut b = other.clone();
384
385        while b != Mpz::zero() {
386            let t = b.clone();
387
388            b = a.ref_euclidean(&b).1;
389
390            b.normalize();
391            a = t;
392        }
393        a
394    }
395    
396    // Z variant
397    fn extended_gcd(&self, other: Self) -> (Self, Self, Self) {
398        let mut gcd = self.clone();
399        let mut new_r = other.clone();
400        let mut bezout_1 = Mpz::one();
401        let mut new_s = Mpz::zero();
402        let mut bezout_2 = Mpz::zero();
403        let mut new_t = Mpz::one();
404
405        while new_r != Mpz::zero() {
406            let (quo, _rem) = gcd.euclidean_div(new_r.clone());
407            let mut temp = new_r.clone();
408            new_r = gcd.ref_subtraction(&quo.ref_product(&temp));
409            gcd = temp.clone();
410
411            temp = new_s.clone();
412            new_s = bezout_1.ref_subtraction(&quo.ref_product(&temp));
413            bezout_1 = temp.clone();
414
415            temp = new_t.clone();
416            new_t = bezout_2.ref_subtraction(&quo.ref_product(&temp));
417            bezout_2 = temp.clone();
418        }
419        (gcd, bezout_1, bezout_2)
420    }
421
422    fn lcm(&self, other: Self) -> NTResult<Self> {
423        if self == &Mpz::zero() && other  == Mpz::zero(){
424          return NTResult::Eval(Mpz::zero())
425        }
426        
427        let cf = self.gcd(other.clone());
428        NTResult::Eval(self.euclidean_div(cf).0.ref_product(&other))
429    }
430
431// FIXME 
432    fn factor(&self) -> NTResult<Factorization<Self>> {
433        let mut n = self.clone();
434        let mut factors: Factorization<Self> = Factorization::new();
435        let twofactor = n.trailing_zeros();
436
437        if twofactor > 0 {
438            n.mut_shr(twofactor as usize);
439            factors.add_factor(Mpz::from(2u64));
440            factors.add_power(twofactor);
441        }
442
443        for i in PRIMELIST[1..].iter() {
444            // skips two as it has already been eliminated
445            let (quo, rem) = n.word_div(*i as u64);
446
447            if rem == 0 {
448                let mut count = 1u64;
449                factors.add_factor(Mpz::from(*i as u64));
450                n = quo;
451
452                'inner: loop {
453                    let (inner_quo, inner_rem) = n.word_div(*i as u64);
454
455                    if inner_rem != 0 {
456                        break 'inner;
457                    }
458                    n = inner_quo;
459                    count += 1;
460                }
461
462                factors.add_power(count);
463            }
464        }
465        n.normalize();
466
467        if n.is_unit(){
468             //       println!("Final {}",factors);
469            return NTResult::Eval(factors);
470        }
471
472        if n.sprp_check(5) {
473            // stops if prime
474            factors.add_factor(n.clone());
475            factors.add_power(1);
476            return NTResult::Eval(factors);
477        }
478
479        'outer: while !n.is_unit() {
480            let k = n.rho_mpz();
481
482            factors.add_factor(k.clone());
483            let mut count = 0u64;
484
485            'secinner: loop {
486                let (inner_quo, inner_rem) = n.ref_euclidean(&k);
487
488                if !inner_rem.is_zero() {
489                    break 'secinner;
490                }
491                n = inner_quo;
492                n.normalize(); // remove ?
493                count += 1;
494            }
495            factors.add_power(count);
496
497            if n.sprp_check(5) {
498                // stops if  n is prime
499                factors.add_factor(n);
500                factors.add_power(1);
501                break 'outer;
502            }
503        }
504
505        NTResult::Eval(factors)
506    }
507    
508
509    /// Returns the integer component of sqrt(x)
510    fn sqrt(&self) -> (Self, Self) {
511        if self.len() < 3 {
512            let k = Mpz::from(self.to_u128().unwrap().sqrt().0 as u64);
513            if !self.is_positive() {
514                return (k, Mpz::one());
515            }
516            return (k, Mpz::zero());
517        }
518        let mut est = self.shr(((self.bit_length() / 2) - 1) as usize).abs();
519
520        loop {
521            let s = est.clone();
522            let t = s.ref_addition(&self.euclidean_div(s.clone()).0);
523            est = t.shr(1);
524            remove_lead_zeros(&mut est.limbs);
525            if est.u_cmp(&s) == Ordering::Greater || est.u_cmp(&s) == Ordering::Equal {
526                if self.sign == Sign::Negative {
527                    return (s, Mpz::one());
528                }
529                return (s, Mpz::zero());
530            }
531        }
532    }
533
534    fn nth_root(&self, n: Self) -> (Self, Self) {
535        let y = u64::try_from(n).unwrap();
536        let shift = ((self.bit_length() / y) - 1) * (y - 1);
537        let mut est = self.shr(shift as usize).abs();
538        let scalar = Mpz::from(y - 1);
539        let ymp = Mpz::from(y);
540
541        loop {
542            let s = est.clone();
543            let t = s
544                .ref_product(&scalar)
545                .ref_addition(&self.euclidean_div(s.pow(y - 1)).0);
546            est = t.euclidean_div(ymp.clone()).0;
547            if est.u_cmp(&s) == Ordering::Greater || est.u_cmp(&s) == Ordering::Equal {
548                if self.sign == Sign::Negative && self.is_even() {
549                    return (s, Mpz::one());
550                }
551                return (s, Mpz::zero());
552            }
553        }
554    }
555    // FIXME : Use prime exponents and fix for x < 0
556   fn max_exp(&self) -> (Self,Self){
557     let mut expo = Mpz::from(self.bit_length());
558     let mut flag = true;
559     if  self.is_positive(){
560       flag = false
561     }
562      loop{
563        let mut base = self.abs().nth_root(expo.clone()).0;
564         if base.pow(u64::try_from(expo.clone()).unwrap()) == self.abs(){
565         if flag{
566           base.neg();
567         }
568           return(base,expo.clone())
569         }
570         expo.predecessor();
571         if u64::try_from(expo.clone()).unwrap() == 1{
572         let mut val = self.clone();
573          if flag {
574            val.neg()
575          }
576           return (val,Mpz::one())
577         }
578      }
579    }
580
581
582    
583    fn radical(&self) -> NTResult<Self>{
584       match u64::try_from(self.clone()){
585        Ok(x) => x.radical().map(Mpz::from),
586        Err(_) => {
587          let mut rad = Mpz::one();
588           for i in self.factor().unwrap().factor_iter(){
589              rad = rad.ref_product(i)
590           }
591          NTResult::Eval(rad)
592         },
593       }
594    }
595
596    fn k_free(&self, x: Self) -> bool {
597       let x = u64::try_from(x).unwrap();
598        for i in self.factor().unwrap().power_iter(){
599            if *i == x {
600                return false;
601            }
602        }
603        true
604    }
605
606    fn euler_totient(&self) -> Self {
607    
608       match self.to_u128(){
609      
610         Some(x)=> return Mpz::from(x.euler_totient()),
611         None => (),
612       }
613       
614        let factors = self.factor().unwrap();
615
616        let mut denominator = Mpz::one();
617        let mut numerator = Mpz::one();
618
619        for i in factors.factor_iter() {
620            denominator = denominator.ref_product(i)
621        }
622        for i in factors.factor_iter() {
623            let mut fctr = i.clone();
624            fctr.predecessor();
625            numerator = numerator.ref_product(&fctr);
626        }
627
628        (self.ref_euclidean(&denominator).0).ref_product(&numerator)
629    }
630
631    fn jordan_totient(&self, k: Self) -> NTResult<Self> {
632        if self.clone() == Mpz::zero() || self.is_unit(){
633          return NTResult::Eval(self.clone())
634        }
635        let fctr = self.factor().unwrap();
636        let mut denom = Mpz::one();
637        let mut numer = Mpz::one();
638        let negone = Mpz::unchecked_new(Sign::Negative, vec![1]);
639        let value = u64::try_from(k).unwrap();
640        for i in fctr.factor_iter() {
641            let pow = i.pow(value);
642
643            denom = denom.ref_product(&pow);
644
645            numer = numer.ref_product(&pow.ref_addition(&negone));
646        }
647
648        NTResult::Eval(
649            numer
650                .ref_product(&self.pow(value))
651                .ref_euclidean(&denom)
652                .0,
653        )
654    }
655    
656    fn exponent(&self) -> NTResult<Self>{
657       
658       match self.to_u128(){
659         Some(x)=> return x.exponent().map(Mpz::from),
660         None => (),
661       }
662       
663       let fctr = self.factor().unwrap();
664       //let base = fctr.iter()..map(|z| z.clone()).collect::<Vec<Self>>();
665       let power = fctr.power_iter().cloned().collect::<Vec<u64>>();
666       let two = Mpz::two();
667       let mut result = Mpz::one();
668      for (idx,el) in fctr.factor_iter().enumerate(){
669        if el == &two && power[0] > 2{
670         let phi = ((el.pow(power[idx]).ref_euclidean(el).0).ref_product(&el.ref_subtraction(&Mpz::one()))).shr(1);
671          result = result.lcm(phi).unwrap();
672        }
673       else{
674         let phi = (el.pow(power[idx]).ref_euclidean(el).0).ref_product(&el.ref_subtraction(&Mpz::one()));
675         result = result.lcm(phi).unwrap();
676       } 
677      }
678     NTResult::Eval(result)
679    }
680
681    fn dedekind_psi(&self, k: Self) -> NTResult<Self> {
682     if *self == Mpz::zero(){
683         return NTResult::Infinite
684       }
685        let k2 = k.shl(1);
686     self.jordan_totient(k2).map(|y| y.ref_euclidean(&self.jordan_totient(k).unwrap()).0)
687    }
688
689    fn legendre(&self, p: Self) -> i8 {
690        let mut p_minus = p.clone();
691        p_minus.set_sign(Sign::Positive);
692        p_minus.predecessor();
693
694        let pow = p_minus.ref_euclidean(&Mpz::two()).0;
695        let k = self.exp_residue(pow, p);
696
697        if k == Mpz::one() {
698            return 1;
699        };
700        if k == p_minus {
701            return -1;
702        };
703        0i8
704    }
705
706    fn checked_legendre(&self, p: Self) -> NTResult<i8> {
707        let mut plus = p.clone();
708        plus.set_sign(Sign::Positive);
709        if plus == Mpz::two() {
710            return NTResult::Undefined;
711        }
712        match plus.is_prime() {
713            true => NTResult::Eval(self.legendre(plus)),
714            false => NTResult::Undefined,
715        }
716    }
717
718    fn liouville(&self) -> i8 {
719        match self.to_u128(){
720          Some(an) => return an.liouville(),
721          None => (),
722        }
723        
724        let mut primeomega = 0;
725
726        for i in self.factor().unwrap().power_iter() {
727            primeomega+=i;
728        }
729        if primeomega&1==0 {
730            return 1;
731        }
732        -1
733    }
734    
735    fn derivative(&self) -> NTResult<Self> {
736       if *self == Mpz::zero(){
737         return NTResult::Eval(Mpz::zero())
738       }
739       if self.is_unit(){
740        return NTResult::Eval(Mpz::zero())
741       }
742       let fctr = self.factor().unwrap();
743       let mut sum = Mpz::zero();
744       
745     for i in 0..fctr.base.len() / 2 {
746        sum.mut_addition(Mpz::from(fctr.power[i]).ref_product(&self.ref_euclidean(&fctr.base[i]).0))
747      }
748    NTResult::Eval(sum)
749    }
750
751    fn mangoldt(&self) -> f64 {
752     match self.to_u128(){
753          Some(eburum) => return eburum.mangoldt(),
754          None => (),
755        }
756        let fctr = self.factor().unwrap();
757        if fctr.base.len() != 2 {
758            return 0f64;
759        }
760        fctr.base[0].ln()
761    }
762    
763    
764    fn mobius(&self) -> i8 {
765     match self.to_u128(){
766          Some(eburum) => return eburum.mobius(),
767          None => (),
768        }
769      let fctr = self.factor().unwrap();
770      /*
771      if fctr.len() == 1{ // if only one factor then return -1
772         return -1
773      }
774      */
775      for i in fctr.power_iter(){
776        if *i  > 1{
777         return 0
778        }
779      }
780      
781      let fctrsum = fctr.power_iter().sum::<u64>();
782      if fctrsum&1 == 1{// if odd number of factors and square free
783        return -1
784      }
785      1
786    }
787
788    fn jacobi(&self, k: Self) -> NTResult<i8> {
789      if k.sign == Sign::Positive && k != Mpz::zero() && !k.is_even() {
790        let mut n = self.clone();
791        let mut p = k.clone();
792        let mut t = 1i8;
793        n = n.euclidean_div(k).1;
794
795        while n != Mpz::zero() {
796            let zeros = n.trailing_zeros();
797            n.mut_shr(zeros as usize);
798
799            if (p.congruence_u64(8, 3) || p.congruence_u64(8, 5)) && (zeros % 2 == 1) {
800                t = -t
801            }
802
803            std::mem::swap(&mut n, &mut p);
804
805            if n.congruence_u64(4, 3) && p.congruence_u64(4, 3) {
806                t = -t;
807            }
808
809            n = n.euclidean_div(p.clone()).1;
810        }
811
812        if p == Mpz::one() {
813            NTResult::Eval(t)
814        } else {
815           NTResult::Eval(0)
816        }
817        }
818        else{
819          NTResult::Undefined
820        }
821    }
822    
823     fn kronecker(&self, k: Self) -> i8{
824     let mut multiplier = 1;
825    
826    if !k.is_positive() && !self.is_positive(){
827      multiplier = -1;
828    }
829     let x = self.clone();
830     if k.is_zero(){
831      if x.is_unit(){
832         return 1
833      }
834     return 0
835    }
836   if k.is_unit(){
837      return 1
838   }
839   let fctr = k.factor().unwrap();
840   let mut start = 0;
841   let mut res = 1;
842   let two = Mpz::two();
843   if fctr.base[0] == two{
844     start = 1;
845     if x.is_even(){
846     res = 0;
847     }
848     else if x.congruence_u64(8,1) || x.congruence_u64(8,7){
849      res=1
850     }
851     else{
852       res = (-1i8).pow(fctr.power[1] as u32)
853     }
854   }
855   if fctr.base[0] == two && fctr.base.len() == 2{
856     return res*multiplier
857   }
858   for i in start..fctr.base.len(){
859     res*=self.legendre(fctr.base[i].clone()).pow(fctr.power[i] as u32);
860   }
861   res*multiplier
862}
863    
864    fn smooth(&self) -> NTResult<Self> {
865       if *self == Mpz::zero(){
866         return NTResult::Infinite
867       }
868       if self.is_unit(){
869        return NTResult::DNE
870       }
871      NTResult::Eval(self.factor().unwrap().max())
872    }
873    
874    
875
876    fn is_smooth(&self, b: Self) -> bool {
877     match self.smooth(){
878      NTResult::Infinite => false,
879      NTResult::Eval(x) => {
880      let mut k = b.clone();
881        k.set_sign(Sign::Positive);
882        matches!(x.u_cmp(&k), Ordering::Less)
883      }, 
884      _=> false,
885     }
886   }
887   
888   fn ord(&self, n: Self) -> NTResult<Self>{
889
890        let ord_2 = |a: Mpz, p: u64| -> Mpz{
891        let modulo = Mpz::one().shl(p as usize);
892        let mut b = a.ref_euclidean(&modulo).1;
893   
894         if b == Mpz::one(){
895            return Mpz::one();
896         }
897      for i in 1..p{
898         b = b.quadratic_residue(modulo.clone());
899         if b == Mpz::one(){
900           return Mpz::one().shl(i as usize);
901         }
902    }
903    return p.into();
904    };
905
906// Given ord(a,p)  calculate ord(a,p^n)
907    let  pp_ord = |a: Mpz, b: Mpz, p: Mpz, e: u64| -> Mpz{
908        let fullpow = p.pow(e);
909
910     	for i in 0..e+1{
911          let partial = b.ref_product(&p.pow(i));
912          if a.exp_residue(partial.clone(),fullpow.clone()) == Mpz::one(){
913             return partial;
914           }
915        }
916    return fullpow;
917    };
918
919    let p_ord = |a: Mpz, p: Mpz| -> Mpz{
920
921 	  let pminus = p.ref_subtraction(&Mpz::one());  
922  	  let fctr = pminus.factor().unwrap();
923   	  let mut m = pminus.clone();
924
925 	for i in fctr.pair_iter(){
926
927 	    for _ in 0..*i.1{
928              let reduced = m.ref_euclidean(&i.0).0;
929
930	        if a.exp_residue(reduced.clone(),p.clone()) == Mpz::one(){
931                   m = reduced;
932                }
933          	else{
934                   break;
935          	}
936            }
937        }
938       m
939     };
940
941
942    if !self.coprime(n.clone()){
943       return NTResult::DNE;
944    }
945
946    let fctr = n.clone().factor().unwrap();
947    let mut fullord = Mpz::one();
948
949    for i in fctr.pair_iter(){
950     let mut ord : Self;
951      if *i.0 == Mpz::two(){
952         ord = ord_2(self.clone(),*i.1);
953      }
954      else{
955        ord = p_ord(self.clone(),i.0.clone());
956        if *i.1 > 1{
957           ord=pp_ord(self.clone(),ord,i.0.clone(),*i.1);
958        }
959      }
960       fullord = fullord.lcm(ord).unwrap(); 
961    }
962    NTResult::Eval(fullord)
963       
964   }
965   
966   fn mul_inverse(&self, n: Self) -> NTResult<Self>{
967       let (xinv,_yinv,gcd) = self.extended_gcd(n.clone());
968       if gcd != Mpz::one(){
969          return NTResult::DNE;
970       }
971       NTResult::Eval(xinv.residue(n))
972   }
973   
974}
number_theory/arithmetic/mpz_ent.rs

number_theory/arithmetic/
mpz_ent.rs
