Add Jacobi check for LL

Args.cpp

@@ -52,6 +52,7 @@ void Args::printHelp() {
 -fft <spec>        : specify FFT e.g.: 1152K, 5M, 5.5M, 256:10:1K
 -block <value>     : PRP GEC block size, or LL iteration-block size. Must divide 10'000.
 -log <step>        : log every <step> iterations. Multiple of 10'000.
+-jacobi <step>     : (LL-only): do Jacobi check every <step> iterations. Default 1'000'000.
 -carry long|short  : force carry type. Short carry may be faster, but requires high bits/word.
 -B1                : P-1 B1 bound, default %u
 -B2                : P-1 B2 bound, default B1 * 30
@@ -157,6 +158,11 @@ void Args::parse(string line) {
         log("BlockSize %u must divide 10'000\n", blockSize);
         throw "invalid block size";
       }
+    } else if (key == "-jacobi") {
+      jacobiStep = stoi(s);
+      if (jacobiStep <= 0 || jacobiStep % 10000) {
+        throw "invalid jacobi step";
+      }
     } else if (key == "-use") {
       string ss = s;
       std::replace(ss.begin(), ss.end(), ',', ' ');
@@ -174,6 +180,11 @@ void Args::parse(string line) {
     log("log step (%u) must be a multiple of 10'000\n", logStep);
     throw "invalid log step";
   }
+
+  if (jacobiStep < 100'000 && !logStep) { logStep = jacobiStep; }
+  if (jacobiStep && ((logStep && jacobiStep % logStep) || (!logStep && jacobiStep % 100'000))) {
+    log("jacobi step (%u) must be a multiple of log step (%u)\n", jacobiStep, logStep);
+  }
 
   if (!masterDir.empty()) {
     if (resultsFile.find_first_of('/') == std::string::npos) {

Args.h

@@ -40,6 +40,7 @@ class Args {
   int carry = CARRY_AUTO;
   u32 blockSize = 0;
   u32 logStep   = 0;
+  u32 jacobiStep = 0;
   string fftSpec;
 
   u32 B1 = 1000000;

GmpUtil.cpp

@@ -52,5 +52,11 @@ std::string GCD(u32 exp, const std::vector<u32>& words, u32 sub) {
   return (resultGcd == 1) ? ""s : resultGcd.get_str();
 }
 
-// "Rev" means: most significant bit first (at index 0).
+// MSB: Most Significant Bit first (at index 0).
 vector<bool> powerSmoothMSB(u32 exp, u32 B1) { return bitsMSB(powerSmooth(exp, B1)); }
+
+int jacobi(u32 exp, const std::vector<u32>& words) {
+  mpz_class w = mpz(words) - 2;
+  mpz_class m = (mpz_class{1} << exp) - 1;
+  return mpz_jacobi(w.get_mpz_t(), m.get_mpz_t());
+}

GmpUtil.h

@@ -15,3 +15,6 @@ std::string GCD(u32 exp, const std::vector<u32>& words, u32 sub = 0);
 vector<bool> bitsMSB(mpz_class a);
 
 vector<bool> powerSmoothMSB(u32 exp, u32 B1);
+
+// Returns jacobi-symbol(words - 2, 2**exp - 1)
+int jacobi(u32 exp, const std::vector<u32>& words);

Gpu.cpp

@@ -903,10 +903,18 @@ void Gpu::printRoundoff(u32 E) {
 #endif
 }
 
+namespace {
+  struct JacobiData {
+    vector<u32> words;
+    u32 k;
+  };
+}
+
 tuple<bool, u64> Gpu::isPrimeLL(u32 E, const Args &args) {
   Buffer<double> buf1{queue, "buf1", N};
   Buffer<double> buf2{queue, "buf2", N};
 
+ reload:
   LLState loaded(E);
   writeData(loaded.data);
   u64 res64 = dataResidue();
@@ -925,37 +933,78 @@ tuple<bool, u64> Gpu::isPrimeLL(u32 E, const Args &args) {
   bool displayRoundoff = args.flags.count("STATS");
 
   u32 blockSize = args.blockSize ? args.blockSize : 1000; // default LL block size
-  u32 logStep = args.logStep ? args.logStep : 100000;
+  u32 logStep = args.logStep ? args.logStep : 100'000;
+  u32 jacobiStep = args.jacobiStep ? args.jacobiStep : 1'000'000;
+
+  future<int> jacobiFuture;
+  JacobiData jacobiData;
 
   while (true) {
     assert(k < kEnd);    
     u32 nextK = std::min(kEnd, k + blockSize);
     modSqLoopLL(nextK - k, buf1, buf2, bufData);
     k = nextK;
 
-    if (!args.noSpin) { spin(); }
-
     queue->finish();
 
+    if (!args.noSpin) { spin(); }
+
     bool doStop = (args.iters && k - startK == args.iters);
     if (signal.stopRequested()) {
       doStop = true;
       log("Stopping, please wait..\n");
       signal.release();
     }
 
-    if (doStop || (k % logStep == 0) || (k >= kEnd)) {
+    bool atEnd = k >= kEnd;
+
+    if (doStop || (k % logStep == 0) || atEnd) {
       if (displayRoundoff) { printRoundoff(E); }
+
       vector<u32> data = readData();
       u64 res64 = residue(data);
       double time = itTimer.reset(k);
       logLL(E, k, res64, time, kEnd);
-      if (k >= kEnd) {
+
+      if (!jacobiFuture.valid() && doStop) {
+        // log("starting Jacobi check\n");
+        jacobiData.words = data; // copy the vector
+        jacobiData.k = k;
+        jacobiFuture = async(launch::async, jacobi, E, std::cref(jacobiData.words));
+      }
+
+      if (jacobiFuture.valid()) {
+        if (atEnd || doStop) {
+          log("waiting for the Jacobi check to finish..\n");
+          jacobiFuture.wait();
+        }
+
+        if (jacobiFuture.wait_for(chrono::steady_clock::duration::zero()) == future_status::ready) {
+          int jacobi = jacobiFuture.get();
+          bool jacobiOK = (jacobi == -1);
+          log("%u %2s %8d (jacobi == %d)\n", E, jacobiOK ? "OK" : "EE", jacobiData.k, jacobi);
+
+          if (jacobiOK) {          
+            LLState{E, jacobiData.k, std::move(jacobiData.words)}.save();
+          }
+
+          if (doStop) { throw "stop requested"; }
+
+          if (!jacobiOK) { goto reload; }
+        }
+      }
+
+      if (!jacobiFuture.valid() && (k % jacobiStep == 0)) {
+        // log("starting Jacobi check\n");
+        jacobiData.words = data; // copy the vector
+        jacobiData.k = k;
+        jacobiFuture = async(launch::async, jacobi, E, std::cref(jacobiData.words));
+      }
+
+      if (atEnd) {
         assert(k == kEnd);
         bool allZero = std::all_of(data.begin(), data.end(), [](u32 x) { return (x == 0); } );
         return {allZero, res64};
-      } else {
-        LLState{E, k, std::move(data)}.save();
       }
       if (doStop) { throw "stop requested"; }
     }