ECM with trial division and primality test. Floating point exception …

.vscode/settings.json

@@ -59,6 +59,20 @@
     "typeinfo": "cpp",
     "unordered_map": "cpp",
     "utility": "cpp",
-    "vector": "cpp"
+    "vector": "cpp",
+    "clocale": "cpp",
+    "compare": "cpp",
+    "condition_variable": "cpp",
+    "deque": "cpp",
+    "numeric": "cpp",
+    "queue": "cpp",
+    "random": "cpp",
+    "thread": "cpp",
+    "*.tcc": "cpp",
+    "concepts": "cpp",
+    "memory_resource": "cpp",
+    "ranges": "cpp",
+    "stop_token": "cpp",
+    "cinttypes": "cpp"
   }
 }

FibPseudoprime.code-workspace

@@ -0,0 +1,71 @@
+{
+  "folders": [
+    {
+      "path": "."
+    }
+  ],
+  "settings": {
+    "files.associations": {
+      "__bit_reference": "cpp",
+      "__config": "cpp",
+      "__debug": "cpp",
+      "__errc": "cpp",
+      "__functional_base": "cpp",
+      "__hash_table": "cpp",
+      "__locale": "cpp",
+      "__mutex_base": "cpp",
+      "__node_handle": "cpp",
+      "__nullptr": "cpp",
+      "__split_buffer": "cpp",
+      "__string": "cpp",
+      "__threading_support": "cpp",
+      "__tuple": "cpp",
+      "algorithm": "cpp",
+      "array": "cpp",
+      "atomic": "cpp",
+      "bit": "cpp",
+      "bitset": "cpp",
+      "cctype": "cpp",
+      "chrono": "cpp",
+      "cmath": "cpp",
+      "complex": "cpp",
+      "cstdarg": "cpp",
+      "cstddef": "cpp",
+      "cstdint": "cpp",
+      "cstdio": "cpp",
+      "cstdlib": "cpp",
+      "cstring": "cpp",
+      "ctime": "cpp",
+      "cwchar": "cpp",
+      "cwctype": "cpp",
+      "exception": "cpp",
+      "functional": "cpp",
+      "initializer_list": "cpp",
+      "ios": "cpp",
+      "iosfwd": "cpp",
+      "iostream": "cpp",
+      "istream": "cpp",
+      "iterator": "cpp",
+      "limits": "cpp",
+      "locale": "cpp",
+      "memory": "cpp",
+      "mutex": "cpp",
+      "new": "cpp",
+      "optional": "cpp",
+      "ostream": "cpp",
+      "ratio": "cpp",
+      "sstream": "cpp",
+      "stdexcept": "cpp",
+      "streambuf": "cpp",
+      "string": "cpp",
+      "string_view": "cpp",
+      "system_error": "cpp",
+      "tuple": "cpp",
+      "type_traits": "cpp",
+      "typeinfo": "cpp",
+      "unordered_map": "cpp",
+      "utility": "cpp",
+      "vector": "cpp"
+    }
+  }
+}

FibPseudoprime/main.cpp

@@ -1,18 +1,18 @@
 #include "nfunctions.hpp"
 #include <iostream>
 
-int main() {
-  long n = 5563 * 2 * 3 * 5 * 7;
-
-  while (true) {
-    long factor = ECFactor(n, 3000, 200);
+void printFactors(std::vector<std::pair<long, int>> factors) {
+  std::cout << "[";
+  for (std::pair<long, int> f : factors) {
+    std::cout << "(" << f.first << ", " << f.second << "), ";
+  }
+  std::cout << "]";
+}
 
-    std::cout << "n: " << n << ", factor: " << factor << std::endl;
+int main() {
+  long n = 5563 * 2 * 3 * 5 * 49 * 686989L;
 
-    if (factor == 1)
-      break;
+  printFactors(factorize(n));
 
-    n = n / factor;
-  }
   return 0;
 }

FibPseudoprime/nfunctions/nfunctions.cpp

@@ -1,4 +1,7 @@
 #include "nfunctions.hpp"
+#include "ptests.hpp"
+#include "smallprimes.hpp"
+#include <cmath>
 
 long gcd(long a, long b) {
   if (b == 0)
@@ -248,6 +251,72 @@ long ECFactor(long n, int r, int tries) {
   return 1;
 }
 
+// returns factorization using primes < 65521
+std::vector<std::pair<long, int>> smallPrimeFactors(long &n) {
+  std::vector<long>::iterator primePtr = smallPrimes.begin();
+  std::vector<std::pair<long, int>> res;
+
+  while (n > 1 || primePtr != smallPrimes.end()) {
+    if (n % (*primePtr) == 0) {
+      std::pair<long, int> factor(*primePtr, 1);
+      n /= *primePtr;
+      while (n % (*primePtr) == 0) {
+        factor.second++;
+        n /= *primePtr;
+      }
+
+      res.push_back(factor);
+    }
+    primePtr++;
+  }
+
+  return res;
+}
+
+std::pair<long, int> highPowers(long &n, long p) {
+  std::pair<long, int> pfactor(p, 1);
+
+  n /= p;
+  while (n % p == 0) {
+    n /= p;
+    pfactor.second++;
+  }
+
+  return pfactor;
+}
+
+// returns array of pairs of (factor, power)
+std::vector<std::pair<long, int>> factorize(long n) {
+  if (n == 0)
+    return std::vector<std::pair<long, int>>();
+  n = std::abs(n);
+
+  // remove small prime factors from n
+  std::vector<std::pair<long, int>> factors = smallPrimeFactors(n);
+
+  if (n == 1)
+    return factors;
+
+  while (n > 1) {
+    long probablePrimeFactor = ECFactor(n, 3000, 200);
+
+    if (probablePrimeFactor == 1)
+      break;
+
+    if (!millerRabin(probablePrimeFactor, 100)) {
+      std::vector<std::pair<long, int>> primeFactors =
+          factorize(probablePrimeFactor);
+      for (std::pair<long, int> f : primeFactors) {
+        factors.push_back(highPowers(f.first, n));
+      }
+    } else {
+      factors.push_back(highPowers(n, probablePrimeFactor));
+    }
+  }
+
+  return factors;
+}
+
 void FibProduct(std::vector<size_t> one, std::vector<size_t> two) {
   size_t x = 0;
   size_t y = 1;

FibPseudoprime/nfunctions/nfunctions.hpp

@@ -29,15 +29,17 @@ std::vector<bool> baseTwo(int n);
 
 std::vector<Mod> chooseEAndP(long n);
 
-long ECAdd(Mod a, Mod b, Mod &x1, Mod &y1, Mod x2, Mod y2);
+// long ECAdd(Mod a, Mod b, Mod &x1, Mod &y1, Mod x2, Mod y2);
 
-long ECDouble(Mod a, Mod b, Mod &x1, Mod &y1);
+// long ECDouble(Mod a, Mod b, Mod &x1, Mod &y1);
 
-long ECMultiply(Mod a, Mod b, Mod &x, Mod &y, int t);
+// long ECMultiply(Mod a, Mod b, Mod &x, Mod &y, int t);
 
-long findRFactorialPOnE(int r, Mod a, Mod b, Mod x, Mod y);
+// long findRFactorialPOnE(int r, Mod a, Mod b, Mod x, Mod y);
 
-long ECFactor(long n, int r, int tries);
+// long ECFactor(long n, int r, int tries);
+
+std::vector<std::pair<long, int>> factorize(long n);
 
 // forgot to make you do this one in the code review
 void FibProduct(std::vector<size_t> one, std::vector<size_t> two);