Proof v2

Args.cpp

@@ -62,8 +62,7 @@ void Args::printHelp() {
 -cleanup           : delete save files at end of run
 -prp <exponent>    : run a single PRP test and exit, ignoring worktodo.txt
 -verify <file>     : verify PRP-proof contained in <file>
--proof <power>     : Valid <power> values are 6 to 9.
-                     By default a proof of power 8 is generated, using 3GB of temporary disk space for a 100M exponent.
+-proof <power>     : By default a proof of power 8 is generated, using 3GB of temporary disk space for a 100M exponent.
                      A lower power reduces disk space requirements but increases the verification cost.
                      A proof of power 9 uses 6GB of disk space for a 100M exponent and enables faster verification.
 -tmpDir <dir>      : specify a folder with plenty of disk space where temporary proof checkpoints will be stored.
@@ -130,12 +129,12 @@ void Args::parse(string line) {
     if (key == "-h" || key == "--help") { printHelp(); throw "help"; }
     else if (key == "-proof") {
       int power = 0;
-      if (s.empty() || (power = stoi(s)) < 6 || power > 10) {
-        log("-proof expects <power> 6 - 9 (found '%s')\n", s.c_str());
-        throw "-proof expects <power> 6 - 9";
+      if (s.empty() || (power = stoi(s)) < 1 || power > 10) {
+        log("-proof expects <power> 1 - 10 (found '%s')\n", s.c_str());
+        throw "-proof <power>";
       }
       proofPow = power;
-      assert(proofPow >= 6 && proofPow <= 10);
+      assert(proofPow > 0);
     } else if (key == "-tmpDir" || key == "-tmpdir") {
       if (s.empty()) {
         log("-tmpDir needs <dir>\n");

Gpu.cpp

@@ -720,6 +720,12 @@ Words Gpu::expMul(const Words& A, u64 h, const Words& B) {
   return readData();
 }
 
+Words Gpu::expMul2(const Words& A, u64 h, const Words& B) {
+  expMul(A, h, B);
+  modMul(bufData, bufData, bufCheck, buf1, buf2, buf3);
+  return readData();
+}
+
 void Gpu::exponentiate(Buffer<int>& bufInOut, u64 exp, Buffer<double>& buf1, Buffer<double>& buf2, Buffer<double>& buf3) {
   if (exp == 0) {
     bufInOut.set(1);
@@ -1489,17 +1495,21 @@ PRPResult Gpu::isPrimePRP(const Args &args, const Task& task) {
   u64 finalRes64 = 0;
 
   // We extract the res64 at kEnd.
-  const u32 kEnd = E; // Type-1 per http://www.mersenneforum.org/showpost.php?p=468378&postcount=209
+  // For M=2^E-1, residue "type-3" == 3^(M+1), and residue "type-1" == type-3 / 9,
+  // See http://www.mersenneforum.org/showpost.php?p=468378&postcount=209
+  // For both type-1 and type-3 we need to do E squarings (as M+1==2^E).
+  const u32 kEnd = E;
   assert(k < kEnd);
 
   // We continue beyound kEnd: up to the next multiple of 1024 if proof is enabled (kProofEnd), and up to the next blockSize
-  u32 kEndEnd = roundUp(proofSet.kProofEnd(kEnd), blockSize);
+  // u32 kEndEnd = roundUp(proofSet.kProofEnd(kEnd), blockSize);
+  u32 kEndEnd = roundUp(kEnd, blockSize);
 
   bool printStats = args.flags.count("STATS");
 
   bool skipNextCheckUpdate = false;
 
-  u32 persistK = proofSet.firstPersistAt(k + 1);
+  u32 persistK = proofSet.next(k);
   bool leadIn = true;
 
   assert(k % blockSize == 0);
@@ -1556,7 +1566,7 @@ PRPResult Gpu::isPrimePRP(const Args &args, const Task& task) {
         goto reload;
       }
       proofSet.save(k, data);
-      persistK = proofSet.firstPersistAt(k + 1);
+      persistK = proofSet.next(k);
     }
 
     if (k == kEnd) {

Gpu.h

@@ -229,6 +229,9 @@ class Gpu {
   // return A^h * B
   Words expMul(const Words& A, u64 h, const Words& B);
 
+  // return A^h * B^2
+  Words expMul2(const Words& A, u64 h, const Words& B);
+
   // A:= A^h * B
   void expMul(Buffer<i32>& A, u64 h, Buffer<i32>& B);
 

Pm1Plan.cpp

@@ -226,7 +226,7 @@ pair<u32, vector<Pm1Plan::BitBlock>> Pm1Plan::makePlan() {
 
   // All primes <= cutValue can be transposed by mulitplying with firstMissingFactor.
   u32 cutValue = lastCovered / firstMissingFactor(D);
-  u32 startValue = max(primeAfter(B1), primeAfter(cutValue));
+  u32 startValue = max(primeAfter(B1), cutValue + 1);
 
   u32 beginBlock = lowerBlock(startValue);
   u32 endBlock = lastBlock + 1;

Proof.cpp

@@ -42,7 +42,7 @@ ProofInfo getInfo(const fs::path& proofFile) {
   File fi = File::openRead(proofFile, true);
   u32 E = 0, power = 0;
   char c = 0;
-  if (fi.scanf(Proof::HEADER, &power, &E, &c) != 3 || c != '\n') {
+  if (fi.scanf(Proof::HEADER_v2, &power, &E, &c) != 3 || c != '\n') {
     log("Proof file '%s' has invalid header\n", proofFile.string().c_str());
     throw "Invalid proof header";
   }
@@ -58,7 +58,7 @@ fs::path Proof::save(const fs::path& proofResultDir) {
   u32 power = middles.size();
   fs::path fileName = proofResultDir / (strE + '-' + to_string(power) + ".proof");
   File fo = File::openWrite(fileName);
-  fo.printf(HEADER, power, E, '\n');
+  fo.printf(HEADER_v2, power, E, '\n');
   fo.write(B.data(), (E-1)/8+1);
   for (const Words& w : middles) { fo.write(w.data(), (E-1)/8+1); }
   return fileName;
@@ -68,7 +68,7 @@ Proof Proof::load(const fs::path& path) {
   File fi = File::openRead(path, true);
   u32 E = 0, power = 0;
   char c = 0;
-  if (fi.scanf(HEADER, &power, &E, &c) != 3 || c != '\n') {
+  if (fi.scanf(HEADER_v2, &power, &E, &c) != 3 || c != '\n') {
     log("Proof file '%s' has invalid header\n", path.string().c_str());
     throw "Invalid proof header";
   }
@@ -80,38 +80,38 @@ Proof Proof::load(const fs::path& path) {
 }
 
 bool Proof::verify(Gpu *gpu) {
+  log("B         %016" PRIx64 "\n", res64(B));
+  for (u32 i = 0; i < middles.size(); ++i) {
+    log("Middle[%u] %016" PRIx64 "\n", i, res64(middles[i]));
+  }
+
   u32 power = middles.size();
   assert(power > 0);
 
-  u32 topK = roundUp(E, (1 << power));
-  assert(topK % (1 << power) == 0);
-  assert(topK > E);
-  u32 step = topK / (1 << power);
-
-  bool isPrime = false;
-  {
-    Words A{makeWords(E, 3)};
-    log("proof: doing %d iterations\n", topK - E + 1);
-    A = gpu->expExp2(A, topK - E + 1);
-    isPrime = (A == B);
-    // log("the proof indicates %s (%016" PRIx64 " vs. %016" PRIx64 " for a PRP)\n",
-    //    isPrime ? "probable prime" : "composite", res64(B), res64(A));
-  }
+  bool isPrime = (B == makeWords(E, 9));
 
   Words A{makeWords(E, 3)};
 
   auto hash = proof::hashWords(E, B);
-
-  for (u32 i = 0; i < power; ++i) {
+
+  u32 span = E;
+  for (u32 i = 0; i < power; ++i, span = (span + 1) / 2) {
     Words& M = middles[i];
     hash = proof::hashWords(E, hash, M);
-    u64 h = hash[0];
+    u64 h = hash[0];    
     A = gpu->expMul(A, h, M);
-    B = gpu->expMul(M, h, B);
+
+    if (span % 2) {
+      B = gpu->expMul2(M, h, B);
+    } else {
+      B = gpu->expMul(M, h, B);
+    }
+
+    log("%u : A %016" PRIx64 ", M %016" PRIx64 ", B %016" PRIx64 ", h %016" PRIx64 "\n", i, res64(A), res64(M), res64(B), h); 
   }
 
-  log("proof verification: doing %d iterations\n", step);
-  A = gpu->expExp2(A, step);
+  log("proof verification: doing %d iterations\n", span);
+  A = gpu->expExp2(A, span);
 
   bool ok = (A == B);
   if (ok) {
@@ -124,10 +124,76 @@ bool Proof::verify(Gpu *gpu) {
 
 // ---- ProofSet ----
 
-Proof ProofSet::computeProof(Gpu *gpu) {
+ProofSet::ProofSet(const fs::path& tmpDir, u32 E, u32 power)
+  : E{E}, power{power}, exponentDir(tmpDir / to_string(E)) {
+
+  assert(E & 1); // E is supposed to be prime
   assert(power > 0);
 
-  Words B = load(topK);
+  fs::create_directory(exponentDir);
+  fs::create_directory(proofPath);
+
+  for (u32 span = (E + 1) / 2; spans.size() < power; span = (span + 1) / 2) { spans.push_back(span); }
+
+  points.push_back(0);
+  for (u32 p = 0; p < power; ++p) {
+    u32 s = spans[p];
+    for (u32 i = 0, end = points.size(); i < end; ++i) {
+      points.push_back(points[i] + s);
+    }
+  }
+
+  assert(points.size() == (1u << power));
+
+  sorted = points;
+  std::sort(sorted.begin(), sorted.end());
+
+  sorted.erase(sorted.begin());
+  assert(sorted.back() < E);
+  sorted.push_back(E);
+  assert(sorted.size() == (1u << power));
+}
+
+bool ProofSet::canDo(const fs::path& tmpDir, u32 E, u32 power, u32 currentK) {
+  assert(power > 0 && power <= 10);
+  return ProofSet{tmpDir, E, power}.isValidTo(currentK);
+}
+
+u32 ProofSet::effectivePower(const fs::path& tmpDir, u32 E, u32 power, u32 currentK) {
+  for (u32 p = power; p > 0; --p) {
+    log("validating proof residues for power %u\n", p);
+    if (canDo(tmpDir, E, p, currentK)) { return p; }      
+  }
+  assert(false);
+}
+
+bool ProofSet::isValidTo(u32 limitK) const {
+  for (u32 k : sorted) {
+    if (k > limitK) { break; }
+    try { load(k); } catch (...) { return false; }
+  }
+  return true;
+}
+
+u32 ProofSet::next(u32 k) const {
+  auto it = upper_bound(sorted.begin(), sorted.end(), k);
+  return (it == sorted.end()) ? u32(-1) : *it;
+}
+
+void ProofSet::save(u32 k, const Words& words) {
+  assert(k > 0 && k <= E);
+  assert(k == *lower_bound(sorted.begin(), sorted.end(), k));
+  cache.save(k, words);
+}
+
+Words ProofSet::load(u32 k) const {
+  assert(k > 0 && k <= E);
+  assert(k == *lower_bound(sorted.begin(), sorted.end(), k));
+  return cache.load(k);
+}
+
+Proof ProofSet::computeProof(Gpu *gpu) {
+  Words B = load(E);
   Words A = makeWords(E, 3);
 
   vector<Words> middles;
@@ -136,16 +202,19 @@ Proof ProofSet::computeProof(Gpu *gpu) {
   auto hash = proof::hashWords(E, B);
 
   vector<Buffer<i32>> bufVect = gpu->makeBufVector(power);
-    
+
   for (u32 p = 0; p < power; ++p) {
     auto bufIt = bufVect.begin();
     assert(p == hashes.size());
-    log("proof: building level %d, hash %016" PRIx64 "\n", (p + 1), hash[0]);
-    u32 s = topK / (1 << (p + 1));
-    for (int i = 0; i < (1 << p); ++i) {
-      Words w = load(s * (2 * i + 1));
+    u32 s = (1u << (power - p));
+    for (u32 i = 0; i < (1u << p); ++i) {
+      u32 idx = s/2 - 1 + i * s;
+      Words w = load(sorted[idx]);
+      // log("loaded %u %u\n", idx, sorted[idx]);
+
       gpu->writeIn(*bufIt++, w);
-      for (int k = 0; i & (1 << k); ++k) {
+      for (u32 k = 0; i & (1u << k); ++k) {
+        assert(k <= p - 1);
         --bufIt;
         u64 h = hashes[p - 1 - k];
         gpu->expMul(*(bufIt - 1), h, *bufIt);
@@ -155,6 +224,8 @@ Proof ProofSet::computeProof(Gpu *gpu) {
     middles.push_back(gpu->readAndCompress(bufVect.front()));
     hash = proof::hashWords(E, hash, middles.back());
     hashes.push_back(hash[0]);
+
+    log("proof level %u : M %016" PRIx64 ", h %016" PRIx64 "\n", p, res64(middles.back()), hashes.back()); 
   }
   return Proof{E, std::move(B), std::move(middles)};
 }