1 files changed, 112 insertions, 0 deletions
diff --git a/comm/third_party/botan/src/lib/pubkey/mce/mce_workfactor.cpp b/comm/third_party/botan/src/lib/pubkey/mce/mce_workfactor.cpp
new file mode 100644
index 0000000000..ce38781c8b
--- /dev/null
+++ b/comm/third_party/botan/src/lib/pubkey/mce/mce_workfactor.cpp
@@ -0,0 +1,112 @@
+/*
+ * (C) Copyright Projet SECRET, INRIA, Rocquencourt
+ * (C) Bhaskar Biswas and Nicolas Sendrier
+ * (C) 2014 Jack Lloyd
+ *
+ * Botan is released under the Simplified BSD License (see license.txt)
+ *
+ */
+
+#include <botan/mceliece.h>
+#include <botan/internal/bit_ops.h>
+#include <cmath>
+
+namespace Botan {
+
+namespace {
+
+double binomial(size_t n, size_t k)
+   {
+   double x = 1;
+
+   for(size_t i = 0; i != k; ++i)
+      {
+      x *= n - i;
+      x /= k -i;
+      }
+
+   return x;
+   }
+
+double log_binomial(size_t n, size_t k)
+   {
+   double x = 0;
+
+   for(size_t i = 0; i != k; ++i)
+      {
+      x += std::log(n - i);
+      x -= std::log(k - i);
+      }
+
+   return x / std::log(2);
+   }
+
+double nb_iter(size_t n, size_t k, size_t w, size_t p, size_t l)
+   {
+   double x = 2 * log_binomial(k / 2, p);
+   x += log_binomial(n - k - l, w - 2 * p);
+   x = log_binomial(n, w) - x;
+   return x;
+   }
+
+double cout_iter(size_t n, size_t k, size_t p, size_t l)
+   {
+   double x = binomial(k / 2, p);
+   const size_t i = static_cast<size_t>(std::log(x) / std::log(2));
+   double res = 2 * p * (n - k - l) * std::ldexp(x * x, -static_cast<int>(l));
+
+   // x <- binomial(k/2,p)*2*(2*l+log[2](binomial(k/2,p)))
+   x *= 2 * (2 * l + i);
+
+   // res <- k*(n-k)/2 +
+   // binomial(k/2,p)*2*(2*l+log[2](binomial(k/2,p))) +
+   // 2*p*(n-k-l)*binomial(k/2,p)^2/2^l
+   res += x + k * ((n - k) / 2.0);
+
+   return std::log(res) / std::log(2); // convert to bits
+   }
+
+double cout_total(size_t n, size_t k, size_t w, size_t p, size_t l)
+   {
+   return nb_iter(n, k, w, p, l) + cout_iter(n, k, p, l);
+   }
+
+double best_wf(size_t n, size_t k, size_t w, size_t p)
+   {
+   if(p >= k / 2)
+      return -1;
+
+   double min = cout_total(n, k, w, p, 0);
+
+   for(size_t l = 1; l < n - k; ++l)
+      {
+      const double lwf = cout_total(n, k, w, p, l);
+      if(lwf < min)
+         min = lwf;
+      else
+         break;
+      }
+
+   return min;
+   }
+
+}
+
+size_t mceliece_work_factor(size_t n, size_t t)
+   {
+   const size_t k = n - ceil_log2(n) * t;
+
+   double min = cout_total(n, k, t, 0, 0); // correspond a p=1
+   for(size_t p = 0; p != t / 2; ++p)
+      {
+      double lwf = best_wf(n, k + 1, t, p);
+      if(lwf < 0)
+         break;
+
+      min = std::min(min, lwf);
+      }
+
+   return static_cast<size_t>(min);
+   }
+
+}