From cc1116c2cacd2cf0511d9f2f08de81e448a58568 Mon Sep 17 00:00:00 2001
From: PlexSheep <PlexSheep@protonmail.com>
Date: Sat, 13 May 2023 16:33:49 +0200
Subject: [PATCH] p minus one working

---
 Cargo.toml      |   2 +
 src/main.rs     |  32 +++++++++++++++
 src/math/mod.rs |   1 +
 src/math/pm1.rs | 105 ++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 140 insertions(+)
 create mode 100644 src/math/pm1.rs

diff --git a/Cargo.toml b/Cargo.toml
index c9789ba..2429e32 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,6 +16,8 @@ path = "src/main.rs"
 [dependencies]
 clap = { version = "4.2.7", features = ["derive"]}
 clap-num = "1.0.2"
+num = "0.4.0"
 num-bigint = "0.4.3"
 num-traits = "0.2.15"
+primes = "0.3.0"
 pyo3 = "0.18.1"
diff --git a/src/main.rs b/src/main.rs
index a5124c5..9d5feab 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -78,6 +78,7 @@ struct AlgoCommand {
 enum MathActions {
     #[command(name="modexp")]
     Modexp(ModexpArgs),
+    Pm1(PM1Args),
 }
 
 #[derive(Args, Clone, Debug, PartialEq, Eq)]
@@ -87,6 +88,12 @@ struct ModexpArgs {
     field: String
 }
 
+#[derive(Args, Clone, Debug, PartialEq, Eq)]
+struct PM1Args {
+    n: u128,
+    max_prime: u128,
+}
+
 #[derive(Subcommand, Clone, Debug, PartialEq, Eq)]
 enum BinaryActions {
     /// bit rotation/circular shifting (only 32bit)
@@ -165,6 +172,31 @@ pub fn main() {
                         println!("result is {}", result)
                     }
                 }
+                MathActions::Pm1(pm1_args) => {
+                    let result: Result<Vec<u128>, String> = math::pm1::p_minus_one(
+                        pm1_args.n, 
+                        pm1_args.max_prime,
+                        args.verbose
+                        );
+                    match result {
+                        Ok(vec) => {
+                            if args.machine {
+                                println!("{:?}", vec)
+                            }
+                            else {
+                                println!("result is {:?}", vec)
+                            }
+                        }
+                        Err(e) => {
+                            if args.machine {
+                                println!("{:?}", e)
+                            }
+                            else {
+                                println!("could not compute: {:?}", e)
+                            }
+                        }
+                    }
+                }
             }
         }
         Commands::Binary(action) => {
diff --git a/src/math/mod.rs b/src/math/mod.rs
index 0dfe59f..65c4dca 100644
--- a/src/math/mod.rs
+++ b/src/math/mod.rs
@@ -7,3 +7,4 @@
 /// License:    MIT
 /// Source:     <https://git.cscherr.de/PlexSheep/plexcryptool/>
 pub mod modexp;
+pub mod pm1;
diff --git a/src/math/pm1.rs b/src/math/pm1.rs
new file mode 100644
index 0000000..5c3b565
--- /dev/null
+++ b/src/math/pm1.rs
@@ -0,0 +1,105 @@
+#![allow(dead_code)]
+/// P minus 1 method
+///
+/// Determine the prime factors of a number with the p minus 1 method.
+/// Effecient for numbers with low ranged prime factors.
+///
+/// Author:     Christoph J. Scherr <software@cscherr.de>
+/// License:    MIT
+/// Source:     <https://git.cscherr.de/PlexSheep/plexcryptool/>
+
+use pyo3::prelude::*;
+
+use num::integer::gcd;
+
+use num_bigint::BigInt;
+
+use primes::{Sieve, PrimeSet};
+
+use crate::math::modexp;
+
+/// excecute the p minus one calculation
+pub fn p_minus_one(n: u128, max_prime: u128, verbose: bool) -> Result<Vec<u128>, String> {
+    assert!(n > 2);
+    let m1: u128 = n -1;
+    let mut k_parts: Vec<(u128, u32)> = Vec::new();
+    let mut prime_parts: Vec<u128> = Vec::new();
+    //
+    // get a list of the early primes
+    let mut pset = Sieve::new();
+    if verbose {
+        println!("getting list of first {max_prime} primes");
+}
+    for (_i_prime, prime) in pset.iter().enumerate().take(max_prime as usize) {
+        let num: u128 = prime as u128;
+        if num > max_prime {
+            break;
+        }
+        let mut exp: u32 = 1;
+        if verbose {
+            println!("current prime: {num}");
+        }
+        loop {
+            if num.pow(exp + 1) < max_prime {
+                exp += 1;
+            }
+            else {
+                break;
+            }
+        }
+        if verbose {
+            println!("exponented prime: {}", num.pow(exp));
+        }
+        k_parts.push((num, exp));
+    }
+    let mut k = 1u128;
+    for (num, exp) in k_parts.clone() {
+        k = num.pow(exp) * k;
+        if verbose {
+            println!("k at step: {k}");
+        }
+    }
+    if verbose {
+        println!("k: {k}\nk parts: {:?}", k_parts);
+    }
+
+    let a = 2u128;
+    let akn1: u128 = ((modexp::modular_exponentiation(
+            BigInt::from(a), 
+            BigInt::from(k), 
+            BigInt::from(n), 
+            false)
+        ) - BigInt::from(1)).try_into().expect("Number too big");
+    if verbose {
+        println!("a: {a}\na**k-1 {akn1}");
+    }
+    let mut next_gcd = gcd(akn1, n);
+
+    if next_gcd == 1 {
+        return Err(format!("P minus one does not offer divisor for {n} with max_prime: {max_prime}"));
+    }
+    let mut q: u128;
+    while next_gcd > 1 {
+        prime_parts.push(next_gcd);
+        q = n / next_gcd;
+        next_gcd = gcd(q, n);
+        if verbose {
+            println!("nextgcd: {next_gcd}|q: {q}");
+        }
+        if prime_parts.contains(&next_gcd) {
+            break;
+        }
+    }
+    return Ok(prime_parts);
+}
+
+/// alternative simple implementation for gcd
+pub fn alt_gcd(mut a: u128, mut b: u128) -> u128 {
+    let mut tmp: u128;
+    while b > 0 {
+        tmp = b;
+        b = a % b;
+        a = tmp;
+    }
+    return a;
+}