Go fountain code library
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..d645695
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,202 @@
+
+ Apache License
+ Version 2.0, January 2004
+ http://www.apache.org/licenses/
+
+ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+ 1. Definitions.
+
+ "License" shall mean the terms and conditions for use, reproduction,
+ and distribution as defined by Sections 1 through 9 of this document.
+
+ "Licensor" shall mean the copyright owner or entity authorized by
+ the copyright owner that is granting the License.
+
+ "Legal Entity" shall mean the union of the acting entity and all
+ other entities that control, are controlled by, or are under common
+ control with that entity. For the purposes of this definition,
+ "control" means (i) the power, direct or indirect, to cause the
+ direction or management of such entity, whether by contract or
+ otherwise, or (ii) ownership of fifty percent (50%) or more of the
+ outstanding shares, or (iii) beneficial ownership of such entity.
+
+ "You" (or "Your") shall mean an individual or Legal Entity
+ exercising permissions granted by this License.
+
+ "Source" form shall mean the preferred form for making modifications,
+ including but not limited to software source code, documentation
+ source, and configuration files.
+
+ "Object" form shall mean any form resulting from mechanical
+ transformation or translation of a Source form, including but
+ not limited to compiled object code, generated documentation,
+ and conversions to other media types.
+
+ "Work" shall mean the work of authorship, whether in Source or
+ Object form, made available under the License, as indicated by a
+ copyright notice that is included in or attached to the work
+ (an example is provided in the Appendix below).
+
+ "Derivative Works" shall mean any work, whether in Source or Object
+ form, that is based on (or derived from) the Work and for which the
+ editorial revisions, annotations, elaborations, or other modifications
+ represent, as a whole, an original work of authorship. For the purposes
+ of this License, Derivative Works shall not include works that remain
+ separable from, or merely link (or bind by name) to the interfaces of,
+ the Work and Derivative Works thereof.
+
+ "Contribution" shall mean any work of authorship, including
+ the original version of the Work and any modifications or additions
+ to that Work or Derivative Works thereof, that is intentionally
+ submitted to Licensor for inclusion in the Work by the copyright owner
+ or by an individual or Legal Entity authorized to submit on behalf of
+ the copyright owner. For the purposes of this definition, "submitted"
+ means any form of electronic, verbal, or written communication sent
+ to the Licensor or its representatives, including but not limited to
+ communication on electronic mailing lists, source code control systems,
+ and issue tracking systems that are managed by, or on behalf of, the
+ Licensor for the purpose of discussing and improving the Work, but
+ excluding communication that is conspicuously marked or otherwise
+ designated in writing by the copyright owner as "Not a Contribution."
+
+ "Contributor" shall mean Licensor and any individual or Legal Entity
+ on behalf of whom a Contribution has been received by Licensor and
+ subsequently incorporated within the Work.
+
+ 2. Grant of Copyright License. Subject to the terms and conditions of
+ this License, each Contributor hereby grants to You a perpetual,
+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+ copyright license to reproduce, prepare Derivative Works of,
+ publicly display, publicly perform, sublicense, and distribute the
+ Work and such Derivative Works in Source or Object form.
+
+ 3. Grant of Patent License. Subject to the terms and conditions of
+ this License, each Contributor hereby grants to You a perpetual,
+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+ (except as stated in this section) patent license to make, have made,
+ use, offer to sell, sell, import, and otherwise transfer the Work,
+ where such license applies only to those patent claims licensable
+ by such Contributor that are necessarily infringed by their
+ Contribution(s) alone or by combination of their Contribution(s)
+ with the Work to which such Contribution(s) was submitted. If You
+ institute patent litigation against any entity (including a
+ cross-claim or counterclaim in a lawsuit) alleging that the Work
+ or a Contribution incorporated within the Work constitutes direct
+ or contributory patent infringement, then any patent licenses
+ granted to You under this License for that Work shall terminate
+ as of the date such litigation is filed.
+
+ 4. Redistribution. You may reproduce and distribute copies of the
+ Work or Derivative Works thereof in any medium, with or without
+ modifications, and in Source or Object form, provided that You
+ meet the following conditions:
+
+ (a) You must give any other recipients of the Work or
+ Derivative Works a copy of this License; and
+
+ (b) You must cause any modified files to carry prominent notices
+ stating that You changed the files; and
+
+ (c) You must retain, in the Source form of any Derivative Works
+ that You distribute, all copyright, patent, trademark, and
+ attribution notices from the Source form of the Work,
+ excluding those notices that do not pertain to any part of
+ the Derivative Works; and
+
+ (d) If the Work includes a "NOTICE" text file as part of its
+ distribution, then any Derivative Works that You distribute must
+ include a readable copy of the attribution notices contained
+ within such NOTICE file, excluding those notices that do not
+ pertain to any part of the Derivative Works, in at least one
+ of the following places: within a NOTICE text file distributed
+ as part of the Derivative Works; within the Source form or
+ documentation, if provided along with the Derivative Works; or,
+ within a display generated by the Derivative Works, if and
+ wherever such third-party notices normally appear. The contents
+ of the NOTICE file are for informational purposes only and
+ do not modify the License. You may add Your own attribution
+ notices within Derivative Works that You distribute, alongside
+ or as an addendum to the NOTICE text from the Work, provided
+ that such additional attribution notices cannot be construed
+ as modifying the License.
+
+ You may add Your own copyright statement to Your modifications and
+ may provide additional or different license terms and conditions
+ for use, reproduction, or distribution of Your modifications, or
+ for any such Derivative Works as a whole, provided Your use,
+ reproduction, and distribution of the Work otherwise complies with
+ the conditions stated in this License.
+
+ 5. Submission of Contributions. Unless You explicitly state otherwise,
+ any Contribution intentionally submitted for inclusion in the Work
+ by You to the Licensor shall be under the terms and conditions of
+ this License, without any additional terms or conditions.
+ Notwithstanding the above, nothing herein shall supersede or modify
+ the terms of any separate license agreement you may have executed
+ with Licensor regarding such Contributions.
+
+ 6. Trademarks. This License does not grant permission to use the trade
+ names, trademarks, service marks, or product names of the Licensor,
+ except as required for reasonable and customary use in describing the
+ origin of the Work and reproducing the content of the NOTICE file.
+
+ 7. Disclaimer of Warranty. Unless required by applicable law or
+ agreed to in writing, Licensor provides the Work (and each
+ Contributor provides its Contributions) on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+ implied, including, without limitation, any warranties or conditions
+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+ PARTICULAR PURPOSE. You are solely responsible for determining the
+ appropriateness of using or redistributing the Work and assume any
+ risks associated with Your exercise of permissions under this License.
+
+ 8. Limitation of Liability. In no event and under no legal theory,
+ whether in tort (including negligence), contract, or otherwise,
+ unless required by applicable law (such as deliberate and grossly
+ negligent acts) or agreed to in writing, shall any Contributor be
+ liable to You for damages, including any direct, indirect, special,
+ incidental, or consequential damages of any character arising as a
+ result of this License or out of the use or inability to use the
+ Work (including but not limited to damages for loss of goodwill,
+ work stoppage, computer failure or malfunction, or any and all
+ other commercial damages or losses), even if such Contributor
+ has been advised of the possibility of such damages.
+
+ 9. Accepting Warranty or Additional Liability. While redistributing
+ the Work or Derivative Works thereof, You may choose to offer,
+ and charge a fee for, acceptance of support, warranty, indemnity,
+ or other liability obligations and/or rights consistent with this
+ License. However, in accepting such obligations, You may act only
+ on Your own behalf and on Your sole responsibility, not on behalf
+ of any other Contributor, and only if You agree to indemnify,
+ defend, and hold each Contributor harmless for any liability
+ incurred by, or claims asserted against, such Contributor by reason
+ of your accepting any such warranty or additional liability.
+
+ END OF TERMS AND CONDITIONS
+
+ APPENDIX: How to apply the Apache License to your work.
+
+ To apply the Apache License to your work, attach the following
+ boilerplate notice, with the fields enclosed by brackets "[]"
+ replaced with your own identifying information. (Don't include
+ the brackets!) The text should be enclosed in the appropriate
+ comment syntax for the file format. We also recommend that a
+ file or class name and description of purpose be included on the
+ same "printed page" as the copyright notice for easier
+ identification within third-party archives.
+
+ Copyright [yyyy] [name of copyright owner]
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
diff --git a/binary.go b/binary.go
new file mode 100644
index 0000000..c752837
--- /dev/null
+++ b/binary.go
@@ -0,0 +1,116 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math/rand"
+)
+
+// Random binary fountain code. In this code, the constituent source blocks in
+// a code block are selected randomly and independently.
+
+// BinaryCodec contains the codec information for the random binary fountain
+// encoder and decoder.
+// Implements fountain.Codec
+type binaryCodec struct {
+ // numSourceBlocks is the number of source blocks (N) the source message is split into.
+ numSourceBlocks int
+}
+
+// NewBinaryCodec returns a codec implementing the binary fountain code,
+// where source blocks composing each LT block are chosen randomly and independently.
+func NewBinaryCodec(numSourceBlocks int) Codec {
+ return &binaryCodec{numSourceBlocks: numSourceBlocks}
+}
+
+// SourceBlocks returns the number of source blocks used in the codec.
+func (c *binaryCodec) SourceBlocks() int {
+ return c.numSourceBlocks
+}
+
+// PickIndices finds the source indices for a code block given an ID and
+// a random seed. Uses the Mersenne Twister internally.
+func (c *binaryCodec) PickIndices(codeBlockIndex int64) []int {
+ random := rand.New(NewMersenneTwister(codeBlockIndex))
+
+ var indices []int
+ for b := 0; b < c.SourceBlocks(); b++ {
+ if random.Intn(2) == 1 {
+ indices = append(indices, b)
+ }
+ }
+
+ return indices
+}
+
+// GenerateIntermediateBlocks simply returns the partition of the input message
+// into source blocks. It does not perform any additional precoding.
+func (c *binaryCodec) GenerateIntermediateBlocks(message []byte, numBlocks int) []block {
+ long, short := partitionBytes(message, c.numSourceBlocks)
+ source := equalizeBlockLengths(long, short)
+
+ return source
+}
+
+// NewDecoder creates a new binary fountain code decoder
+func (c *binaryCodec) NewDecoder(messageLength int) Decoder {
+ return newBinaryDecoder(c, messageLength)
+}
+
+// binaryDecoder is the state required to decode a combinatoric fountain
+// code message.
+type binaryDecoder struct {
+ codec binaryCodec
+ messageLength int
+
+ // The sparse equation matrix used for decoding.
+ matrix sparseMatrix
+}
+
+// newBinaryDecoder creates a new decoder for a particular message.
+// The codec parameters used to create the original encoding blocks must be provided.
+// The decoder is only valid for decoding code blocks for a particular message.
+func newBinaryDecoder(c *binaryCodec, length int) *binaryDecoder {
+ return &binaryDecoder{
+ codec: *c,
+ messageLength: length,
+ matrix: sparseMatrix{
+ coeff: make([][]int, c.numSourceBlocks),
+ v: make([]block, c.numSourceBlocks),
+ }}
+}
+
+// AddBlocks adds a set of encoded blocks to the decoder. Returns true if the
+// message can be fully decoded. False if there is insufficient information.
+func (d *binaryDecoder) AddBlocks(blocks []LTBlock) bool {
+ for i := range blocks {
+ d.matrix.addEquation(d.codec.PickIndices(blocks[i].BlockCode),
+ block{data: blocks[i].Data})
+ }
+ return d.matrix.determined()
+}
+
+// Decode extracts the decoded message from the decoder. If the decoder does
+// not have sufficient information to produce an output, returns a nil slice.
+func (d *binaryDecoder) Decode() []byte {
+ if !d.matrix.determined() {
+ return nil
+ }
+
+ d.matrix.reduce()
+
+ lenLong, lenShort, numLong, numShort := partition(d.messageLength, d.codec.numSourceBlocks)
+ return d.matrix.reconstruct(d.messageLength, lenLong, lenShort, numLong, numShort)
+}
diff --git a/binary_test.go b/binary_test.go
new file mode 100644
index 0000000..1f60b90
--- /dev/null
+++ b/binary_test.go
@@ -0,0 +1,122 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+// TestBinaryDecoder ensures that the binary fountain produces a round-trip
+// decodeable message using the corresponding decoder.
+func TestBinaryDecoder(t *testing.T) {
+ c := NewBinaryCodec(13)
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ ids := make([]int64, 45)
+ random := rand.New(rand.NewSource(8923489))
+ for i := range ids {
+ ids[i] = int64(random.Intn(100000))
+ }
+
+ blocks := EncodeLTBlocks(message, ids, c)
+ t.Log("blocks =", blocks)
+
+ d := newBinaryDecoder(c.(*binaryCodec), len(message))
+
+ for i := 0; i < 16; i++ {
+ d.AddBlocks([]LTBlock{blocks[i]})
+ if testing.Verbose() {
+ printMatrix(d.matrix, t)
+ }
+ }
+
+ d.matrix.reduce()
+ t.Log("REDUCE")
+ printMatrix(d.matrix, t)
+
+ decoded := d.Decode()
+ printMatrix(d.matrix, t)
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Decoded message doesn't match original. Got %v, want %v", decoded, message)
+ }
+}
+
+// TestbinaryDecoderBlockTable tests many combinations of fountain block ID
+// combinations to ensure that the codec has the expected reconstruction
+// properties.
+func TestBinaryDecoderBlockTable(t *testing.T) {
+ c := NewBinaryCodec(13)
+
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ random := rand.New(rand.NewSource(8234923))
+
+ moreBlocksNeeded := 0
+ for i := 0; i < 100; i++ {
+ r := rand.New(rand.NewSource(random.Int63()))
+ ids := make([]int64, 45)
+ for i := range ids {
+ ids[i] = int64(r.Intn(100000))
+ }
+ blocks := EncodeLTBlocks(message, ids, c)
+
+ d := newBinaryDecoder(c.(*binaryCodec), len(message))
+ d.AddBlocks(blocks[0:30])
+ if !d.matrix.determined() {
+ moreBlocksNeeded++
+ d.AddBlocks(blocks[31:46])
+ }
+ decoded := d.Decode()
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Decoded message doesn't match original. Got %v, want %v", decoded, message)
+ }
+ }
+
+ if moreBlocksNeeded > 2 {
+ t.Errorf("Needed too many high-block-count decoding sequences: %d", moreBlocksNeeded)
+ }
+}
+
+// TestBinaryDecodeMessageTable tests a large number of source messages to make
+// sure they are all reconstructed accurately. This provides assurance that the
+// decoder is functioning accurately.
+func TestBinaryDecodeMessageTable(t *testing.T) {
+ c := NewBinaryCodec(10)
+ random := rand.New(rand.NewSource(8234982))
+ for i := 0; i < 100; i++ {
+ r := rand.New(rand.NewSource(random.Int63()))
+ messageLen := r.Intn(1000) + 1000
+ message := make([]byte, messageLen)
+ for j := 0; j < len(message); j++ {
+ message[j] = byte(r.Intn(200))
+ }
+ ids := make([]int64, 50)
+ for i := range ids {
+ ids[i] = int64(r.Intn(100000))
+ }
+ blocks := EncodeLTBlocks(message, ids, c)
+
+ d := newBinaryDecoder(c.(*binaryCodec), len(message))
+ d.AddBlocks(blocks[0:25])
+ if !d.matrix.determined() {
+ t.Errorf("Message should be determined after 25 blocks")
+ } else {
+ decoded := d.Decode()
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Incorrect message decode. Length=%d, message=%v", len(message), message)
+ }
+ }
+ }
+}
diff --git a/block.go b/block.go
new file mode 100644
index 0000000..e8af62c
--- /dev/null
+++ b/block.go
@@ -0,0 +1,248 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+// A block represents a contiguous range of data being encoded or decoded,
+// or a block of coded data. Details of how the source text is split into blocks
+// is governed by the particular fountain code used.
+type block struct {
+ // Data content of this source or code block.
+ data []byte
+
+ // How many padding bytes this block has at the end.
+ padding int
+}
+
+// newBlock creates a new block with a given length. The block will initially be
+// all padding.
+func newBlock(len int) *block {
+ return &block{padding: len}
+}
+
+// length returns the length of the block in bytes. Counts data bytes as well
+// as any padding.
+func (b *block) length() int {
+ return len(b.data) + b.padding
+}
+
+func (b *block) empty() bool {
+ return b.length() == 0
+}
+
+// A common operation is to XOR entire code blocks together with other blocks.
+// When this is done, padding bytes count as 0 (that is XOR identity), and the
+// destination block will be modified so that its data is large enough to
+// contain the result of the XOR.
+func (b *block) xor(a block) {
+ if len(b.data) < len(a.data) {
+ var inc = len(a.data) - len(b.data)
+ b.data = append(b.data, make([]byte, inc)...)
+ if b.padding > inc {
+ b.padding -= inc
+ } else {
+ b.padding = 0
+ }
+ }
+
+ for i := 0; i < len(a.data); i++ {
+ b.data[i] ^= a.data[i]
+ }
+}
+
+// partitionBytes partitions an input text into a sequence of p blocks. The
+// sizes of the blocks will be given by the partition() function. The last
+// block may have padding.
+// Return values: the slice of longer blocks, the slice of shorter blocks.
+// Within each block slice, all will have uniform lengths.
+func partitionBytes(in []byte, p int) ([]block, []block) {
+ sliceIntoBlocks := func(in []byte, num, length int) ([]block, []byte) {
+ blocks := make([]block, num)
+ for i := range blocks {
+ if len(in) > length {
+ blocks[i].data, in = in[:length], in[length:]
+ } else {
+ blocks[i].data, in = in, []byte{}
+ }
+ if len(blocks[i].data) < length {
+ blocks[i].padding = length - len(blocks[i].data)
+ }
+ }
+ return blocks, in
+ }
+
+ lenLong, lenShort, numLong, numShort := partition(len(in), p)
+ long, in := sliceIntoBlocks(in, numLong, lenLong)
+ short, _ := sliceIntoBlocks(in, numShort, lenShort)
+ return long, short
+}
+
+// equalizeBlockLengths adds padding to all short blocks to make them equal in
+// size to the long blocks. The caller should ensure that all the longBlocks
+// have the same length.
+// Returns a block slice containing all the long and short blocks.
+func equalizeBlockLengths(longBlocks, shortBlocks []block) []block {
+ if len(longBlocks) == 0 {
+ return shortBlocks
+ }
+ if len(shortBlocks) == 0 {
+ return longBlocks
+ }
+
+ for i := range shortBlocks {
+ shortBlocks[i].padding += longBlocks[0].length() - shortBlocks[i].length()
+ }
+
+ blocks := make([]block, len(longBlocks)+len(shortBlocks))
+ copy(blocks, longBlocks)
+ copy(blocks[len(longBlocks):], shortBlocks)
+ return blocks
+}
+
+// sparseMatrix is the block decoding data structure. It is a sparse matrix of
+// XOR equations. The coefficients are the indices of the source blocks which
+// are XORed together to produce the values. So if equation _i_ of the matrix is
+// block_0 ^ block_2 ^ block_3 ^ block_9 = [0xD2, 0x38]
+// that would be represented as coeff[i] = [0, 2, 3, 9], v[i].data = [0xD2, 0x38]
+// Example: The sparse coefficient matrix
+// | 0 1 1 0 |
+// | 0 1 0 0 |
+// | 0 0 1 1 |
+// | 0 0 0 1 |
+// would be represented as
+// [ [ 1, 2],
+// [ 1 ],
+// [ 2, 3],
+// [ 3 ]]
+// Every row has M[i][0] >= i. If we added components [2] to this matrix, it
+// would replace the M[2] row ([2, 3]), and then the resulting component vector
+// after cancellation against that row ([3]) would be used instead of the
+// original equation. In this case, M[3] is already populated, so the new
+// equation is redundant (and could be used for ECC, theoretically), but if we
+// didn't have an entry in M[3], it would be placed there.
+// The values were omitted from this discussion, but they follow along by doing
+// XOR operations as the components are reduced during insertion.
+type sparseMatrix struct {
+ coeff [][]int
+ v []block
+}
+
+// xorRow performs a reduction of the given candidate equation (indices, b)
+// with the specified matrix row (index s). It does so by XORing the values,
+// and then taking the symmetric difference of the coefficients of that matrix
+// row and the provided indices. (That is, the "set XOR".) Assumes both
+// coefficient slices are sorted.
+func (m *sparseMatrix) xorRow(s int, indices []int, b block) ([]int, block) {
+ b.xor(m.v[s])
+
+ var newIndices []int
+ coeffs := m.coeff[s]
+ var i, j int
+ for i < len(coeffs) && j < len(indices) {
+ index := indices[j]
+ if coeffs[i] == index {
+ i++
+ j++
+ } else if coeffs[i] < index {
+ newIndices = append(newIndices, coeffs[i])
+ i++
+ } else {
+ newIndices = append(newIndices, index)
+ j++
+ }
+ }
+
+ newIndices = append(newIndices, coeffs[i:]...)
+ newIndices = append(newIndices, indices[j:]...)
+ return newIndices, b
+}
+
+// addEquation adds an XOR equation to the decode matrix. The online decode
+// strategy is a variant of that of Bioglio, Grangetto, and Gaeta
+// (http://www.di.unito.it/~bioglio/Papers/CL2009-lt.pdf) It maintains the
+// invariant that either coeff[i][0] == i or len(coeff[i]) == 0. That is, while
+// adding an equation to the matrix, it ensures that the decode matrix remains
+// triangular.
+func (m *sparseMatrix) addEquation(components []int, b block) {
+ // This loop reduces the incoming equation by XOR until it either fits into
+ // an empty row in the decode matrix or is discarded as redundant.
+ for len(components) > 0 && len(m.coeff[components[0]]) > 0 {
+ s := components[0]
+ if len(components) >= len(m.coeff[s]) {
+ components, b = m.xorRow(s, components, b)
+ } else {
+ // Swap the existing row for the new one, reduce the existing one and
+ // see if it fits elsewhere.
+ components, m.coeff[s] = m.coeff[s], components
+ b, m.v[s] = m.v[s], b
+ }
+ }
+
+ if len(components) > 0 {
+ m.coeff[components[0]] = components
+ m.v[components[0]] = b
+ }
+}
+
+// Check to see if the decode matrix is fully specified. This is true when
+// all rows have non-empty coefficient slices.
+// TODO(gbillock): is there a weakness here if an auxiliary block is unpopulated?
+func (m *sparseMatrix) determined() bool {
+ for _, r := range m.coeff {
+ if len(r) == 0 {
+ return false
+ }
+ }
+ return true
+}
+
+// reduce performs Gaussian Elimination over the whole matrix. Presumes
+// the matrix is triangular, and that the method is not called unless there is
+// enough data for a solution.
+// TODO(gbillock): Could profitably do this online as well?
+func (m *sparseMatrix) reduce() {
+ for i := len(m.coeff) - 1; i >= 0; i-- {
+ for j := 0; j < i; j++ {
+ ci, cj := m.coeff[i], m.coeff[j]
+ for k := 1; k < len(cj); k++ {
+ if cj[k] == ci[0] {
+ m.v[j].xor(m.v[i])
+ continue
+ }
+ }
+ }
+ // All but the leading coefficient in the rows have been reduced out.
+ m.coeff[i] = m.coeff[i][0:1]
+ }
+}
+
+// reconstruct pastes the fully reduced values in the sparse matrix result column
+// into a new byte array and returns it. The length/number parameters are typically
+// those given by partition().
+// lenLong is how many long blocks there are.
+// lenShort is how many short blocks there are (following the long blocks).
+// numLong is how many bytes are in the long blocks.
+// numShort is how many bytes the short blocks are.
+func (m *sparseMatrix) reconstruct(totalLength, lenLong, lenShort, numLong, numShort int) []byte {
+ out := make([]byte, totalLength)
+ out = out[0:0]
+ for i := 0; i < numLong; i++ {
+ out = append(out, m.v[i].data[0:lenLong]...)
+ }
+ for i := numLong; i < numLong+numShort; i++ {
+ out = append(out, m.v[i].data[0:lenShort]...)
+ }
+
+ return out
+}
diff --git a/block_test.go b/block_test.go
new file mode 100644
index 0000000..0b2f7de
--- /dev/null
+++ b/block_test.go
@@ -0,0 +1,280 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "bytes"
+ "reflect"
+ "testing"
+)
+
+func TestBlockLength(t *testing.T) {
+ var lengthTests = []struct {
+ b block
+ len int
+ }{
+ {block{}, 0},
+ {block{[]byte{1, 0, 1}, 0}, 3},
+ {block{[]byte{1, 0, 1}, 1}, 4},
+ }
+
+ for _, i := range lengthTests {
+ if i.b.length() != i.len {
+ t.Errorf("Length of b is %d, should be %d", i.b.length(), i.len)
+ }
+ if (i.len == 0) != i.b.empty() {
+ t.Errorf("Emptiness check error. Got %v, want %v", i.b.empty(), i.len == 0)
+ }
+ }
+}
+
+func TestBlockXor(t *testing.T) {
+ var xorTests = []struct {
+ a block
+ b block
+ out block
+ }{
+ {block{[]byte{1, 0, 1}, 0}, block{[]byte{1, 1, 1}, 0}, block{[]byte{0, 1, 0}, 0}},
+ {block{[]byte{1}, 0}, block{[]byte{0, 14, 6}, 0}, block{[]byte{1, 14, 6}, 0}},
+ {block{}, block{[]byte{100, 200}, 0}, block{[]byte{100, 200}, 0}},
+ {block{[]byte{}, 5}, block{[]byte{0, 1, 0}, 0}, block{[]byte{0, 1, 0}, 2}},
+ {block{[]byte{}, 5}, block{[]byte{0, 1, 0, 2, 3}, 0}, block{[]byte{0, 1, 0, 2, 3}, 0}},
+ {block{[]byte{}, 5}, block{[]byte{0, 1, 0, 2, 3, 7}, 0}, block{[]byte{0, 1, 0, 2, 3, 7}, 0}},
+ {block{[]byte{1}, 4}, block{[]byte{0, 1, 0, 2, 3, 7}, 0}, block{[]byte{1, 1, 0, 2, 3, 7}, 0}},
+ }
+
+ for _, i := range xorTests {
+ t.Logf("...Testing %v XOR %v", i.a, i.b)
+ originalLength := i.a.length()
+ i.a.xor(i.b)
+ if i.a.length() < originalLength {
+ t.Errorf("Length shrunk. Got %d, want length >= %d", i.a.length(), originalLength)
+ }
+ if len(i.a.data) != len(i.b.data) {
+ t.Errorf("a and b data should be same length after xor. a len=%d, b len=%d", len(i.a.data), len(i.b.data))
+ }
+
+ if !bytes.Equal(i.a.data, i.out.data) {
+ t.Errorf("XOR value is %v : should be %v", i.a.data, i.out.data)
+ }
+ }
+}
+
+func TestPartitionBytes(t *testing.T) {
+ a := make([]byte, 100)
+ for i := 0; i < len(a); i++ {
+ a[i] = byte(i)
+ }
+
+ var partitionTests = []struct {
+ numPartitions int
+ lenLong, lenShort int
+ }{
+ {11, 1, 10},
+ {3, 1, 2},
+ }
+
+ for _, i := range partitionTests {
+ t.Logf("Partitioning %v into %d", a, i.numPartitions)
+ long, short := partitionBytes(a, i.numPartitions)
+ if len(long) != i.lenLong {
+ t.Errorf("Got %d long blocks, should have %d", len(long), i.lenLong)
+ }
+ if len(short) != i.lenShort {
+ t.Errorf("Got %d short blocks, should have %d", len(short), i.lenShort)
+ }
+ if short[len(short)-1].padding != 0 {
+ t.Errorf("Should fit blocks exactly, have last padding %d", short[len(short)-1].padding)
+ }
+ if long[0].data[0] != 0 {
+ t.Errorf("Long block should be first. First value is %v", long[0].data)
+ }
+ }
+}
+
+func TestEqualizeBlockLengths(t *testing.T) {
+ b := []byte("abcdefghijklmnopq")
+ var equalizeTests = []struct {
+ numPartitions int
+ length int
+ padding int
+ }{
+ {1, 17, 0},
+ {2, 9, 1},
+ {3, 6, 1},
+ {4, 5, 1},
+ {5, 4, 1},
+ {6, 3, 1},
+ {7, 3, 1},
+ {8, 3, 1},
+ {9, 2, 1},
+ {10, 2, 1},
+ {16, 2, 1},
+ {17, 1, 0},
+ }
+
+ for _, i := range equalizeTests {
+ long, short := partitionBytes(b, i.numPartitions)
+ blocks := equalizeBlockLengths(long, short)
+ if len(blocks) != i.numPartitions {
+ t.Errorf("Got %d blocks, should have %d", len(blocks), i.numPartitions)
+ }
+ for k := range blocks {
+ if blocks[k].length() != i.length {
+ t.Errorf("Got block length %d for block %d, should be %d",
+ blocks[0].length(), k, i.length)
+ }
+ }
+ if blocks[len(blocks)-1].padding != i.padding {
+ t.Errorf("Padding of last block is %d, should be %d",
+ blocks[len(blocks)-1].padding, i.padding)
+ }
+ }
+}
+
+func printMatrix(m sparseMatrix, t *testing.T) {
+ t.Log("------- matrix -----------")
+ for i := range m.coeff {
+ t.Logf("%v = %v\n", m.coeff[i], m.v[i].data)
+ }
+}
+
+func TestMatrixXorRow(t *testing.T) {
+ var xorRowTests = []struct {
+ arow []int
+ r []int
+ result []int
+ }{
+ {[]int{0, 1}, []int{2, 3}, []int{0, 1, 2, 3}},
+ {[]int{0, 1}, []int{1, 2, 3}, []int{0, 2, 3}},
+ {[]int{}, []int{1, 2, 3}, []int{1, 2, 3}},
+ {[]int{1, 2, 3}, []int{}, []int{1, 2, 3}},
+ {[]int{1}, []int{2}, []int{1, 2}},
+ {[]int{1}, []int{1}, []int{}},
+ {[]int{1, 2}, []int{1, 2, 3, 4}, []int{3, 4}},
+ {[]int{3, 4}, []int{1, 2, 3, 4}, []int{1, 2}},
+ {[]int{1, 2, 3, 4}, []int{1, 2}, []int{3, 4}},
+ {[]int{0, 1, 2, 3, 4}, []int{1, 2}, []int{0, 3, 4}},
+ {[]int{3, 4}, []int{1, 2, 3, 4, 5}, []int{1, 2, 5}},
+ {[]int{3, 4, 8}, []int{1, 2, 3, 4, 5}, []int{1, 2, 5, 8}},
+ }
+
+ for _, test := range xorRowTests {
+ m := sparseMatrix{coeff: [][]int{test.arow}, v: []block{block{[]byte{1}, 0}}}
+
+ testb := block{[]byte{2}, 0}
+ test.r, testb = m.xorRow(0, test.r, testb)
+
+ // Needed since under DeepEqual the nil and the empty slice are not equal.
+ if test.r == nil {
+ test.r = make([]int, 0)
+ }
+ if !reflect.DeepEqual(test.r, test.result) {
+ t.Errorf("XOR row result got %v, should be %v", test.r, test.result)
+ }
+ if !reflect.DeepEqual(testb, block{[]byte{3}, 0}) {
+ t.Errorf("XOR row block got %v, should be %v", testb, block{[]byte{3}, 0})
+ }
+ }
+}
+
+func TestMatrixBasic(t *testing.T) {
+ m := sparseMatrix{coeff: [][]int{{}, {}}, v: []block{block{}, block{}}}
+
+ m.addEquation([]int{0}, block{data: []byte{1}})
+ if m.determined() {
+ t.Errorf("2-row matrix should not be determined after 1 equation")
+ printMatrix(m, t)
+ }
+
+ m.addEquation([]int{0, 1}, block{data: []byte{2}})
+ if !m.determined() {
+ t.Errorf("2-row matrix should be determined after 2 equations")
+ printMatrix(m, t)
+ }
+
+ printMatrix(m, t)
+
+ if !reflect.DeepEqual(m.coeff[0], []int{0}) ||
+ !reflect.DeepEqual(m.v[0].data, []byte{1}) {
+ t.Errorf("Equation 0 got (%v = %v), want ([0] = [1])", m.coeff[0], m.v[0].data)
+ }
+ if !reflect.DeepEqual(m.coeff[1], []int{1}) ||
+ !reflect.DeepEqual(m.v[1].data, []byte{3}) {
+ t.Errorf("Equation 1 got (%v = %v), want ([1] = [3])", m.coeff[0], m.v[0].data)
+ }
+
+ m.reduce()
+ if !reflect.DeepEqual(m.coeff[0], []int{0}) ||
+ !reflect.DeepEqual(m.v[0].data, []byte{1}) {
+ t.Errorf("Equation 0 got (%v = %v), want ([0] = [1])", m.coeff[0], m.v[0].data)
+ }
+ if !reflect.DeepEqual(m.coeff[1], []int{1}) ||
+ !reflect.DeepEqual(m.v[1].data, []byte{3}) {
+ t.Errorf("Equation 1 got (%v = %v), want ([1] = [3])", m.coeff[0], m.v[0].data)
+ }
+}
+
+func TestMatrixLarge(t *testing.T) {
+ m := sparseMatrix{coeff: make([][]int, 4), v: make([]block, 4)}
+
+ m.addEquation([]int{2, 3}, block{data: []byte{1}})
+ m.addEquation([]int{2}, block{data: []byte{2}})
+ if m.determined() {
+ t.Errorf("4-row matrix should not be determined after 2 equations")
+ printMatrix(m, t)
+ }
+ printMatrix(m, t)
+
+ // Should have triangular entries in {2} and {3} now.
+ if len(m.coeff[2]) != 1 || m.v[2].data[0] != 2 {
+ t.Errorf("Equation 2 got %v = %v, should be [2] = [2]", m.coeff[2], m.v[2])
+ }
+ if len(m.coeff[3]) != 1 || m.v[3].data[0] != 3 {
+ t.Errorf("Equation 3 got %v = %v, should be [3] = [3]", m.coeff[3], m.v[3])
+ }
+ if len(m.coeff[0]) != 0 || len(m.coeff[1]) != 0 {
+ t.Errorf("Equations 0 and 1 should be empty")
+ printMatrix(m, t)
+ }
+
+ m.addEquation([]int{0, 1, 2, 3}, block{data: []byte{4}})
+ if m.determined() {
+ t.Errorf("4-row matrix should not be determined after 3 equations")
+ printMatrix(m, t)
+ }
+
+ m.addEquation([]int{3}, block{data: []byte{3}})
+ if m.determined() {
+ t.Errorf("4-row matrix should not be determined after redundant equation")
+ printMatrix(m, t)
+ }
+
+ m.addEquation([]int{0, 2}, block{data: []byte{8}})
+ if !m.determined() {
+ t.Errorf("4-row matrix should be determined after 4 equations")
+ printMatrix(m, t)
+ }
+
+ // The matrix should now have entries in rows 0 and 1, but not equal to the
+ // original equations.
+ printMatrix(m, t)
+ if !reflect.DeepEqual(m.coeff[0], []int{0, 2}) {
+ t.Errorf("Got %v for coeff[0], expect [0, 2]", m.coeff[0])
+ }
+ if !reflect.DeepEqual(m.coeff[1], []int{1, 3}) {
+ t.Errorf("Got %v for coeff[1], expect [1, 3]", m.coeff[1])
+ }
+}
diff --git a/constants.go b/constants.go
new file mode 100644
index 0000000..01e2049
--- /dev/null
+++ b/constants.go
@@ -0,0 +1,722 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+var smallPrimes = []int{
+ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
+ 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109,
+ 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
+ 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
+ 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
+ 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487,
+ 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
+ 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
+ 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
+ 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919,
+ 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021,
+ 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103,
+ 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
+ 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297,
+ 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
+ 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489,
+ 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
+ 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669,
+ 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783,
+ 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879,
+ 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999}
+
+// The V0 table from RFC 5053 Section 5.6.1
+var v0table = [256]uint32{
+ 251291136, 3952231631, 3370958628, 4070167936, 123631495, 3351110283,
+ 3218676425, 2011642291, 774603218, 2402805061, 1004366930,
+ 1843948209, 428891132, 3746331984, 1591258008, 3067016507,
+ 1433388735, 504005498, 2032657933, 3419319784, 2805686246,
+ 3102436986, 3808671154, 2501582075, 3978944421, 246043949,
+ 4016898363, 649743608, 1974987508, 2651273766, 2357956801, 689605112,
+ 715807172, 2722736134, 191939188, 3535520147, 3277019569, 1470435941,
+ 3763101702, 3232409631, 122701163, 3920852693, 782246947, 372121310,
+ 2995604341, 2045698575, 2332962102, 4005368743, 218596347,
+ 3415381967, 4207612806, 861117671, 3676575285, 2581671944,
+ 3312220480, 681232419, 307306866, 4112503940, 1158111502, 709227802,
+ 2724140433, 4201101115, 4215970289, 4048876515, 3031661061,
+ 1909085522, 510985033, 1361682810, 129243379, 3142379587, 2569842483,
+ 3033268270, 1658118006, 932109358, 1982290045, 2983082771,
+ 3007670818, 3448104768, 683749698, 778296777, 1399125101, 1939403708,
+ 1692176003, 3868299200, 1422476658, 593093658, 1878973865,
+ 2526292949, 1591602827, 3986158854, 3964389521, 2695031039,
+ 1942050155, 424618399, 1347204291, 2669179716, 2434425874,
+ 2540801947, 1384069776, 4123580443, 1523670218, 2708475297,
+ 1046771089, 2229796016, 1255426612, 4213663089, 1521339547,
+ 3041843489, 420130494, 10677091, 515623176, 3457502702, 2115821274,
+ 2720124766, 3242576090, 854310108, 425973987, 325832382, 1796851292,
+ 2462744411, 1976681690, 1408671665, 1228817808, 3917210003,
+ 263976645, 2593736473, 2471651269, 4291353919, 650792940, 1191583883,
+ 3046561335, 2466530435, 2545983082, 969168436, 2019348792,
+ 2268075521, 1169345068, 3250240009, 3963499681, 2560755113,
+ 911182396, 760842409, 3569308693, 2687243553, 381854665, 2613828404,
+ 2761078866, 1456668111, 883760091, 3294951678, 1604598575,
+ 1985308198, 1014570543, 2724959607, 3062518035, 3115293053,
+ 138853680, 4160398285, 3322241130, 2068983570, 2247491078,
+ 3669524410, 1575146607, 828029864, 3732001371, 3422026452,
+ 3370954177, 4006626915, 543812220, 1243116171, 3928372514,
+ 2791443445, 4081325272, 2280435605, 885616073, 616452097, 3188863436,
+ 2780382310, 2340014831, 1208439576, 258356309, 3837963200,
+ 2075009450, 3214181212, 3303882142, 880813252, 1355575717, 207231484,
+ 2420803184, 358923368, 1617557768, 3272161958, 1771154147,
+ 2842106362, 1751209208, 1421030790, 658316681, 194065839, 3241510581,
+ 38625260, 301875395, 4176141739, 297312930, 2137802113, 1502984205,
+ 3669376622, 3728477036, 234652930, 2213589897, 2734638932,
+ 1129721478, 3187422815, 2859178611, 3284308411, 3819792700,
+ 3557526733, 451874476, 1740576081, 3592838701, 1709429513,
+ 3702918379, 3533351328, 1641660745, 179350258, 2380520112,
+ 3936163904, 3685256204, 3156252216, 1854258901, 2861641019,
+ 3176611298, 834787554, 331353807, 517858103, 3010168884, 4012642001,
+ 2217188075, 3756943137, 3077882590, 2054995199, 3081443129,
+ 3895398812, 1141097543, 2376261053, 2626898255, 2554703076,
+ 401233789, 1460049922, 678083952, 1064990737, 940909784, 1673396780,
+ 528881783, 1712547446, 3629685652, 1358307511}
+
+// The V1 table from RFC 5053 Section 5.6.2
+var v1table = [256]uint32{
+ 807385413, 2043073223, 3336749796, 1302105833, 2278607931, 541015020,
+ 1684564270, 372709334, 3508252125, 1768346005, 1270451292,
+ 2603029534, 2049387273, 3891424859, 2152948345, 4114760273,
+ 915180310, 3754787998, 700503826, 2131559305, 1308908630, 224437350,
+ 4065424007, 3638665944, 1679385496, 3431345226, 1779595665,
+ 3068494238, 1424062773, 1033448464, 4050396853, 3302235057,
+ 420600373, 2868446243, 311689386, 259047959, 4057180909, 1575367248,
+ 4151214153, 110249784, 3006865921, 4293710613, 3501256572, 998007483,
+ 499288295, 1205710710, 2997199489, 640417429, 3044194711, 486690751,
+ 2686640734, 2394526209, 2521660077, 49993987, 3843885867, 4201106668,
+ 415906198, 19296841, 2402488407, 2137119134, 1744097284, 579965637,
+ 2037662632, 852173610, 2681403713, 1047144830, 2982173936, 910285038,
+ 4187576520, 2589870048, 989448887, 3292758024, 506322719, 176010738,
+ 1865471968, 2619324712, 564829442, 1996870325, 339697593, 4071072948,
+ 3618966336, 2111320126, 1093955153, 957978696, 892010560, 1854601078,
+ 1873407527, 2498544695, 2694156259, 1927339682, 1650555729,
+ 183933047, 3061444337, 2067387204, 228962564, 3904109414, 1595995433,
+ 1780701372, 2463145963, 307281463, 3237929991, 3852995239,
+ 2398693510, 3754138664, 522074127, 146352474, 4104915256, 3029415884,
+ 3545667983, 332038910, 976628269, 3123492423, 3041418372, 2258059298,
+ 2139377204, 3243642973, 3226247917, 3674004636, 2698992189,
+ 3453843574, 1963216666, 3509855005, 2358481858, 747331248,
+ 1957348676, 1097574450, 2435697214, 3870972145, 1888833893,
+ 2914085525, 4161315584, 1273113343, 3269644828, 3681293816,
+ 412536684, 1156034077, 3823026442, 1066971017, 3598330293,
+ 1979273937, 2079029895, 1195045909, 1071986421, 2712821515,
+ 3377754595, 2184151095, 750918864, 2585729879, 4249895712,
+ 1832579367, 1192240192, 946734366, 31230688, 3174399083, 3549375728,
+ 1642430184, 1904857554, 861877404, 3277825584, 4267074718,
+ 3122860549, 666423581, 644189126, 226475395, 307789415, 1196105631,
+ 3191691839, 782852669, 1608507813, 1847685900, 4069766876,
+ 3931548641, 2526471011, 766865139, 2115084288, 4259411376,
+ 3323683436, 568512177, 3736601419, 1800276898, 4012458395, 1823982,
+ 27980198, 2023839966, 869505096, 431161506, 1024804023, 1853869307,
+ 3393537983, 1500703614, 3019471560, 1351086955, 3096933631,
+ 3034634988, 2544598006, 1230942551, 3362230798, 159984793, 491590373,
+ 3993872886, 3681855622, 903593547, 3535062472, 1799803217, 772984149,
+ 895863112, 1899036275, 4187322100, 101856048, 234650315, 3183125617,
+ 3190039692, 525584357, 1286834489, 455810374, 1869181575, 922673938,
+ 3877430102, 3422391938, 1414347295, 1971054608, 3061798054,
+ 830555096, 2822905141, 167033190, 1079139428, 4210126723, 3593797804,
+ 429192890, 372093950, 1779187770, 3312189287, 204349348, 452421568,
+ 2800540462, 3733109044, 1235082423, 1765319556, 3174729780,
+ 3762994475, 3171962488, 442160826, 198349622, 45942637, 1324086311,
+ 2901868599, 678860040, 3812229107, 19936821, 1119590141, 3640121682,
+ 3545931032, 2102949142, 2828208598, 3603378023, 4135048896}
+
+// The systematic indices table from RFC 5053 Section 5.7
+var systematicIndextable = [8193]uint16{0, 0, 0, 0,
+ 18, 14, 61, 46, 14, 22, 20, 40, 48, 1, 29, 40, 43, 46, 18, 8, 20, 2,
+ 61, 26, 13, 29, 36, 19, 58, 5, 58, 0, 54, 56, 24, 14, 5, 67, 39, 31,
+ 25, 29, 24, 19, 14, 56, 49, 49, 63, 30, 4, 39, 2, 1, 20, 19, 61, 4,
+ 54, 70, 25, 52, 9, 26, 55, 69, 27, 68, 75, 19, 64, 57, 45, 3, 37, 31,
+ 100, 41, 25, 41, 53, 23, 9, 31, 26, 30, 30, 46, 90, 50, 13, 90, 77,
+ 61, 31, 54, 54, 3, 21, 66, 21, 11, 23, 11, 29, 21, 7, 1, 27, 4, 34,
+ 17, 85, 69, 17, 75, 93, 57, 0, 53, 71, 88, 119, 88, 90, 22, 0, 58,
+ 41, 22, 96, 26, 79, 118, 19, 3, 81, 72, 50, 0, 32, 79, 28, 25, 12,
+ 25, 29, 3, 37, 30, 30, 41, 84, 32, 31, 61, 32, 61, 7, 56, 54, 39, 33,
+ 66, 29, 3, 14, 75, 75, 78, 84, 75, 84, 25, 54, 25, 25, 107, 78, 27,
+ 73, 0, 49, 96, 53, 50, 21, 10, 73, 58, 65, 27, 3, 27, 18, 54, 45, 69,
+ 29, 3, 65, 31, 71, 76, 56, 54, 76, 54, 13, 5, 18, 142, 17, 3, 37,
+ 114, 41, 25, 56, 0, 23, 3, 41, 22, 22, 31, 18, 48, 31, 58, 37, 75,
+ 88, 3, 56, 1, 95, 19, 73, 52, 52, 4, 75, 26, 1, 25, 10, 1, 70, 31,
+ 31, 12, 10, 54, 46, 11, 74, 84, 74, 8, 58, 23, 74, 8, 36, 11, 16, 94,
+ 76, 14, 57, 65, 8, 22, 10, 36, 36, 96, 62, 103, 6, 75, 103, 58, 10,
+ 15, 41, 75, 125, 58, 15, 10, 34, 29, 34, 4, 16, 29, 18, 18, 28, 71,
+ 28, 43, 77, 18, 41, 41, 41, 62, 29, 96, 15, 106, 43, 15, 3, 43, 61,
+ 3, 18, 103, 77, 29, 103, 19, 58, 84, 58, 1, 146, 32, 3, 70, 52, 54,
+ 29, 70, 69, 124, 62, 1, 26, 38, 26, 3, 16, 26, 5, 51, 120, 41, 16, 1,
+ 43, 34, 34, 29, 37, 56, 29, 96, 86, 54, 25, 84, 50, 34, 34, 93, 84,
+ 96, 29, 29, 50, 50, 6, 1, 105, 78, 15, 37, 19, 50, 71, 36, 6, 54, 8,
+ 28, 54, 75, 75, 16, 75, 131, 5, 25, 16, 69, 17, 69, 6, 96, 53, 96,
+ 41, 119, 6, 6, 88, 50, 88, 52, 37, 0, 124, 73, 73, 7, 14, 36, 69, 79,
+ 6, 114, 40, 79, 17, 77, 24, 44, 37, 69, 27, 37, 29, 33, 37, 50, 31,
+ 69, 29, 101, 7, 61, 45, 17, 73, 37, 34, 18, 94, 22, 22, 63, 3, 25,
+ 25, 17, 3, 90, 34, 34, 41, 34, 41, 54, 41, 54, 41, 41, 41, 163, 143,
+ 96, 18, 32, 39, 86, 104, 11, 17, 17, 11, 86, 104, 78, 70, 52, 78, 17,
+ 73, 91, 62, 7, 128, 50, 124, 18, 101, 46, 10, 75, 104, 73, 58, 132,
+ 34, 13, 4, 95, 88, 33, 76, 74, 54, 62, 113, 114, 103, 32, 103, 69,
+ 54, 53, 3, 11, 72, 31, 53, 102, 37, 53, 11, 81, 41, 10, 164, 10, 41,
+ 31, 36, 113, 82, 3, 125, 62, 16, 4, 41, 41, 4, 128, 49, 138, 128, 74,
+ 103, 0, 6, 101, 41, 142, 171, 39, 105, 121, 81, 62, 41, 81, 37, 3,
+ 81, 69, 62, 3, 69, 70, 21, 29, 4, 91, 87, 37, 79, 36, 21, 71, 37, 41,
+ 75, 128, 128, 15, 25, 3, 108, 73, 91, 62, 114, 62, 62, 36, 36, 15,
+ 58, 114, 61, 114, 58, 105, 114, 41, 61, 176, 145, 46, 37, 30, 220,
+ 77, 138, 15, 1, 128, 53, 50, 50, 58, 8, 91, 114, 105, 63, 91, 37, 37,
+ 13, 169, 51, 102, 6, 102, 23, 105, 23, 58, 6, 29, 29, 19, 82, 29, 13,
+ 36, 27, 29, 61, 12, 18, 127, 127, 12, 44, 102, 18, 4, 15, 206, 53,
+ 127, 53, 17, 69, 69, 69, 29, 29, 109, 25, 102, 25, 53, 62, 99, 62,
+ 62, 29, 62, 62, 45, 91, 125, 29, 29, 29, 4, 117, 72, 4, 30, 71, 71,
+ 95, 79, 179, 71, 30, 53, 32, 32, 49, 25, 91, 25, 26, 26, 103, 123,
+ 26, 41, 162, 78, 52, 103, 25, 6, 142, 94, 45, 45, 94, 127, 94, 94,
+ 94, 47, 209, 138, 39, 39, 19, 154, 73, 67, 91, 27, 91, 84, 4, 84, 91,
+ 12, 14, 165, 142, 54, 69, 192, 157, 185, 8, 95, 25, 62, 103, 103, 95,
+ 71, 97, 62, 128, 0, 29, 51, 16, 94, 16, 16, 51, 0, 29, 85, 10, 105,
+ 16, 29, 29, 13, 29, 4, 4, 132, 23, 95, 25, 54, 41, 29, 50, 70, 58,
+ 142, 72, 70, 15, 72, 54, 29, 22, 145, 29, 127, 29, 85, 58, 101, 34,
+ 165, 91, 46, 46, 25, 185, 25, 77, 128, 46, 128, 46, 188, 114, 46, 25,
+ 45, 45, 114, 145, 114, 15, 102, 142, 8, 73, 31, 139, 157, 13, 79, 13,
+ 114, 150, 8, 90, 91, 123, 69, 82, 132, 8, 18, 10, 102, 103, 114, 103,
+ 8, 103, 13, 115, 55, 62, 3, 8, 154, 114, 99, 19, 8, 31, 73, 19, 99,
+ 10, 6, 121, 32, 13, 32, 119, 32, 29, 145, 30, 13, 13, 114, 145, 32,
+ 1, 123, 39, 29, 31, 69, 31, 140, 72, 72, 25, 25, 123, 25, 123, 8, 4,
+ 85, 8, 25, 39, 25, 39, 85, 138, 25, 138, 25, 33, 102, 70, 25, 25, 31,
+ 25, 25, 192, 69, 69, 114, 145, 120, 120, 8, 33, 98, 15, 212, 155, 8,
+ 101, 8, 8, 98, 68, 155, 102, 132, 120, 30, 25, 123, 123, 101, 25,
+ 123, 32, 24, 94, 145, 32, 24, 94, 118, 145, 101, 53, 53, 25, 128,
+ 173, 142, 81, 81, 69, 33, 33, 125, 4, 1, 17, 27, 4, 17, 102, 27, 13,
+ 25, 128, 71, 13, 39, 53, 13, 53, 47, 39, 23, 128, 53, 39, 47, 39,
+ 135, 158, 136, 36, 36, 27, 157, 47, 76, 213, 47, 156, 25, 25, 53, 25,
+ 53, 25, 86, 27, 159, 25, 62, 79, 39, 79, 25, 145, 49, 25, 143, 13,
+ 114, 150, 130, 94, 102, 39, 4, 39, 61, 77, 228, 22, 25, 47, 119, 205,
+ 122, 119, 205, 119, 22, 119, 258, 143, 22, 81, 179, 22, 22, 143, 25,
+ 65, 53, 168, 36, 79, 175, 37, 79, 70, 79, 103, 70, 25, 175, 4, 96,
+ 96, 49, 128, 138, 96, 22, 62, 47, 95, 105, 95, 62, 95, 62, 142, 103,
+ 69, 103, 30, 103, 34, 173, 127, 70, 127, 132, 18, 85, 22, 71, 18,
+ 206, 206, 18, 128, 145, 70, 193, 188, 8, 125, 114, 70, 128, 114, 145,
+ 102, 25, 12, 108, 102, 94, 10, 102, 1, 102, 124, 22, 22, 118, 132,
+ 22, 116, 75, 41, 63, 41, 189, 208, 55, 85, 69, 8, 71, 53, 71, 69,
+ 102, 165, 41, 99, 69, 33, 33, 29, 156, 102, 13, 251, 102, 25, 13,
+ 109, 102, 164, 102, 164, 102, 25, 29, 228, 29, 259, 179, 222, 95, 94,
+ 30, 30, 30, 142, 55, 142, 72, 55, 102, 128, 17, 69, 164, 165, 3, 164,
+ 36, 165, 27, 27, 45, 21, 21, 237, 113, 83, 231, 106, 13, 154, 13,
+ 154, 128, 154, 148, 258, 25, 154, 128, 3, 27, 10, 145, 145, 21, 146,
+ 25, 1, 185, 121, 0, 1, 95, 55, 95, 95, 30, 0, 27, 95, 0, 95, 8, 222,
+ 27, 121, 30, 95, 121, 0, 98, 94, 131, 55, 95, 95, 30, 98, 30, 0, 91,
+ 145, 66, 179, 66, 58, 175, 29, 0, 31, 173, 146, 160, 39, 53, 28, 123,
+ 199, 123, 175, 146, 156, 54, 54, 149, 25, 70, 178, 128, 25, 70, 70,
+ 94, 224, 54, 4, 54, 54, 25, 228, 160, 206, 165, 143, 206, 108, 220,
+ 234, 160, 13, 169, 103, 103, 103, 91, 213, 222, 91, 103, 91, 103, 31,
+ 30, 123, 13, 62, 103, 50, 106, 42, 13, 145, 114, 220, 65, 8, 8, 175,
+ 11, 104, 94, 118, 132, 27, 118, 193, 27, 128, 127, 127, 183, 33, 30,
+ 29, 103, 128, 61, 234, 165, 41, 29, 193, 33, 207, 41, 165, 165, 55,
+ 81, 157, 157, 8, 81, 11, 27, 8, 8, 98, 96, 142, 145, 41, 179, 112,
+ 62, 180, 206, 206, 165, 39, 241, 45, 151, 26, 197, 102, 192, 125,
+ 128, 67, 128, 69, 128, 197, 33, 125, 102, 13, 103, 25, 30, 12, 30,
+ 12, 30, 25, 77, 12, 25, 180, 27, 10, 69, 235, 228, 343, 118, 69, 41,
+ 8, 69, 175, 25, 69, 25, 125, 41, 25, 41, 8, 155, 146, 155, 146, 155,
+ 206, 168, 128, 157, 27, 273, 211, 211, 168, 11, 173, 154, 77, 173,
+ 77, 102, 102, 102, 8, 85, 95, 102, 157, 28, 122, 234, 122, 157, 235,
+ 222, 241, 10, 91, 179, 25, 13, 25, 41, 25, 206, 41, 6, 41, 158, 206,
+ 206, 33, 296, 296, 33, 228, 69, 8, 114, 148, 33, 29, 66, 27, 27, 30,
+ 233, 54, 173, 108, 106, 108, 108, 53, 103, 33, 33, 33, 176, 27, 27,
+ 205, 164, 105, 237, 41, 27, 72, 165, 29, 29, 259, 132, 132, 132, 364,
+ 71, 71, 27, 94, 160, 127, 51, 234, 55, 27, 95, 94, 165, 55, 55, 41,
+ 0, 41, 128, 4, 123, 173, 6, 164, 157, 121, 121, 154, 86, 164, 164,
+ 25, 93, 164, 25, 164, 210, 284, 62, 93, 30, 25, 25, 30, 30, 260, 130,
+ 25, 125, 57, 53, 166, 166, 166, 185, 166, 158, 94, 113, 215, 159, 62,
+ 99, 21, 172, 99, 184, 62, 259, 4, 21, 21, 77, 62, 173, 41, 146, 6,
+ 41, 128, 121, 41, 11, 121, 103, 159, 164, 175, 206, 91, 103, 164, 72,
+ 25, 129, 72, 206, 129, 33, 103, 102, 102, 29, 13, 11, 251, 234, 135,
+ 31, 8, 123, 65, 91, 121, 129, 65, 243, 10, 91, 8, 65, 70, 228, 220,
+ 243, 91, 10, 10, 30, 178, 91, 178, 33, 21, 25, 235, 165, 11, 161,
+ 158, 27, 27, 30, 128, 75, 36, 30, 36, 36, 173, 25, 33, 178, 112, 162,
+ 112, 112, 112, 162, 33, 33, 178, 123, 123, 39, 106, 91, 106, 106,
+ 158, 106, 106, 284, 39, 230, 21, 228, 11, 21, 228, 159, 241, 62, 10,
+ 62, 10, 68, 234, 39, 39, 138, 62, 22, 27, 183, 22, 215, 10, 175, 175,
+ 353, 228, 42, 193, 175, 175, 27, 98, 27, 193, 150, 27, 173, 17, 233,
+ 233, 25, 102, 123, 152, 242, 108, 4, 94, 176, 13, 41, 219, 17, 151,
+ 22, 103, 103, 53, 128, 233, 284, 25, 265, 128, 39, 39, 138, 42, 39,
+ 21, 86, 95, 127, 29, 91, 46, 103, 103, 215, 25, 123, 123, 230, 25,
+ 193, 180, 30, 60, 30, 242, 136, 180, 193, 30, 206, 180, 60, 165, 206,
+ 193, 165, 123, 164, 103, 68, 25, 70, 91, 25, 82, 53, 82, 186, 53, 82,
+ 53, 25, 30, 282, 91, 13, 234, 160, 160, 126, 149, 36, 36, 160, 149,
+ 178, 160, 39, 294, 149, 149, 160, 39, 95, 221, 186, 106, 178, 316,
+ 267, 53, 53, 164, 159, 164, 165, 94, 228, 53, 52, 178, 183, 53, 294,
+ 128, 55, 140, 294, 25, 95, 366, 15, 304, 13, 183, 77, 230, 6, 136,
+ 235, 121, 311, 273, 36, 158, 235, 230, 98, 201, 165, 165, 165, 91,
+ 175, 248, 39, 185, 128, 39, 39, 128, 313, 91, 36, 219, 130, 25, 130,
+ 234, 234, 130, 234, 121, 205, 304, 94, 77, 64, 259, 60, 60, 60, 77,
+ 242, 60, 145, 95, 270, 18, 91, 199, 159, 91, 235, 58, 249, 26, 123,
+ 114, 29, 15, 191, 15, 30, 55, 55, 347, 4, 29, 15, 4, 341, 93, 7, 30,
+ 23, 7, 121, 266, 178, 261, 70, 169, 25, 25, 158, 169, 25, 169, 270,
+ 270, 13, 128, 327, 103, 55, 128, 103, 136, 159, 103, 327, 41, 32,
+ 111, 111, 114, 173, 215, 173, 25, 173, 180, 114, 173, 173, 98, 93,
+ 25, 160, 157, 159, 160, 159, 159, 160, 320, 35, 193, 221, 33, 36,
+ 136, 248, 91, 215, 125, 215, 156, 68, 125, 125, 1, 287, 123, 94, 30,
+ 184, 13, 30, 94, 123, 206, 12, 206, 289, 128, 122, 184, 128, 289,
+ 178, 29, 26, 206, 178, 65, 206, 128, 192, 102, 197, 36, 94, 94, 155,
+ 10, 36, 121, 280, 121, 368, 192, 121, 121, 179, 121, 36, 54, 192,
+ 121, 192, 197, 118, 123, 224, 118, 10, 192, 10, 91, 269, 91, 49, 206,
+ 184, 185, 62, 8, 49, 289, 30, 5, 55, 30, 42, 39, 220, 298, 42, 347,
+ 42, 234, 42, 70, 42, 55, 321, 129, 172, 173, 172, 13, 98, 129, 325,
+ 235, 284, 362, 129, 233, 345, 175, 261, 175, 60, 261, 58, 289, 99,
+ 99, 99, 206, 99, 36, 175, 29, 25, 432, 125, 264, 168, 173, 69, 158,
+ 273, 179, 164, 69, 158, 69, 8, 95, 192, 30, 164, 101, 44, 53, 273,
+ 335, 273, 53, 45, 128, 45, 234, 123, 105, 103, 103, 224, 36, 90, 211,
+ 282, 264, 91, 228, 91, 166, 264, 228, 398, 50, 101, 91, 264, 73, 36,
+ 25, 73, 50, 50, 242, 36, 36, 58, 165, 204, 353, 165, 125, 320, 128,
+ 298, 298, 180, 128, 60, 102, 30, 30, 53, 179, 234, 325, 234, 175, 21,
+ 250, 215, 103, 21, 21, 250, 91, 211, 91, 313, 301, 323, 215, 228,
+ 160, 29, 29, 81, 53, 180, 146, 248, 66, 159, 39, 98, 323, 98, 36, 95,
+ 218, 234, 39, 82, 82, 230, 62, 13, 62, 230, 13, 30, 98, 0, 8, 98, 8,
+ 98, 91, 267, 121, 197, 30, 78, 27, 78, 102, 27, 298, 160, 103, 264,
+ 264, 264, 175, 17, 273, 273, 165, 31, 160, 17, 99, 17, 99, 234, 31,
+ 17, 99, 36, 26, 128, 29, 214, 353, 264, 102, 36, 102, 264, 264, 273,
+ 273, 4, 16, 138, 138, 264, 128, 313, 25, 420, 60, 10, 280, 264, 60,
+ 60, 103, 178, 125, 178, 29, 327, 29, 36, 30, 36, 4, 52, 183, 183,
+ 173, 52, 31, 173, 31, 158, 31, 158, 31, 9, 31, 31, 353, 31, 353, 173,
+ 415, 9, 17, 222, 31, 103, 31, 165, 27, 31, 31, 165, 27, 27, 206, 31,
+ 31, 4, 4, 30, 4, 4, 264, 185, 159, 310, 273, 310, 173, 40, 4, 173, 4,
+ 173, 4, 250, 250, 62, 188, 119, 250, 233, 62, 121, 105, 105, 54, 103,
+ 111, 291, 236, 236, 103, 297, 36, 26, 316, 69, 183, 158, 206, 129,
+ 160, 129, 184, 55, 179, 279, 11, 179, 347, 160, 184, 129, 179, 351,
+ 179, 353, 179, 129, 129, 351, 11, 111, 93, 93, 235, 103, 173, 53, 93,
+ 50, 111, 86, 123, 94, 36, 183, 60, 55, 55, 178, 219, 253, 321, 178,
+ 235, 235, 183, 183, 204, 321, 219, 160, 193, 335, 121, 70, 69, 295,
+ 159, 297, 231, 121, 231, 136, 353, 136, 121, 279, 215, 366, 215, 353,
+ 159, 353, 353, 103, 31, 31, 298, 298, 30, 30, 165, 273, 25, 219, 35,
+ 165, 259, 54, 36, 54, 54, 165, 71, 250, 327, 13, 289, 165, 196, 165,
+ 165, 94, 233, 165, 94, 60, 165, 96, 220, 166, 271, 158, 397, 122, 53,
+ 53, 137, 280, 272, 62, 30, 30, 30, 105, 102, 67, 140, 8, 67, 21, 270,
+ 298, 69, 173, 298, 91, 179, 327, 86, 179, 88, 179, 179, 55, 123, 220,
+ 233, 94, 94, 175, 13, 53, 13, 154, 191, 74, 83, 83, 325, 207, 83, 74,
+ 83, 325, 74, 316, 388, 55, 55, 364, 55, 183, 434, 273, 273, 273, 164,
+ 213, 11, 213, 327, 321, 21, 352, 185, 103, 13, 13, 55, 30, 323, 123,
+ 178, 435, 178, 30, 175, 175, 30, 481, 527, 175, 125, 232, 306, 232,
+ 206, 306, 364, 206, 270, 206, 232, 10, 30, 130, 160, 130, 347, 240,
+ 30, 136, 130, 347, 136, 279, 298, 206, 30, 103, 273, 241, 70, 206,
+ 306, 434, 206, 94, 94, 156, 161, 321, 321, 64, 161, 13, 183, 183, 83,
+ 161, 13, 169, 13, 159, 36, 173, 159, 36, 36, 230, 235, 235, 159, 159,
+ 335, 312, 42, 342, 264, 39, 39, 39, 34, 298, 36, 36, 252, 164, 29,
+ 493, 29, 387, 387, 435, 493, 132, 273, 105, 132, 74, 73, 206, 234,
+ 273, 206, 95, 15, 280, 280, 280, 280, 397, 273, 273, 242, 397, 280,
+ 397, 397, 397, 273, 397, 280, 230, 137, 353, 67, 81, 137, 137, 353,
+ 259, 312, 114, 164, 164, 25, 77, 21, 77, 165, 30, 30, 231, 234, 121,
+ 234, 312, 121, 364, 136, 123, 123, 136, 123, 136, 150, 264, 285, 30,
+ 166, 93, 30, 39, 224, 136, 39, 355, 355, 397, 67, 67, 25, 67, 25,
+ 298, 11, 67, 264, 374, 99, 150, 321, 67, 70, 67, 295, 150, 29, 321,
+ 150, 70, 29, 142, 355, 311, 173, 13, 253, 103, 114, 114, 70, 192, 22,
+ 128, 128, 183, 184, 70, 77, 215, 102, 292, 30, 123, 279, 292, 142,
+ 33, 215, 102, 468, 123, 468, 473, 30, 292, 215, 30, 213, 443, 473,
+ 215, 234, 279, 279, 279, 279, 265, 443, 206, 66, 313, 34, 30, 206,
+ 30, 51, 15, 206, 41, 434, 41, 398, 67, 30, 301, 67, 36, 3, 285, 437,
+ 136, 136, 22, 136, 145, 365, 323, 323, 145, 136, 22, 453, 99, 323,
+ 353, 9, 258, 323, 231, 128, 231, 382, 150, 420, 39, 94, 29, 29, 353,
+ 22, 22, 347, 353, 39, 29, 22, 183, 8, 284, 355, 388, 284, 60, 64, 99,
+ 60, 64, 150, 95, 150, 364, 150, 95, 150, 6, 236, 383, 544, 81, 206,
+ 388, 206, 58, 159, 99, 231, 228, 363, 363, 121, 99, 121, 121, 99,
+ 422, 544, 273, 173, 121, 427, 102, 121, 235, 284, 179, 25, 197, 25,
+ 179, 511, 70, 368, 70, 25, 388, 123, 368, 159, 213, 410, 159, 236,
+ 127, 159, 21, 373, 184, 424, 327, 250, 176, 176, 175, 284, 316, 176,
+ 284, 327, 111, 250, 284, 175, 175, 264, 111, 176, 219, 111, 427, 427,
+ 176, 284, 427, 353, 428, 55, 184, 493, 158, 136, 99, 287, 264, 334,
+ 264, 213, 213, 292, 481, 93, 264, 292, 295, 295, 6, 367, 279, 173,
+ 308, 285, 158, 308, 335, 299, 137, 137, 572, 41, 137, 137, 41, 94,
+ 335, 220, 36, 224, 420, 36, 265, 265, 91, 91, 71, 123, 264, 91, 91,
+ 123, 107, 30, 22, 292, 35, 241, 356, 298, 14, 298, 441, 35, 121, 71,
+ 63, 130, 63, 488, 363, 71, 63, 307, 194, 71, 71, 220, 121, 125, 71,
+ 220, 71, 71, 71, 71, 235, 265, 353, 128, 155, 128, 420, 400, 130,
+ 173, 183, 183, 184, 130, 173, 183, 13, 183, 130, 130, 183, 183, 353,
+ 353, 183, 242, 183, 183, 306, 324, 324, 321, 306, 321, 6, 6, 128,
+ 306, 242, 242, 306, 183, 183, 6, 183, 321, 486, 183, 164, 30, 78,
+ 138, 158, 138, 34, 206, 362, 55, 70, 67, 21, 375, 136, 298, 81, 298,
+ 298, 298, 230, 121, 30, 230, 311, 240, 311, 311, 158, 204, 136, 136,
+ 184, 136, 264, 311, 311, 312, 312, 72, 311, 175, 264, 91, 175, 264,
+ 121, 461, 312, 312, 238, 475, 350, 512, 350, 312, 313, 350, 312, 366,
+ 294, 30, 253, 253, 253, 388, 158, 388, 22, 388, 22, 388, 103, 321,
+ 321, 253, 7, 437, 103, 114, 242, 114, 114, 242, 114, 114, 242, 242,
+ 242, 306, 242, 114, 7, 353, 335, 27, 241, 299, 312, 364, 506, 409,
+ 94, 462, 230, 462, 243, 230, 175, 175, 462, 461, 230, 428, 426, 175,
+ 175, 165, 175, 175, 372, 183, 572, 102, 85, 102, 538, 206, 376, 85,
+ 85, 284, 85, 85, 284, 398, 83, 160, 265, 308, 398, 310, 583, 289,
+ 279, 273, 285, 490, 490, 211, 292, 292, 158, 398, 30, 220, 169, 368,
+ 368, 368, 169, 159, 368, 93, 368, 368, 93, 169, 368, 368, 443, 368,
+ 298, 443, 368, 298, 538, 345, 345, 311, 178, 54, 311, 215, 178, 175,
+ 222, 264, 475, 264, 264, 475, 478, 289, 63, 236, 63, 299, 231, 296,
+ 397, 299, 158, 36, 164, 164, 21, 492, 21, 164, 21, 164, 403, 26, 26,
+ 588, 179, 234, 169, 465, 295, 67, 41, 353, 295, 538, 161, 185, 306,
+ 323, 68, 420, 323, 82, 241, 241, 36, 53, 493, 301, 292, 241, 250, 63,
+ 63, 103, 442, 353, 185, 353, 321, 353, 185, 353, 353, 185, 409, 353,
+ 589, 34, 271, 271, 34, 86, 34, 34, 353, 353, 39, 414, 4, 95, 95, 4,
+ 225, 95, 4, 121, 30, 552, 136, 159, 159, 514, 159, 159, 54, 514, 206,
+ 136, 206, 159, 74, 235, 235, 312, 54, 312, 42, 156, 422, 629, 54,
+ 465, 265, 165, 250, 35, 165, 175, 659, 175, 175, 8, 8, 8, 8, 206,
+ 206, 206, 50, 435, 206, 432, 230, 230, 234, 230, 94, 299, 299, 285,
+ 184, 41, 93, 299, 299, 285, 41, 285, 158, 285, 206, 299, 41, 36, 396,
+ 364, 364, 120, 396, 514, 91, 382, 538, 807, 717, 22, 93, 412, 54,
+ 215, 54, 298, 308, 148, 298, 148, 298, 308, 102, 656, 6, 148, 745,
+ 128, 298, 64, 407, 273, 41, 172, 64, 234, 250, 398, 181, 445, 95,
+ 236, 441, 477, 504, 102, 196, 137, 364, 60, 453, 137, 364, 367, 334,
+ 364, 299, 196, 397, 630, 589, 589, 196, 646, 337, 235, 128, 128, 343,
+ 289, 235, 324, 427, 324, 58, 215, 215, 461, 425, 461, 387, 440, 285,
+ 440, 440, 285, 387, 632, 325, 325, 440, 461, 425, 425, 387, 627, 191,
+ 285, 440, 308, 55, 219, 280, 308, 265, 538, 183, 121, 30, 236, 206,
+ 30, 455, 236, 30, 30, 705, 83, 228, 280, 468, 132, 8, 132, 132, 128,
+ 409, 173, 353, 132, 409, 35, 128, 450, 137, 398, 67, 432, 423, 235,
+ 235, 388, 306, 93, 93, 452, 300, 190, 13, 452, 388, 30, 452, 13, 30,
+ 13, 30, 306, 362, 234, 721, 635, 809, 784, 67, 498, 498, 67, 353,
+ 635, 67, 183, 159, 445, 285, 183, 53, 183, 445, 265, 432, 57, 420,
+ 432, 420, 477, 327, 55, 60, 105, 183, 218, 104, 104, 475, 239, 582,
+ 151, 239, 104, 732, 41, 26, 784, 86, 300, 215, 36, 64, 86, 86, 675,
+ 294, 64, 86, 528, 550, 493, 565, 298, 230, 312, 295, 538, 298, 295,
+ 230, 54, 374, 516, 441, 54, 54, 323, 401, 401, 382, 159, 837, 159,
+ 54, 401, 592, 159, 401, 417, 610, 264, 150, 323, 452, 185, 323, 323,
+ 185, 403, 185, 423, 165, 425, 219, 407, 270, 231, 99, 93, 231, 631,
+ 756, 71, 364, 434, 213, 86, 102, 434, 102, 86, 23, 71, 335, 164, 323,
+ 409, 381, 4, 124, 41, 424, 206, 41, 124, 41, 41, 703, 635, 124, 493,
+ 41, 41, 487, 492, 124, 175, 124, 261, 600, 488, 261, 488, 261, 206,
+ 677, 261, 308, 723, 908, 704, 691, 723, 488, 488, 441, 136, 476, 312,
+ 136, 550, 572, 728, 550, 22, 312, 312, 22, 55, 413, 183, 280, 593,
+ 191, 36, 36, 427, 36, 695, 592, 19, 544, 13, 468, 13, 544, 72, 437,
+ 321, 266, 461, 266, 441, 230, 409, 93, 521, 521, 345, 235, 22, 142,
+ 150, 102, 569, 235, 264, 91, 521, 264, 7, 102, 7, 498, 521, 235, 537,
+ 235, 6, 241, 420, 420, 631, 41, 527, 103, 67, 337, 62, 264, 527, 131,
+ 67, 174, 263, 264, 36, 36, 263, 581, 253, 465, 160, 286, 91, 160, 55,
+ 4, 4, 631, 631, 608, 365, 465, 294, 427, 427, 335, 669, 669, 129, 93,
+ 93, 93, 93, 74, 66, 758, 504, 347, 130, 505, 504, 143, 505, 550, 222,
+ 13, 352, 529, 291, 538, 50, 68, 269, 130, 295, 130, 511, 295, 295,
+ 130, 486, 132, 61, 206, 185, 368, 669, 22, 175, 492, 207, 373, 452,
+ 432, 327, 89, 550, 496, 611, 527, 89, 527, 496, 550, 516, 516, 91,
+ 136, 538, 264, 264, 124, 264, 264, 264, 264, 264, 535, 264, 150, 285,
+ 398, 285, 582, 398, 475, 81, 694, 694, 64, 81, 694, 234, 607, 723,
+ 513, 234, 64, 581, 64, 124, 64, 607, 234, 723, 717, 367, 64, 513,
+ 607, 488, 183, 488, 450, 183, 550, 286, 183, 363, 286, 414, 67, 449,
+ 449, 366, 215, 235, 95, 295, 295, 41, 335, 21, 445, 225, 21, 295,
+ 372, 749, 461, 53, 481, 397, 427, 427, 427, 714, 481, 714, 427, 717,
+ 165, 245, 486, 415, 245, 415, 486, 274, 415, 441, 456, 300, 548, 300,
+ 422, 422, 757, 11, 74, 430, 430, 136, 409, 430, 749, 191, 819, 592,
+ 136, 364, 465, 231, 231, 918, 160, 589, 160, 160, 465, 465, 231, 157,
+ 538, 538, 259, 538, 326, 22, 22, 22, 179, 22, 22, 550, 179, 287, 287,
+ 417, 327, 498, 498, 287, 488, 327, 538, 488, 583, 488, 287, 335, 287,
+ 335, 287, 41, 287, 335, 287, 327, 441, 335, 287, 488, 538, 327, 498,
+ 8, 8, 374, 8, 64, 427, 8, 374, 417, 760, 409, 373, 160, 423, 206,
+ 160, 106, 499, 160, 271, 235, 160, 590, 353, 695, 478, 619, 590, 353,
+ 13, 63, 189, 420, 605, 427, 643, 121, 280, 415, 121, 415, 595, 417,
+ 121, 398, 55, 330, 463, 463, 123, 353, 330, 582, 309, 582, 582, 405,
+ 330, 550, 405, 582, 353, 309, 308, 60, 353, 7, 60, 71, 353, 189, 183,
+ 183, 183, 582, 755, 189, 437, 287, 189, 183, 668, 481, 384, 384, 481,
+ 481, 481, 477, 582, 582, 499, 650, 481, 121, 461, 231, 36, 235, 36,
+ 413, 235, 209, 36, 689, 114, 353, 353, 235, 592, 36, 353, 413, 209,
+ 70, 308, 70, 699, 308, 70, 213, 292, 86, 689, 465, 55, 508, 128, 452,
+ 29, 41, 681, 573, 352, 21, 21, 648, 648, 69, 509, 409, 21, 264, 21,
+ 509, 514, 514, 409, 21, 264, 443, 443, 427, 160, 433, 663, 433, 231,
+ 646, 185, 482, 646, 433, 13, 398, 172, 234, 42, 491, 172, 234, 234,
+ 832, 775, 172, 196, 335, 822, 461, 298, 461, 364, 1120, 537, 169,
+ 169, 364, 694, 219, 612, 231, 740, 42, 235, 321, 279, 960, 279, 353,
+ 492, 159, 572, 321, 159, 287, 353, 287, 287, 206, 206, 321, 287, 159,
+ 321, 492, 159, 55, 572, 600, 270, 492, 784, 173, 91, 91, 443, 443,
+ 582, 261, 497, 572, 91, 555, 352, 206, 261, 555, 285, 91, 555, 497,
+ 83, 91, 619, 353, 488, 112, 4, 592, 295, 295, 488, 235, 231, 769,
+ 568, 581, 671, 451, 451, 483, 299, 1011, 432, 422, 207, 106, 701,
+ 508, 555, 508, 555, 125, 870, 555, 589, 508, 125, 749, 482, 125, 125,
+ 130, 544, 643, 643, 544, 488, 22, 643, 130, 335, 544, 22, 130, 544,
+ 544, 488, 426, 426, 4, 180, 4, 695, 35, 54, 433, 500, 592, 433, 262,
+ 94, 401, 401, 106, 216, 216, 106, 521, 102, 462, 518, 271, 475, 365,
+ 193, 648, 206, 424, 206, 193, 206, 206, 424, 299, 590, 590, 364, 621,
+ 67, 538, 488, 567, 51, 51, 513, 194, 81, 488, 486, 289, 567, 563,
+ 749, 563, 338, 338, 502, 563, 822, 338, 563, 338, 502, 201, 230, 201,
+ 533, 445, 175, 201, 175, 13, 85, 960, 103, 85, 175, 30, 445, 445,
+ 175, 573, 196, 877, 287, 356, 678, 235, 489, 312, 572, 264, 717, 138,
+ 295, 6, 295, 523, 55, 165, 165, 295, 138, 663, 6, 295, 6, 353, 138,
+ 6, 138, 169, 129, 784, 12, 129, 194, 605, 784, 445, 234, 627, 563,
+ 689, 627, 647, 570, 627, 570, 647, 206, 234, 215, 234, 816, 627, 816,
+ 234, 627, 215, 234, 627, 264, 427, 427, 30, 424, 161, 161, 916, 740,
+ 180, 616, 481, 514, 383, 265, 481, 164, 650, 121, 582, 689, 420, 669,
+ 589, 420, 788, 549, 165, 734, 280, 224, 146, 681, 788, 184, 398, 784,
+ 4, 398, 417, 417, 398, 636, 784, 417, 81, 398, 417, 81, 185, 827,
+ 420, 241, 420, 41, 185, 185, 718, 241, 101, 185, 185, 241, 241, 241,
+ 241, 241, 185, 324, 420, 420, 1011, 420, 827, 241, 184, 563, 241,
+ 183, 285, 529, 285, 808, 822, 891, 822, 488, 285, 486, 619, 55, 869,
+ 39, 567, 39, 289, 203, 158, 289, 710, 818, 158, 818, 355, 29, 409,
+ 203, 308, 648, 792, 308, 308, 91, 308, 6, 592, 792, 106, 106, 308,
+ 41, 178, 91, 751, 91, 259, 734, 166, 36, 327, 166, 230, 205, 205,
+ 172, 128, 230, 432, 623, 838, 623, 432, 278, 432, 42, 916, 432, 694,
+ 623, 352, 452, 93, 314, 93, 93, 641, 88, 970, 914, 230, 61, 159, 270,
+ 159, 493, 159, 755, 159, 409, 30, 30, 836, 128, 241, 99, 102, 984,
+ 538, 102, 102, 273, 639, 838, 102, 102, 136, 637, 508, 627, 285, 465,
+ 327, 327, 21, 749, 327, 749, 21, 845, 21, 21, 409, 749, 1367, 806,
+ 616, 714, 253, 616, 714, 714, 112, 375, 21, 112, 375, 375, 51, 51,
+ 51, 51, 393, 206, 870, 713, 193, 802, 21, 1061, 42, 382, 42, 543,
+ 876, 42, 876, 382, 696, 543, 635, 490, 353, 353, 417, 64, 1257, 271,
+ 64, 377, 127, 127, 537, 417, 905, 353, 538, 465, 605, 876, 427, 324,
+ 514, 852, 427, 53, 427, 557, 173, 173, 7, 1274, 563, 31, 31, 31, 745,
+ 392, 289, 230, 230, 230, 91, 218, 327, 420, 420, 128, 901, 552, 420,
+ 230, 608, 552, 476, 347, 476, 231, 159, 137, 716, 648, 716, 627, 740,
+ 718, 679, 679, 6, 718, 740, 6, 189, 679, 125, 159, 757, 1191, 409,
+ 175, 250, 409, 67, 324, 681, 605, 550, 398, 550, 931, 478, 174, 21,
+ 316, 91, 316, 654, 409, 425, 425, 699, 61, 699, 321, 698, 321, 698,
+ 61, 425, 699, 321, 409, 699, 299, 335, 321, 335, 61, 698, 699, 654,
+ 698, 299, 425, 231, 14, 121, 515, 121, 14, 165, 81, 409, 189, 81,
+ 373, 465, 463, 1055, 507, 81, 81, 189, 1246, 321, 409, 886, 104, 842,
+ 689, 300, 740, 380, 656, 656, 832, 656, 380, 300, 300, 206, 187, 175,
+ 142, 465, 206, 271, 468, 215, 560, 83, 215, 83, 215, 215, 83, 175,
+ 215, 83, 83, 111, 206, 756, 559, 756, 1367, 206, 559, 1015, 559, 559,
+ 946, 1015, 548, 559, 756, 1043, 756, 698, 159, 414, 308, 458, 997,
+ 663, 663, 347, 39, 755, 838, 323, 755, 323, 159, 159, 717, 159, 21,
+ 41, 128, 516, 159, 717, 71, 870, 755, 159, 740, 717, 374, 516, 740,
+ 51, 148, 335, 148, 335, 791, 120, 364, 335, 335, 51, 120, 251, 538,
+ 251, 971, 1395, 538, 78, 178, 538, 538, 918, 129, 918, 129, 538, 538,
+ 656, 129, 538, 538, 129, 538, 1051, 538, 128, 838, 931, 998, 823,
+ 1095, 334, 870, 334, 367, 550, 1061, 498, 745, 832, 498, 745, 716,
+ 498, 498, 128, 997, 832, 716, 832, 130, 642, 616, 497, 432, 432, 432,
+ 432, 642, 159, 432, 46, 230, 788, 160, 230, 478, 46, 693, 103, 920,
+ 230, 589, 643, 160, 616, 432, 165, 165, 583, 592, 838, 784, 583, 710,
+ 6, 583, 583, 6, 35, 230, 838, 592, 710, 6, 589, 230, 838, 30, 592,
+ 583, 6, 583, 6, 6, 583, 30, 30, 6, 375, 375, 99, 36, 1158, 425, 662,
+ 417, 681, 364, 375, 1025, 538, 822, 669, 893, 538, 538, 450, 409,
+ 632, 527, 632, 563, 632, 527, 550, 71, 698, 550, 39, 550, 514, 537,
+ 514, 537, 111, 41, 173, 592, 173, 648, 173, 173, 173, 1011, 514, 173,
+ 173, 514, 166, 648, 355, 161, 166, 648, 497, 327, 327, 550, 650, 21,
+ 425, 605, 555, 103, 425, 605, 842, 836, 1011, 636, 138, 756, 836,
+ 756, 756, 353, 1011, 636, 636, 1158, 741, 741, 842, 756, 741, 1011,
+ 677, 1011, 770, 366, 306, 488, 920, 920, 665, 775, 502, 500, 775,
+ 775, 648, 364, 833, 207, 13, 93, 500, 364, 500, 665, 500, 93, 295,
+ 183, 1293, 313, 272, 313, 279, 303, 93, 516, 93, 1013, 381, 6, 93,
+ 93, 303, 259, 643, 168, 673, 230, 1261, 230, 230, 673, 1060, 1079,
+ 1079, 550, 741, 741, 590, 527, 741, 741, 442, 741, 442, 848, 741,
+ 590, 925, 219, 527, 925, 335, 442, 590, 239, 590, 590, 590, 239, 527,
+ 239, 1033, 230, 734, 241, 741, 230, 549, 548, 1015, 1015, 32, 36,
+ 433, 465, 724, 465, 73, 73, 73, 465, 808, 73, 592, 1430, 250, 154,
+ 154, 250, 538, 353, 353, 353, 353, 353, 175, 194, 206, 538, 632,
+ 1163, 960, 175, 175, 538, 452, 632, 1163, 175, 538, 960, 194, 175,
+ 194, 632, 960, 632, 94, 632, 461, 960, 1163, 1163, 461, 632, 960,
+ 755, 707, 105, 382, 625, 382, 382, 784, 707, 871, 559, 387, 387, 871,
+ 784, 559, 784, 88, 36, 570, 314, 1028, 975, 335, 335, 398, 573, 573,
+ 573, 21, 215, 562, 738, 612, 424, 21, 103, 788, 870, 912, 23, 186,
+ 757, 73, 818, 23, 73, 563, 952, 262, 563, 137, 262, 1022, 952, 137,
+ 1273, 442, 952, 604, 137, 308, 384, 913, 235, 325, 695, 398, 95, 668,
+ 776, 713, 309, 691, 22, 10, 364, 682, 682, 578, 481, 1252, 1072,
+ 1252, 825, 578, 825, 1072, 1149, 592, 273, 387, 273, 427, 155, 1204,
+ 50, 452, 50, 1142, 50, 367, 452, 1142, 611, 367, 50, 50, 367, 50,
+ 1675, 99, 367, 50, 1501, 1099, 830, 681, 689, 917, 1089, 453, 425,
+ 235, 918, 538, 550, 335, 161, 387, 859, 324, 21, 838, 859, 1123, 21,
+ 723, 21, 335, 335, 206, 21, 364, 1426, 21, 838, 838, 335, 364, 21,
+ 21, 859, 920, 838, 838, 397, 81, 639, 397, 397, 588, 933, 933, 784,
+ 222, 830, 36, 36, 222, 1251, 266, 36, 146, 266, 366, 581, 605, 366,
+ 22, 966, 681, 681, 433, 730, 1013, 550, 21, 21, 938, 488, 516, 21,
+ 21, 656, 420, 323, 323, 323, 327, 323, 918, 581, 581, 830, 361, 830,
+ 364, 259, 364, 496, 496, 364, 691, 705, 691, 475, 427, 1145, 600,
+ 179, 427, 527, 749, 869, 689, 335, 347, 220, 298, 689, 1426, 183,
+ 554, 55, 832, 550, 550, 165, 770, 957, 67, 1386, 219, 683, 683, 355,
+ 683, 355, 355, 738, 355, 842, 931, 266, 325, 349, 256, 1113, 256,
+ 423, 960, 554, 554, 325, 554, 508, 22, 142, 22, 508, 916, 767, 55,
+ 1529, 767, 55, 1286, 93, 972, 550, 931, 1286, 1286, 972, 93, 1286,
+ 1392, 890, 93, 1286, 93, 1286, 972, 374, 931, 890, 808, 779, 975,
+ 975, 175, 173, 4, 681, 383, 1367, 173, 383, 1367, 383, 173, 175, 69,
+ 238, 146, 238, 36, 148, 888, 238, 173, 238, 148, 238, 888, 185, 925,
+ 925, 797, 925, 815, 925, 469, 784, 289, 784, 925, 797, 925, 925,
+ 1093, 925, 925, 925, 1163, 797, 797, 815, 925, 1093, 784, 636, 663,
+ 925, 187, 922, 316, 1380, 709, 916, 916, 187, 355, 948, 916, 187,
+ 916, 916, 948, 948, 916, 355, 316, 316, 334, 300, 1461, 36, 583,
+ 1179, 699, 235, 858, 583, 699, 858, 699, 1189, 1256, 1189, 699, 797,
+ 699, 699, 699, 699, 427, 488, 427, 488, 175, 815, 656, 656, 150, 322,
+ 465, 322, 870, 465, 1099, 582, 665, 767, 749, 635, 749, 600, 1448,
+ 36, 502, 235, 502, 355, 502, 355, 355, 355, 172, 355, 355, 95, 866,
+ 425, 393, 1165, 42, 42, 42, 393, 939, 909, 909, 836, 552, 424, 1333,
+ 852, 897, 1426, 1333, 1446, 1426, 997, 1011, 852, 1198, 55, 32, 239,
+ 588, 681, 681, 239, 1401, 32, 588, 239, 462, 286, 1260, 984, 1160,
+ 960, 960, 486, 828, 462, 960, 1199, 581, 850, 663, 581, 751, 581,
+ 581, 1571, 252, 252, 1283, 264, 430, 264, 430, 430, 842, 252, 745,
+ 21, 307, 681, 1592, 488, 857, 857, 1161, 857, 857, 857, 138, 374,
+ 374, 1196, 374, 1903, 1782, 1626, 414, 112, 1477, 1040, 356, 775,
+ 414, 414, 112, 356, 775, 435, 338, 1066, 689, 689, 1501, 689, 1249,
+ 205, 689, 765, 220, 308, 917, 308, 308, 220, 327, 387, 838, 917, 917,
+ 917, 220, 662, 308, 220, 387, 387, 220, 220, 308, 308, 308, 387,
+ 1009, 1745, 822, 279, 554, 1129, 543, 383, 870, 1425, 241, 870, 241,
+ 383, 716, 592, 21, 21, 592, 425, 550, 550, 550, 427, 230, 57, 483,
+ 784, 860, 57, 308, 57, 486, 870, 447, 486, 433, 433, 870, 433, 997,
+ 486, 443, 433, 433, 997, 486, 1292, 47, 708, 81, 895, 394, 81, 935,
+ 81, 81, 81, 374, 986, 916, 1103, 1095, 465, 495, 916, 667, 1745, 518,
+ 220, 1338, 220, 734, 1294, 741, 166, 828, 741, 741, 1165, 1371, 1371,
+ 471, 1371, 647, 1142, 1878, 1878, 1371, 1371, 822, 66, 327, 158, 427,
+ 427, 465, 465, 676, 676, 30, 30, 676, 676, 893, 1592, 93, 455, 308,
+ 582, 695, 582, 629, 582, 85, 1179, 85, 85, 1592, 1179, 280, 1027,
+ 681, 398, 1027, 398, 295, 784, 740, 509, 425, 968, 509, 46, 833, 842,
+ 401, 184, 401, 464, 6, 1501, 1501, 550, 538, 883, 538, 883, 883, 883,
+ 1129, 550, 550, 333, 689, 948, 21, 21, 241, 2557, 2094, 273, 308, 58,
+ 863, 893, 1086, 409, 136, 1086, 592, 592, 830, 830, 883, 830, 277,
+ 68, 689, 902, 277, 453, 507, 129, 689, 630, 664, 550, 128, 1626,
+ 1626, 128, 902, 312, 589, 755, 755, 589, 755, 407, 1782, 589, 784,
+ 1516, 1118, 407, 407, 1447, 589, 235, 755, 1191, 235, 235, 407, 128,
+ 589, 1118, 21, 383, 1331, 691, 481, 383, 1129, 1129, 1261, 1104,
+ 1378, 1129, 784, 1129, 1261, 1129, 947, 1129, 784, 784, 1129, 1129,
+ 35, 1104, 35, 866, 1129, 1129, 64, 481, 730, 1260, 481, 970, 481,
+ 481, 481, 481, 863, 481, 681, 699, 863, 486, 681, 481, 481, 55, 55,
+ 235, 1364, 944, 632, 822, 401, 822, 952, 822, 822, 99, 550, 2240,
+ 550, 70, 891, 860, 860, 550, 550, 916, 1176, 1530, 425, 1530, 916,
+ 628, 1583, 916, 628, 916, 916, 628, 628, 425, 916, 1062, 1265, 916,
+ 916, 916, 280, 461, 916, 916, 1583, 628, 1062, 916, 916, 677, 1297,
+ 924, 1260, 83, 1260, 482, 433, 234, 462, 323, 1656, 997, 323, 323,
+ 931, 838, 931, 1933, 1391, 367, 323, 931, 1391, 1391, 103, 1116,
+ 1116, 1116, 769, 1195, 1218, 312, 791, 312, 741, 791, 997, 312, 334,
+ 334, 312, 287, 287, 633, 1397, 1426, 605, 1431, 327, 592, 705, 1194,
+ 592, 1097, 1118, 1503, 1267, 1267, 1267, 618, 1229, 734, 1089, 785,
+ 1089, 1129, 1148, 1148, 1089, 915, 1148, 1129, 1148, 1011, 1011,
+ 1229, 871, 1560, 1560, 1560, 563, 1537, 1009, 1560, 632, 985, 592,
+ 1308, 592, 882, 145, 145, 397, 837, 383, 592, 592, 832, 36, 2714,
+ 2107, 1588, 1347, 36, 36, 1443, 1453, 334, 2230, 1588, 1169, 650,
+ 1169, 2107, 425, 425, 891, 891, 425, 2532, 679, 274, 274, 274, 325,
+ 274, 1297, 194, 1297, 627, 314, 917, 314, 314, 1501, 414, 1490, 1036,
+ 592, 1036, 1025, 901, 1218, 1025, 901, 280, 592, 592, 901, 1461, 159,
+ 159, 159, 2076, 1066, 1176, 1176, 516, 327, 516, 1179, 1176, 899,
+ 1176, 1176, 323, 1187, 1229, 663, 1229, 504, 1229, 916, 1229, 916,
+ 1661, 41, 36, 278, 1027, 648, 648, 648, 1626, 648, 646, 1179, 1580,
+ 1061, 1514, 1008, 1741, 2076, 1514, 1008, 952, 1089, 427, 952, 427,
+ 1083, 425, 427, 1089, 1083, 425, 427, 425, 230, 920, 1678, 920, 1678,
+ 189, 189, 953, 189, 133, 189, 1075, 189, 189, 133, 1264, 725, 189,
+ 1629, 189, 808, 230, 230, 2179, 770, 230, 770, 230, 21, 21, 784,
+ 1118, 230, 230, 230, 770, 1118, 986, 808, 916, 30, 327, 918, 679,
+ 414, 916, 1165, 1355, 916, 755, 733, 433, 1490, 433, 433, 433, 605,
+ 433, 433, 433, 1446, 679, 206, 433, 21, 2452, 206, 206, 433, 1894,
+ 206, 822, 206, 2073, 206, 206, 21, 822, 21, 206, 206, 21, 383, 1513,
+ 375, 1347, 432, 1589, 172, 954, 242, 1256, 1256, 1248, 1256, 1256,
+ 1248, 1248, 1256, 842, 13, 592, 13, 842, 1291, 592, 21, 175, 13, 592,
+ 13, 13, 1426, 13, 1541, 445, 808, 808, 863, 647, 219, 1592, 1029,
+ 1225, 917, 1963, 1129, 555, 1313, 550, 660, 550, 220, 660, 552, 663,
+ 220, 533, 220, 383, 550, 1278, 1495, 636, 842, 1036, 425, 842, 425,
+ 1537, 1278, 842, 554, 1508, 636, 554, 301, 842, 792, 1392, 1021, 284,
+ 1172, 997, 1021, 103, 1316, 308, 1210, 848, 848, 1089, 1089, 848,
+ 848, 67, 1029, 827, 1029, 2078, 827, 1312, 1029, 827, 590, 872, 1312,
+ 427, 67, 67, 67, 67, 872, 827, 872, 2126, 1436, 26, 2126, 67, 1072,
+ 2126, 1610, 872, 1620, 883, 883, 1397, 1189, 555, 555, 563, 1189,
+ 555, 640, 555, 640, 1089, 1089, 610, 610, 1585, 610, 1355, 610, 1015,
+ 616, 925, 1015, 482, 230, 707, 231, 888, 1355, 589, 1379, 151, 931,
+ 1486, 1486, 393, 235, 960, 590, 235, 960, 422, 142, 285, 285, 327,
+ 327, 442, 2009, 822, 445, 822, 567, 888, 2611, 1537, 323, 55, 1537,
+ 323, 888, 2611, 323, 1537, 323, 58, 445, 593, 2045, 593, 58, 47, 770,
+ 842, 47, 47, 842, 842, 648, 2557, 173, 689, 2291, 1446, 2085, 2557,
+ 2557, 2291, 1780, 1535, 2291, 2391, 808, 691, 1295, 1165, 983, 948,
+ 2000, 948, 983, 983, 2225, 2000, 983, 983, 705, 948, 2000, 1795,
+ 1592, 478, 592, 1795, 1795, 663, 478, 1790, 478, 592, 1592, 173, 901,
+ 312, 4, 1606, 173, 838, 754, 754, 128, 550, 1166, 551, 1480, 550,
+ 550, 1875, 1957, 1166, 902, 1875, 550, 550, 551, 2632, 551, 1875,
+ 1875, 551, 2891, 2159, 2632, 3231, 551, 815, 150, 1654, 1059, 1059,
+ 734, 770, 555, 1592, 555, 2059, 770, 770, 1803, 627, 627, 627, 2059,
+ 931, 1272, 427, 1606, 1272, 1606, 1187, 1204, 397, 822, 21, 1645,
+ 263, 263, 822, 263, 1645, 280, 263, 605, 1645, 2014, 21, 21, 1029,
+ 263, 1916, 2291, 397, 397, 496, 270, 270, 1319, 264, 1638, 264, 986,
+ 1278, 1397, 1278, 1191, 409, 1191, 740, 1191, 754, 754, 387, 63, 948,
+ 666, 666, 1198, 548, 63, 1248, 285, 1248, 169, 1248, 1248, 285, 918,
+ 224, 285, 1426, 1671, 514, 514, 717, 514, 51, 1521, 1745, 51, 605,
+ 1191, 51, 128, 1191, 51, 51, 1521, 267, 513, 952, 966, 1671, 897, 51,
+ 71, 592, 986, 986, 1121, 592, 280, 2000, 2000, 1165, 1165, 1165,
+ 1818, 222, 1818, 1165, 1252, 506, 327, 443, 432, 1291, 1291, 2755,
+ 1413, 520, 1318, 227, 1047, 828, 520, 347, 1364, 136, 136, 452, 457,
+ 457, 132, 457, 488, 1087, 1013, 2225, 32, 1571, 2009, 483, 67, 483,
+ 740, 740, 1013, 2854, 866, 32, 2861, 866, 887, 32, 2444, 740, 32, 32,
+ 866, 2225, 866, 32, 1571, 2627, 32, 850, 1675, 569, 1158, 32, 1158,
+ 1797, 2641, 1565, 1158, 569, 1797, 1158, 1797, 55, 1703, 42, 55,
+ 2562, 675, 1703, 42, 55, 749, 488, 488, 347, 1206, 1286, 1286, 488,
+ 488, 1206, 1286, 1206, 1286, 550, 550, 1790, 860, 550, 2452, 550,
+ 550, 2765, 1089, 1633, 797, 2244, 1313, 194, 2129, 194, 194, 194,
+ 818, 32, 194, 450, 1313, 2387, 194, 1227, 2387, 308, 2232, 526, 476,
+ 278, 830, 830, 194, 830, 194, 278, 194, 714, 476, 830, 714, 830, 278,
+ 830, 2532, 1218, 1759, 1446, 960, 1747, 187, 1446, 1759, 960, 105,
+ 1446, 1446, 1271, 1446, 960, 960, 1218, 1446, 1446, 105, 1446, 960,
+ 488, 1446, 427, 534, 842, 1969, 2460, 1969, 842, 842, 1969, 427, 941,
+ 2160, 427, 230, 938, 2075, 1675, 1675, 895, 1675, 34, 129, 1811, 239,
+ 749, 1957, 2271, 749, 1908, 129, 239, 239, 129, 129, 2271, 2426,
+ 1355, 1756, 194, 1583, 194, 194, 1583, 194, 1355, 194, 1628, 2221,
+ 1269, 2425, 1756, 1355, 1355, 1583, 1033, 427, 582, 30, 582, 582,
+ 935, 1444, 1962, 915, 733, 915, 938, 1962, 767, 353, 1630, 1962,
+ 1962, 563, 733, 563, 733, 353, 822, 1630, 740, 2076, 2076, 2076, 589,
+ 589, 2636, 866, 589, 947, 1528, 125, 273, 1058, 1058, 1161, 1635,
+ 1355, 1161, 1161, 1355, 1355, 650, 1206, 1206, 784, 784, 784, 784,
+ 784, 412, 461, 412, 2240, 412, 679, 891, 461, 679, 679, 189, 189,
+ 1933, 1651, 2515, 189, 1386, 538, 1386, 1386, 1187, 1386, 2423, 2601,
+ 2285, 175, 175, 2331, 194, 3079, 384, 538, 2365, 2294, 538, 2166,
+ 1841, 3326, 1256, 3923, 976, 85, 550, 550, 1295, 863, 863, 550, 1249,
+ 550, 1759, 146, 1069, 920, 2633, 885, 885, 1514, 1489, 166, 1514,
+ 2041, 885, 2456, 885, 2041, 1081, 1948, 362, 550, 94, 324, 2308, 94,
+ 2386, 94, 550, 874, 1329, 1759, 2280, 1487, 493, 493, 2099, 2599,
+ 1431, 1086, 1514, 1086, 2099, 1858, 368, 1330, 2599, 1858, 2846,
+ 2846, 2907, 2846, 713, 713, 1854, 1123, 713, 713, 3010, 1123, 3010,
+ 538, 713, 1123, 447, 822, 555, 2011, 493, 508, 2292, 555, 1736, 2135,
+ 2704, 555, 2814, 555, 2000, 555, 555, 822, 914, 327, 679, 327, 648,
+ 537, 2263, 931, 1496, 537, 1296, 1745, 1592, 1658, 1795, 650, 1592,
+ 1745, 1745, 1658, 1592, 1745, 1592, 1745, 1658, 1338, 2124, 1592,
+ 1745, 1745, 1745, 837, 1726, 2897, 1118, 1118, 230, 1118, 1118, 1118,
+ 1388, 1748, 514, 128, 1165, 931, 514, 2974, 2041, 2387, 2041, 979,
+ 185, 36, 1269, 550, 173, 812, 36, 1165, 2676, 2562, 1473, 2885, 1982,
+ 1578, 1578, 383, 383, 2360, 383, 1578, 2360, 1584, 1982, 1578, 1578,
+ 1578, 2019, 1036, 355, 724, 2023, 205, 303, 355, 1036, 1966, 355,
+ 1036, 401, 401, 401, 830, 401, 849, 578, 401, 849, 849, 578, 1776,
+ 1123, 552, 2632, 808, 1446, 1120, 373, 1529, 1483, 1057, 893, 1284,
+ 1430, 1529, 1529, 2632, 1352, 2063, 1606, 1352, 1606, 2291, 3079,
+ 2291, 1529, 506, 838, 1606, 1606, 1352, 1529, 1529, 1483, 1529, 1606,
+ 1529, 259, 902, 259, 902, 612, 612, 284, 398, 2991, 1534, 1118, 1118,
+ 1118, 1118, 1118, 734, 284, 2224, 398, 734, 284, 734, 398, 3031, 398,
+ 734, 1707, 2643, 1344, 1477, 475, 1818, 194, 1894, 691, 1528, 1184,
+ 1207, 1501, 6, 2069, 871, 2069, 3548, 1443, 2069, 2685, 3265, 1350,
+ 3265, 2069, 2069, 128, 1313, 128, 663, 414, 1313, 414, 2000, 128,
+ 2000, 663, 1313, 699, 1797, 550, 327, 550, 1526, 699, 327, 1797,
+ 1526, 550, 550, 327, 550, 1426, 1426, 1426, 2285, 1123, 890, 728,
+ 1707, 728, 728, 327, 253, 1187, 1281, 1364, 1571, 2170, 755, 3232,
+ 925, 1496, 2170, 2170, 1125, 443, 902, 902, 925, 755, 2078, 2457,
+ 902, 2059, 2170, 1643, 1129, 902, 902, 1643, 1129, 606, 36, 103, 338,
+ 338, 1089, 338, 338, 338, 1089, 338, 36, 340, 1206, 1176, 2041, 833,
+ 1854, 1916, 1916, 1501, 2132, 1736, 3065, 367, 1934, 833, 833, 833,
+ 2041, 3017, 2147, 818, 1397, 828, 2147, 398, 828, 818, 1158, 818,
+ 689, 327, 36, 1745, 2132, 582, 1475, 189, 582, 2132, 1191, 582, 2132,
+ 1176, 1176, 516, 2610, 2230, 2230, 64, 1501, 537, 1501, 173, 2230,
+ 2988, 1501, 2694, 2694, 537, 537, 173, 173, 1501, 537, 64, 173, 173,
+ 64, 2230, 537, 2230, 537, 2230, 2230, 2069, 3142, 1645, 689, 1165,
+ 1165, 1963, 514, 488, 1963, 1145, 235, 1145, 1078, 1145, 231, 2405,
+ 552, 21, 57, 57, 57, 1297, 1455, 1988, 2310, 1885, 2854, 2014, 734,
+ 1705, 734, 2854, 734, 677, 1988, 1660, 734, 677, 734, 677, 677, 734,
+ 2854, 1355, 677, 1397, 2947, 2386, 1698, 128, 1698, 3028, 2386, 2437,
+ 2947, 2386, 2643, 2386, 2804, 1188, 335, 746, 1187, 1187, 861, 2519,
+ 1917, 2842, 1917, 675, 1308, 234, 1917, 314, 314, 2339, 2339, 2592,
+ 2576, 902, 916, 2339, 916, 2339, 916, 2339, 916, 1089, 1089, 2644,
+ 1221, 1221, 2446, 308, 308, 2225, 2225, 3192, 2225, 555, 1592, 1592,
+ 555, 893, 555, 550, 770, 3622, 2291, 2291, 3419, 465, 250, 2842,
+ 2291, 2291, 2291, 935, 160, 1271, 308, 325, 935, 1799, 1799, 1891,
+ 2227, 1799, 1598, 112, 1415, 1840, 2014, 1822, 2014, 677, 1822, 1415,
+ 1415, 1822, 2014, 2386, 2159, 1822, 1415, 1822, 179, 1976, 1033, 179,
+ 1840, 2014, 1415, 1970, 1970, 1501, 563, 563, 563, 462, 563, 1970,
+ 1158, 563, 563, 1541, 1238, 383, 235, 1158, 383, 1278, 383, 1898,
+ 2938, 21, 2938, 1313, 2201, 2059, 423, 2059, 1313, 872, 1313, 2044,
+ 89, 173, 3327, 1660, 2044, 1623, 173, 1114, 1114, 1592, 1868, 1651,
+ 1811, 383, 3469, 1811, 1651, 869, 383, 383, 1651, 1651, 3223, 2166,
+ 3469, 767, 383, 1811, 767, 2323, 3355, 1457, 3341, 2640, 2976, 2323,
+ 3341, 2323, 2640, 103, 103, 1161, 1080, 2429, 370, 2018, 2854, 2429,
+ 2166, 2429, 2094, 2207, 871, 1963, 1963, 2023, 2023, 2336, 663, 2893,
+ 1580, 691, 663, 705, 2046, 2599, 409, 2295, 1118, 2494, 1118, 1950,
+ 549, 2494, 2453, 2046, 2494, 2453, 2046, 2453, 2046, 409, 1118, 4952,
+ 2291, 2225, 1894, 1423, 2498, 567, 4129, 1475, 1501, 795, 463, 2084,
+ 828, 828, 232, 828, 232, 232, 1818, 1818, 666, 463, 232, 220, 220,
+ 2162, 2162, 833, 4336, 913, 35, 913, 21, 2927, 886, 3037, 383, 886,
+ 876, 1747, 383, 916, 916, 916, 2927, 916, 1747, 837, 1894, 717, 423,
+ 481, 1894, 1059, 2262, 3206, 4700, 1059, 3304, 2262, 871, 1831, 871,
+ 3304, 1059, 1158, 1934, 1158, 756, 1511, 41, 978, 1934, 2603, 720,
+ 41, 756, 41, 325, 2611, 1158, 173, 1123, 1934, 1934, 1511, 2045,
+ 2045, 2045, 1423, 3206, 3691, 2512, 3206, 2512, 2000, 1811, 2504,
+ 2504, 2611, 2437, 2437, 2437, 1455, 893, 150, 2665, 1966, 605, 398,
+ 2331, 1177, 516, 1962, 4241, 94, 1252, 760, 1292, 1962, 1373, 2000,
+ 1990, 3684, 42, 1868, 3779, 1811, 1811, 2041, 3010, 5436, 1780, 2041,
+ 1868, 1811, 1780, 1811, 1868, 1811, 2041, 1868, 1811, 5627, 4274,
+ 1811, 1868, 4602, 1811, 1811, 1474, 2665, 235, 1474, 2665}
diff --git a/luby.go b/luby.go
new file mode 100644
index 0000000..cefdc23
--- /dev/null
+++ b/luby.go
@@ -0,0 +1,231 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+/*
+Package fountain includes implementations of several fountain codes.
+
+Fountain codes have the property that a very large (more or less unlimited)
+number of code blocks can be generated from a fixed number of source
+blocks. The original message can be recovered from any subset of sufficient
+size of these code blocks, so even if some code blocks are lost, the message
+can still be reconstructed once a sufficient number have been received.
+So in a transmission system, the receiver need not notify the
+transmitter about every code block, it only need notify the transmitter
+when the source message has been fully reconstructed.
+
+The overall approach used by this package is that there are various codec
+implementations which follow the same overall algorithm -- splitting a source
+message into source blocks, manipulating those source blocks to produce a
+set of precode blocks, then for each code block to be produced, picking
+constituent precode blocks to use to create the code block, and then using
+an LT (Luby Transform) process to produce the code blocks.
+*/
+package fountain
+
+import (
+ "math/rand"
+)
+
+// Codec is an interface for fountain codes which follow the general
+// scheme of preparing intermediate encoding representations based on the input
+// message and picking LT composition indices given an integer code block number.
+type Codec interface {
+ // SourceBlocks returns the number of source blocks to be used in the
+ // codec. Note that this may not be the same number as the number of intermediate
+ // code blocks. It is the minimum number of encoded blocks necessary to
+ // reconstruct the original message.
+ SourceBlocks() int
+
+ // GenerateIntermediateBlocks prepares a set of precode blocks given the input
+ // message. The precode blocks may just be identically the blocks in the
+ // original message, or it may be a transformation on those source blocks
+ // derived from a codec-specific relationship.
+ // For example, in Online codes, this consists of adding auxiliary blocks.
+ // In a Raptor code, the entire set of source blocks is transformed into a
+ // different set of precode blocks.
+ GenerateIntermediateBlocks(message []byte, numBlocks int) []block
+
+ // PickIndices asks the codec to select the (non-strict subset of the) precode
+ // blocks to be used in the LT composition of a particular code block. These
+ // blocks will then be XORed to produce the code block.
+ PickIndices(codeBlockIndex int64) []int
+
+ // NewDecoder creates a decoder suitable for use with blocks encoded using this
+ // codec for a known message size (in bytes). The decoder will be initialized
+ // and ready to receive incoming blocks for decoding.
+ NewDecoder(messageLength int) Decoder
+}
+
+// LTBlock is an encoded block structure representing a block created using
+// the LT transform.
+type LTBlock struct {
+ // BlockCode is the ID used to construct the encoded block.
+ // The way in which the ID is mapped to the choice of blocks will vary by
+ // codec.
+ BlockCode int64
+
+ // Data is the contents of the encoded block.
+ Data []byte
+}
+
+// Decoder is an interface allowing decoding of fountain-code-encoded messages
+// as the blocks are received.
+type Decoder interface {
+ // AddBlocks adds a set of encoded blocks to the decoder. Returns true if the
+ // message can be fully decoded. False if there is insufficient information.
+ AddBlocks(blocks []LTBlock) bool
+
+ // Decode extracts the decoded message from the decoder. If the decoder does
+ // not have sufficient information to produce an output, returns a nil slice.
+ Decode() []byte
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Implementation of Luby Transform codes.
+// The Luby Transform (LT) converts a source text split into a number of source
+// blocks into an unbounded set of code blocks, each of which is formed of an
+// XOR operation of some subset of the source blocks.
+// See "LT Codes" -- M. Luby (2002)
+
+// lubyCodec contains the codec information for the Luby Transform encoder
+// and decoder.
+// Implements fountain.Codec.
+type lubyCodec struct {
+ // sourceBlocks is the number of source blocks (N) the source message is split into.
+ sourceBlocks int
+
+ // random is a source of randomness for sampling the degree distribution
+ // and the source blocks when composing a code block.
+ random *rand.Rand
+
+ // degreeCDF is the degree distribution function from which encoding block
+ // compositions are chosen.
+ degreeCDF []float64
+}
+
+// NewLubyCodec creates a new Codec using the provided number of source blocks,
+// PRNG, and degree distribution function.
+// The intermediate blocks will be a roughly-equal-sized partition of the source
+// message padded so that all blocks have equal size. The indices will be picked
+// using the provided PRNG seeded with the BlockCode ID of the LTBlock
+// to be created, according to the degree CDF provided.
+func NewLubyCodec(sourceBlocks int, random *rand.Rand, degreeCDF []float64) Codec {
+ return &lubyCodec{
+ sourceBlocks: sourceBlocks,
+ random: random,
+ degreeCDF: degreeCDF}
+}
+
+// SourceBlocks retrieves the number of source blocks the codec is configured to use.
+func (c *lubyCodec) SourceBlocks() int {
+ return c.sourceBlocks
+}
+
+// PickIndices uses the provided PRNG to select a random number of source
+// blocks with degree d, given by a random selection in the degreeCDF parameter.
+// The degree distribution is how likely the encoder is to pick code blocks composed
+// of d source blocks.
+func (c *lubyCodec) PickIndices(codeBlockIndex int64) []int {
+ c.random.Seed(codeBlockIndex)
+ d := pickDegree(c.random, c.degreeCDF)
+ return sampleUniform(c.random, d, c.sourceBlocks)
+}
+
+// GenerateIntermediateEncoding for the LubyCodec simply splits the source message
+// into numBlocks blocks of roughly equal size, padding shorter ones so that all
+// blocks are the same length.
+func (c *lubyCodec) GenerateIntermediateBlocks(message []byte, numBlocks int) []block {
+ long, short := partitionBytes(message, c.sourceBlocks)
+ return equalizeBlockLengths(long, short)
+}
+
+// generateLubyTransformBlock generates a single code block from the set of
+// source blocks, given the composition indices, by XORing the source blocks
+// together.
+func generateLubyTransformBlock(source []block, indices []int) block {
+ var symbol block
+
+ for _, i := range indices {
+ if i < len(source) {
+ symbol.xor(source[i])
+ }
+ }
+
+ return symbol
+}
+
+// EncodeLTBlocks encodes a sequence of LT-encoded code blocks from the given message
+// and the block IDs. Suitable for use with any fountain.Codec.
+// Note: This method is destructive to the message array.
+func EncodeLTBlocks(message []byte, encodedBlockIDs []int64, c Codec) []LTBlock {
+ source := c.GenerateIntermediateBlocks(message, c.SourceBlocks())
+
+ ltBlocks := make([]LTBlock, len(encodedBlockIDs))
+ for i := range encodedBlockIDs {
+ indices := c.PickIndices(encodedBlockIDs[i])
+ ltBlocks[i].BlockCode = encodedBlockIDs[i]
+ b := generateLubyTransformBlock(source, indices)
+ ltBlocks[i].Data = make([]byte, b.length())
+ copy(ltBlocks[i].Data, b.data)
+ }
+ return ltBlocks
+}
+
+// NewDecoder creates a luby transform decoder
+func (c *lubyCodec) NewDecoder(messageLength int) Decoder {
+ return newLubyDecoder(c, messageLength)
+}
+
+// lubyDecoder is the state required to decode a Luby Transform message.
+type lubyDecoder struct {
+ codec *lubyCodec
+ messageLength int
+
+ // The sparse equation matrix used for decoding.
+ matrix sparseMatrix
+}
+
+// newLubyDecoder creates a new decoder for a particular Luby Transform message.
+// The codec parameters used to create the original encoding blocks must be provided.
+// The decoder is only valid for decoding code blocks for a particular message.
+func newLubyDecoder(c *lubyCodec, length int) *lubyDecoder {
+ d := &lubyDecoder{codec: c, messageLength: length}
+ d.matrix.coeff = make([][]int, c.SourceBlocks())
+ d.matrix.v = make([]block, c.SourceBlocks())
+
+ return d
+}
+
+// AddBlocks adds a set of encoded blocks to the decoder. Returns true if the
+// message can be fully decoded. False if there is insufficient information.
+func (d *lubyDecoder) AddBlocks(blocks []LTBlock) bool {
+ for i := range blocks {
+ indices := d.codec.PickIndices(blocks[i].BlockCode)
+ d.matrix.addEquation(indices, block{data: blocks[i].Data})
+ }
+ return d.matrix.determined()
+}
+
+// Decode extracts the decoded message from the decoder. If the decoder does
+// not have sufficient information to produce an output, returns a nil slice.
+func (d *lubyDecoder) Decode() []byte {
+ if !d.matrix.determined() {
+ return nil
+ }
+
+ d.matrix.reduce()
+
+ lenLong, lenShort, numLong, numShort := partition(d.messageLength, d.codec.SourceBlocks())
+ return d.matrix.reconstruct(d.messageLength, lenLong, lenShort, numLong, numShort)
+}
diff --git a/luby_test.go b/luby_test.go
new file mode 100644
index 0000000..44bcba1
--- /dev/null
+++ b/luby_test.go
@@ -0,0 +1,95 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+func TestLubyTransformBlockGenerator(t *testing.T) {
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ codec := NewLubyCodec(4, rand.New(NewMersenneTwister(200)), solitonDistribution(4))
+
+ wantIndices := [][]int{
+ {0},
+ {1},
+ {3},
+ {0, 1},
+ {1, 2, 3},
+ }
+
+ // These magic seeds are chosen such that they will generate the block compositions
+ // in wantIndices given the PRNG with which we initialized the codec above (the
+ // Mersenne Twister with seed=200).
+ encodeBlocks := []int64{7, 34, 5, 31, 25}
+ for i := range wantIndices {
+ indices := codec.PickIndices(encodeBlocks[i])
+ if !reflect.DeepEqual(indices, wantIndices[i]) {
+ t.Logf("Got %v indices for %d, want %v", indices, encodeBlocks[i], wantIndices[i])
+ }
+ }
+
+ source := codec.GenerateIntermediateBlocks(message, codec.SourceBlocks())
+ lubyBlocks := make([]LTBlock, len(encodeBlocks))
+ for i := range encodeBlocks {
+ b := generateLubyTransformBlock(source, wantIndices[i])
+ lubyBlocks[i].BlockCode = encodeBlocks[i]
+ lubyBlocks[i].Data = make([]byte, b.length())
+ copy(lubyBlocks[i].Data, b.data)
+ }
+
+ if len(source) != codec.SourceBlocks() {
+ t.Logf("Got %d encoded blocks, want %d", len(source), codec.SourceBlocks())
+ }
+
+ if string(lubyBlocks[0].Data) != "abcdefg" {
+ t.Errorf("Data for {0} block is %v, should be 'abcdefg'", string(lubyBlocks[0].Data))
+ }
+ if string(lubyBlocks[1].Data) != "hijklmn" {
+ t.Errorf("Data for {1} block is %v, should be 'hijklmn'", string(lubyBlocks[1].Data))
+ }
+ if string(lubyBlocks[2].Data) != "uvwxyz" {
+ t.Errorf("Data for {1} block is %v, should be 'uvwxyz'", string(lubyBlocks[2].Data))
+ }
+ if lubyBlocks[3].Data[0] != 'a'^'h' {
+ t.Errorf("Data[0] for {0, 1} block is %d, should be 'a'^'h' (%d)",
+ lubyBlocks[3].Data[0], 'a'^'h')
+ }
+ if lubyBlocks[4].Data[0] != 'h'^'o'^'u' {
+ t.Errorf("Data[0] for {1,2,3} block is %d, should be 'h'^'o'^'u' (%d)",
+ lubyBlocks[3].Data[0], 'h'^'o'^'u')
+ }
+}
+
+func TestLubyDecoder(t *testing.T) {
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ codec := NewLubyCodec(4, rand.New(NewMersenneTwister(200)), solitonDistribution(4))
+
+ encodeBlocks := []int64{7, 34, 5, 31, 25}
+ lubyBlocks := EncodeLTBlocks(message, encodeBlocks, codec)
+
+ decoder := codec.NewDecoder(len(message))
+ determined := decoder.AddBlocks(lubyBlocks)
+ if !determined {
+ t.Errorf("After adding code blocks, decoder is still undetermined.")
+ }
+ decoded := decoder.Decode()
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Decoded luby transform message is %v, expected %v", decoded, message)
+ t.Logf("String value = %v", string(decoded))
+ }
+}
diff --git a/mersenne.go b/mersenne.go
new file mode 100644
index 0000000..24da7fc
--- /dev/null
+++ b/mersenne.go
@@ -0,0 +1,200 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/rand"
+)
+
+// MersenneTwister is an implementation of the MT19937 PRNG of Matsumoto and Nishimura.
+// Following http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
+// Uses the 32-bit version of the algorithm.
+// Satisfies math/rand.Source
+type MersenneTwister struct {
+ mt [624]uint32
+ index int
+ initialized bool
+}
+
+// NewMersenneTwister creates a new MT19937 PRNG with the given seed. The seed
+// is converted to a 32-bit seed by XORing the high and low halves.
+func NewMersenneTwister(seed int64) rand.Source {
+ t := &MersenneTwister{}
+ t.Seed(seed)
+
+ return t
+}
+
+func (t *MersenneTwister) Seed(seed int64) {
+ t.initialize(uint32(((seed >> 32) ^ seed) & math.MaxUint32))
+}
+
+// Int63 produces a new int64 value between 0 and 2^63-1 by combining the bits
+// of two Uint32 values.
+func (t *MersenneTwister) Int63() int64 {
+ a := t.Uint32()
+ b := t.Uint32()
+ return (int64(a) << 31) ^ int64(b)
+}
+
+func (t *MersenneTwister) Uint32() uint32 {
+ if !t.initialized {
+ t.initialize(4357) // value from original paper; lets the twister do something reasonable when empty
+ }
+
+ // Every 624 calls, revolve the untempered seed matrix
+ if t.index == 0 {
+ t.generateUntempered()
+ }
+
+ y := t.mt[t.index]
+ t.index++
+ if t.index >= len(t.mt) {
+ t.index = 0
+ }
+ y ^= y >> 11
+ y ^= (y << 7) & 0x9d2c5680
+ y ^= (y << 15) & 0xefc60000
+ y ^= y >> 18
+
+ return y
+}
+
+func (t *MersenneTwister) initialize(seed uint32) {
+ t.index = 0
+ t.mt[0] = seed
+
+ for i := 1; i < len(t.mt); i++ {
+ // Improved initialization: not subject to the same runs of correlated numbers.
+ t.mt[i] = (1812433253*(t.mt[i-1]^(t.mt[i-1]>>30)) + uint32(i)) & math.MaxUint32
+ // Original paper did this
+ // t.mt[i] = uint32(69069 * int64(t.mt[i-1]) & math.MaxUint32)
+ }
+ t.initialized = true
+}
+
+func (t *MersenneTwister) generateUntempered() {
+ mag01 := [2]uint32{0x0, 0x9908b0df}
+ for i := 0; i < len(t.mt); i++ {
+ y := (t.mt[i] & 0x80000000) | (t.mt[(i+1)%len(t.mt)] & 0x7fffffff)
+ t.mt[i] = (t.mt[(i+397)%len(t.mt)] ^ (y >> 1)) ^ mag01[y&0x01]
+ }
+}
+
+// MersenneTwister64 is a 64-bit MT19937 PRNG after Nishimura.
+// See http://dl.acm.org/citation.cfm?id=369540&dl=ACM&coll=DL&CFID=261426526&CFTOKEN=25107569
+// Satisfies math/rand.Source
+type MersenneTwister64 struct {
+ mt [312]uint64
+ index int
+ initialized bool
+}
+
+// NewMersenneTwister64 creates a new 64-bit version of the MT19937 PRNG.
+func NewMersenneTwister64(seed int64) rand.Source {
+ t := &MersenneTwister64{}
+ t.Seed(seed)
+
+ return t
+}
+
+// Seed initializes the state of the PRNG with the given seed value.
+func (t *MersenneTwister64) Seed(seed int64) {
+ t.initialize(uint64(seed))
+}
+
+// Int63 returns the next value from the PRNG. This value is the low 63 bits
+// of the Uint64 value.
+func (t *MersenneTwister64) Int63() int64 {
+ return int64(t.Uint64() & math.MaxInt64)
+}
+
+func (t *MersenneTwister64) initialize(seed uint64) {
+ t.index = 0
+ t.mt[0] = seed
+
+ for i := 1; i < len(t.mt); i++ {
+ t.mt[i] = 6364136223846793005*(t.mt[i-1]^(t.mt[i-1]>>62)) + uint64(i)
+ }
+ t.initialized = true
+}
+
+// SeedSlice allows the twister to be initialized with a slice of seed values.
+func (t *MersenneTwister64) SeedSlice(seed []uint64) {
+ t.initialize(19650218)
+
+ length := len(seed)
+ if len(t.mt) > length {
+ length = len(t.mt)
+ }
+
+ i, j := 1, 0
+ for k := 0; k < length; k++ {
+ t.mt[i] = (t.mt[i] ^ ((t.mt[i-1] ^ (t.mt[i-1] >> 62)) * 3935559000370003845)) + seed[j] + uint64(j)
+ i++
+ j++
+ if i >= len(t.mt) {
+ t.mt[0] = t.mt[len(t.mt)-1]
+ i = 1
+ }
+ if j >= len(seed) {
+ j = 0
+ }
+ }
+
+ for k := 0; k < len(t.mt)-1; k++ {
+ t.mt[i] = (t.mt[i] ^ ((t.mt[i-1] ^ (t.mt[i-1] >> 62)) * 2862933555777941757)) - uint64(i)
+ i++
+ if i >= len(t.mt) {
+ t.mt[0] = t.mt[len(t.mt)-1]
+ i = 1
+ }
+ }
+
+ t.mt[0] = 1 << 63
+}
+
+// Uint64 returns the next pseudo-random value from the twister.
+func (t *MersenneTwister64) Uint64() uint64 {
+ if !t.initialized {
+ t.initialize(5489)
+ }
+
+ // Every 312 calls, revolve the untempered seed matrix
+ if t.index == 0 {
+ t.generateUntempered()
+ }
+
+ y := t.mt[t.index]
+ t.index++
+ if t.index >= len(t.mt) {
+ t.index = 0
+ }
+ y ^= (y >> 29) & 0x5555555555555555
+ y ^= (y << 17) & 0x71d67fffeda60000
+ y ^= (y << 37) & 0xfff7eee000000000
+ y ^= y >> 43
+
+ return y
+}
+
+func (t *MersenneTwister64) generateUntempered() {
+ mag01 := [2]uint64{0x0, 0xb5026f5aa96619e9}
+ for i := 0; i < len(t.mt); i++ {
+ y := (t.mt[i] & 0xffffffff80000000) | (t.mt[(i+1)%len(t.mt)] & 0x7fffffff)
+ t.mt[i] = (t.mt[(i+156)%len(t.mt)] ^ (y >> 1)) ^ mag01[y&0x01]
+ }
+}
diff --git a/mersenne_test.go b/mersenne_test.go
new file mode 100644
index 0000000..7ba15f9
--- /dev/null
+++ b/mersenne_test.go
@@ -0,0 +1,596 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/big"
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+func NewTwister(seed int64) *MersenneTwister {
+ t := &MersenneTwister{}
+ t.Seed(seed)
+
+ return t
+}
+
+func TestMersenne(t *testing.T) {
+ var mt MersenneTwister
+
+ mt.Uint32()
+ if !mt.initialized {
+ t.Error("MersenneTwister should be initialized after Uint32 call.")
+ }
+
+ var mt2 MersenneTwister
+ mt2.Int63()
+ if !mt2.initialized {
+ t.Error("MersenneTwister should be initialized after Int63 call.")
+ }
+
+ var mt3 MersenneTwister
+ mt3.Seed(1)
+ if !mt3.initialized {
+ t.Error("MersenneTwister should be initialized after Seed call.")
+ }
+
+}
+
+func TestMersenneGen(t *testing.T) {
+ mt := NewTwister(1)
+
+ // See implementation at github.com/cslarsen/mersenne-twister
+ seed1Expected := []uint32{
+ 1791095845, 4282876139, 3093770124, 4005303368, 491263,
+ 550290313, 1298508491, 4290846341, 630311759, 1013994432,
+ 396591248, 1703301249, 799981516, 1666063943, 1484172013,
+ 2876537340, 1704103302, 4018109721, 2314200242, 3634877716,
+ 1800426750, 1345499493, 2942995346, 2252917204, 878115723,
+ 1904615676, 3771485674, 986026652, 117628829, 2295290254,
+ 2879636018, 3925436996, 1792310487, 1963679703, 2399554537,
+ 1849836273, 602957303, 4033523166, 850839392, 3343156310,
+ 3439171725, 3075069929, 4158651785, 3447817223, 1346146623,
+ 398576445, 2973502998, 2225448249, 3764062721, 3715233664,
+ 3842306364, 3561158865, 365262088, 3563119320, 167739021,
+ 1172740723, 729416111, 254447594, 3771593337, 2879896008,
+ 422396446, 2547196999, 1808643459, 2884732358, 4114104213,
+ 1768615473, 2289927481, 848474627, 2971589572, 1243949848,
+ 1355129329, 610401323, 2948499020, 3364310042, 3584689972,
+ 1771840848, 78547565, 146764659, 3221845289, 2680188370,
+ 4247126031, 2837408832, 3213347012, 1282027545, 1204497775,
+ 1916133090, 3389928919, 954017671, 443352346, 315096729,
+ 1923688040, 2015364118, 3902387977, 413056707, 1261063143,
+ 3879945342, 1235985687, 513207677, 558468452, 2253996187,
+ 83180453, 359158073, 2915576403, 3937889446, 908935816,
+ 3910346016, 1140514210, 1283895050, 2111290647, 2509932175,
+ 229190383, 2430573655, 2465816345, 2636844999, 630194419,
+ 4108289372, 2531048010, 1120896190, 3005439278, 992203680,
+ 439523032, 2291143831, 1778356919, 4079953217, 2982425969,
+ 2117674829, 1778886403, 2321861504, 214548472, 3287733501,
+ 2301657549, 194758406, 2850976308, 601149909, 2211431878,
+ 3403347458, 4057003596, 127995867, 2519234709, 3792995019,
+ 3880081671, 2322667597, 590449352, 1924060235, 598187340,
+ 3831694379, 3467719188, 1621712414, 1708008996, 2312516455,
+ 710190855, 2801602349, 3983619012, 1551604281, 1493642992,
+ 2452463100, 3224713426, 2739486816, 3118137613, 542518282,
+ 3793770775, 2964406140, 2678651729, 2782062471, 3225273209,
+ 1520156824, 1498506954, 3278061020, 1159331476, 1531292064,
+ 3847801996, 3233201345, 1838637662, 3785334332, 4143956457,
+ 50118808, 2849459538, 2139362163, 2670162785, 316934274,
+ 492830188, 3379930844, 4078025319, 275167074, 1932357898,
+ 1526046390, 2484164448, 4045158889, 1752934226, 1631242710,
+ 1018023110, 3276716738, 3879985479, 3313975271, 2463934640,
+ 1294333494, 12327951, 3318889349, 2650617233, 656828586}
+
+ got := make([]uint32, len(seed1Expected))
+
+ for i := range got {
+ got[i] = mt.Uint32()
+ }
+
+ if !reflect.DeepEqual(got, seed1Expected) {
+ t.Errorf("Got bad sequence for seed=1. Got %v, want %v", got, seed1Expected)
+ }
+
+ mt = NewTwister(1)
+
+ u := mt.Uint32()
+ if u != 1791095845 {
+ t.Errorf("Seed=1 value is %d, should be 1791095845", u)
+ }
+}
+
+func countBits(v uint32) int {
+ c := 0
+ for i := 0; i < 32; i++ {
+ if v&0x1 == 0x1 {
+ c++
+ }
+ v >>= 1
+ }
+ return c
+}
+
+func TestMersenneDifference(t *testing.T) {
+ // Checks that very similar seeds produce different outputs.
+ mt1 := NewTwister(123456)
+ mt2 := NewTwister(123457)
+ v1 := mt1.Uint32()
+ v2 := mt2.Uint32()
+ if math.Abs(float64(v1-v2)) < 0xF000 {
+ t.Errorf("PRNG values unseemly close for two close initial seeds. Got %d, %d", v1, v2)
+ }
+ if countBits(v1^v2) != 16 {
+ t.Errorf("Got common bits %d between random values from two close seeds. Want 16.", countBits(v1^v2))
+ }
+}
+
+// calculateBitStatistics performs about the simplest DieHard randomness test --
+// checks the bit count statistics of random numbers to make sure they match
+// a binomial distribution.
+// Returns the sum-squared-difference metric of the observed distribution and the
+// binomial distribution, normalized to the number of trials.
+func calculateBitStatistics(mt *MersenneTwister, trials int) float64 {
+ var bitsSet [33]int64
+ for i := 0; i < trials; i++ {
+ bitsSet[countBits(mt.Uint32())]++
+ }
+
+ // bitsSet should follow a binomial distribution pretty closely
+ var ssd int64
+ sum := 4294967296 // Sum of the binomial distribution from 0-32.
+ for i := 0; i < 33; i++ {
+ z := big.Int{}
+ binomial := int64(z.Binomial(32, int64(i)).Uint64()) * int64(trials) / int64(sum)
+ ssd += (bitsSet[i] - binomial) * (bitsSet[i] - binomial)
+ }
+ return math.Sqrt(float64(ssd)) / float64(trials)
+}
+
+func TestMersenneStatistics(t *testing.T) {
+ mt := NewTwister(892349823)
+ for i := 0; i < 25; i++ {
+ seed := mt.Uint32()
+ mt = NewTwister(int64(seed))
+ if calculateBitStatistics(mt, 10000) > 0.03 {
+ t.Errorf("Trial bits-set distribution not binomial within 3%% for seed %d", seed)
+ }
+ }
+}
+
+func TestMersenneSource(t *testing.T) {
+ r := rand.New(NewMersenneTwister(8189023))
+ count := 0
+ for i := 0; i < 100; i++ {
+ if r.Intn(2) == 1 {
+ count++
+ }
+ }
+ if math.Abs(float64(count)-50) > 3 {
+ t.Errorf("Expect coin toss about 50-50, was %d", count)
+ }
+}
+
+// Initial run on x86 hardware: 17.8ns/op
+func BenchmarkMersenne(b *testing.B) {
+ mt := NewTwister(903490345)
+ b.ResetTimer()
+ for i := 0; i < b.N; i++ {
+ mt.Uint32()
+ }
+}
+
+func NewTwister64(seed uint64) *MersenneTwister64 {
+ t := &MersenneTwister64{}
+ t.initialize(seed)
+
+ return t
+}
+
+func TestMersenne64(t *testing.T) {
+ mt := &MersenneTwister64{}
+ seed := []uint64{0x12345, 0x23456, 0x34567, 0x45678}
+ mt.SeedSlice(seed)
+ if !mt.initialized {
+ t.Errorf("Twister should be initialized after SeedSlice")
+ }
+ r := mt.Uint64()
+ if r != 7266447313870364031 {
+ t.Errorf("MersenneTwister64(initialized with %v).Uint64() = %d, want %d",
+ seed, r, 7266447313870364031)
+ }
+}
+
+func TestMersenne64Gen(t *testing.T) {
+ // See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html
+ var expected64 = []uint64{
+ 7266447313870364031, 4946485549665804864, 16945909448695747420,
+ 16394063075524226720, 4873882236456199058, 14877448043947020171,
+ 6740343660852211943, 13857871200353263164, 5249110015610582907,
+ 10205081126064480383, 1235879089597390050, 17320312680810499042,
+ 16489141110565194782, 8942268601720066061, 13520575722002588570,
+ 14226945236717732373, 9383926873555417063, 15690281668532552105,
+ 11510704754157191257, 15864264574919463609, 6489677788245343319,
+ 5112602299894754389, 10828930062652518694, 15942305434158995996,
+ 15445717675088218264, 4764500002345775851, 14673753115101942098,
+ 236502320419669032, 13670483975188204088, 14931360615268175698,
+ 8904234204977263924, 12836915408046564963, 12120302420213647524,
+ 15755110976537356441, 5405758943702519480, 10951858968426898805,
+ 17251681303478610375, 4144140664012008120, 18286145806977825275,
+ 13075804672185204371, 10831805955733617705, 6172975950399619139,
+ 12837097014497293886, 12903857913610213846, 560691676108914154,
+ 1074659097419704618, 14266121283820281686, 11696403736022963346,
+ 13383246710985227247, 7132746073714321322, 10608108217231874211,
+ 9027884570906061560, 12893913769120703138, 15675160838921962454,
+ 2511068401785704737, 14483183001716371453, 3774730664208216065,
+ 5083371700846102796, 9583498264570933637, 17119870085051257224,
+ 5217910858257235075, 10612176809475689857, 1924700483125896976,
+ 7171619684536160599, 10949279256701751503, 15596196964072664893,
+ 14097948002655599357, 615821766635933047, 5636498760852923045,
+ 17618792803942051220, 580805356741162327, 425267967796817241,
+ 8381470634608387938, 13212228678420887626, 16993060308636741960,
+ 957923366004347591, 6210242862396777185, 1012818702180800310,
+ 15299383925974515757, 17501832009465945633, 17453794942891241229,
+ 15807805462076484491, 8407189590930420827, 974125122787311712,
+ 1861591264068118966, 997568339582634050, 18046771844467391493,
+ 17981867688435687790, 3809841506498447207, 9460108917638135678,
+ 16172980638639374310, 958022432077424298, 4393365126459778813,
+ 13408683141069553686, 13900005529547645957, 15773550354402817866,
+ 16475327524349230602, 6260298154874769264, 12224576659776460914,
+ 6405294864092763507, 7585484664713203306, 5187641382818981381,
+ 12435998400285353380, 13554353441017344755, 646091557254529188,
+ 11393747116974949255, 16797249248413342857, 15713519023537495495,
+ 12823504709579858843, 4738086532119935073, 4429068783387643752,
+ 585582692562183870, 1048280754023674130, 6788940719869959076,
+ 11670856244972073775, 2488756775360218862, 2061695363573180185,
+ 6884655301895085032, 3566345954323888697, 12784319933059041817,
+ 4772468691551857254, 6864898938209826895, 7198730565322227090,
+ 2452224231472687253, 13424792606032445807, 10827695224855383989,
+ 11016608897122070904, 14683280565151378358, 7077866519618824360,
+ 17487079941198422333, 3956319990205097495, 5804870313319323478,
+ 8017203611194497730, 3310931575584983808, 5009341981771541845,
+ 6930001938490791874, 14415278059151389495, 11001114762641844083,
+ 6715939435439735925, 411419160297131328, 4522402260441335284,
+ 3381955501804126859, 15935778656111987797, 4345051260540166684,
+ 13978444093099579683, 9219789505504949817, 9245142924137529075,
+ 11628184459157386459, 7242398879359936370, 8511401943157540109,
+ 11948130810477009827, 6865450671488705049, 13965005347172621081,
+ 15956599226522058336, 7737868921014130584, 2107342503741411693,
+ 15818996300425101108, 16399939197527488760, 13971145494081508107,
+ 3910681448359868691, 4249175367970221090, 9735751321242454020,
+ 12418107929362160460, 241792245481991138, 5806488997649497146,
+ 10724207982663648949, 1121862814449214435, 1326996977123564236,
+ 4902706567834759475, 12782714623891689967, 7306216312942796257,
+ 15681656478863766664, 957364844878149318, 5651946387216554503,
+ 8197027112357634782, 6302075516351125977, 13454588464089597862,
+ 15638309200463515550, 10116604639722073476, 12052913535387714920,
+ 2889379661594013754, 15383926144832314187, 7841953313015471731,
+ 17310575136995821873, 9820021961316981626, 15319619724109527290,
+ 15349724127275899898, 10511508162402504492, 6289553862380300393,
+ 15046218882019267110, 11772020174577005930, 3537640779967351792,
+ 6801855569284252424, 17687268231192623388, 12968358613633237218,
+ 1429775571144180123, 10427377732172208413, 12155566091986788996,
+ 16465954421598296115, 12710429690464359999, 9547226351541565595,
+ 12156624891403410342, 2985938688676214686, 18066917785985010959,
+ 5975570403614438776, 11541343163022500560, 11115388652389704592,
+ 9499328389494710074, 9247163036769651820, 3688303938005101774,
+ 2210483654336887556, 15458161910089693228, 6558785204455557683,
+ 1288373156735958118, 18433986059948829624, 3435082195390932486,
+ 16822351800343061990, 3120532877336962310, 16681785111062885568,
+ 7835551710041302304, 2612798015018627203, 15083279177152657491,
+ 6591467229462292195, 10592706450534565444, 7438147750787157163,
+ 323186165595851698, 7444710627467609883, 8473714411329896576,
+ 2782675857700189492, 3383567662400128329, 3200233909833521327,
+ 12897601280285604448, 3612068790453735040, 8324209243736219497,
+ 15789570356497723463, 1083312926512215996, 4797349136059339390,
+ 5556729349871544986, 18266943104929747076, 1620389818516182276,
+ 172225355691600141, 3034352936522087096, 1266779576738385285,
+ 3906668377244742888, 6961783143042492788, 17159706887321247572,
+ 4676208075243319061, 10315634697142985816, 13435140047933251189,
+ 716076639492622016, 13847954035438697558, 7195811275139178570,
+ 10815312636510328870, 6214164734784158515, 16412194511839921544,
+ 3862249798930641332, 1005482699535576005, 4644542796609371301,
+ 17600091057367987283, 4209958422564632034, 5419285945389823940,
+ 11453701547564354601, 9951588026679380114, 7425168333159839689,
+ 8436306210125134906, 11216615872596820107, 3681345096403933680,
+ 5770016989916553752, 11102855936150871733, 11187980892339693935,
+ 396336430216428875, 6384853777489155236, 7551613839184151117,
+ 16527062023276943109, 13429850429024956898, 9901753960477271766,
+ 9731501992702612259, 5217575797614661659, 10311708346636548706,
+ 15111747519735330483, 4353415295139137513, 1845293119018433391,
+ 11952006873430493561, 3531972641585683893, 16852246477648409827,
+ 15956854822143321380, 12314609993579474774, 16763911684844598963,
+ 16392145690385382634, 1545507136970403756, 17771199061862790062,
+ 12121348462972638971, 12613068545148305776, 954203144844315208,
+ 1257976447679270605, 3664184785462160180, 2747964788443845091,
+ 15895917007470512307, 15552935765724302120, 16366915862261682626,
+ 8385468783684865323, 10745343827145102946, 2485742734157099909,
+ 916246281077683950, 15214206653637466707, 12895483149474345798,
+ 1079510114301747843, 10718876134480663664, 1259990987526807294,
+ 8326303777037206221, 14104661172014248293, 15531278677382192198,
+ 3874303698666230242, 3611366553819264523, 1358753803061653874,
+ 1552102816982246938, 14492630642488100979, 15001394966632908727,
+ 2273140352787320862, 17843678642369606172, 2903980458593894032,
+ 16971437123015263604, 12969653681729206264, 3593636458822318001,
+ 9719758956915223015, 7437601263394568346, 3327758049015164431,
+ 17851524109089292731, 14769614194455139039, 8017093497335662337,
+ 12026985381690317404, 739616144640253634, 15535375191850690266,
+ 2418267053891303448, 15314073759564095878, 10333316143274529509,
+ 16565481511572123421, 16317667579273275294, 13991958187675987741,
+ 3753596784796798785, 9078249094693663275, 8459506356724650587,
+ 12579909555010529099, 7827737296967050903, 5489801927693999341,
+ 10995988997350541459, 14721747867313883304, 7915884580303296560,
+ 4105766302083365910, 12455549072515054554, 13602111324515032467,
+ 5205971628932290989, 5034622965420036444, 9134927878875794005,
+ 11319873529597990213, 14815445109496752058, 2266601052460299470,
+ 5696993487088103383, 6540200741841280242, 6631495948031875490,
+ 5328340585170897740, 17897267040961463930, 9030000260502624168,
+ 14285709137129830926, 12854071997824681544, 15408328651008978682,
+ 1063314403033437073, 13765209628446252802, 242013711116865605,
+ 4772374239432528212, 2515855479965038648, 5872624715703151235,
+ 14237704570091006662, 678604024776645862, 12329607334079533339,
+ 17570877682732917020, 2695443415284373666, 4312672841405514468,
+ 6454343485137106900, 8425658828390111343, 16335501385875554899,
+ 5551095603809016713, 11781094401885925035, 9395557946368382509,
+ 9765123360948816956, 18107191819981188154, 16049267500594757404,
+ 16349966108299794199, 1040405303135858246, 2366386386131378192,
+ 223761048139910454, 15375217587047847934, 15231693398695187454,
+ 12916726640254571028, 8878036960829635584, 1626201782473074365,
+ 5758998126998248293, 18077917959300292758, 10585588923088536745,
+ 15072345664541731497, 3559348759319842667, 12744591691872202375,
+ 2388494115860283059, 6414691845696331748, 3069528498807764495,
+ 8737958486926519702, 18059264986425101074, 3139684427605102737,
+ 12378931902986734693, 410666675039477949, 12139894855769838924,
+ 5780722552400398675, 7039346665375142557, 3020733445712569008,
+ 2612305843503943561, 13651771214166527665, 16478681918975800939,
+ 566088527565499576, 4715785502295754870, 6957318344287196220,
+ 11645756868405128885, 13139951104358618000, 17650948583490040612,
+ 18168787973649736637, 5486282999836125542, 6122201977153895166,
+ 17324241605502052782, 10063523107521105867, 17537430712468011382,
+ 10828407533637104262, 10294139354198325113, 12557151830240236401,
+ 16673044307512640231, 10918020421896090419, 11077531235278014145,
+ 5499571814940871256, 2334252435740638702, 18177461912527387031,
+ 2000007376901262542, 7968425560071444214, 1472650787501520648,
+ 3115849849651526279, 7980970700139577536, 12153253535907642097,
+ 8109716914843248719, 3154976533165008908, 5553369513523832559,
+ 10345792701798576501, 3677445364544507875, 10637177623943913351,
+ 7380255087060498096, 14479400372337014801, 15381362583330700960,
+ 204531043189704802, 13699106540959723942, 3817903465872254783,
+ 10972364467110284934, 2701394334530963810, 2931625600749229147,
+ 16428252083632828910, 11873166501966812913, 5566810080537233762,
+ 7840617383807795056, 10699413880206684652, 18259119259617231436,
+ 10332714341486317526, 10137911902863059694, 669146221352346842,
+ 8373571610024623455, 10620002450820868661, 12220730820779815970,
+ 5902974968095412898, 7931010481705150841, 16413777368097063650,
+ 11273457888324769727, 13719113891065284171, 8327795098009702553,
+ 10333342364827584837, 6202832891413866653, 9137034567886143162,
+ 14514450826524340059, 473610156015331016, 813689571029117640,
+ 13776316799690285717, 10429708855338427756, 8995290140880620858,
+ 2320123852041754384, 8082864073645003641, 6961777411740398590,
+ 10008644283003991179, 3239064015890722333, 16762634970725218787,
+ 16467281536733948427, 10563290046315192938, 5108560603794851559,
+ 15121667220761532906, 14155440077372845941, 10050536352394623377,
+ 15474881667376037792, 3448088038819200619, 3692020001240358871,
+ 6444847992258394902, 8687650838094264665, 3028124591188972359,
+ 16945232313401161629, 15547830510283682816, 3982930188609442149,
+ 14270781928849894661, 13768475593433447867, 13815150225221307677,
+ 8502397232429564693, 718377350715476994, 7459266877697905475,
+ 8353375565171101521, 7807281661994435472, 16924127046922196149,
+ 10157812396471387805, 2519858716882670232, 7384148884750265792,
+ 8077153156180046901, 3499231286164597752, 2700106282881469611,
+ 14679824700835879737, 14188324938219126828, 3016120398601032793,
+ 10858152824243889420, 9412371965669250534, 4857522662584941069,
+ 984331743838900386, 4094160040294753142, 2368635764350388458,
+ 15101240511397838657, 15584415763303953578, 7831857200208015446,
+ 1952643641639729063, 4184323302594028609, 16795120381104846695,
+ 3541559381538365280, 15408472870896842474, 5628362450757896366,
+ 16277348886873708846, 12437047172652330846, 10172715019035948149,
+ 1999700669649752791, 6217957085626135027, 11220551167830336823,
+ 16478747645632411810, 5437280487207382147, 11382378739613087836,
+ 15866932785489521505, 5502694314775516684, 16440179278067648435,
+ 15510104554374162846, 15722061259110909195, 10760687291786964354,
+ 10736868329920212671, 4166148127664495614, 14303518358120527892,
+ 9122250801678898571, 10028508179936801946, 216630713752669403,
+ 10655207865433859491, 4041437116174699233, 6280982262534375348,
+ 297501356638818866, 13976146806363377485, 13752396481560145603,
+ 11472199956603637419, 16393728429143900496, 14752844047515986640,
+ 1524477318846038424, 6596889774254235440, 1591982099532234960,
+ 8065146456116391065, 3964696017750868345, 17040425970526664920,
+ 11511165586176539991, 3443401252003315103, 16314977947073778249,
+ 16860120454903458341, 5370503221561340846, 15362920279125264094,
+ 2822458124714999779, 14575378304387898337, 9689406052675046032,
+ 2872149351415175149, 13019620945255883050, 14929026760148695825,
+ 8503417349692327218, 9677798905341573754, 828949921821462483,
+ 16110482368362750196, 15794218816553655671, 14942910774764855088,
+ 12026350906243760195, 13610867176871462505, 18324536557697872582,
+ 2658962269666727629, 327225403251576027, 9207535177029277544,
+ 8744129291351887858, 6129603385168921503, 18385497655031085907,
+ 13024478718952333892, 14547683159720717167, 5932119629366981711,
+ 325385464632594563, 3559879386019806291, 6629264948665231298,
+ 14358245326238118181, 15662449672706340765, 13975503159145803297,
+ 3609534220891499022, 4224273587485638227, 9274084767162416370,
+ 13156843921244091998, 18284750575626858789, 14664767920489118779,
+ 11292057742031803221, 13919998707305829132, 14473305049457001422,
+ 9696877879685767807, 1406758246007973837, 2429517644459056881,
+ 14361215588101587430, 11386164476149757528, 10474116023593331839,
+ 2921165656527786564, 15604610369733358953, 12955027028676000544,
+ 10314281035410779907, 3167047178514709947, 1088721329408346700,
+ 17930425515478182741, 7466411836095405617, 15534027454610690575,
+ 10879629128927506091, 11502219301371200635, 13915106894453889418,
+ 4226784327815861027, 12335222183627106346, 3648499746356007767,
+ 18441388887898023393, 18117929843327093625, 4237736098094830438,
+ 14229123019768296655, 3930112058127932690, 12663879236019645778,
+ 9281161952002617309, 4978473890680876319, 845759387067546611,
+ 1386164484606776333, 8008554770639925512, 11159581016793288971,
+ 18065390393740782906, 17647985458967631018, 9092379465737744314,
+ 2914678236848656327, 4376066698447630270, 16057186499919087528,
+ 3031333261848790078, 2926746602873431597, 7931945763526885287,
+ 147649915388326849, 15801792398814946230, 5265900391686545347,
+ 16173686275871890830, 7562781050481886043, 5853506575839330404,
+ 14957980734704564792, 10944286556353523404, 1783009880614150597,
+ 9529762028588888983, 822992871011696119, 2130074274744257510,
+ 8000279549284809219, 3514744284158856431, 128770032569293263,
+ 3737367602618100572, 16364836605077998543, 783266423471782696,
+ 4569418252658970391, 11093950688157406886, 14888808512267628166,
+ 4217786261273670948, 17047486076688645713, 14133826721458860485,
+ 17539744882220127106, 12394675039129853905, 5757634999463277090,
+ 9621947619435861331, 1182210208559436772, 14603391040490913939,
+ 17481976703660945893, 14063388816234683976, 2046622692581829572,
+ 8294969799792017441, 5293778434844788058, 17976364049306763808,
+ 399482430848083948, 16495545010129798933, 15241340958282367519,
+ 989828753826900814, 17616558773874893537, 2471817920909589004,
+ 11764082277667899978, 9618755269550400950, 1240014743757147125,
+ 1887649378641563002, 1842982574728131416, 13243531042427194002,
+ 7688268125537013927, 3080422097287486736, 2562894809975407783,
+ 12428984115620094788, 1355581933694478148, 9895969242586224966,
+ 8628445623963160889, 4298916726468199239, 12773165416305557280,
+ 5240726258301567487, 4975412836403427561, 1842172398579595303,
+ 7812151462958058676, 17974510987263071769, 14980707022065991200,
+ 18294903201142729875, 12911672684850242753, 8979482998667235743,
+ 16808468362384462073, 5981317232108359798, 12373702800369335100,
+ 16119707581920094765, 2782738549717633602, 15454155188515389391,
+ 16495638000603654629, 16348757069342790497, 7769562861984504567,
+ 17504300515449231559, 5557710032938318996, 11846125204788401203,
+ 13957316349928882624, 2738350683717432043, 15738068448047700954,
+ 6224714837294524999, 6081930777706411111, 11366312928059597928,
+ 4355315799925031482, 12393324728734964015, 15277140291994338591,
+ 1406052433297386355, 15859448364509213398, 1672805458341158435,
+ 2926095111610982994, 11056431822276774455, 12083767323511977430,
+ 3296968762229741153, 12312076899982286460, 17769284994682227273,
+ 15349428916826953443, 1056147296359223910, 18305757538706977431,
+ 6214378374180465222, 14279648441175008454, 17791306410319136644,
+ 956593013486324072, 2921235772936241950, 10002890515925652606,
+ 10399654693663712506, 6446247931049971441, 6380465770144534958,
+ 11439178472613251620, 10131486500045494660, 3692642123868351947,
+ 10972816599561388940, 4931112976348785580, 8213967169213816566,
+ 15336469859637867841, 15026830342847689383, 7524668622380765825,
+ 17309937346758783807, 372780684412666438, 5642417144539399955,
+ 18303842993081194577, 11085303253831702827, 15658163165983586950,
+ 8517521928922081563, 16091186344159989860, 17614656488010863910,
+ 4736067146481515156, 13449945221374241354, 17755469346196579408,
+ 13300502638545717375, 6611828134763118043, 14177591906740276597,
+ 9340430243077460347, 7499765399826404087, 3409518087967832469,
+ 9013253864026602045, 4444307427984430192, 3729283608700519712,
+ 13642048880719588383, 16486557958022946240, 2996465014991157904,
+ 10020049344596426576, 12302485648009883778, 8492591321344423126,
+ 17407986443716172520, 10530482934957373052, 15740662350540828750,
+ 1790629986901049436, 6305948377669917188, 15092985352503125323,
+ 928505047232899787, 14404651977039851607, 7564177565277805597,
+ 3411236815351677870, 7752718145953236134, 12315979971311483798,
+ 12477729506691004724, 14654956300924793305, 6689803038918974388,
+ 1540738812233000153, 13508351811701989957, 15864432023192136053,
+ 7990997967273843917, 7424300239290765161, 39585249496300263,
+ 3877436595063283319, 10710642254398044448, 4653804418844456375,
+ 1232267496410380283, 3690525514009038824, 15459770765077428485,
+ 13240346522153894145, 5674964360688390624, 16973644653010587289,
+ 15924280764204855206, 15196708627253442662, 17596174821341373274,
+ 16196745023027393691, 6980050627399795351, 17582264380857746637,
+ 18170372407506856324, 12108126025631005514, 15687749089493373169,
+ 5814107289258228434, 9381977959648494876, 15895601183088112734,
+ 16267869075651604263, 15228381979765852785, 11949618678312581999,
+ 4545324791131029438, 582725409406225185, 15282520250746126790,
+ 14758446535973412711, 7605613563088071833, 1111140641057375915,
+ 5364843095234852245, 218335432181198977, 4891472444796201742,
+ 4564628942836375772, 15500501278323817088, 4913946328556108657,
+ 2684786251736694229, 12090498456116310122, 5310885782157038567,
+ 5032788439854011923, 12627401038822728242, 11869662610126430929,
+ 17650156853043540226, 12126672500118808436, 10437658933435653256,
+ 13133995470637873311, 4601324715591152820, 1874350460376708372,
+ 5808688626286061164, 13777088437302430376, 5018451954762213522,
+ 2588296738534474754, 5503414509154170711, 5230497186769951796,
+ 13261090710400573914, 8515217303152165705, 11074538219737365303,
+ 15481562385740613213, 12705484409881007350, 14221931471178549498,
+ 12905633420087112297, 17337759164357146506, 14081997515778175224,
+ 17384320185513122939, 7131793076779216692, 17483217190312403109,
+ 900692047897995877, 14723287313048560400, 6132094372965340305,
+ 7572797575350925726, 12725160700431903514, 380860122911632449,
+ 1900504978569024571, 8423729759529914138, 7305587201606052334,
+ 12446871355267313320, 4615812356515386206, 3361817115406652303,
+ 17690418922000878428, 14632214537567910559, 2709702289926174775,
+ 3459675155951086144, 7788364399926538150, 16043992474431955950,
+ 15830963823784930267, 4216893617835797954, 538159724689093771,
+ 16029152738918251363, 14444848757576686696, 12941757045272633696,
+ 10900480525147953314, 12547307449905859302, 16001571796892398181,
+ 407942194622690676, 13873235372903944444, 18071603799493008777,
+ 1015646077646778622, 9387605808959554815, 11566702442022019410,
+ 7061722181092883183, 2629032108249254109, 5271820053177594520,
+ 12640880742139693547, 10098688629735675775, 5716304472850923064,
+ 3312674502353063071, 7295926377425759633, 833281439103466115,
+ 16316743519466861667, 9912050326606348167, 11651133878100804242,
+ 18026798122431692459, 6157758321723692663, 4856021830695749349,
+ 7074321707293278978, 10748097797809573561, 2949954440753264783,
+ 9813922580940661152, 9949237950172138336, 15643982711269455885,
+ 16078663425810239127, 12508044395364228880, 12920301578340189344,
+ 15368071871011048915, 1610400750626363239, 11994736084146033126,
+ 6042574085746186088, 4154587549267685807, 15915752367312946034,
+ 1191196620621769193, 467437822242538360, 2836463788873877488,
+ 10476401302029164984, 1716169985450737419, 5327734953288310341,
+ 3994170067185955262, 884431883768190063, 11019001754831208284,
+ 14322807384384895215, 161011537360955545, 1466223959660131656,
+ 5227048585229497539, 12410731857504225031, 2142243279080761103,
+ 17682826799106851430, 1792612570704179953, 14727410295243056025,
+ 1459567192481221274, 5669760721687603135, 17507918443756456845,
+ 10354471145847018200, 10362475129248202288, 13143844410150939443,
+ 6861184673150072028, 18396524361124732580, 543906666394301875,
+ 12476817828199026728, 11853496871128122868, 12747674713108891748,
+ 7986179867749890282, 9158195177777627533, 2217320706811118570,
+ 8631389005200569973, 5538133061362648855, 3369942850878700758,
+ 7813559982698427184, 509051590411815948, 10197035660403006684,
+ 13004818533162292132, 9831652587047067687, 7619315254749630976,
+ 994412663058993407}
+
+ mt := &MersenneTwister64{}
+ mt.SeedSlice([]uint64{0x12345, 0x23456, 0x34567, 0x45678})
+ out := []uint64{}
+ for _ = range expected64 {
+ out = append(out, mt.Uint64())
+ }
+
+ if !reflect.DeepEqual(out, expected64) {
+ t.Errorf("Mersenne Twister 64-bit expected output incorrect")
+ for i := range out {
+ if out[i] != expected64[i] {
+ t.Errorf("Bad at %d: Got %d, want %d", i, out[i], expected64[i])
+ }
+ }
+ }
+}
+
+func TestMersenne64Source(t *testing.T) {
+ r := rand.New(NewMersenneTwister64(8189023))
+ count := 0
+ for i := 0; i < 100; i++ {
+ if r.Intn(2) == 1 {
+ count++
+ }
+ }
+ if math.Abs(float64(count)-50) > 3 {
+ t.Errorf("Expect coin toss about 50-50, was %d", count)
+ }
+}
+
+// Initial run on x86 hardware: 16ns/op
+func BenchmarkMersenne64(b *testing.B) {
+ mt := &MersenneTwister64{}
+ mt.SeedSlice([]uint64{0x12345, 0x23456, 0x34567, 0x45678})
+ b.ResetTimer()
+
+ for i := 0; i < b.N; i++ {
+ mt.Uint64()
+ }
+}
diff --git a/online.go b/online.go
new file mode 100644
index 0000000..d78fa93
--- /dev/null
+++ b/online.go
@@ -0,0 +1,281 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/rand"
+)
+
+// Implemention of Online Codes. See
+// http://cs.nyu.edu/web/Research/TechReports/TR2002-833/TR2002-833.pdf
+// After Maymounkov and Mazieres
+//
+// The input message is a sequence of bytes which is split into N blocks.
+//
+// The constraints on the codec parameters are that 0.55*e*q*N >= q, probably
+// like q*2, which implies that 0.55*e*N >= 1 or e*N >= 1.82.
+//
+// Taking epsilon small (i.e. 0.01), this means N > 200 or so. A more
+// reasonable approach gives N > 400 or so for e=0.01.
+//
+// What this means for small texts is that e should be fairly big, like say
+// 0.3 or something. That means 0.3*N > 4 or so, meaning N is like 12-15.
+// This still allows us to get good probability of recovery. It means source
+// symbols are small, but (e/2)^(q+1) can still be 10^-12 with e=0.3 if
+// q is large enough (like 15).
+//
+// The number of code blocks expected to provide almost certain message
+// recovery is (1+epsilon)(1+0.55*Q*epsilon)N where N is the number of source blocks.
+//
+// For large texts, the primary concern is picking N such that the packet
+// size is convenient. Ideally epsilon is quite small, meaning not much extra
+// data transmission required, and since there are many blocks, that can be
+// done without increasing q very much.
+
+// onlineCodec contains the parameters for application of an online code, typically
+// for a particular known message.
+// Recommended parameters for large NumSourceBlocks: Epsilon=0.01, Quality=3.
+// Note that it requires quite large numbers (that is, thousands) of input blocks
+// to approach optimality. For example, NumSourceBlocks=1000 requires about 3%
+// overhead at these settings to achieve recovery error rate of 1e-8.
+// Implements fountain.Codec
+type onlineCodec struct {
+ // epsilon is the suboptimality parameter. ("Efficiency" or "e")
+ // A message of N blocks can be decoded with high probability
+ // from (1+3*epsilon)*numSourceBlocks received blocks.
+ epsilon float64
+
+ // quality is the decoder quality factor ("q"). This parameter influences the
+ // failure rate of the decoder.
+ // Given (1+3*epsilon)*N blocks, the algorithm will fail with probability
+ // (epsilon/2)^(quality+1)
+ quality int
+
+ // numSourceBlocks is the number of source blocks ("N") to construct from the
+ // input message. This parameter interacts with the message length to set the
+ // packet size, so should be picked with that in mind.
+ numSourceBlocks int
+
+ // randomSeed is a source of randomness for selecting auxiliary encoding blocks.
+ // This seeds a psuedorandom source identically for both encoding and decoding.
+ randomSeed int64
+
+ // cdf is the cumulative distribution function of the degree distribution.
+ cdf []float64
+}
+
+// NewOnlineCodec creates a new encoder for an Online code.
+// epsilon is the suboptimality parameter. ("Efficiency" or "e")
+// A message of N blocks can be decoded with high probability
+// from (1+3*epsilon)*numSourceBlocks received blocks.
+// quality is the decoder quality factor ("q"). This parameter influences the
+// failure rate of the decoder.
+// Given (1+3*epsilon)*N blocks, the algorithm will fail with probability
+// (epsilon/2)^(quality+1)
+// seed is the random seed used to pick auxiliary encoding blocks.
+func NewOnlineCodec(sourceBlocks int, epsilon float64, quality int, seed int64) Codec {
+ return &onlineCodec{
+ epsilon: epsilon,
+ quality: quality,
+ numSourceBlocks: sourceBlocks,
+ randomSeed: seed,
+ cdf: onlineSolitonDistribution(epsilon)}
+}
+
+// SourceBlocks returns the number of source blocks into which the codec will
+// partition an input message.
+func (c *onlineCodec) SourceBlocks() int {
+ return c.numSourceBlocks
+}
+
+// numAuxBlocks returns the number of auxiliary blocks to create for the outer
+// encoding.
+func (c onlineCodec) numAuxBlocks() int {
+ // Note: equation is from the paper.
+ return int(math.Ceil(0.55 * float64(c.quality) * c.epsilon * float64(c.numSourceBlocks)))
+}
+
+// estimateDecodeBlocksNeeded returns a rough lower bound on the number of decode
+// blocks likely needed to successfully decode a message. This number is about
+// (1+epsilon)(NumSourceBlocks + numAuxBlocks)
+func (c onlineCodec) estimateDecodeBlocksNeeded() int {
+ return int(math.Ceil((1 + c.epsilon) * float64(c.numSourceBlocks+c.numAuxBlocks())))
+}
+
+// GenerateIntermediateBlocks finds a set of auxiliary encoding blocks using an
+// LT process, which it then appends to the original set of message blocks.
+func (c *onlineCodec) GenerateIntermediateBlocks(message []byte, numBlocks int) []block {
+ src, aux := generateOuterEncoding(message, *c)
+ intermediate := make([]block, len(src), len(src)+len(aux))
+ copy(intermediate, src)
+ intermediate = append(intermediate, aux...)
+ return intermediate
+}
+
+// generateOuterEncoding creates the source and auxiliary blocks after section
+// 3.1 of http://pdos.csail.mit.edu/~petar/papers/maymounkov-bigdown-lncs.ps
+// Returns two slices of blocks: the original source blocks and the auxiliary
+// blocks.
+// Basic idea: the auxiliary blocks are randomly composed of the source blocks
+// and then used to generate code blocks. This makes recovery of the full
+// original message from code blocks more robust.
+func generateOuterEncoding(message []byte, codec onlineCodec) ([]block, []block) {
+ numAuxBlocks := codec.numAuxBlocks()
+ long, short := partitionBytes(message, codec.numSourceBlocks)
+ source := equalizeBlockLengths(long, short)
+
+ aux := make([]block, numAuxBlocks)
+ // Ensure all aux blocks have the same length as the source blocks,
+ // even if they don't happen to get loaded with data.
+ for i := range aux {
+ aux[i].padding = source[0].length()
+ }
+
+ random := rand.New(NewMersenneTwister(codec.randomSeed))
+ for i := 0; i < codec.numSourceBlocks; i++ {
+ touchAuxBlocks := sampleUniform(random, codec.quality, numAuxBlocks)
+ for _, j := range touchAuxBlocks {
+ aux[j].xor(source[i])
+ }
+ }
+
+ return source, aux
+}
+
+// generateCodeBlock creates a new code symbol, which is the XOR of
+// outer blocks [b_k1, b_k2, b_k3, ... b_kd]
+// Where the sequence k1, k2, k3, ..., kd is provided in the indices.
+func generateCodeBlock(source []block, aux []block, indices []int) block {
+ var symbol block
+
+ for _, i := range indices {
+ if i < len(source) {
+ symbol.xor(source[i])
+ } else {
+ symbol.xor(aux[i-len(source)])
+ }
+ }
+
+ return symbol
+}
+
+// PickIndices finds the source indices for a code block given an ID using
+// the CDF for the online degree distribution.
+func (c *onlineCodec) PickIndices(codeBlockIndex int64) []int {
+ random := rand.New(NewMersenneTwister(codeBlockIndex))
+
+ degree := pickDegree(random, c.cdf)
+ // Pick blocks from the augmented set of original+aux blocks produced
+ // by GenerateIntermediateBlocks.
+ s := sampleUniform(random, degree, c.SourceBlocks()+c.numAuxBlocks())
+ return s
+}
+
+// encodeOnlineBlocks creates a set of online code blocks given the ids provided.
+// An easy way to generate the ids is to pick a pseudo-random sequence and then
+// just grab the first M members of the sequence.
+// The characteristic of an online code is that this method may be called
+// repeatedly with different ids, generating different code blocks for the same
+// message, all of which can then be used interchangeably by the decoder.
+// For each code block, we pick a random set of outer-encoding blocks to XOR
+// to compose it.
+func encodeOnlineBlocks(message []byte, ids []int64, codec onlineCodec) []LTBlock {
+ source, aux := generateOuterEncoding(message, codec)
+ blocks := make([]LTBlock, len(ids))
+ for i := range blocks {
+ indices := codec.PickIndices(ids[i])
+ block := generateCodeBlock(source, aux, indices)
+ blocks[i].BlockCode = ids[i]
+ blocks[i].Data = make([]byte, source[0].length())
+ copy(blocks[i].Data, block.data)
+ }
+ return blocks
+}
+
+// onlineDecoder is the state required for decoding a particular message prepared
+// with the onlineCodec. It must be initialized with the same parameters
+// used for encoding, as well as the expected message length.
+// Implements fountain.Decoder
+type onlineDecoder struct {
+ codec *onlineCodec
+ messageLength int
+
+ // The sparse equation matrix used for decoding.
+ matrix sparseMatrix
+}
+
+// NewDecoder creates an online transform decoder
+func (c *onlineCodec) NewDecoder(messageLength int) Decoder {
+ return newOnlineDecoder(c, messageLength)
+}
+
+// newOnlineDecoder creates a new decoder for a particular message. The codec
+// parameters as well as the original message length must be provided. The
+// decoder is only valid for decoding blocks for a particular source message.
+func newOnlineDecoder(c *onlineCodec, length int) *onlineDecoder {
+ d := &onlineDecoder{codec: c, messageLength: length}
+
+ numAuxBlocks := c.numAuxBlocks()
+ d.matrix.coeff = make([][]int, c.numSourceBlocks+numAuxBlocks)
+ d.matrix.v = make([]block, c.numSourceBlocks+numAuxBlocks)
+
+ // Now we add the initial auxiliary equations into the decode matrix.
+ // These come in as synthetic decode blocks, which have value 0 and
+ // coefficient bits set indicating their constituent outer blocks.
+ auxBlockComposition := make([][]int, numAuxBlocks)
+ random := rand.New(NewMersenneTwister(c.randomSeed))
+ for i := 0; i < c.numSourceBlocks; i++ {
+ touchAuxBlocks := sampleUniform(random, c.quality, numAuxBlocks)
+ for _, j := range touchAuxBlocks {
+ auxBlockComposition[j] = append(auxBlockComposition[j], i)
+ }
+ }
+ for i := range auxBlockComposition {
+ auxBlockComposition[i] = append(auxBlockComposition[i], i+c.numSourceBlocks)
+ }
+ // Note: these composition slices are guaranteed sorted since we added constituent
+ // source blocks in order, followed by the aux block index. So we can now just
+ // add them to the equation matrix as if they were received.
+
+ for i := range auxBlockComposition {
+ d.matrix.addEquation(auxBlockComposition[i], block{})
+ }
+
+ return d
+}
+
+// AddBlocks adds a set of encoded blocks to the decoder. Returns true if the
+// message can be fully decoded. False if there is insufficient information.
+func (d *onlineDecoder) AddBlocks(blocks []LTBlock) bool {
+ for i := range blocks {
+ indices := d.codec.PickIndices(blocks[i].BlockCode)
+ d.matrix.addEquation(indices, block{data: blocks[i].Data})
+ }
+ return d.matrix.determined()
+}
+
+// Decode extracts the decoded message from the decoder. If the decoder does
+// not have sufficient information to produce an output, returns a nil slice.
+func (d *onlineDecoder) Decode() []byte {
+ if !d.matrix.determined() {
+ return nil
+ }
+
+ d.matrix.reduce()
+
+ lenLong, lenShort, numLong, numShort := partition(d.messageLength, d.codec.numSourceBlocks)
+ return d.matrix.reconstruct(d.messageLength, lenLong, lenShort, numLong, numShort)
+}
diff --git a/online_test.go b/online_test.go
new file mode 100644
index 0000000..d2ecded
--- /dev/null
+++ b/online_test.go
@@ -0,0 +1,175 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+func printEncoding(source []block, aux []block, t *testing.T) {
+ t.Log("Outer Encoding Blocks")
+ t.Log("---------------------")
+ kb := 0
+ for s := range source {
+ t.Log("src", kb, source[s].data)
+ kb++
+ }
+ for a := range aux {
+ t.Log("aux", kb, aux[a].data)
+ kb++
+ }
+}
+
+func TestOnlineBlocks(t *testing.T) {
+ c := onlineCodec{epsilon: 0.01, quality: 5, numSourceBlocks: 6, randomSeed: 200}
+ if c.numAuxBlocks() != 1 {
+ t.Errorf("Got %d aux blocks, want 1", c.numAuxBlocks())
+ }
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+
+ source, aux := generateOuterEncoding(message, c)
+ printEncoding(source, aux, t)
+
+ if source[0].data[0] != 'a' {
+ t.Errorf("Source data should start with message beginning. Got %s", source[0].data)
+ }
+
+ block := generateCodeBlock(source, aux, []int{4})
+ if !reflect.DeepEqual(block.data, source[4].data) {
+ t.Errorf("Single data block is %v, should be %v", block.data, source[4].data)
+ }
+ block = generateCodeBlock(source, aux, []int{2, 5, 6})
+ if block.data[0] != 107^119^7 {
+ t.Errorf("XOR data block got %v, should be %v", block.data[0], 107^119^7)
+ }
+ t.Log("block =", block)
+
+ codec := NewOnlineCodec(6, 0.01, 5, 200)
+ ltblocks := EncodeLTBlocks(message, []int64{252}, codec)
+ indices := codec.PickIndices(252)
+ if !reflect.DeepEqual(indices, []int{4}) {
+ t.Errorf("Indices for 252 are %v, should be [4]", indices)
+ }
+ if !reflect.DeepEqual(ltblocks[0].Data, source[4].data) {
+ t.Errorf("Single data block is %v, should be %v", ltblocks[0].Data, source[4].data)
+ }
+ t.Log("block =", ltblocks[0])
+}
+
+func TestDecoder(t *testing.T) {
+ c := NewOnlineCodec(13, 0.3, 10, 200).(*onlineCodec)
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ ids := make([]int64, 45)
+ random := rand.New(rand.NewSource(8923489))
+ for i := range ids {
+ ids[i] = int64(random.Intn(100000))
+ }
+ source, aux := generateOuterEncoding(message, *c)
+ printEncoding(source, aux, t)
+
+ blocks := encodeOnlineBlocks(message, ids, *c)
+ t.Log("blocks =", blocks)
+
+ d := newOnlineDecoder(c, len(message))
+
+ for i := 0; i < 16; i++ {
+ d.AddBlocks([]LTBlock{blocks[i]})
+ if testing.Verbose() {
+ printMatrix(d.matrix, t)
+ }
+ }
+
+ d.matrix.reduce()
+ t.Log("REDUCE")
+ printMatrix(d.matrix, t)
+
+ decoded := d.Decode()
+ printMatrix(d.matrix, t)
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Got %v, want %v", decoded, message)
+ }
+}
+
+func TestDecoderBlockTable(t *testing.T) {
+ c := NewOnlineCodec(13, 0.3, 10, 0).(*onlineCodec)
+ if c.numAuxBlocks() != 22 {
+ t.Errorf("Got %d aux blocks, want 22", c.numAuxBlocks())
+ }
+ needed := c.estimateDecodeBlocksNeeded()
+ if needed != 46 {
+ t.Errorf("Got %d blocks expected to be needed, want 17", needed)
+ }
+
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ random := rand.New(rand.NewSource(8234923))
+
+ moreBlocksNeeded := 0
+ for i := 0; i < 100; i++ {
+ c.randomSeed = random.Int63()
+ r := rand.New(rand.NewSource(random.Int63()))
+ ids := make([]int64, 45)
+ for i := range ids {
+ ids[i] = int64(r.Intn(100000))
+ }
+ blocks := encodeOnlineBlocks(message, ids, *c)
+
+ d := newOnlineDecoder(c, len(message))
+ d.AddBlocks(blocks[0:30])
+ if !d.matrix.determined() {
+ moreBlocksNeeded++
+ d.AddBlocks(blocks[31:46])
+ }
+ decoded := d.Decode()
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Got %v, want %v", decoded, message)
+ }
+ }
+
+ if moreBlocksNeeded > 2 {
+ t.Errorf("Needed too many high-block-count decoding sequences: %d", moreBlocksNeeded)
+ }
+}
+
+func TestDecodeMessageTable(t *testing.T) {
+ c := NewOnlineCodec(10, 0.2, 7, 0).(*onlineCodec)
+ random := rand.New(rand.NewSource(8234982))
+ for i := 0; i < 100; i++ {
+ c.randomSeed = random.Int63()
+ r := rand.New(rand.NewSource(random.Int63()))
+ messageLen := r.Intn(1000) + 1000
+ message := make([]byte, messageLen)
+ for j := 0; j < len(message); j++ {
+ message[j] = byte(r.Intn(200))
+ }
+ ids := make([]int64, 50)
+ for i := range ids {
+ ids[i] = int64(r.Intn(100000))
+ }
+ blocks := encodeOnlineBlocks(message, ids, *c)
+
+ d := newOnlineDecoder(c, len(message))
+ d.AddBlocks(blocks[0:25])
+ if !d.matrix.determined() {
+ t.Errorf("Message should be determined after 25 blocks")
+ } else {
+ decoded := d.Decode()
+ if !reflect.DeepEqual(decoded, message) {
+ t.Errorf("Incorrect message decode. Length=%d", len(message))
+ }
+ }
+ }
+}
diff --git a/raptor.go b/raptor.go
new file mode 100644
index 0000000..3dfa3a4
--- /dev/null
+++ b/raptor.go
@@ -0,0 +1,370 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "sort"
+)
+
+// The Raptor fountain code (also called the R10 code) from RFC 5053.
+// This code nearly matches the performance of the random binary fountain, but
+// does so with a hard cap on the degree distribution on code blocks, meaning
+// that the decode process becomes linear instead of quadratic.
+//
+// In addition, it is a systematic code, meaning that the original source blocks
+// are part of the encoding. This enables the source blocks to be sent simply,
+// and then repair blocks constructed as needed using the code.
+//
+// A limitation is that the code supports a maximum of 8192 source blocks,
+// thus requiring very large source messages to be split up into sub-messages if
+// smaller packet sizes are a goal. Performance varies from the random fountain
+// the most for higher loss rates and smaller numbers of source blocks. A reasonable
+// expectation is that the encoding overhead due to using the code is a few percent.
+//
+// When encoding the message, the codec takes a list of block IDs and returns a
+// set of raptor coded blocks generated according to the RFC. The contract in the
+// RFC is that ids from 0 to K (codec.NumSourceSymbols) are identical to the source
+// symbols. (That's what makes it a systematic code.) Such symbols should not be in
+// the same packet as repair symbols with id >= K.
+//
+// A typical usage in a transmission system might be to just split the message and
+// send the first K symbols normally. Then create the codec and generate repair
+// symbols using random ESI values >= K until the message is reconstructed by
+// the receiver.
+//
+// The BlockCode in the resulting LTBlocks will be a uint16-compatible value.
+//
+// IMPORTANT NOTE: encoding is destructive to the input message.
+
+// raptorCodec describes the parameters needed to construct a raptor code. The codec
+// governs the production of an unbounded set of LTBlocks from a given source message.
+// If the total transfer size is such that the application wants to split the message
+// into sub-blocks before transfer, that should be done independently of this codec
+// in accordance with the instructions in RFC 5053. For many uses, there will be a
+// value of NumSourceSymbols and SymbolAlignmentSize such that the resulting LTBlock
+// size is of an acceptable length. Note: the RFC provides a way to pack multiple
+// symbols per packet. If packet erasure is the loss model for the code, this
+// ends up behaving like choosing a smaller NumSourceSymbols, in which case it's
+// simpler to pass just one code symbol per packet.
+// Implements fountain.Codec
+type raptorCodec struct {
+ // SymbolAlignmentSize = Al is the size of each symbol in the source message in bytes.
+ // Usually 4. This is the XOR granularity in bytes. On 32-byte machines 4-byte XORs
+ // will be most efficient. On the other hand, the code will perform with less overhead
+ // with larger numbers of source blocks.
+ SymbolAlignmentSize int
+
+ // NumSourceSymbols = K. Must be in the range [4, 8192] (inclusive). This is
+ // how many source symbols the input message will be divided into. If NumSourceSymbols
+ // doesn't evenly divide the length of the message in units of SymbolAlignmentSize,
+ // there will be null padding applied to the block.
+ NumSourceSymbols int
+}
+
+// NewRaptorCodec creates a new R10 raptor codec using the provided number of
+// source blocks and alignment size.
+func NewRaptorCodec(sourceBlocks int, alignmentSize int) Codec {
+ return &raptorCodec{
+ NumSourceSymbols: sourceBlocks,
+ SymbolAlignmentSize: alignmentSize}
+}
+
+// SourceBlocks returns the number of source symbols used by the codec.
+func (c *raptorCodec) SourceBlocks() int {
+ return c.NumSourceSymbols
+}
+
+// RAND function from section 5.4.4.1
+// x, i should be non-negative, m positive.
+// Produces a pseudo-random value in the range [0, m-1]
+func raptorRand(x, i, m uint32) uint32 {
+ v0 := v0table[(x+i)%256]
+ v1 := v1table[((x/256)+i)%256]
+ return (v0 ^ v1) % m
+}
+
+// Deg function from section 5.4.4.2
+// deg calculates the degree to be used in code block generation.
+func deg(v uint32) int {
+ f := [...]uint32{0, 10241, 491582, 712794, 831695, 948446, 1032189, 1048576}
+ d := [...]int{0, 1, 2, 3, 4, 10, 11, 40}
+
+ for j := 1; j < len(f)-1; j++ {
+ if v < f[j] {
+ return d[j]
+ }
+ }
+
+ return d[len(d)-1]
+}
+
+// From RFC section 5.4.2.3 This function computes L, S, and H from K.
+// K is the number of source symbols (limited to 2**16).
+// The return values are:
+// L is the number of intermediate symbols desired (K+S+H), followed by
+// S, the number of LDPC symbols, followed by
+// H, the number of half-symbols.
+func intermediateSymbols(k int) (int, int, int) {
+ // X is the smallest positive integer such that X*(X-1) >= 2*K
+ x := int(math.Floor(math.Sqrt(2 * float64(k))))
+ if x < 1 {
+ x = 1
+ }
+
+ for (x * (x - 1)) < (2 * k) {
+ x++
+ }
+
+ // S is the smallest prime such that S >= ceil(0.01*K) + X
+ s := int(math.Ceil(0.01*float64(k))) + x
+ s = smallestPrimeGreaterOrEqual(s)
+
+ // H is the smallest integer such that choose(H, ceil(H/2)) >= K + S
+ // choose(h, h/2) <= 4^(h/2), so begin with
+ // h/2 ln(4) = ln K+S
+ // h = ln(K+S)/ln(4)
+ h := int(math.Floor(math.Log(float64(s)+float64(k)) / math.Log(4)))
+ for centerBinomial(h) < k+s {
+ h++
+ }
+
+ return k + s + h, s, h
+}
+
+// Triple generator from RFC section 5.4.4.4
+// k is the number of source symbols.
+// x is the (random) code symbol ID.
+// The generator creates values (d, a, b) to be used in constructing intermediate blocks.
+func tripleGenerator(k int, x uint16) (int, uint32, uint32) {
+ l, _, _ := intermediateSymbols(k)
+ lprime := smallestPrimeGreaterOrEqual(l)
+ q := uint32(65521) // largest prime < 2^16
+ jk := uint32(systematicIndextable[k])
+
+ a := uint32((53591 + (uint64(jk) * 997)) % uint64(q))
+ b := (10267 * (jk + 1)) % q
+ y := uint32((uint64(b) + (uint64(x) * uint64(a))) % uint64(q))
+ v := raptorRand(y, 0, 1048576) // 1048576 == 2^20
+ d := deg(v)
+ a = 1 + raptorRand(y, 1, uint32(lprime-1))
+ b = raptorRand(y, 2, uint32(lprime))
+
+ return d, a, b
+}
+
+// findLTIndices discovers the composition of the ESI=x LT code block for a
+// raptor code. k is the number of source blocks.
+func findLTIndices(k int, x uint16) []int {
+ l, _, _ := intermediateSymbols(k)
+ lprime := uint32(smallestPrimeGreaterOrEqual(l))
+ d, a, b := tripleGenerator(k, x)
+
+ if d > l {
+ d = l
+ }
+
+ indices := make([]int, 0)
+ for b >= uint32(l) {
+ b = (b + a) % lprime
+ }
+ indices = append(indices, int(b))
+
+ for j := 1; j < d; j++ {
+ b = (b + a) % lprime
+ for b >= uint32(l) {
+ b = (b + a) % lprime
+ }
+ indices = append(indices, int(b))
+ }
+
+ sort.Ints(indices)
+ return indices
+}
+
+// ltEncode is the LT encoding function. RFC section 5.4.4.3
+// c is the intermediate symbol vector, k is the number of source symbols.
+// x is the symbol ID we are generating.
+// The output is an code block containing the bytes of that symbol.
+func ltEncode(k int, x uint16, c []block) block {
+ indices := findLTIndices(k, x)
+
+ result := block{}
+ for _, i := range indices {
+ result.xor(c[i])
+ }
+
+ return result
+}
+
+// raptorIntermediateBlocks takes the source blocks and follows the algorithm
+// described in the RFC to generate the intermediate encoding which is then used
+// as the LT source encoding for the systematic code. The properties of the intermediate
+// encoding is that we will create L resulting blocks from K source blocks.
+// The first K intermediate blocks have an LT relationship with the K source blocks --
+// we basically use an unsystematic LT code using the tripleGenerator values where
+// the X "ESI" is simply the source block index. We then generate S+H additional
+// intermediate blocks.
+//
+// The S blocks are a low density parity check group of blocks which have contributions
+// from three of the first K intermediate blocks arranged in successive clusters.
+// This low density arrangement cycles column-wise in the decode matrix, meaning
+// that there is extra decode-friendly information in the decode matrix to make the
+// expected decoding operations use a sparse enough matrix so that it is efficient.
+//
+// The H blocks are composed in the decode matrix of many (half) of the first K
+// intermediate blocks. The coding is taken from a gray code, and is intended to
+// be more-or-less equivalent to a random binary fountain from the coding
+// perspective, so when these blocks are chosen in the coded blocks, it ensures
+// good performance by making sure there are no unrepresented source blocks.
+//
+// We then decode this -- the J(K) systematic index values are chosen such that
+// the first L=K+S+H members of this set will produce an invertible decode matrix.
+// The guarantees that the re-encoding of the first K code symbols are equal to
+// the source symbols (ensuring a systematic code), and that the analysis of the
+// overall code will follow the performance of the random binary fountain closely.
+//
+// This method is destructive to the source blocks.
+func raptorIntermediateBlocks(source []block) []block {
+ ltdecoder := newRaptorDecoder(&raptorCodec{SymbolAlignmentSize: 1,
+ NumSourceSymbols: len(source)}, 1)
+ for i := 0; i < len(source); i++ {
+ indices := findLTIndices(len(source), uint16(i))
+ ltdecoder.matrix.addEquation(indices, source[i])
+ }
+
+ ltdecoder.matrix.reduce()
+
+ // panics if ~ltdecoder.determined. The J(K) selection should ensure that
+ // never happens.
+ intermediate := ltdecoder.matrix.v
+ return intermediate
+}
+
+// GenerateIntermediateBlocks creates the pre-code representation given the
+// message argument blocks. For the raptor code, this pre-code is generated by
+// a reverse-coding process which ensures that for BlockCode=0, the 0th block of
+// the incoming message is produced, and so on up to the 'len(message)-1'th BlockCode.
+func (c *raptorCodec) GenerateIntermediateBlocks(message []byte, numBlocks int) []block {
+ sourceLong, sourceShort := partitionBytes(message, numBlocks)
+ source := equalizeBlockLengths(sourceLong, sourceShort)
+ return raptorIntermediateBlocks(source)
+}
+
+// PickIndices chooses a set of indices for the provided CodeBlock index value
+// which are used to compose an LTBlock. It functions by
+func (c *raptorCodec) PickIndices(codeBlockIndex int64) []int {
+ return findLTIndices(int(c.SourceBlocks()), uint16(codeBlockIndex))
+}
+
+// NewDecoder creates a new raptor decoder
+func (c *raptorCodec) NewDecoder(messageLength int) Decoder {
+ return newRaptorDecoder(c, messageLength)
+}
+
+// raptorDecoder is the state required for decoding a particular message prepared
+// with the Raptor code. It must be initialized with the same raptorCodec parameters
+// used for encoding, as well as the expected message length.
+type raptorDecoder struct {
+ codec raptorCodec
+ messageLength int
+
+ // The sparse equation matrix used for decoding.
+ matrix sparseMatrix
+}
+
+// newRaptorDecoder creates a new raptor decoder for a given message. The
+// codec supplied must be the same one as the message was encoded with.
+func newRaptorDecoder(c *raptorCodec, length int) *raptorDecoder {
+ d := &raptorDecoder{codec: *c, messageLength: length}
+
+ l, s, h := intermediateSymbols(c.NumSourceSymbols)
+
+ // Add the S + H intermediate symbol composition equations.
+ d.matrix.coeff = make([][]int, l)
+ d.matrix.v = make([]block, l)
+
+ k := c.NumSourceSymbols
+ compositions := make([][]int, s)
+
+ for i := 0; i < k; i++ {
+ a := 1 + (int(math.Floor(float64(i)/float64(s))) % (s - 1))
+ b := i % s
+ compositions[b] = append(compositions[b], i)
+ b = (b + a) % s
+ compositions[b] = append(compositions[b], i)
+ b = (b + a) % s
+ compositions[b] = append(compositions[b], i)
+ }
+ for i := 0; i < s; i++ {
+ compositions[i] = append(compositions[i], k+i)
+ d.matrix.addEquation(compositions[i], block{})
+ }
+
+ compositions = make([][]int, h)
+
+ hprime := int(math.Ceil(float64(h) / 2))
+ m := buildGraySequence(k+s, hprime)
+ for i := 0; i < h; i++ {
+ for j := 0; j < k+s; j++ {
+ if bitSet(uint(m[j]), uint(i)) {
+ compositions[i] = append(compositions[i], j)
+ }
+ }
+ compositions[i] = append(compositions[i], k+s+i)
+ d.matrix.addEquation(compositions[i], block{})
+ }
+
+ return d
+}
+
+// AddBlocks adds a set of encoded blocks to the decoder. Returns true if the
+// message can be fully decoded. False if there is insufficient information.
+func (d *raptorDecoder) AddBlocks(blocks []LTBlock) bool {
+ for i := range blocks {
+ indices := findLTIndices(d.codec.NumSourceSymbols, uint16(blocks[i].BlockCode))
+ d.matrix.addEquation(indices, block{data: blocks[i].Data})
+ }
+ return d.matrix.determined()
+}
+
+// Decode extracts the decoded message from the decoder. If the decoder does
+// not have sufficient information to produce an output, returns a nil slice.
+func (d *raptorDecoder) Decode() []byte {
+ if !d.matrix.determined() {
+ return nil
+ }
+
+ d.matrix.reduce()
+
+ // Now the intermediate blocks are held in d.matrix.v. Use the encoder function
+ // to recover the source blocks.
+ intermediate := d.matrix.v
+ source := make([]block, d.codec.NumSourceSymbols)
+ for i := 0; i < d.codec.NumSourceSymbols; i++ {
+ source[i] = ltEncode(d.codec.NumSourceSymbols, uint16(i), intermediate)
+ }
+
+ lenLong, lenShort, numLong, numShort := partition(d.messageLength, d.codec.NumSourceSymbols)
+ out := make([]byte, d.messageLength)
+ out = out[0:0]
+ for i := 0; i < numLong; i++ {
+ out = append(out, source[i].data[0:lenLong]...)
+ }
+ for i := numLong; i < numLong+numShort; i++ {
+ out = append(out, source[i].data[0:lenShort]...)
+ }
+ return out
+}
diff --git a/raptor_test.go b/raptor_test.go
new file mode 100644
index 0000000..5858f4f
--- /dev/null
+++ b/raptor_test.go
@@ -0,0 +1,323 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+func TestRaptorRand(t *testing.T) {
+ var randTests = []struct {
+ x uint32
+ i uint32
+ m uint32
+ r uint32
+ }{
+ {1, 4, 150, 50},
+ {20005, 19, 25, 6},
+ {2180, 11, 1383483, 1166141},
+ }
+
+ for _, test := range randTests {
+ if test.r != raptorRand(test.x, test.i, test.m) {
+ t.Errorf("raptorRand(%d, %d, %d) = %d, should be %d",
+ test.x, test.i, test.m, raptorRand(test.x, test.i, test.m), test.r)
+ }
+ }
+}
+
+func TestDeg(t *testing.T) {
+ var degreeTests = []struct {
+ x uint32
+ d int
+ }{
+ {0, 1},
+ {10000, 1},
+ {10240, 1},
+ {10241, 2},
+ {10242, 2},
+ {715000, 4},
+ {1000000, 11},
+ {1034300, 40},
+ {1048575, 40},
+ {1048576, 40},
+ }
+
+ for _, test := range degreeTests {
+ if test.d != deg(test.x) {
+ t.Errorf("deg(%d) = %d, should be %d", test.x, deg(test.x), test.d)
+ }
+ }
+}
+
+func TestIntermediateSymbols(t *testing.T) {
+ var intermediateTests = []struct {
+ k int
+ l int
+ s int
+ h int
+ }{
+ {0, 4, 2, 2},
+ {1, 8, 3, 4},
+ {10, 23, 7, 6}, // from a Luby, Shokrollahi paper
+ {13, 26, 7, 6},
+ {14, 28, 7, 7},
+ {500, 553, 41, 12},
+ {5000, 5166, 151, 15},
+ }
+
+ for _, test := range intermediateTests {
+ l, s, h := intermediateSymbols(test.k)
+ if l != test.l || s != test.s || h != test.h {
+ t.Errorf("intermediateSymbols(%d) = (%d, %d, %d), should be %d, %d, %d",
+ test.k, l, s, h, test.l, test.s, test.h)
+ }
+ }
+}
+
+func TestTripleGenerator(t *testing.T) {
+ var tripleTests = []struct {
+ k int
+ x uint16
+ d int
+ a uint32
+ b uint32
+ }{
+ {0, 3, 2, 4, 3},
+ {1, 4, 4, 2, 5},
+ {4, 0, 10, 13, 1},
+ {4, 4, 4, 6, 2},
+ {500, 514, 2, 107, 279},
+ {1000, 52918, 3, 1070, 121},
+ }
+
+ for _, test := range tripleTests {
+ d, a, b := tripleGenerator(test.k, test.x)
+ if d != test.d || a != test.a || b != test.b {
+ t.Errorf("tripleGenerator(%d, %d) = (%d, %d, %d), should be %d, %d, %d",
+ test.k, test.x, d, a, b, test.d, test.a, test.b)
+ }
+ }
+}
+
+func TestSystematicIndices(t *testing.T) {
+ if systematicIndextable[4] != 18 {
+ t.Errorf("Systematic index for 4 was %d, must be 18", systematicIndextable[4])
+ }
+ if systematicIndextable[21] != 2 {
+ t.Errorf("Systematic index for 4 was %d, must be 2", systematicIndextable[4])
+ }
+ if systematicIndextable[8192] != 2665 {
+ t.Errorf("Systematic index for 4 was %d, must be 2665", systematicIndextable[4])
+ }
+}
+
+func TestLTIndices(t *testing.T) {
+ var ltIndexTests = []struct {
+ k int
+ x uint16
+ indices []int
+ }{
+ {4, 0, []int{1, 2, 3, 4, 6, 7, 8, 10, 11, 12}},
+ {4, 4, []int{2, 3, 8, 9}},
+ {100, 1, []int{51, 104}},
+ {1000, 727, []int{306, 687, 1040}},
+ {10, 57279, []int{19, 20, 21, 22}},
+ }
+
+ for _, test := range ltIndexTests {
+ indices := findLTIndices(test.k, test.x)
+ if !reflect.DeepEqual(indices, test.indices) {
+ t.Errorf("findLTIndices(%d, %d) = %v, should be %v",
+ test.k, test.x, indices, test.indices)
+ }
+ }
+}
+
+func TestRaptorDecoderConstruction(t *testing.T) {
+ decoder := newRaptorDecoder(&raptorCodec{SymbolAlignmentSize: 1,
+ NumSourceSymbols: 10}, 1)
+ printMatrix(decoder.matrix, t)
+ // From the first row of the constraint matrix. Test vectors from a paper by
+ // Luby and Shokrollahi.
+ if !reflect.DeepEqual(decoder.matrix.coeff[0], []int{0, 5, 6, 7, 10}) {
+ t.Errorf("First matrix equation was %v, should be {0, 5, 6, 7, 10}",
+ decoder.matrix.coeff[0])
+ }
+ // Fourth row
+ if !reflect.DeepEqual(decoder.matrix.coeff[1], []int{1, 2, 3, 8, 13}) {
+ t.Errorf("Second matrix equation was %v, should be {1, 2, 3, 8, 13}",
+ decoder.matrix.coeff[0])
+ }
+ // Fifth row
+ if !reflect.DeepEqual(decoder.matrix.coeff[2], []int{2, 3, 4, 7, 9, 14}) {
+ t.Errorf("Third matrix equation was %v, should be {2, 3, 4, 7, 9, 14}",
+ decoder.matrix.coeff[0])
+ }
+}
+
+func printIntermediateEncoding(intermediate []block, t *testing.T) {
+ t.Log("Intermediate Encoding Blocks")
+ t.Log("----------------------------")
+ kb := 0
+ for s := range intermediate {
+ t.Log("intermediate", kb, intermediate[s].data)
+ kb++
+ }
+}
+
+func TestIntermediateBlocks(t *testing.T) {
+ blocks := [4]block{
+ {data: []byte{0, 0, 0, 1}},
+ {data: []byte{0, 0, 1, 0}},
+ {data: []byte{0, 1, 0, 0}},
+ {data: []byte{1, 0, 0, 0}},
+ }
+
+ srcBlocks := make([]block, len(blocks))
+ for i := range blocks {
+ srcBlocks[i].xor(blocks[i])
+ }
+ intermediate := raptorIntermediateBlocks(srcBlocks) // destructive to srcBlocks
+ if len(intermediate) != 14 {
+ t.Errorf("Length of intermediate blocks is %d, should be 14", len(intermediate))
+ }
+ printIntermediateEncoding(intermediate, t)
+
+ t.Log("Finding intermediate equations")
+ // test that intermediate ltEnc equations hold
+ for i := 0; i < 4; i++ {
+ block := ltEncode(4, uint16(i), intermediate)
+ if !reflect.DeepEqual(block, blocks[i]) {
+ t.Errorf("The result of LT encoding on the intermediate blocks for "+
+ "block %d is %v, should be the source blocks %v", i,
+ block.data, blocks[i].data)
+ }
+ }
+}
+
+func TestSystematicRaptorCode(t *testing.T) {
+ c := NewRaptorCodec(13, 2)
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ blocks := c.GenerateIntermediateBlocks(message, c.SourceBlocks())
+
+ messageCopy := []byte("abcdefghijklmnopqrstuvwxyz")
+ sourceLong, sourceShort := partitionBytes(messageCopy, c.SourceBlocks())
+ sourceCopy := equalizeBlockLengths(sourceLong, sourceShort)
+
+ for _, testIndex := range []int64{0, 1, 2, 3, 4, 5} {
+ t.Logf("Testing %d", testIndex)
+ b := ltEncode(13, uint16(testIndex), blocks)
+ if !reflect.DeepEqual(b.data, sourceCopy[testIndex].data) {
+ t.Errorf("LT encoding of CodeBlock=%d was (%v), should be the %d'th source block (%v)",
+ testIndex, b.data, testIndex, sourceCopy[testIndex].data)
+ }
+ }
+}
+
+func TestIntermediateBlocks13(t *testing.T) {
+ blocks := make([]block, 13)
+ for i := range blocks {
+ blocks[i].data = make([]byte, 13)
+ blocks[i].data[i] = 1
+ }
+
+ srcBlocks := make([]block, 13)
+ for i := range blocks {
+ srcBlocks[i].xor(blocks[i])
+ }
+ intermediate := raptorIntermediateBlocks(srcBlocks) // destructive to srcBlocks
+ if len(intermediate) != 26 {
+ t.Errorf("Length of intermediate blocks is %d, should be 26", len(intermediate))
+ }
+ printIntermediateEncoding(intermediate, t)
+
+ t.Log("Finding intermediate equations")
+ // test that intermediate ltEnc equations hold
+ for i := 0; i < 13; i++ {
+ block := ltEncode(13, uint16(i), intermediate)
+ if !reflect.DeepEqual(block, blocks[i]) {
+ t.Errorf("The result of LT encoding on the intermediate blocks for "+
+ "block %d is %v, should be the source blocks %v", i,
+ block.data, blocks[i].data)
+ }
+ }
+
+ ids := make([]int64, 45)
+ random := rand.New(rand.NewSource(8923489))
+ for i := range ids {
+ ids[i] = int64(random.Intn(60000))
+ }
+ c := raptorCodec{NumSourceSymbols: 13, SymbolAlignmentSize: 13}
+ srcBlocks = make([]block, 13)
+ for i := range blocks {
+ srcBlocks[i].xor(blocks[i])
+ }
+
+ message := make([]byte, 0)
+ for i := range blocks {
+ message = append(message, blocks[i].data...)
+ }
+ codeBlocks := EncodeLTBlocks(message, ids, &c) // destructive to srcBlocks
+
+ t.Log("DECODE--------")
+ decoder := newRaptorDecoder(&c, len(message))
+ for i := 0; i < 17; i++ {
+ decoder.AddBlocks([]LTBlock{codeBlocks[i]})
+ }
+
+ message = make([]byte, 0)
+ for i := range blocks {
+ message = append(message, blocks[i].data...)
+ }
+
+ if decoder.matrix.determined() {
+ t.Log("Recovered:\n", decoder.matrix.v)
+ out := decoder.Decode()
+ if !reflect.DeepEqual(message, out) {
+ t.Errorf("Decoding result must equal %v, got %v", message, out)
+ }
+ }
+}
+
+func TestRaptorCodec(t *testing.T) {
+ c := NewRaptorCodec(13, 2)
+ message := []byte("abcdefghijklmnopqrstuvwxyz")
+ ids := make([]int64, 45)
+ random := rand.New(rand.NewSource(8923489))
+ for i := range ids {
+ ids[i] = int64(random.Intn(60000))
+ }
+
+ messageCopy := make([]byte, len(message))
+ copy(messageCopy, message)
+
+ codeBlocks := EncodeLTBlocks(messageCopy, ids, c)
+
+ t.Log("DECODE--------")
+ decoder := newRaptorDecoder(c.(*raptorCodec), len(message))
+ for i := 0; i < 17; i++ {
+ decoder.AddBlocks([]LTBlock{codeBlocks[i]})
+ }
+ if decoder.matrix.determined() {
+ t.Log("Recovered:\n", decoder.matrix.v)
+ out := decoder.Decode()
+ if !reflect.DeepEqual(message, out) {
+ t.Errorf("Decoding result must equal %s, got %s", string(message), string(out))
+ }
+ }
+}
diff --git a/ru10.go b/ru10.go
new file mode 100644
index 0000000..84e13c8
--- /dev/null
+++ b/ru10.go
@@ -0,0 +1,210 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/rand"
+ "sort"
+)
+
+// The RU10 fountain is an unsystematic(*) fountain code which uses a degree
+// distribution and intermediate block generation scheme similar to the
+// R10 (Raptor) code. The unsystematic code has much lower independent block
+// generation setup cost -- it requires only generating the auxiliary blocks,
+// and does not require the full decode run of the algorithm to generate the
+// intermediate encoding. This makes it cheaper in two ways for parallel
+// execution. If all receivers can be expected to have the decoder in place,
+// there is no decoding advantage to requiring the source blocks to appear in the
+// code block set. But more to the point, if multiple encoders are working from
+// the same source blocks, they can generate code blocks independently without
+// having to worry if they are re-transmitting source blocks. Since the encoder
+// state is cheaper and code blocks interchangeable, we sacrifice a little
+// performance for that parallelizability and statelessness. We also gain a little
+// in that this encoder variant needn't rely on the systematic number table
+// of the R10 codec, and without that constraint, the number of possible code block
+// ESI IDs is basically unlimited.
+//
+// (*) Well, not by design at least.
+
+// This triple generator uses the Mersenne Twister to generate random seeds.
+// k is the number of source symbols.
+// x is the (random) code symbol ID.
+// The generator creates values (d, a, b) to be used in constructing intermediate blocks.
+func ru10TripleGenerator(k int, x int64) (int, uint32, uint32) {
+ l, _, _ := intermediateSymbols(k)
+ lprime := smallestPrimeGreaterOrEqual(l)
+
+ // TODO(gbillock): nudge x as a function of k to get better overhead-failure curve?
+ rand := rand.New(NewMersenneTwister64(x))
+
+ v := uint32(rand.Int63() % 1048576)
+ a := uint32(1 + (rand.Int63() % int64(lprime - 1)))
+ b := uint32(rand.Int63() % int64(lprime))
+ d := deg(v)
+
+ return d, a, b
+}
+
+// ru10Codec implements the Raptor-alike fountain code.
+// Implements fountain.Codec.
+type ru10Codec struct {
+ numSourceSymbols int
+
+ symbolAlignmentSize int
+}
+
+// NewRU10Codec creates an unsystematic raptor-like fountain codec which uses an
+// intermediate block generation algorithm similar to the Raptor R10 codec.
+func NewRU10Codec(numSourceSymbols int, symbolAlignmentSize int) Codec {
+ return &ru10Codec{
+ numSourceSymbols: numSourceSymbols,
+ symbolAlignmentSize: symbolAlignmentSize}
+}
+
+// SourceBlocks returns the number of source blocks the codec uses in the
+// source message plus intermediate blocks added.
+func (c *ru10Codec) SourceBlocks() int {
+ return c.numSourceSymbols
+}
+
+// PickIndices uses the R10 distribution function to pick indices. It gets
+// numbers from the triple generator.
+func (c *ru10Codec) PickIndices(codeBlockIndex int64) []int {
+ d, a, b := ru10TripleGenerator(c.numSourceSymbols, codeBlockIndex)
+ l, _, _ := intermediateSymbols(c.numSourceSymbols)
+ lprime := uint32(smallestPrimeGreaterOrEqual(l))
+
+ if d > l {
+ d = l
+ }
+
+ indices := make([]int, 0)
+ for b >= uint32(l) {
+ b = (b + a) % lprime
+ }
+ indices = append(indices, int(b))
+
+ for j := 1; j < d; j++ {
+ b = (b + a) % lprime
+ for b >= uint32(l) {
+ b = (b + a) % lprime
+ }
+ indices = append(indices, int(b))
+ }
+
+ sort.Ints(indices)
+ return indices
+}
+
+// RU10 intermediate encoding consists of the source symbols plus additional
+// intermediate symbols consisting of exactly the S and H blocks the R10 code
+// uses. The difference is that the code is unsystematic -- the source blocks
+// aren't necessarily going to be represented at the output -- so when we do
+// the decode we don't need to translate from the intermediate symbols back to
+// the source symbols: the source symbols are just the first K intermediate symbols.
+func (c *ru10Codec) GenerateIntermediateBlocks(message []byte, numBlocks int) []block {
+ sourceLong, sourceShort := partitionBytes(message, c.numSourceSymbols)
+ source := equalizeBlockLengths(sourceLong, sourceShort)
+
+ _, s, h := intermediateSymbols(c.numSourceSymbols)
+
+ k := c.numSourceSymbols
+ compositions := make([][]int, s)
+
+ for i := 0; i < k; i++ {
+ a := 1 + (int(math.Floor(float64(i)/float64(s))) % (s - 1))
+ b := i % s
+ compositions[b] = append(compositions[b], i)
+ b = (b + a) % s
+ compositions[b] = append(compositions[b], i)
+ b = (b + a) % s
+ compositions[b] = append(compositions[b], i)
+ }
+ for i := 0; i < s; i++ {
+ b := generateLubyTransformBlock(source, compositions[i])
+ source = append(source, b)
+ }
+
+ hprime := int(math.Ceil(float64(h) / 2))
+ m := buildGraySequence(k+s, hprime)
+ for i := 0; i < h; i++ {
+ hcomposition := make([]int, 0)
+ for j := 0; j < k+s; j++ {
+ if bitSet(uint(m[j]), uint(i)) {
+ hcomposition = append(hcomposition, j)
+ }
+ }
+ b := generateLubyTransformBlock(source, hcomposition)
+ source = append(source, b)
+ }
+ return source
+}
+
+// NewDecoder creates a new RU10 decoder
+func (c *ru10Codec) NewDecoder(messageLength int) Decoder {
+ return newRU10Decoder(c, messageLength)
+}
+
+// ru10Decoder is the corresponding decoder for fountain codes using the RU10 encoder.
+type ru10Decoder struct {
+ decoder *raptorDecoder
+}
+
+// newRU10Decoder creates a new raptor decoder for a given message. The
+// codec supplied must be the same one as the message was encoded with.
+func newRU10Decoder(c *ru10Codec, length int) *ru10Decoder {
+ return &ru10Decoder{
+ decoder: newRaptorDecoder(&raptorCodec{
+ SymbolAlignmentSize: c.symbolAlignmentSize,
+ NumSourceSymbols: c.numSourceSymbols},
+ length),
+ }
+}
+
+func (d *ru10Decoder) AddBlocks(blocks []LTBlock) bool {
+ c := ru10Codec{
+ symbolAlignmentSize: d.decoder.codec.SymbolAlignmentSize,
+ numSourceSymbols: d.decoder.codec.NumSourceSymbols}
+ for i := range blocks {
+ indices := c.PickIndices(blocks[i].BlockCode)
+ d.decoder.matrix.addEquation(indices, block{data: blocks[i].Data})
+ }
+ return d.decoder.matrix.determined()
+}
+
+func (d *ru10Decoder) Decode() []byte {
+ if !d.decoder.matrix.determined() {
+ return nil
+ }
+
+ d.decoder.matrix.reduce()
+
+ // Now the intermediate blocks are held in d.decoder.matrix.v. The source
+ // blocks are the first K intermediate blocks.
+ intermediate := d.decoder.matrix.v
+
+ lenLong, lenShort, numLong, numShort :=
+ partition(d.decoder.messageLength, d.decoder.codec.NumSourceSymbols)
+ out := make([]byte, d.decoder.messageLength)
+ out = out[0:0]
+ for i := 0; i < numLong; i++ {
+ out = append(out, intermediate[i].data[0:lenLong]...)
+ }
+ for i := numLong; i < numLong+numShort; i++ {
+ out = append(out, intermediate[i].data[0:lenShort]...)
+ }
+ return out
+}
diff --git a/util.go b/util.go
new file mode 100644
index 0000000..196173d
--- /dev/null
+++ b/util.go
@@ -0,0 +1,291 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/rand"
+ "sort"
+)
+
+// Note that these CDFs (cumulative distribution function) will be used for
+// selecting source blocks for code block generation. A typical algorithm
+// is to choose a number from a distribution, then pick uniformly that many
+// source blocks to XOR into a code block. To use CDF mapping values, pick a
+// random number r (0 <= r < 1) and then find the smallest i such that
+// CDF[i] >= r.
+
+// solitonDistribution returns a CDF mapping for the soliton distribution.
+// N (the number of elements in the CDF) cannot be less than 1
+// The CDF is one-based: the probability of picking 1 from the distribution
+// is CDF[1].
+func solitonDistribution(n int) []float64 {
+ cdf := make([]float64, n+1)
+ cdf[1] = 1 / float64(n)
+ for i := 2; i < len(cdf); i++ {
+ cdf[i] = cdf[i-1] + (1 / (float64(i) * float64(i-1)))
+ }
+ return cdf
+}
+
+// robustSolitonDistribution returns a CDF mapping for the robust solition
+// distribution.
+// This is an addition to the soliton distribution with three parameters,
+// N, M, and delta.
+// Before normalization, the correction pdf(i) = 1/i*M, for i=1..M-1,
+// pdf(M) = ln(N/(M*delta))/M
+// pdf(i) = 0 for i = M+1..N
+// These values are added to the ideal soliton distribution, and then the
+// result normalized.
+// The CDF is one-based: the probability of picking 1 from the distribution
+// is CDF[1].
+func robustSolitonDistribution(n int, m int, delta float64) []float64 {
+ pdf := make([]float64, n+1)
+
+ pdf[1] = 1/float64(n) + 1/float64(m)
+ total := pdf[1]
+ for i := 2; i < len(pdf); i++ {
+ pdf[i] = (1 / (float64(i) * float64(i-1)))
+ if i < m {
+ pdf[i] = 1 / (float64(i) * float64(m))
+ }
+ if i == m {
+ pdf[i] += math.Log(float64(n)/(float64(m)*delta)) / float64(m)
+ }
+ total += pdf[i]
+ }
+
+ cdf := make([]float64, n+1)
+ for i := 1; i < len(pdf); i++ {
+ pdf[i] /= total
+ cdf[i] = cdf[i-1] + pdf[i]
+ }
+ return cdf
+}
+
+// onlineSolitionDistribution returns a soliton-like distribution for
+// Online Codes
+// See http://pdos.csail.mit.edu/~petar/papers/maymounkov-bigdown-lncs.ps
+// 'Rateless Codes and Big Downloads' by Maymounkov and Mazieres
+// The distribution is described by a parameter epsilon, which is the
+// the "overage factor" required to reconstruct the source message in an
+// Online Code.
+// F = ciel(ln(eps^2/4 / ln(1 - eps/2))
+// and the pdf is pdf[1] = 1 - (1 + 1/F)/(1 + eps)
+// pdf[i] = ((1 - pdf[1])F) / ((F-1)i(i-1)) for 2 <= i <= F
+func onlineSolitonDistribution(eps float64) []float64 {
+ f := math.Ceil(math.Log(eps*eps/4) / math.Log(1-(eps/2)))
+
+ cdf := make([]float64, int(f+1))
+
+ // First coefficient is 1 - ( (1 + 1/f) / (1+e) )
+ rho := 1 - ((1 + (1 / f)) / (1 + eps))
+ cdf[1] = rho
+
+ // Subsequent i'th coefficient is (1-rho)*F / (F-1)i*(i-1)
+ for i := 2; i <= int(f); i++ {
+ rhoI := ((1 - rho) * f) / ((f - 1) * float64(i-1) * float64(i))
+ cdf[i] = cdf[i-1] + rhoI
+ }
+
+ return cdf
+}
+
+// pickDegree returns the smallest index i such that cdf[i] > r
+// (r a random number from the random generator)
+// cdf must be sorted in ascending order.
+func pickDegree(random *rand.Rand, cdf []float64) int {
+ r := random.Float64()
+ d := sort.SearchFloat64s(cdf, r)
+ if cdf[d] > r {
+ return d
+ }
+
+ if d < len(cdf)-1 {
+ return d + 1
+ } else {
+ return len(cdf) - 1
+ }
+}
+
+// sampleUniform picks num numbers from [0,max) uniformly.
+// There will be no duplicates.
+// If num >= max, simply returns a slice with all indices from 0 to max-1
+// without touching the random number generator.
+// The returned slice is sorted.
+func sampleUniform(random *rand.Rand, num, max int) []int {
+ if num >= max {
+ picks := make([]int, max)
+ for i := 0; i < max; i++ {
+ picks[i] = i
+ }
+ return picks
+ }
+
+ picks := make([]int, num)
+ seen := make(map[int]bool)
+ for i := 0; i < num; i++ {
+ p := random.Intn(max)
+ for seen[p] {
+ p = random.Intn(max)
+ }
+ picks[i] = p
+ seen[p] = true
+ }
+ sort.Ints(picks)
+ return picks
+}
+
+// partition is the block partitioning function from RFC 5053 S.5.3.1.2
+// See http://tools.ietf.org/html/rfc5053
+// Partitions a number i (a size) into j semi-equal pieces. The details are
+// in the return values: there are jl longer pieces of size il, and js shorter
+// pieces of size is.
+func partition(i, j int) (il int, is int, jl int, js int) {
+ il = int(math.Ceil(float64(i) / float64(j)))
+ is = int(math.Floor(float64(i) / float64(j)))
+ jl = i - (is * j)
+ js = j - jl
+
+ if jl == 0 {
+ il = 0
+ }
+ if js == 0 {
+ is = 0
+ }
+
+ return
+}
+
+// factorial calculates the factorial (x!) of the input argument.
+func factorial(x int) int {
+ f := 1
+ for i := 1; i <= x; i++ {
+ f *= i
+ }
+ return f
+}
+
+// centerBinomial calculates choose(x, ceil(x/2)) = x!/(x/2)!(x-(x/2)m!)
+func centerBinomial(x int) int {
+ return choose(x, x/2)
+}
+
+// choose calculates (n k) or n choose k. Tolerant of quite large n/k.
+func choose(n int, k int) int {
+ if k > n/2 {
+ k = n - k
+ }
+ numerator := make([]int, n-k)
+ denominator := make([]int, n-k)
+ for i, j := k+1, 1; i <= n; i, j = i+1, j+1 {
+ numerator[j-1] = int(i)
+ denominator[j-1] = j
+ }
+
+ if len(denominator) > 0 {
+ // find the first member of numerator not in denominator
+ z := sort.SearchInts(numerator, denominator[len(denominator)-1])
+ if z > 0 {
+ numerator = numerator[z+1:]
+ denominator = denominator[0 : len(denominator)-z-1]
+ }
+ }
+
+ for j := len(denominator) - 1; j > 0; j-- {
+ for i := len(numerator) - 1; i >= 0; i-- {
+ if numerator[i]%denominator[j] == 0 {
+ numerator[i] /= denominator[j]
+ denominator[j] = 1
+ break
+ }
+ }
+ }
+ f := 1
+ for _, i := range numerator {
+ f *= i
+ }
+ return f
+}
+
+// bitSet returns true if x has the b'th bit set
+func bitSet(x uint, b uint) bool {
+ return (x>>b)&1 == 1
+}
+
+// bitsSet returns how many bits in x are set.
+// This algorithm basically uses shifts and ANDs to sum up the bits in
+// a tree fashion.
+func bitsSet(x uint64) int {
+ x -= (x >> 1) & 0x5555555555555555
+ x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
+ x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f
+ return int((x * 0x0101010101010101) >> 56)
+}
+
+// grayCode calculates the gray code representation of the input argument
+// The Gray code is a binary representation in which successive values differ
+// by exactly one bit. See http://en.wikipedia.org/wiki/Gray_code
+func grayCode(x uint64) uint64 {
+ return (x >> 1) ^ x
+}
+
+// buildGraySequence returns a sequence (in ascending order) of "length" Gray numbers,
+// all of which have exactly "b" bits set.
+func buildGraySequence(length int, b int) []int {
+ s := make([]int, length)
+ i := 0
+ for x := uint64(0); ; x++ {
+ g := grayCode(x)
+ if bitsSet(g) == b {
+ s[i] = int(g)
+ i++
+ if i >= length {
+ break
+ }
+ }
+ }
+ return s
+}
+
+// isPrime tests x for primality. Works on numbers less than the square of
+// the largest smallPrimes array.
+func isPrime(x int) bool {
+ for _, p := range smallPrimes {
+ if p*p > x {
+ return true
+ }
+ if x%p == 0 {
+ return false
+ }
+ }
+ // Well, not really, but we don't know any higher primes for sure.
+ return true
+}
+
+// smallestPrimeGreaterOrEqual returns the smallest prime greater than or equal to x
+// TODO(gbillock): should handle up to 70000 or so
+func smallestPrimeGreaterOrEqual(x int) int {
+ if x <= smallPrimes[len(smallPrimes)-1] {
+ p := sort.Search(len(smallPrimes), func(i int) bool {
+ return smallPrimes[i] >= x
+ })
+ return smallPrimes[p]
+ }
+ for !isPrime(x) {
+ x++
+ }
+ return x
+}
diff --git a/util_test.go b/util_test.go
new file mode 100644
index 0000000..f979d9b
--- /dev/null
+++ b/util_test.go
@@ -0,0 +1,408 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package fountain
+
+import (
+ "math"
+ "math/rand"
+ "reflect"
+ "testing"
+)
+
+func almostEqual(a, b float64) bool {
+ return math.Abs(a-b) < 1e-5
+}
+
+func TestSolitonDistribution(t *testing.T) {
+ // Test several randomly chosen values of n.
+ tests := make([]int, 100)
+ for i := range tests {
+ tests[i] = rand.Intn(1e3) + 1 // Add 1 so that n will never be 0.
+ }
+ // Plus a couple of specific ones.
+ tests = append(tests, 1, 10)
+
+ for _, n := range tests {
+ cdf := solitonDistribution(n)
+ if len(cdf) != n+1 {
+ t.Errorf("n=%d: Wrong length CDF: %d", n, len(cdf))
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[n], 1) {
+ t.Errorf("n=%d: CDF[max] = %f, should be 1", n, cdf[n])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[0], 0) {
+ t.Errorf("n=%d: CDF[0] = %f, should be 0.0", n, cdf[0])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[1], 1/float64(n)) {
+ t.Errorf("n=%d: CDF[1] = %f, should be 1/n (=%f)", n, cdf[1], 1/float64(n))
+ t.Log("CDF=", cdf)
+ }
+ }
+}
+
+func TestRobustSolitonDistribution(t *testing.T) {
+ cdf := robustSolitonDistribution(10, 8, 0.1)
+ if len(cdf) != 11 {
+ t.Errorf("Wrong length CDF: %d, should be 11", len(cdf))
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[0], 0) {
+ t.Errorf("CDF[0] = %f, should be 0.0", cdf[0])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[1], .287474) {
+ t.Errorf("CDF[1] = %f, should be 0.287474", cdf[1])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[len(cdf)-1], 1) {
+ t.Errorf("CDF[max] = %f, should be very nearly 1", cdf[9])
+ t.Log("CDF=", cdf)
+ }
+
+ if (cdf[8]-cdf[7] < cdf[7]-cdf[6]) ||
+ (cdf[8]-cdf[7] < cdf[9]-cdf[8]) {
+ t.Errorf("CDF must have mode at M position(8): %v", cdf)
+ t.Logf("PDF[7, 8, 9] = %v, %v, %v", cdf[7]-cdf[6], cdf[8]-cdf[7], cdf[9]-cdf[8])
+ t.Log("CDF=", cdf)
+ }
+}
+
+func TestOnlineSolitonDistribution(t *testing.T) {
+ cdf := onlineSolitonDistribution(0.1)
+ if len(cdf) != 118 {
+ t.Errorf("Wrong length CDF: %d. Should be 118", len(cdf))
+ t.Log("CDF=", cdf)
+ }
+
+ if (cdf[2] - cdf[1]) < (cdf[1] - cdf[0]) {
+ t.Errorf("CDF should have mode in second position: %v", cdf[0:5])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[0], 0) {
+ t.Errorf("CDF[0] = %f, should be 0.0", cdf[0])
+ t.Log("CDF=", cdf)
+ }
+
+ if !almostEqual(cdf[len(cdf)-1], 1) {
+ t.Errorf("CDF[max] = %f, should be very nearly 1", cdf[len(cdf)-1])
+ t.Log("CDF=", cdf)
+ }
+
+ cdf = onlineSolitonDistribution(0.01)
+ if len(cdf) != 2116 {
+ t.Errorf("Wrong length CDF for 0.0: %d, should be 2116", len(cdf))
+ }
+}
+
+func TestPickDegree(t *testing.T) {
+ cdf := onlineSolitonDistribution(0.25)
+ random := rand.New(rand.NewSource(25))
+ var numLessThanFive int
+ for i := 0; i < 100; i++ {
+ d := pickDegree(random, cdf)
+ if d < 1 || d > len(cdf)-1 {
+ t.Errorf("Degree out of bounds: %d", d)
+ }
+ if d < 5 {
+ numLessThanFive++
+ }
+ }
+ if numLessThanFive < 80 {
+ t.Errorf("Too many large degrees: %d, should be < 80", numLessThanFive)
+ }
+}
+
+func TestSampleUniform(t *testing.T) {
+ random := rand.New(rand.NewSource(256))
+
+ var sampleTests = []struct {
+ num int
+ length int
+ expectedSamples []int
+ }{
+ {2, 5, []int{0, 4}},
+ {2, 2, []int{0, 1}},
+ {12, 2, []int{0, 1}},
+ }
+
+ for _, i := range sampleTests {
+ out := sampleUniform(random, i.num, i.length)
+ if !reflect.DeepEqual(out, i.expectedSamples) {
+ t.Errorf("Bad sample. Got %v, want %v", out, i.expectedSamples)
+ }
+ }
+}
+
+func TestPartition(t *testing.T) {
+ var partitionTests = []struct {
+ totalSize int
+ numPartitions int
+ numLong, numShort, lenLong, lenShort int
+ }{
+ {100, 10, 0, 10, 0, 10},
+ {100, 9, 1, 8, 12, 11},
+ {100, 11, 1, 10, 10, 9},
+ }
+
+ for _, i := range partitionTests {
+ il, is, jl, js := partition(i.totalSize, i.numPartitions)
+ if jl+js != i.numPartitions {
+ t.Errorf("Total blocks = %d, must be %d", il+jl, i.numPartitions)
+ }
+ if il*jl+is*js != i.totalSize {
+ t.Errorf("Total sampled size = %d, must be %d", il*jl+is*js, i.totalSize)
+ }
+ if jl != i.numLong {
+ t.Errorf("Bad number of long blocks. got %d, want %d", jl, i.numLong)
+ }
+ if js != i.numShort {
+ t.Errorf("Bad number of short blocks. got %d, want %d", js, i.numShort)
+ }
+ if il != i.lenLong {
+ t.Errorf("Bad long block length. got %d, want %d", il, i.lenLong)
+ }
+ if is != i.lenShort {
+ t.Errorf("Bad short block length. got %d, want %d", is, i.lenShort)
+ }
+ }
+}
+
+func TestFactorial(t *testing.T) {
+ var factorialTests = []struct {
+ x int
+ f int
+ }{
+ {0, 1},
+ {1, 1},
+ {5, 120},
+ {10, 3628800},
+ {14, 87178291200},
+ }
+
+ for _, test := range factorialTests {
+ if test.f != factorial(test.x) {
+ t.Errorf("%d! = %d, should be %d", test.x, factorial(test.x), test.f)
+ }
+ }
+}
+
+func TestBinomial(t *testing.T) {
+ var binomialTests = []struct {
+ x int
+ b int
+ }{
+ {2, 2},
+ {6, 20},
+ {7, 35},
+ {11, 462},
+ {12, 924},
+ }
+
+ for _, test := range binomialTests {
+ if test.b != centerBinomial(test.x) {
+ t.Errorf("(%d, %d/2) = %d, should be %d", test.x, test.x, centerBinomial(test.x), test.b)
+ }
+ }
+}
+
+func TestChoose(t *testing.T) {
+ var chooseTests = []struct {
+ n int
+ k int
+ comb int
+ }{
+ {0, 0, 1},
+ {1, 0, 1},
+ {2, 1, 2},
+ {7, 2, factorial(7) / (factorial(5) * factorial(2))},
+ {5, 3, factorial(5) / (factorial(2) * factorial(3))},
+ {12, 7, factorial(12) / (factorial(5) * factorial(7))},
+ {12, 1, 12},
+ {12, 2, 66},
+ {52, 5, 2598960},
+ {52, 1, 52},
+ {52, 52, 1},
+ {52, 0, 1},
+ }
+ for _, test := range chooseTests {
+ if choose(test.n, test.k) != test.comb {
+ t.Errorf("choose(%d, %d) = %d, should be %d",
+ test.n, test.k, choose(test.n, test.k), test.comb)
+ }
+ }
+}
+
+func TestBitSet(t *testing.T) {
+ var bitTests = []struct {
+ x uint
+ b uint
+ equal bool
+ }{
+ {0, 0, false},
+ {0, 1, false},
+ {1, 0, true},
+ {7, 1, true},
+ {16, 3, false},
+ {16, 4, true},
+ {16, 5, false},
+ {0x1000, 12, true},
+ {0x4000, 14, true},
+ {0x4000, 15, false},
+ }
+
+ for _, test := range bitTests {
+ if bitSet(test.x, test.b) != test.equal {
+ t.Errorf("%d bit set in %d = %t, should be %t", test.b, test.x, bitSet(test.x, test.b), test.equal)
+ }
+ }
+}
+
+// simpleBitsSet returns how many bits in x are set.
+func simpleBitsSet(x uint64) int {
+ var count int
+ var i uint64
+ s := uint64(1)
+ for i = 0; i < 64; i++ {
+ if (x & s) != 0 {
+ count++
+ }
+ s <<= 1
+ }
+ return count
+}
+
+func TestBitsSet(t *testing.T) {
+ var bitsSetTests = []struct {
+ x uint64
+ b int
+ }{
+ {0, 0},
+ {1, 1},
+ {2, 1},
+ {4, 1},
+ {6, 2},
+ {7, 3},
+ {11, 3},
+ {12, 2},
+ {1<<63 | 1<<62, 2},
+ {0x66666666, 16},
+ {0x46464646, 4*1 + 4*2},
+ {0xffdd4411, 2*4 + 2*3 + 2*1 + 2*1},
+ }
+
+ for _, test := range bitsSetTests {
+ if test.b != bitsSet(test.x) {
+ t.Errorf("bits set in %d = %d, should be %d", test.x, bitsSet(test.x), test.b)
+ }
+ }
+
+ for i := 0; i < 1000; i++ {
+ x := uint64(rand.Int63())
+ if simpleBitsSet(x) != bitsSet(x) {
+ t.Errorf("bits set in %d = %d, should be %d", x, bitsSet(x), simpleBitsSet(x))
+ }
+ }
+}
+
+func TestGrayCode(t *testing.T) {
+ var grayTests = []struct {
+ x uint64
+ g uint64
+ }{
+ {0, 0},
+ {1, 1},
+ {2, 3},
+ {5, 7},
+ {6, 5},
+ {9, 13},
+ }
+
+ for _, test := range grayTests {
+ if test.g != grayCode(test.x) {
+ t.Errorf("grayCode(%d) = %d, should be %d", test.x, grayCode(test.x), test.g)
+ }
+ }
+}
+
+func TestGraySequence(t *testing.T) {
+ var grayTests = []struct {
+ length int
+ b int
+ seq []int
+ }{
+ {4, 3, []int{7, 13, 14, 11}},
+ {6, 2, []int{3, 6, 5, 12, 10, 9}},
+ }
+
+ for _, test := range grayTests {
+ if !reflect.DeepEqual(buildGraySequence(test.length, test.b), test.seq) {
+ t.Errorf("gray sequence for %d = %v, should be %v", test.b, buildGraySequence(test.length, test.b), test.seq)
+ }
+ }
+}
+
+func TestSmallestPrime(t *testing.T) {
+ var primeTests = []struct {
+ x int
+ p int
+ }{
+ {0, 2},
+ {1, 2},
+ {2, 2},
+ {20, 23},
+ {1426, 1427},
+ {1427, 1427},
+ {1427, 1427},
+ {1997, 1997},
+ {1998, 1999},
+ {1999, 1999},
+ {3301, 3301},
+ {8522, 8527},
+ }
+
+ for _, test := range primeTests {
+ if test.p != smallestPrimeGreaterOrEqual(test.x) {
+ t.Errorf("smallestPrimeGreaterOrEqual(%d) = %d, should be %d", test.x, smallestPrimeGreaterOrEqual(test.x), test.p)
+ }
+ }
+}
+
+func TestIsPrime(t *testing.T) {
+ var primeTests = []struct {
+ x int
+ prime bool
+ }{
+ {2000, false},
+ {2099, true},
+ {2607, false},
+ {9007, true},
+ }
+
+ for _, test := range primeTests {
+ if test.prime != isPrime(test.x) {
+ t.Errorf("isPrime(%d) = %t, should be %t", test.x, isPrime(test.x), test.prime)
+ }
+ }
+}
