1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
/**
* (C) Copyright Projet SECRET, INRIA, Rocquencourt
* (C) Bhaskar Biswas and Nicolas Sendrier
*
* Distributed under the terms of the Botan license
*
*/
#include <botan/mceliece.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)
{
// x <- binomial(k/2,p)
double x = binomial(k / 2, p);
// i <- log[2](binomial(k/2,p))
size_t i = (size_t) (std::log(x) / std::log(2)); // normalement i < 2^31
// res <- 2*p*(n-k-l)*binomial(k/2,p)^2/2^l
double res = 2 * p * (n - k - l) * ldexp(x * x, -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);
}
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;
// On part de l = u, en faisant croitre l.
// On s'arrète dés que le work factor croit.
// Puis on explore les valeurs <u, mais en tenant de la convexite'
double min = cout_total(n, k, w, p, 0);
for(size_t l = 1; l < n - k; ++l)
{
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 k, size_t 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 min;
}
}
|
