number_theory/primitive/
eightbytes.rs

1use crate::structs::{Certificate,Factorization};
2use crate::ntrait::{NumberTheory,MontArith,Reduction};
3use crate::primitive::factorprim::{poly_factor,drbg};
4use crate::data::primes::PRIMELIST;
5use crate::result::NTResult;
6
7impl NumberTheory for u64 {
8
9    fn is_unit(&self) -> bool{
10	if *self == 1{
11	    return true;
12	}
13	false
14    }
15
16    fn rng() -> Self {
17        #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
18    {
19        if is_x86_feature_detected!("rdrand"){// USE RDRAND if possible
20            let mut x: u64 = 0;
21            unsafe { core::arch::x86_64::_rdrand64_step(&mut x) };
22            return x;
23        }// If processor is x86 and does not support RDRAND use xor shift
24       return drbg(std::time::SystemTime::now().duration_since(std::time::UNIX_EPOCH).unwrap().as_nanos() as u64)
25    }
26        #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
27    {// All other architectures use xor shift
28        drbg(std::time::SystemTime::now().duration_since(std::time::UNIX_EPOCH).unwrap().as_nanos() as u64)
29    }
30    }
31
32    fn residue(&self, ring: Self) -> Self{
33      if ring == 0{
34        return *self
35      }
36        *self % ring
37    }
38     
39    fn euclidean_div(&self, other: Self) -> (Self, Self) {
40        (*self / other, *self % other)
41    }
42    
43    fn mul_inverse(&self, n: Self) -> NTResult<Self>{
44       let (gcd,inv,_) = self.extended_gcd(n);
45       
46       if gcd != 1{
47          return NTResult::DNE
48       }
49       NTResult::Eval(inv)
50    }
51    
52    fn fermat(&self,base: Self) -> bool{
53        if *self&1 == 1{
54           let one = self.n_identity();
55           let mbase = base.to_mont(*self);
56           let inv = self.inv_2();
57           return mbase.mont_pow(one,*self-1,inv,*self)==one;
58       }
59       base.exp_residue(*self-1,*self)==1
60       
61    }
62
63    fn strong_fermat(&self, base: Self) -> bool {
64        if *self&1==0{
65            return self.fermat(base);
66        }
67        let inv = self.inv_2();
68        let tzc = self.wrapping_sub(1).trailing_zeros();
69        let one = self.n_identity();
70        let oneinv = self.wrapping_sub(one);
71        let b = machine_prime::to_mont(base,*self);
72    
73        machine_prime::strong_fermat(*self, tzc, b, one, oneinv, inv)
74    }
75
76    fn is_prime(&self) -> bool {
77        machine_prime::is_prime(*self)
78    }
79
80    fn prime_proof(&self) -> Certificate<Self> {
81
82        if *self == 2 {
83            return Certificate::new(*self,3,vec![]);
84        }
85
86        let x_minus = *self - 1;
87        let fctrs = x_minus.factor().unwrap()
88            .factor_iter().map(|x| x_minus/ *x)
89            .collect::<Vec<Self>>();
90
91        loop {
92            // loops until we get a 
93
94            let mut witness = Self::rng() % (*self - 2) + 2;
95            
96            'witness: loop {
97
98                if witness.coprime(*self) {
99
100		  break 'witness;
101		}
102              
103                witness += 1;
104            }
105              
106	   if !self.strong_fermat(witness){
107		return Certificate::new(*self,witness,fctrs)	
108	   }
109
110           
111            'inner: for (idx, i) in fctrs.iter().enumerate() {
112                if witness.exp_residue(*i, *self).is_unit(){
113                    break 'inner;
114                }
115                if idx == fctrs.len() - 1 {
116                    return Certificate::new(*self,witness,fctrs);
117                }
118            }
119	  }
120    }
121
122    fn prime_list(&self, sup: Self) -> Vec<Self> {
123       
124        let inf = std::cmp::min(*self, sup);
125        let mut hi = std::cmp::max(*self, sup);
126
127        hi = hi.saturating_add(1);
128
129        let mut primevector = vec![];
130
131        for i in inf..hi {
132            if i.is_prime() {
133                primevector.push(i)
134            }
135        }
136        primevector
137    }
138
139    fn nth_prime(&self) -> NTResult<Self> {
140        let mut count = 0u64;
141        let mut start = 0u64;
142        
143/*        if *self < 203280221{
144          return (*self as u32).nth_prime().map(|y| y as u64)
145        }
146  */      
147        if *self > 425656284035217743{
148          return NTResult::Overflow
149        }
150    
151        loop {
152            start += 1;
153
154            if start == Self::MAX {
155                return NTResult::Overflow;
156            }
157            if start.is_prime() {
158                count += 1;
159            }
160            if count == *self {
161                return NTResult::Eval(start);
162            }
163        }
164    }
165
166    fn pi(&self) -> Self {
167      if  self.reducible(){
168         return (*self as u32).pi() as u64
169       }
170       
171        let mut count = 203280221u64;
172        for i in (1u64<<32)..*self {
173            if i.is_prime() {
174                count += 1;
175            }
176        }
177        count
178    }
179
180    fn prime_gen(k: u32) -> NTResult<Self> {
181        if k > 64 {
182           return NTResult::Overflow;
183        }
184        if k < 33 {
185            return u32::prime_gen(k).map(|x| x as u64);
186        }
187        let form = (1 << (k - 1)) + 1;
188        let bitlength = form - 2;
189
190        loop {
191            let p = u64::rng();
192            if ((p & bitlength) | form).is_prime() {
193                return NTResult::Eval((p & bitlength) | form);
194            }
195        }
196    }
197
198    fn factor(&self) -> NTResult<Factorization<Self>> {
199       /* if self < &4294967295{
200          return (*self as u32).factor().iter().map(|x| *x as u64).collect::<Vec<u64>>()
201        }
202       */ 
203        if *self == 0{
204           return NTResult::InfiniteSet;
205        }
206        
207        if *self == 1{
208           return NTResult::DNE
209        }
210        let mut n = *self;
211        let twofactors = n.trailing_zeros();
212        n >>= twofactors;
213
214        let mut factors: Factorization<u64> = Factorization::new();
215
216        if twofactors > 0 {
217            factors.add_factor(2);
218            factors.add_power(twofactors as u64);
219        }
220
221        for i in PRIMELIST[1..128].iter() {
222            // strips out small primes
223            if n % *i as u64 == 0 {
224                factors.add_factor(*i as u64);
225                let mut count = 0u64;
226                while n % *i as u64 == 0 {
227                    count += 1;
228                    n /= *i as u64;
229                }
230                factors.add_power(count);
231            }
232        }
233
234        if n == 1 {
235            return  NTResult::Eval(factors);
236        }
237
238        if n.is_prime() {
239            factors.add_factor(n);
240            factors.add_power(1);
241            return  NTResult::Eval(factors);
242        }
243
244        while n != 1 {
245            let k = poly_factor(n);
246            factors.add_factor(k);
247            let mut count = 0u64;
248            while n % k == 0 {
249                count += 1;
250                n /= k;
251            }
252            factors.add_power(count);
253            
254            if n.is_prime() {
255            factors.add_factor(n);
256            factors.add_power(1);
257            return  NTResult::Eval(factors);
258          }
259          
260        }
261        NTResult::Eval(factors)
262    }
263    
264
265
266    fn sqrt(&self) -> (Self, Self) {
267       // if *self < 0x100000000 {
268       //     return ((*self as u32).sqrt().0 as u64, 0);
269       // }
270        let mut est = (*self as f64).sqrt() as Self + 1;
271
272        loop {
273            let s = est;
274            let t = s + *self / s;
275            est = t >> 1;
276            if est >= s {
277                return (s, 0);
278            }
279        }
280    }
281
282    fn nth_root(&self, n: Self) -> (Self, Self) {
283        if n > 63 {
284            return (1, 0);
285        }
286
287        if n == 1 {
288            return (n, 0);
289        }
290
291        if n == 0 {
292            panic!("No integer is a zeroth factor ")
293        }
294
295        let mut est = (*self as f64).powf((n as f64).recip()) as Self + 1;
296
297        loop {
298            let s = est;
299            let t = (n - 1) * s + *self / s.pow(n as u32 - 1);
300            est = t / n;
301            if est >= s {
302                return (s, 0);
303            }
304        }
305    }
306
307    fn max_exp(&self) -> (Self,Self){
308    
309      for i in 1..64{
310      let p = 64-i;
311      let base = self.nth_root(p).0;
312         if base.pow(p as u32) == *self{
313           return(base,p)
314         }
315      }
316     (*self,1)
317    }    
318    
319    fn radical(&self) -> NTResult<Self> {
320      if self.reducible(){
321         return (*self as u32).radical().map(|kishum| kishum as u64)
322      }
323        self.factor().map(|y| y.factor_iter().product::<Self>())
324    }
325
326    fn k_free(&self, k: Self) -> bool {
327       if self.reducible(){
328         return (*self as u32).k_free(k as u32)
329      }
330      // FIXME Remove unwrap
331        let factors = self.factor().unwrap();
332        factors.k_free(k)
333    }
334
335    fn gcd(&self, other: Self) -> Self {
336        let mut a = *self;
337        let mut b = other;
338        
339        if b == 0 {
340            return a;
341        } else if a == 0 {
342            return b;
343        }
344
345        let min_two_factor = (a|b).trailing_zeros();
346        a >>= min_two_factor;
347        b >>= min_two_factor;
348        loop {
349            if b > a {
350                std::mem::swap(&mut b, &mut a);
351            }
352            a -= b;
353
354            if a == 0 {
355                return b << min_two_factor;
356            }
357            a >>= a.trailing_zeros();
358        }
359    }
360
361    fn extended_gcd(&self, other: Self) -> (Self, Self, Self) {
362     
363        let mut gcd: u64 = *self;
364        let mut new_r: u64 = other;
365        let mut bezout_1: u64 = 1;
366        let mut new_s: u64 = 0;
367        let mut bezout_2: u64 = 0;
368        let mut new_t: u64 = 1;
369
370        while new_r != 0 {
371            let quotient = gcd / new_r;
372            let mut temp_r: u64 = new_r;
373            let prod = quotient.wrapping_mul(temp_r);
374
375            new_r = gcd - prod;
376            gcd = temp_r;
377
378            let mut temp = new_s;
379            new_s = bezout_1.wrapping_sub(temp.wrapping_mul(quotient));
380            bezout_1 = temp;
381
382            temp = new_t;
383            new_t = bezout_2.wrapping_sub(temp.wrapping_mul(quotient));
384            bezout_2 = temp;
385        }
386        
387        if bezout_1 > 1<<63{
388           bezout_1=other-(bezout_1.wrapping_neg());
389        }
390        if bezout_2 > 1<<63{
391           bezout_2=*self-bezout_2.wrapping_neg();
392        }
393        (
394            gcd,
395            bezout_1,
396            bezout_2,
397        )
398    }
399    
400    fn lcm(&self, other: Self) -> NTResult<Self> {
401        if self == &0 && other == 0{
402           return NTResult::Eval(0)
403         }
404        let cf = self.gcd(other);
405        let (v, flag) = (*self / cf).overflowing_mul(other);
406        if flag {
407            return NTResult::Overflow;
408        }
409        NTResult::Eval(v)
410    }
411
412    fn euler_totient(&self) -> Self {
413      if self.reducible(){
414        return (*self as u32).euler_totient() as u64
415      }
416        let factors = self.factor().unwrap();
417        let numerator = factors.factor_iter().map(|x| x - 1u64).product::<u64>();
418        let denominator = factors.factor_iter().product::<u64>();
419        (self / denominator) * numerator
420    }
421
422    fn jordan_totient(&self, k: Self) -> NTResult<Self> {
423        if k > u32::MAX as u64 {
424            return NTResult::Overflow;
425        }
426        
427        if *self < 2{
428           return NTResult::Eval(*self)
429        }
430        
431        let (coef, flag) = self.overflowing_pow(k as u32);
432
433        if flag {
434            return NTResult::CompOverflow;
435        }
436
437        let mut denom = 1u64;
438        let mut numer = 1u64;
439
440        for i in self.factor().unwrap().factor_iter(){
441            let pow = i.pow(k as u32);
442
443            denom *= pow;
444
445            numer *= pow - 1;
446        }
447
448        NTResult::Eval(numer * (coef / denom))
449    }
450    
451    fn exponent(&self) -> NTResult<Self>{ 
452    
453       if self.reducible(){
454         return (*self as u32).exponent().map(|x| x as u64)
455       }
456    
457       if *self == 0{
458          return NTResult::Infinite;
459       }
460       let fctr = self.factor().unwrap();
461       let mut result = 1;
462       for i in 0..fctr.base.len(){
463       if fctr.base[0] == 2 && fctr.power[0] > 2{
464          let phi = 2u64<<(fctr.power[0]-1);
465          result=phi;
466       }
467       else{
468         let el = fctr.base[i];
469         let phi =  (el.pow(fctr.power[i] as u32)/el)*(el-1);
470         match result.lcm(phi){
471             NTResult::Overflow => {return NTResult::Overflow;},
472             NTResult::Eval(x) => result=x,
473             _=> (),
474         }
475       }
476       }
477        return NTResult::Eval(result);
478    }
479    
480    fn dedekind_psi(&self, k: Self) -> NTResult<Self> {
481        if *self == 0{
482         return NTResult::Infinite
483       }
484       
485        let (k2, flag) = k.overflowing_shl(1);
486        if flag {
487            return NTResult::Overflow;
488        }
489        self.jordan_totient(k2).map(|y| y/self.jordan_totient(k).unwrap())
490    }
491
492   fn quadratic_residue(&self, n: Self) -> Self {
493        if n == 0 {
494            return self.wrapping_mul(*self)
495        }
496        ((*self as u128 * *self as u128) % n as u128) as Self
497    }
498    
499    fn checked_quadratic_residue(&self, n: Self) -> NTResult<Self> {
500        if n == 0 {
501            return NTResult::from_option(self.checked_mul(*self),NTResult::Overflow)
502        }
503        NTResult::Eval(((*self as u128 * *self as u128) % n as u128) as Self)
504    }
505    
506
507    fn product_residue(&self, other: Self, n: Self) -> Self {
508        if n == 0 {
509            return self.wrapping_mul(other)
510        }
511        ((*self as u128 * other as u128) % n as u128) as Self
512    }
513    
514    fn checked_product_residue(&self, other: Self, n: Self) -> NTResult<Self> {
515        if n == 0 {
516            return NTResult::from_option(self.checked_mul(*self),NTResult::Overflow)
517        }
518        NTResult::Eval(((*self as u128 * other as u128) % n as u128) as Self)
519    }
520
521    fn exp_residue(&self, p: Self, modulus: Self) -> Self {
522        
523        if modulus == 0 {
524            return self.pow(p as u32);
525        }
526        
527        if modulus & 1 == 0 {
528        
529        let k = modulus.trailing_zeros() as u64;
530        let s = modulus >> k;
531
532        let reducer = (1 << k) - 1; // A shorthand for arithmetic over Z[2k]
533
534        let k_rem = self.even_pow(p,reducer);
535
536        let s_rem = self.odd_pow(p,s);
537        
538        let s_inv = s.inv_2()&reducer;
539    
540        let y = k_rem.wrapping_sub(s_rem).wrapping_mul(s_inv) & reducer;
541
542        s_rem + s * y //)%modulus
543    } else {
544        self.odd_pow(p,modulus) //%modulus
545    }
546    }
547
548    fn checked_exp_residue(&self, p: Self, modulus: Self) -> NTResult<Self> {
549        if modulus == 0 {
550            if  p > u32::MAX as u64 {
551                return NTResult::Overflow;
552            }
553            match self.checked_pow(p as u32) {
554                Some(x) => return NTResult::Eval(x),
555                None => return NTResult::Overflow,
556            };
557        }
558
559        NTResult::Eval(self.exp_residue(p,modulus))
560    }
561
562    fn legendre(&self, p: Self) -> i8 {
563        let k = self.exp_residue((p - 1) >> 1, p);
564        if k == 1 {
565            return 1;
566        };
567        if k == p - 1 {
568            return -1;
569        };
570        0i8
571    }
572
573    fn checked_legendre(&self, p: Self) -> NTResult<i8> {
574        if p == 2 || !p.is_prime() {
575            return NTResult::Undefined;
576        }
577        NTResult::Eval(self.legendre(p))
578    }
579
580    fn liouville(&self) -> i8 {
581      if self.reducible(){
582        return (*self as u32).liouville()
583      }
584        let primeomega = self.factor().unwrap().prime_omega();
585        if primeomega & 1 == 0 {
586            return 1;
587        }
588        -1
589    }
590    
591    fn derivative(&self) -> NTResult<Self> {
592       if *self < 94 {
593         return (*self as u32).derivative().map(|y| y as Self)
594       }
595       
596       let fctr = self.factor().unwrap();
597       let mut sum : u64 = 0;
598       
599       for (power,base) in fctr.power_iter().zip(fctr.factor_iter()){
600          match sum.checked_add((*power as u64)*(*self/base)){
601            Some(x) => sum =x,
602            None => return NTResult::Overflow,
603          }
604       }
605       NTResult::Eval(sum)
606       }
607
608    fn mangoldt(&self) -> f64 {
609     if self.reducible(){
610        return (*self as u32).mangoldt()
611     }
612      let base = self.max_exp().0;
613       if base.is_prime(){
614         return (base as f64).ln()
615       }
616       0f64
617    }
618    
619    fn mobius(&self) -> i8 {
620     //if self.reducible(){
621     //   return (*self as u32).mobius()
622     // }
623      let fctr = self.factor().unwrap();
624     // if fctr.len() == 1{ // if only one factor then return -1
625     //    return -1
626     // }
627      if !fctr.k_free(2){
628         return 0;
629      }
630      let fctrsum = fctr.prime_omega();//fctr[1..].iter().step_by(2).sum::<Self>();
631      if fctrsum&1 == 1{// if odd number of factors and square free
632        return -1
633      }
634      1
635    }
636
637    fn jacobi(&self, k: Self) -> NTResult<i8> {
638        if k == 0 && k % 2 == 0 {
639           return NTResult::Undefined;
640        }
641        let mut n = *self;
642        let mut p = k;
643        let mut t = 1i8;
644        n %= p;
645
646        while n != 0 {
647            let zeros = n.trailing_zeros();
648            n >>= zeros;
649
650            if (p % 8 == 3 || p % 8 == 5) && (zeros % 2 == 1) {
651                t = -t
652            }
653
654            std::mem::swap(&mut n, &mut p);
655            if n % 4 == 3 && p % 4 == 3 {
656                t = -t;
657            }
658            n %= p;
659        }
660
661        if p == 1 {
662            NTResult::Eval(t)
663        } else {
664            NTResult::Eval(0)
665        }
666    }
667
668     fn kronecker(&self, k: Self) -> i8{
669
670       if k==0 && *self==1{
671         return 1i8;
672       }
673       if k==0{
674          return 0i8;
675       }
676       if k&1==1{
677         self.jacobi(k).unwrap()
678       }
679       else{
680         let tzc = k.trailing_zeros();
681         let odd = k>>tzc;
682         let res = *self&7;
683         let mut even_component = 1i8;
684         if (res == 3 || res == 5) && tzc&1==1{
685            even_component = -1i8;
686         }
687         if *self&1==0{
688           even_component = 0i8;
689         }
690        even_component*self.jacobi(odd).unwrap()
691       }
692     }
693
694
695  fn smooth(&self) -> NTResult<Self> {
696       if *self == 0{
697         return NTResult::Infinite
698       }
699       if *self == 1{
700        return NTResult::DNE
701       }
702       self.factor().map(|x| x.max())
703  
704    }
705    
706    
707
708    fn is_smooth(&self, b: Self) -> bool {
709     match self.smooth(){
710      NTResult::Infinite => false,
711      NTResult::Eval(x) => x <= b, 
712      _=> false,
713     }
714   }
715   
716    fn ord(&self, n: Self) -> NTResult<Self>{
717      let ord_2 = |a: u64, p: u64| -> u64{
718        let modulo = (1u64<<p)-1;
719        let mut b = a&modulo;
720   
721         if b == 1{
722            return 1;
723         }
724      for i in 1..p{
725         b = b.wrapping_mul(b)&modulo;
726         if b == 1{
727           return 1<<i;
728         }
729    }
730    return p;
731    };
732
733// Given ord(a,p)  calculate ord(a,p^n)
734    let  pp_ord = |a: u64, b: u64, p: u64, e: u32| -> u64{
735     for i in 0..e+1{
736          if a.exp_residue(b*p.pow(i),p.pow(e)) ==1{
737             return b*p.pow(i);
738           }
739   }
740    return b*p.pow(e);
741    };
742
743    let p_ord = |a: u64, p: u64| -> u64{
744   
745   let fctr = (p-1).factor().unwrap();
746   
747   let mut m = p-1;
748   for i in fctr.pair_iter(){
749     for _ in 0..*i.1{
750          if a.exp_residue(m/ *i.0,p) == 1{
751            m = m/ *i.0;
752          }
753          else{
754            break;
755          }
756     }
757  }
758 m
759};
760
761
762    if self.gcd(n) != 1{
763       return NTResult::DNE;
764    }
765    let fctr = n.factor().unwrap();
766    let mut fullord = 1u64;
767    for i in fctr.pair_iter(){
768     let mut ord : Self;
769      if *i.0 == 2{
770         ord = ord_2(*self,*i.1);
771      }
772      else{
773        ord = p_ord(*self,*i.0);
774        if *i.1 > 1{
775           ord=pp_ord(*self,ord,*i.0,*i.1 as u32);
776        }
777      }
778       fullord = fullord.lcm(ord).unwrap(); 
779    }
780    NTResult::Eval(fullord)
781 }
782
783}
number_theory/primitive/eightbytes.rs

number_theory/primitive/
eightbytes.rs
