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();//machine_prime::mul_inv2(x);
68        let tzc = self.wrapping_sub(1).trailing_zeros();
69        let one = self.n_identity();
70        let oneinv = machine_prime::mont_prod(machine_prime::mont_sub(*self, one, *self), one, *self, inv);
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        if b == 0 {
339            return a;
340        } else if a == 0 {
341            return b;
342        }
343
344        let self_two_factor = a.trailing_zeros();
345        let other_two_factor = b.trailing_zeros();
346        let min_two_factor = std::cmp::min(self_two_factor, other_two_factor);
347        a >>= self_two_factor;
348        b >>= other_two_factor;
349        loop {
350            if b > a {
351                std::mem::swap(&mut b, &mut a);
352            }
353            a -= b;
354
355            if a == 0 {
356                return b << min_two_factor;
357            }
358            a >>= a.trailing_zeros();
359        }
360    }
361
362    fn extended_gcd(&self, other: Self) -> (Self, Self, Self) {
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: u64 = new_r;
373            new_r = gcd - quotient * temp;
374            gcd = temp;
375
376            temp = new_s;
377            if bezout_1 < quotient.product_residue(temp, other) {
378                // First bezout coefficient is computed over Z[q]
379                new_s = other - (quotient.product_residue(temp, other) - bezout_1)
380            } else {
381                new_s = bezout_1.wrapping_sub(quotient * temp);
382            }
383
384            bezout_1 = temp;
385
386            temp = new_t;
387            if bezout_2 < quotient.product_residue(temp, *self) {
388                // Second bezout coefficient is computed over Z[p]
389                new_t = *self - (quotient.product_residue(temp, *self) - bezout_2)
390            } else {
391                new_t = bezout_2.wrapping_sub(quotient.product_residue(temp, *self));
392            }
393            bezout_2 = temp
394        }
395        (gcd, bezout_1, bezout_2)
396    }
397/*
398    fn coprime(&self, other: Self) -> bool{
399       self.gcd(other)==1
400    }
401*/
402    fn lcm(&self, other: Self) -> NTResult<Self> {
403        if self == &0 && other == 0{
404           return NTResult::Eval(0)
405         }
406        let cf = self.gcd(other);
407        let (v, flag) = (*self / cf).overflowing_mul(other);
408        if flag {
409            return NTResult::Overflow;
410        }
411        NTResult::Eval(v)
412    }
413
414    fn euler_totient(&self) -> Self {
415      if self.reducible(){
416        return (*self as u32).euler_totient() as u64
417      }
418        let factors = self.factor().unwrap();
419        let numerator = factors.factor_iter().map(|x| x - 1u64).product::<u64>();
420        let denominator = factors.factor_iter().product::<u64>();
421        (self / denominator) * numerator
422    }
423
424    fn jordan_totient(&self, k: Self) -> NTResult<Self> {
425        if k > u32::MAX as u64 {
426            return NTResult::Overflow;
427        }
428        
429        if *self < 2{
430           return NTResult::Eval(*self)
431        }
432        
433        let (coef, flag) = self.overflowing_pow(k as u32);
434
435        if flag {
436            return NTResult::CompOverflow;
437        }
438
439        let mut denom = 1u64;
440        let mut numer = 1u64;
441
442        for i in self.factor().unwrap().factor_iter(){
443            let pow = i.pow(k as u32);
444
445            denom *= pow;
446
447            numer *= pow - 1;
448        }
449
450        NTResult::Eval(numer * (coef / denom))
451    }
452    
453    fn exponent(&self) -> NTResult<Self>{ 
454    
455       if self.reducible(){
456         return (*self as u32).exponent().map(|x| x as u64)
457       }
458    
459       if *self == 0{
460          return NTResult::Infinite;
461       }
462       let fctr = self.factor().unwrap();
463       let mut result = 1;
464       for i in 0..fctr.base.len(){
465       if fctr.base[0] == 2 && fctr.power[0] > 2{
466          let phi = 2u64<<(fctr.power[0]-1);
467          result=phi;
468       }
469       else{
470         let el = fctr.base[i];
471         let phi =  (el.pow(fctr.power[i] as u32)/el)*(el-1);
472         match result.lcm(phi){
473             NTResult::Overflow => {return NTResult::Overflow;},
474             NTResult::Eval(x) => result=x,
475             _=> (),
476         }
477       }
478       }
479        return NTResult::Eval(result);
480    }
481    
482    fn dedekind_psi(&self, k: Self) -> NTResult<Self> {
483        if *self == 0{
484         return NTResult::Infinite
485       }
486       
487        let (k2, flag) = k.overflowing_shl(1);
488        if flag {
489            return NTResult::Overflow;
490        }
491        self.jordan_totient(k2).map(|y| y/self.jordan_totient(k).unwrap())
492    }
493
494   fn quadratic_residue(&self, n: Self) -> Self {
495        if n == 0 {
496            return self.wrapping_mul(*self)
497        }
498        ((*self as u128 * *self as u128) % n as u128) as Self
499    }
500    
501    fn checked_quadratic_residue(&self, n: Self) -> NTResult<Self> {
502        if n == 0 {
503            return NTResult::from_option(self.checked_mul(*self),NTResult::Overflow)
504        }
505        NTResult::Eval(((*self as u128 * *self as u128) % n as u128) as Self)
506    }
507    
508
509    fn product_residue(&self, other: Self, n: Self) -> Self {
510        if n == 0 {
511            return self.wrapping_mul(other)
512        }
513        ((*self as u128 * other as u128) % n as u128) as Self
514    }
515    
516    fn checked_product_residue(&self, other: Self, n: Self) -> NTResult<Self> {
517        if n == 0 {
518            return NTResult::from_option(self.checked_mul(*self),NTResult::Overflow)
519        }
520        NTResult::Eval(((*self as u128 * other as u128) % n as u128) as Self)
521    }
522
523    fn exp_residue(&self, p: Self, modulus: Self) -> Self {
524        
525        if modulus == 0 {
526            return self.pow(p as u32);
527        }
528        
529        if modulus & 1 == 0 {
530        
531        let k = modulus.trailing_zeros() as u64;
532        let s = modulus >> k;
533
534        let reducer = (1 << k) - 1; // A shorthand for arithmetic over Z[2k]
535
536        let k_rem = self.even_pow(p,reducer);
537
538        let s_rem = self.odd_pow(p,s);
539        
540        let s_inv = s.inv_2()&reducer;
541    
542        let y = k_rem.wrapping_sub(s_rem).wrapping_mul(s_inv) & reducer;
543
544        s_rem + s * y //)%modulus
545    } else {
546        self.odd_pow(p,modulus) //%modulus
547    }
548    }
549
550    fn checked_exp_residue(&self, p: Self, modulus: Self) -> NTResult<Self> {
551        if modulus == 0 {
552            if  p > u32::MAX as u64 {
553                return NTResult::Overflow;
554            }
555            match self.checked_pow(p as u32) {
556                Some(x) => return NTResult::Eval(x),
557                None => return NTResult::Overflow,
558            };
559        }
560
561        NTResult::Eval(self.exp_residue(p,modulus))
562    }
563
564    fn legendre(&self, p: Self) -> i8 {
565        let k = self.exp_residue((p - 1) >> 1, p);
566        if k == 1 {
567            return 1;
568        };
569        if k == p - 1 {
570            return -1;
571        };
572        0i8
573    }
574
575    fn checked_legendre(&self, p: Self) -> NTResult<i8> {
576        if p == 2 || !p.is_prime() {
577            return NTResult::Undefined;
578        }
579        NTResult::Eval(self.legendre(p))
580    }
581
582    fn liouville(&self) -> i8 {
583      if self.reducible(){
584        return (*self as u32).liouville()
585      }
586        let primeomega = self.factor().unwrap().prime_omega();
587        if primeomega & 1 == 0 {
588            return 1;
589        }
590        -1
591    }
592    
593    fn derivative(&self) -> NTResult<Self> {
594       if *self < 94 {
595         return (*self as u32).derivative().map(|y| y as Self)
596       }
597       
598       let fctr = self.factor().unwrap();
599       let mut sum : u64 = 0;
600       
601       for (power,base) in fctr.power_iter().zip(fctr.factor_iter()){
602          match sum.checked_add((*power as u64)*(*self/base)){
603            Some(x) => sum =x,
604            None => return NTResult::Overflow,
605          }
606       }
607       NTResult::Eval(sum)
608       }
609
610    fn mangoldt(&self) -> f64 {
611     if self.reducible(){
612        return (*self as u32).mangoldt()
613     }
614      let base = self.max_exp().0;
615       if base.is_prime(){
616         return (base as f64).ln()
617       }
618       0f64
619    }
620    
621    fn mobius(&self) -> i8 {
622     //if self.reducible(){
623     //   return (*self as u32).mobius()
624     // }
625      let fctr = self.factor().unwrap();
626     // if fctr.len() == 1{ // if only one factor then return -1
627     //    return -1
628     // }
629      if !fctr.k_free(2){
630         return 0;
631      }
632      let fctrsum = fctr.prime_omega();//fctr[1..].iter().step_by(2).sum::<Self>();
633      if fctrsum&1 == 1{// if odd number of factors and square free
634        return -1
635      }
636      1
637    }
638
639    fn jacobi(&self, k: Self) -> NTResult<i8> {
640        if k == 0 && k % 2 == 0 {
641           return NTResult::Undefined;
642        }
643        let mut n = *self;
644        let mut p = k;
645        let mut t = 1i8;
646        n %= p;
647
648        while n != 0 {
649            let zeros = n.trailing_zeros();
650            n >>= zeros;
651
652            if (p % 8 == 3 || p % 8 == 5) && (zeros % 2 == 1) {
653                t = -t
654            }
655
656            std::mem::swap(&mut n, &mut p);
657            if n % 4 == 3 && p % 4 == 3 {
658                t = -t;
659            }
660            n %= p;
661        }
662
663        if p == 1 {
664            NTResult::Eval(t)
665        } else {
666            NTResult::Eval(0)
667        }
668    }
669
670     fn kronecker(&self, k: Self) -> i8{
671     let x = *self;
672     if k == 0{
673      if x == 1{
674         return 1
675      }
676     return 0
677    }
678   if k == 1{
679      return 1
680   }
681   let fctr = k.factor().unwrap();
682   let mut start = 0;
683   let mut res = 1;
684   
685   if fctr.base[0] ==  2{
686     start = 1;
687     if x&1 == 0{
688     res = 0;
689     }
690     else if x % 8 == 1 || x % 8 == 7{
691      res=1
692     }
693     else{
694       res = (-1i8).pow(fctr.power[0] as u32)
695     }
696   }
697   if fctr.base[0] == 2 && fctr.base.len() == 1{
698     return res
699   }
700   for i in start..fctr.base.len(){
701     res*=self.legendre(fctr.base[i]).pow(fctr.power[i] as u32);
702   }
703   res
704}
705
706
707
708  fn smooth(&self) -> NTResult<Self> {
709       if *self == 0{
710         return NTResult::Infinite
711       }
712       if *self == 1{
713        return NTResult::DNE
714       }
715       self.factor().map(|x| x.max())
716  
717    }
718    
719    
720
721    fn is_smooth(&self, b: Self) -> bool {
722     match self.smooth(){
723      NTResult::Infinite => false,
724      NTResult::Eval(x) => x <= b, 
725      _=> false,
726     }
727   }
728   
729    fn ord(&self, n: Self) -> NTResult<Self>{
730      let ord_2 = |a: u64, p: u64| -> u64{
731        let modulo = (1u64<<p)-1;
732        let mut b = a&modulo;
733   
734         if b == 1{
735            return 1;
736         }
737      for i in 1..p{
738         b = b.wrapping_mul(b)&modulo;
739         if b == 1{
740           return 1<<i;
741         }
742    }
743    return p;
744    };
745
746// Given ord(a,p)  calculate ord(a,p^n)
747    let  pp_ord = |a: u64, b: u64, p: u64, e: u32| -> u64{
748     for i in 0..e+1{
749          if a.exp_residue(b*p.pow(i),p.pow(e)) ==1{
750             return b*p.pow(i);
751           }
752   }
753    return b*p.pow(e);
754    };
755
756    let p_ord = |a: u64, p: u64| -> u64{
757   
758   let fctr = (p-1).factor().unwrap();
759   
760   let mut m = p-1;
761   for i in fctr.pair_iter(){
762     for _ in 0..*i.1{
763          if a.exp_residue(m/ *i.0,p) == 1{
764            m = m/ *i.0;
765          }
766          else{
767            break;
768          }
769     }
770  }
771 m
772};
773
774
775    if self.gcd(n) != 1{
776       return NTResult::DNE;
777    }
778    let fctr = n.factor().unwrap();
779    let mut fullord = 1u64;
780    for i in fctr.pair_iter(){
781     let mut ord : Self;
782      if *i.0 == 2{
783         ord = ord_2(*self,*i.1);
784      }
785      else{
786        ord = p_ord(*self,*i.0);
787        if *i.1 > 1{
788           ord=pp_ord(*self,ord,*i.0,*i.1 as u32);
789        }
790      }
791       fullord = fullord.lcm(ord).unwrap(); 
792    }
793    NTResult::Eval(fullord)
794 }
795
796}
number_theory/primitive/eightbytes.rs

number_theory/primitive/
eightbytes.rs
