Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/coding/code_bounds.py
4034 views
1
r"""
2
Bounds for Parameters of Codes
3
4
This module provided some upper and lower bounds for the parameters
5
of codes.
6
7
AUTHORS:
8
9
- David Joyner (2006-07): initial implementation.
10
11
- William Stein (2006-07): minor editing of docs and code (fixed bug
12
in elias_bound_asymp)
13
14
- David Joyner (2006-07): fixed dimension_upper_bound to return an
15
integer, added example to elias_bound_asymp.
16
17
- " (2009-05): removed all calls to Guava but left it as an option.
18
19
Let `F` be a finite field (we denote the finite field with
20
`q` elements by `\GF{q}`).
21
A subset `C` of `V=F^n` is called a code of
22
length `n`. A subspace of `V` (with the standard
23
basis) is called a linear code of length `n`. If its
24
dimension is denoted `k` then we typically store a basis of
25
`C` as a `k\times n` matrix (the rows are the basis
26
vectors). If `F=\GF{2}` then `C` is
27
called a binary code. If `F` has `q` elements
28
then `C` is called a `q`-ary code. The elements
29
of a code `C` are called codewords. The information rate
30
of `C` is
31
32
33
.. math::
34
35
R={\frac{\log_q\vert C\vert}{n}},
36
37
38
where `\vert C\vert` denotes the number of elements of
39
`C`. If `{\bf v}=(v_1,v_2,...,v_n)`,
40
`{\bf w}=(w_1,w_2,...,w_n)` are vectors in
41
`V=F^n` then we define
42
43
44
.. math::
45
46
d({\bf v},{\bf w}) =\vert\{i\ \vert\ 1\leq i\leq n,\ v_i\not= w_i\}\vert
47
48
49
to be the Hamming distance between `{\bf v}` and
50
`{\bf w}`. The function
51
`d:V\times V\rightarrow \Bold{N}` is called the Hamming
52
metric. The weight of a vector (in the Hamming metric) is
53
`d({\bf v},{\bf 0})`. The minimum distance of a linear
54
code is the smallest non-zero weight of a codeword in `C`.
55
The relatively minimum distance is denoted
56
57
58
.. math::
59
60
\delta = d/n.
61
62
A linear code with length
63
`n`, dimension `k`, and minimum distance
64
`d` is called an `[n,k,d]_q`-code and
65
`n,k,d` are called its parameters. A (not necessarily
66
linear) code `C` with length `n`, size
67
`M=|C|`, and minimum distance `d` is called an
68
`(n,M,d)_q`-code (using parentheses instead of square
69
brackets). Of course, `k=\log_q(M)` for linear codes.
70
71
What is the "best" code of a given length? Let `F` be a
72
finite field with `q` elements. Let `A_q(n,d)`
73
denote the largest `M` such that there exists a
74
`(n,M,d)` code in `F^n`. Let `B_q(n,d)`
75
(also denoted `A^{lin}_q(n,d)`) denote the largest
76
`k` such that there exists a `[n,k,d]` code in
77
`F^n`. (Of course, `A_q(n,d)\geq B_q(n,d)`.)
78
Determining `A_q(n,d)` and `B_q(n,d)` is one of
79
the main problems in the theory of error-correcting codes.
80
81
These quantities related to solving a generalization of the
82
childhood game of "20 questions".
83
84
GAME: Player 1 secretly chooses a number from `1` to
85
`M` (`M` is large but fixed). Player 2 asks a
86
series of "yes/no questions" in an attempt to determine that
87
number. Player 1 may lie at most `e` times
88
(`e\geq 0` is fixed). What is the minimum number of "yes/no
89
questions" Player 2 must ask to (always) be able to correctly
90
determine the number Player 1 chose?
91
92
If feedback is not allowed (the only situation considered here),
93
call this minimum number `g(M,e)`.
94
95
Lemma: For fixed `e` and `M`, `g(M,e)` is
96
the smallest `n` such that `A_2(n,2e+1)\geq M`.
97
98
Thus, solving the solving a generalization of the game of "20
99
questions" is equivalent to determining `A_2(n,d)`! Using
100
Sage, you can determine the best known estimates for this number in
101
2 ways:
102
103
(1) Indirectly, using minimum_distance_lower_bound(n,k,F) and
104
minimum_distance_upper_bound(n,k,F) (both of which which connect
105
to the internet using Steven Sivek's linear_code_bound(q,n,k))
106
(2) codesize_upper_bound(n,d,q), dimension_upper_bound(n,d,q).
107
108
This module implements:
109
110
- codesize_upper_bound(n,d,q), for the best known (as of May,
111
2006) upper bound A(n,d) for the size of a code of length n,
112
minimum distance d over a field of size q.
113
114
- dimension_upper_bound(n,d,q), an upper bound
115
`B(n,d)=B_q(n,d)` for the dimension of a linear code of
116
length n, minimum distance d over a field of size q.
117
118
- gilbert_lower_bound(n,q,d), a lower bound for number of
119
elements in the largest code of min distance d in
120
`\GF{q}^n`.
121
122
- gv_info_rate(n,delta,q), `log_q(GLB)/n`, where GLB is
123
the Gilbert lower bound and delta = d/n.
124
125
- gv_bound_asymp(delta,q), asymptotic analog of Gilbert lower
126
bound.
127
128
- plotkin_upper_bound(n,q,d)
129
130
- plotkin_bound_asymp(delta,q), asymptotic analog of Plotkin
131
bound.
132
133
- griesmer_upper_bound(n,q,d)
134
135
- elias_upper_bound(n,q,d)
136
137
- elias_bound_asymp(delta,q), asymptotic analog of Elias bound.
138
139
- hamming_upper_bound(n,q,d)
140
141
- hamming_bound_asymp(delta,q), asymptotic analog of Hamming
142
bound.
143
144
- singleton_upper_bound(n,q,d)
145
146
- singleton_bound_asymp(delta,q), asymptotic analog of Singleton
147
bound.
148
149
- mrrw1_bound_asymp(delta,q), "first" asymptotic
150
McEliese-Rumsey-Rodemich-Welsh bound for the information rate.
151
152
153
PROBLEM: In this module we shall typically either (a) seek bounds
154
on k, given n, d, q, (b) seek bounds on R, delta, q (assuming n is
155
"infinity").
156
157
TODO:
158
159
- Johnson bounds for binary codes.
160
161
- mrrw2_bound_asymp(delta,q), "second" asymptotic
162
McEliese-Rumsey-Rodemich-Welsh bound for the information rate.
163
164
REFERENCES:
165
166
- C. Huffman, V. Pless, Fundamentals of error-correcting codes,
167
Cambridge Univ. Press, 2003.
168
"""
169
170
#*****************************************************************************
171
# Copyright (C) 2006 David Joyner <[email protected]>
172
# 2006 William Stein <[email protected]>
173
#
174
# Distributed under the terms of the GNU General Public License (GPL)
175
#
176
# http://www.gnu.org/licenses/
177
#*****************************************************************************
178
179
from sage.interfaces.all import gap
180
from sage.rings.all import QQ, RR, ZZ, RDF
181
from sage.rings.arith import factorial
182
from sage.functions.all import log, sqrt
183
from sage.misc.decorators import rename_keyword
184
185
@rename_keyword(deprecated='Sage version 4.6', method="algorithm")
186
def codesize_upper_bound(n,d,q,algorithm=None):
187
r"""
188
This computes the minimum value of the upper bound using the
189
algorithms of Singleton, Hamming, Plotkin and Elias.
190
191
If algorithm="gap" then this returns the best known upper
192
bound `A(n,d)=A_q(n,d)` for the size of a code of length n,
193
minimum distance d over a field of size q. The function first
194
checks for trivial cases (like d=1 or n=d), and if the value
195
is in the built-in table. Then it calculates the minimum value
196
of the upper bound using the algorithms of Singleton, Hamming,
197
Johnson, Plotkin and Elias. If the code is binary,
198
`A(n, 2\ell-1) = A(n+1,2\ell)`, so the function
199
takes the minimum of the values obtained from all algorithms for the
200
parameters `(n, 2\ell-1)` and `(n+1, 2\ell)`. This
201
wraps GUAVA's UpperBound( n, d, q ).
202
203
EXAMPLES::
204
sage: codesize_upper_bound(10,3,2)
205
93
206
sage: codesize_upper_bound(10,3,2,algorithm="gap") # requires optional GAP package Guava
207
85
208
209
"""
210
if algorithm=="gap":
211
return int(gap.eval("UpperBound(%s,%s,%s)"%( n, d, q )))
212
else:
213
eub = elias_upper_bound(n,q,d)
214
gub = griesmer_upper_bound(n,q,d)
215
hub = hamming_upper_bound(n,q,d)
216
pub = plotkin_upper_bound(n,q,d)
217
sub = singleton_upper_bound(n,q,d)
218
return min([eub,gub,hub,pub,sub])
219
220
def dimension_upper_bound(n,d,q):
221
r"""
222
Returns an upper bound `B(n,d) = B_q(n,d)` for the
223
dimension of a linear code of length n, minimum distance d over a
224
field of size q.
225
226
EXAMPLES::
227
228
sage: dimension_upper_bound(10,3,2)
229
6
230
"""
231
q = ZZ(q)
232
return int(log(codesize_upper_bound(n,d,q),q))
233
234
def volume_hamming(n,q,r):
235
r"""
236
Returns number of elements in a Hamming ball of radius r in `\GF{q}^n`.
237
Agrees with Guava's SphereContent(n,r,GF(q)).
238
239
EXAMPLES::
240
241
sage: volume_hamming(10,2,3)
242
176
243
"""
244
ans=sum([factorial(n)/(factorial(i)*factorial(n-i))*(q-1)**i for i in range(r+1)])
245
return ans
246
247
def gilbert_lower_bound(n,q,d):
248
r"""
249
Returns lower bound for number of elements in the largest code of
250
minimum distance d in `\GF{q}^n`.
251
252
EXAMPLES::
253
254
sage: gilbert_lower_bound(10,2,3)
255
128/7
256
"""
257
ans=q**n/volume_hamming(n,q,d-1)
258
return ans
259
260
@rename_keyword(deprecated='Sage version 4.6', method="algorithm")
261
def plotkin_upper_bound(n,q,d, algorithm=None):
262
r"""
263
Returns Plotkin upper bound for number of elements in the largest
264
code of minimum distance d in `\GF{q}^n`.
265
266
The algorithm="gap" option wraps Guava's UpperBoundPlotkin.
267
268
EXAMPLES::
269
sage: plotkin_upper_bound(10,2,3)
270
192
271
sage: plotkin_upper_bound(10,2,3,algorithm="gap") # requires optional GAP package Guava
272
192
273
"""
274
if algorithm=="gap":
275
ans=gap.eval("UpperBoundPlotkin(%s,%s,%s)"%(n,d,q))
276
#print "calling Guava ..."
277
return QQ(ans)
278
else:
279
t = 1 - 1/q
280
if (q==2) and (n == 2*d) and (d%2 == 0):
281
return 4*d
282
elif (q==2) and (n == 2*d + 1) and (d%2 == 1):
283
return 4*d + 4
284
elif d > t*n:
285
return int(d/( d - t*n))
286
elif d < t*n + 1:
287
fact = (d-1) / t
288
if RR(fact)==RR(int(fact)):
289
fact = int(fact) + 1
290
return int(d/( d - t * fact)) * q**(n - fact)
291
292
@rename_keyword(deprecated='Sage version 4.6', method="algorithm")
293
def griesmer_upper_bound(n,q,d,algorithm=None):
294
r"""
295
Returns the Griesmer upper bound for number of elements in the
296
largest code of minimum distance d in `\GF{q}^n`.
297
Wraps GAP's UpperBoundGriesmer.
298
299
EXAMPLES::
300
sage: griesmer_upper_bound(10,2,3)
301
128
302
sage: griesmer_upper_bound(10,2,3,algorithm="gap") # requires optional GAP package Guava
303
128
304
"""
305
if algorithm=="gap":
306
#print "calling Guava ..."
307
ans=gap.eval("UpperBoundGriesmer(%s,%s,%s)"%(n,d,q))
308
return QQ(ans)
309
else:
310
den = 1
311
s = 0
312
k = 0
313
add = 0
314
while s <= n:
315
if not(add == 1):
316
if d%den==0:
317
add = int(d/den)
318
else:
319
add = int(d/den)+1
320
s = s + add
321
den = den * q
322
k = k + 1
323
return q**(k-1)
324
325
326
@rename_keyword(deprecated='Sage version 4.6', method="algorithm")
327
def elias_upper_bound(n,q,d,algorithm=None):
328
r"""
329
Returns the Elias upper bound for number of elements in the largest
330
code of minimum distance d in `\GF{q}^n`. Wraps
331
GAP's UpperBoundElias.
332
333
EXAMPLES::
334
335
sage: elias_upper_bound(10,2,3)
336
232
337
sage: elias_upper_bound(10,2,3,algorithm="gap") # requires optional GAP package Guava
338
232
339
340
"""
341
r = 1-1/q
342
if algorithm=="gap":
343
#print "calling Guava ..."
344
ans=gap.eval("UpperBoundElias(%s,%s,%s)"%(n,d,q))
345
return QQ(ans)
346
else:
347
def ff(n,d,w,q):
348
return r*n*d*q**n/((w**2-2*r*n*w+r*n*d)*volume_hamming(n,q,w));
349
def get_list(n,d,q):
350
I = []
351
for i in range(1,int(r*n)+1):
352
if i**2-2*r*n*i+r*n*d>0:
353
I.append(i)
354
return I
355
I = get_list(n,d,q)
356
bnd = min([ff(n,d,w,q) for w in I])
357
return int(bnd)
358
359
def hamming_upper_bound(n,q,d):
360
r"""
361
Returns the Hamming upper bound for number of elements in the
362
largest code of minimum distance d in `\GF{q}^n`.
363
Wraps GAP's UpperBoundHamming.
364
365
The Hamming bound (also known as the sphere packing bound) returns
366
an upper bound on the size of a code of length n, minimum distance
367
d, over a field of size q. The Hamming bound is obtained by
368
dividing the contents of the entire space
369
`\GF{q}^n` by the contents of a ball with radius
370
floor((d-1)/2). As all these balls are disjoint, they can never
371
contain more than the whole vector space.
372
373
374
.. math::
375
376
M \leq {q^n \over V(n,e)},
377
378
379
380
where M is the maximum number of codewords and `V(n,e)` is
381
equal to the contents of a ball of radius e. This bound is useful
382
for small values of d. Codes for which equality holds are called
383
perfect.
384
385
EXAMPLES::
386
387
sage: hamming_upper_bound(10,2,3)
388
93
389
"""
390
return int((q**n)/(volume_hamming(n, q, int((d-1)/2))))
391
392
def singleton_upper_bound(n,q,d):
393
r"""
394
Returns the Singleton upper bound for number of elements in the
395
largest code of minimum distance d in `\GF{q}^n`.
396
Wraps GAP's UpperBoundSingleton.
397
398
This bound is based on the shortening of codes. By shortening an
399
`(n, M, d)` code d-1 times, an `(n-d+1,M,1)` code
400
results, with `M \leq q^n-d+1`. Thus
401
402
403
.. math::
404
405
M \leq q^{n-d+1}.
406
407
408
409
Codes that meet this bound are called maximum distance separable
410
(MDS).
411
412
EXAMPLES::
413
sage: singleton_upper_bound(10,2,3)
414
256
415
"""
416
return q**(n - d + 1)
417
418
def gv_info_rate(n,delta,q):
419
"""
420
GV lower bound for information rate of a q-ary code of length n
421
minimum distance delta\*n
422
423
EXAMPLES::
424
425
sage: RDF(gv_info_rate(100,1/4,3))
426
0.367049926083
427
"""
428
q = ZZ(q)
429
ans=log(gilbert_lower_bound(n,q,int(n*delta)),q)/n
430
return ans
431
432
def entropy(x,q):
433
"""
434
Computes the entropy at `x` on the `q`-ary symmetric channel.
435
436
INPUT:
437
438
- ``x`` - real number in the interval `[0, 1]`.
439
440
- ``q`` - integer greater than 1. This is the base of the logarithm.
441
442
EXAMPLES::
443
444
sage: entropy(0, 2)
445
0
446
sage: entropy(1/5,4)
447
-1/5*log(1/5)/log(4) - 4/5*log(4/5)/log(4) + 1/5*log(3)/log(4)
448
sage: entropy(1, 3)
449
log(2)/log(3)
450
451
Check that values not within the limits are properly handled::
452
453
sage: entropy(1.1, 2)
454
Traceback (most recent call last):
455
...
456
ValueError: The entropy function is defined only for x in the interval [0, 1]
457
sage: entropy(1, 1)
458
Traceback (most recent call last):
459
...
460
ValueError: The value q must be an integer greater than 1
461
"""
462
if x < 0 or x > 1:
463
raise ValueError("The entropy function is defined only for x in the"
464
" interval [0, 1]")
465
q = ZZ(q) # This will error out if q is not an integer
466
if q < 2: # Here we check that q is actually at least 2
467
raise ValueError("The value q must be an integer greater than 1")
468
if x == 0:
469
return 0
470
if x == 1:
471
return log(q-1,q)
472
H = x*log(q-1,q)-x*log(x,q)-(1-x)*log(1-x,q)
473
return H
474
475
def gv_bound_asymp(delta,q):
476
"""
477
Computes the asymptotic GV bound for the information rate, R.
478
479
EXAMPLES::
480
sage: RDF(gv_bound_asymp(1/4,2))
481
0.188721875541
482
sage: f = lambda x: gv_bound_asymp(x,2)
483
sage: plot(f,0,1)
484
"""
485
return (1-entropy(delta,q))
486
487
488
def hamming_bound_asymp(delta,q):
489
"""
490
Computes the asymptotic Hamming bound for the information rate.
491
492
EXAMPLES::
493
sage: RDF(hamming_bound_asymp(1/4,2))
494
0.4564355568
495
sage: f = lambda x: hamming_bound_asymp(x,2)
496
sage: plot(f,0,1)
497
"""
498
return (1-entropy(delta/2,q))
499
500
def singleton_bound_asymp(delta,q):
501
"""
502
Computes the asymptotic Singleton bound for the information rate.
503
504
EXAMPLES::
505
sage: singleton_bound_asymp(1/4,2)
506
3/4
507
sage: f = lambda x: singleton_bound_asymp(x,2)
508
sage: plot(f,0,1)
509
"""
510
return (1-delta)
511
512
def plotkin_bound_asymp(delta,q):
513
"""
514
Computes the asymptotic Plotkin bound for the information rate,
515
provided `0 < \delta < 1-1/q`.
516
517
EXAMPLES::
518
519
sage: plotkin_bound_asymp(1/4,2)
520
1/2
521
"""
522
r = 1-1/q
523
return (1-delta/r)
524
525
def elias_bound_asymp(delta,q):
526
"""
527
Computes the asymptotic Elias bound for the information rate,
528
provided `0 < \delta 1-1/q`.
529
530
EXAMPLES::
531
532
sage: elias_bound_asymp(1/4,2)
533
0.39912396330...
534
"""
535
r = 1-1/q
536
return RDF((1-entropy(r-sqrt(r*(r-delta)), q)))
537
538
def mrrw1_bound_asymp(delta,q):
539
"""
540
Computes the first asymptotic McEliese-Rumsey-Rodemich-Welsh bound
541
for the information rate, provided `0 < \delta < 1-1/q`.
542
543
EXAMPLES::
544
545
sage: mrrw1_bound_asymp(1/4,2)
546
0.354578902665
547
"""
548
return RDF(entropy((q-1-delta*(q-2)-2*sqrt((q-1)*delta*(1-delta)))/q,q))
549
550
551
552
553
554
555