diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..7cd1091
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,21 @@
+# If you prefer the allow list template instead of the deny list, see community template:
+# https://github.com/github/gitignore/blob/main/community/Golang/Go.AllowList.gitignore
+#
+# Binaries for programs and plugins
+*.exe
+*.exe~
+*.dll
+*.so
+*.dylib
+
+# Test binary, built with `go test -c`
+*.test
+
+# Output of the go coverage tool, specifically when used with LiteIDE
+*.out
+
+# Dependency directories (remove the comment below to include it)
+# vendor/
+
+# Go workspace file
+go.work
\ No newline at end of file
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000..9e74a8d
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2022 Matthew Willer
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
\ No newline at end of file
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..abce880
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,4 @@
+.PHONY: test
+
+test:
+	go test -bench . -benchmem ./...
\ No newline at end of file
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..411e3c4
--- /dev/null
+++ b/README.md
@@ -0,0 +1,70 @@
+Quadratic Residue Pseudo-Random Number Generator
+================================================
+
+A pseudo-random sequence generator for Go using [Preshing's method][preshing], based on computing quadratic residues
+modulo some prime number _p_:
+
+n ≡ x² (mod p)
+
+The PRNG is reasonably fast: only about 50% slower than the default source in `math/rand`.  However, it also has the
+added advantage that it produces a permuted sequence with no repeated values, which may be desirable in some contexts.
+
+[preshing]: https://preshing.com/20121224/how-to-generate-a-sequence-of-unique-random-integers/
+
+Installation
+------------
+
+```bash
+go get -u github.com/mattwiller/qrprng
+```
+
+Usage
+-----
+
+### As `rand.Source`
+
+For general-purpose use, the PRNG can be used as a source for all the functionality offered by [`math.Rand`][rand]:
+```go
+rng := rand.New(qrprng.Default())
+rng.Seed(0xb0a710ad)
+normDistDecimal := rng.NormFloat64()
+```
+
+Although the generator only produces `uint64` values directly, this allows it to be used in many different ways.
+
+[rand]: https://pkg.go.dev/math/rand#Rand
+
+### Random sequence access
+
+The generator can calculate any term of the sequence in constant time as a `uint64`:
+
+```go
+rng := qrprng.Default()
+n, _ := rng.Index(7_648_235_472)
+```
+
+### Custom generator
+
+The parameters of the PRNG are fully customizable, and can use any prime _p_ where p ≡ 3 (mod 4):
+
+```go
+p := unit64(11)
+rng, _ := qrprng.New(p, 0, 0)
+
+permutation := make([]uint64, p)
+for i := uint64(0); i < p; i++ {
+	permutation[i], _ = rng.Index(i)
+}
+
+fmt.Printf(permutation)
+// [8 6 4 7 9 3 2 0 5 1 10]
+```
+
+License
+-------
+
+Copyright 2022 Matthew Willer
+
+Licensed under the [MIT License](./LICENSE.txt).  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/benchmark_test.go b/benchmark_test.go
new file mode 100644
index 0000000..a40bdf3
--- /dev/null
+++ b/benchmark_test.go
@@ -0,0 +1,42 @@
+package qrprng_test
+
+import (
+	"math/rand"
+	"testing"
+
+	"github.com/mattwiller/qrprng"
+)
+
+var BenchmarkResult uint64
+
+func BenchmarkUInt64(b *testing.B) {
+	prime := uint64(9_021_057_379)
+	offset := uint64(1_000_014_012)
+	intermediateOffset := uint64(2_947_624_585)
+	rng, _ := qrprng.New(prime, intermediateOffset, offset)
+	for i := 0; i < b.N; i++ {
+		n := rng.Uint64()
+		BenchmarkResult = n
+	}
+}
+
+func BenchmarkDefault(b *testing.B) {
+	rng := qrprng.Default()
+	for i := 0; i < b.N; i++ {
+		n := rng.Uint64()
+		BenchmarkResult = n
+	}
+}
+
+func BenchmarkStdLib(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		BenchmarkResult = rand.Uint64()
+	}
+}
+
+func BenchmarkStdLibWithQRPRNGSource(b *testing.B) {
+	rng := rand.New(qrprng.Default())
+	for i := 0; i < b.N; i++ {
+		BenchmarkResult = rng.Uint64()
+	}
+}
diff --git a/go.mod b/go.mod
new file mode 100644
index 0000000..7ee7037
--- /dev/null
+++ b/go.mod
@@ -0,0 +1,11 @@
+module github.com/mattwiller/qrprng
+
+go 1.17
+
+require github.com/stretchr/testify v1.7.0
+
+require (
+	github.com/davecgh/go-spew v1.1.0 // indirect
+	github.com/pmezard/go-difflib v1.0.0 // indirect
+	gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c // indirect
+)
diff --git a/go.sum b/go.sum
new file mode 100644
index 0000000..acb88a4
--- /dev/null
+++ b/go.sum
@@ -0,0 +1,11 @@
+github.com/davecgh/go-spew v1.1.0 h1:ZDRjVQ15GmhC3fiQ8ni8+OwkZQO4DARzQgrnXU1Liz8=
+github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=
+github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM=
+github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4=
+github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME=
+github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY=
+github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg=
+gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM=
+gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0=
+gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo=
+gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM=
diff --git a/prng_test.go b/prng_test.go
new file mode 100644
index 0000000..1ce7dbf
--- /dev/null
+++ b/prng_test.go
@@ -0,0 +1,56 @@
+package qrprng_test
+
+import (
+	"testing"
+
+	"github.com/mattwiller/qrprng"
+	"github.com/stretchr/testify/assert"
+)
+
+func TestPermutationTrivial(t *testing.T) {
+	testPermutation(t, 11, 0, 0)
+}
+
+func TestPermutationTiny(t *testing.T) {
+	testPermutation(t, 83, 72, 0)
+}
+
+func TestPermutationSmall(t *testing.T) {
+	testPermutation(t, 727, 253, 0)
+}
+
+func TestFullExample(t *testing.T) {
+	prime := uint64(727)
+	rng, err := qrprng.New(prime, 253, 0)
+	// expected mask: 471, maxMask: 511
+	assert.NoError(t, err)
+
+	expected := []uint64{110, 242, 722, 61, 269, 172, 175, 206, 231, 145, 409, 476, 103, 410, 484, 81, 310, 671, 37, 317, 340, 365, 22, 59, 77, 620, 528, 376, 462, 414, 358, 123, 659, 542, 156, 96, 460, 165, 205, 264, 474, 117, 3, 302, 291, 84, 182, 270, 91, 196, 693, 719, 79, 385, 480, 326, 8, 486, 403, 378, 388, 703, 625, 29, 558, 373, 66, 161, 130, 191, 107, 281, 419, 456, 428, 62, 407, 56, 457, 507, 355, 128, 23, 337, 437, 104, 504, 636, 621, 212, 53, 250, 678, 224, 427, 605, 532, 716, 590, 618, 374, 221, 64, 332, 489, 72, 634, 593, 92, 362, 691, 106, 685, 266, 473, 301, 89, 644, 73, 118, 31, 46, 548, 261, 98, 357, 211, 346, 214, 335, 431, 298, 580, 341, 541, 436, 582, 351, 467, 610, 30, 662, 537, 252, 615, 568, 220, 78, 698, 381, 14, 247, 342, 651, 243, 15, 80, 533, 569, 82, 425, 320, 233, 614, 453, 633, 314, 553, 148, 676, 176, 240, 527, 155, 702, 93, 430, 690, 307, 97, 518, 255, 137, 479, 35, 186, 699, 669, 51, 321, 530, 299, 599, 609, 39, 169, 392, 115, 557, 287, 706, 271, 535, 675, 101, 640, 550, 306, 296, 190, 421, 109, 20, 674, 207, 408, 547, 125, 141, 708, 499, 126, 664, 681, 132, 583, 492, 503, 183, 594, 622, 313, 514, 352, 420, 244, 554, 32, 194, 154, 139, 363, 189, 95, 389, 393, 589, 394, 10, 57, 87, 217, 395, 416, 277, 160, 442, 134, 600, 441, 63, 360, 129, 423, 577, 500, 86, 322, 116, 601, 324, 418, 672, 377, 529, 561, 402, 143, 144, 556, 364, 613, 555, 397, 197, 552, 501, 0, 26, 203, 333, 369, 114, 700, 202, 171, 574, 228, 448, 71, 282, 286, 58, 2, 440, 433, 446, 663, 140, 717, 181, 177, 238, 564, 198, 626, 227, 348, 526, 163, 210, 283, 151, 396, 300, 576, 513, 665, 650, 273, 215, 630, 725, 76, 336, 44, 49, 113, 68, 639, 174, 619, 239, 602, 12, 323, 187, 327, 204, 546, 259, 406, 572, 318, 635, 482, 692, 100, 36, 256, 638, 468, 153, 549, 667, 517, 193, 27, 565, 166, 339, 405, 726, 366, 660, 641, 45, 38, 308, 519, 311, 573, 234, 603, 122, 50, 162, 643, 375, 368, 661, 458, 531, 632, 24, 168, 4, 262, 567, 575, 222, 697, 21, 229, 588, 449, 120, 131, 343, 292, 150, 711, 178, 54, 344, 380, 631, 483, 521, 704, 432, 679, 497, 515, 152, 40, 712, 538, 673, 562, 248, 138, 9, 658, 188, 560, 235, 319, 508, 545, 88, 379, 361, 648, 157, 670, 213, 657, 372, 713, 475, 536, 353, 496, 383, 185, 653, 316, 677, 399, 510, 170, 315, 142, 450, 99, 1, 293, 701, 275, 477, 288, 328, 689, 43, 723, 429, 612, 434, 680, 415, 656, 481, 345, 216, 463, 34, 83, 655, 444, 312, 596, 683, 505, 209, 354, 85, 135, 511, 90, 179, 359, 666, 512, 581, 721, 485, 401, 278, 471, 370, 274, 566, 295, 417, 598, 19, 498, 52, 616, 540, 438, 652, 470, 523, 435, 330, 551, 624, 367, 102, 263, 714, 637, 570, 25, 617, 297, 447, 506, 254, 258, 195, 251, 445, 382, 611, 164, 218, 180, 539, 629, 455, 146, 7, 121, 280, 226, 595, 208, 627, 94, 487, 424, 347, 272, 28, 232, 111, 623, 391, 260, 578, 384, 41, 387, 606, 411, 585, 715, 289, 159, 13, 592, 413, 331, 108, 586, 70, 509, 682, 465, 686, 225, 472, 559, 654, 276, 452, 133, 591, 724, 607, 60, 584, 443, 267, 112, 426, 524, 469, 124, 147, 390, 69, 439, 649, 303, 11, 237, 350, 459, 105, 647, 304, 684, 696, 587, 42, 522, 167, 334, 604, 488, 563, 5, 493, 534, 246, 571, 65, 520, 33, 16, 75, 74, 257, 284, 223, 502, 454, 404, 544, 219, 694, 645, 705, 249, 494, 646, 398, 597, 18, 356, 709, 309, 525, 668, 461, 230, 253, 294, 338, 495, 47, 173, 464, 6, 265, 688, 136, 491, 268, 642, 192, 386, 490, 608, 543, 371, 241, 687, 628, 285, 478, 245, 710, 516, 422, 17, 579, 290, 412, 127, 400, 67, 279, 201, 695, 48, 158, 718, 149, 236, 720, 199, 55, 325, 305, 184, 119, 200, 466, 451, 349, 329, 707}
+	for i := uint64(0); i < prime; i++ {
+		n, err := rng.Index(i)
+		assert.NoError(t, err)
+		assert.Equal(t, expected[i], n, "Sequence permutation diverged at index %d: expected %d / got %d", i, expected[i], n)
+	}
+}
+
+func TestPermutationMedium(t *testing.T) {
+	testPermutation(t, 11_587, 1234, 0)
+}
+
+func testPermutation(t *testing.T, prime, intermediateOffset, offset uint64) {
+	t.Helper()
+
+	rng, err := qrprng.New(prime, intermediateOffset, offset)
+	assert.NoError(t, err)
+
+	seen := make(map[uint64]struct{}, prime-1)
+	for i := uint64(0); i < prime; i++ {
+		n, err := rng.Index(i)
+		assert.NoError(t, err)
+		seen[n] = struct{}{}
+	}
+
+	for i := uint64(0); i < prime; i++ {
+		assert.Contains(t, seen, i, "Expected sequence to contain permuted element: %d", i)
+	}
+}
diff --git a/qrprng.go b/qrprng.go
new file mode 100644
index 0000000..b407cab
--- /dev/null
+++ b/qrprng.go
@@ -0,0 +1,144 @@
+package qrprng
+
+import (
+	"fmt"
+	"math"
+	"math/big"
+	"math/bits"
+)
+
+const (
+	INT63_MASK = (1 << 63) - 1
+	// Largest prime (3 mod 4) less than 2^64, permutes [0, 2^64-189)
+	DEFAULT_PRIME               = uint64(math.MaxUint64 - 188)
+	DEFAULT_INTERMEDIATE_OFFSET = 5_577_006_791_947_779_410
+)
+
+// QuadraticResiduePRNG is a thread-unsafe PRNG based on Preshing's method using quadratic residues.
+// The PRNG has the unique advantage of generating a permutation: when `offset` is 0, the output will cycle through
+// all numbers less than `prime` without repeats (until all have been output and the cycle restarts).
+// It implements both rand.Source and rand.Source64, and can be used via rand.New() to generate various random data.
+type QuadraticResiduePRNG struct {
+	prime              uint64
+	intermediateOffset uint64
+	offset             uint64
+	maxMask            uint64
+	mask               uint64
+	idx                uint64
+}
+
+// New creates a new PRNG instance with the given parameters, which are validated for correctness before creation.
+// The chosen prime must be 3 mod 4, and the intermediate offset (seed) can be any number less than the prime.  The
+// offset will be added to all output, effectively placing a floor on the output values.
+func New(prime, intermediateOffset, offset uint64) (*QuadraticResiduePRNG, error) {
+	if err := validate(prime, offset, intermediateOffset); err != nil {
+		return &QuadraticResiduePRNG{}, err
+	}
+
+	maxMask := calculateMaxMask(prime)
+	return &QuadraticResiduePRNG{
+		prime:              prime,
+		intermediateOffset: intermediateOffset,
+		offset:             offset,
+		maxMask:            maxMask,
+		mask:               calculateMask(prime, intermediateOffset, maxMask),
+	}, nil
+}
+
+// Default returns a new PRNG instance suitable for general-purpose use.
+// It uses the largest possible prime to permute 99.999999999999999% of possible uint64 values.
+func Default() *QuadraticResiduePRNG {
+	prng, err := New(DEFAULT_PRIME, DEFAULT_INTERMEDIATE_OFFSET, 0)
+	if err != nil {
+		panic(err)
+	}
+	return prng
+}
+
+// Index generates the ith element of the permutation described by the generator.
+// If i >= prime, then an error is returned.  However, it can be ignored if desired;
+// the sequence will simply cycle.
+func (prng *QuadraticResiduePRNG) Index(i uint64) (uint64, error) {
+	if i >= prng.prime {
+		return i, fmt.Errorf("invalid index %d: must be less than chosen prime", i)
+	}
+
+	intermediate := prng.permuteQPR(i) + prng.intermediateOffset
+	masked := prng.applyMask(intermediate % prng.prime)
+	return prng.offset + prng.permuteQPR(masked), nil
+}
+
+func (prng *QuadraticResiduePRNG) applyMask(i uint64) uint64 {
+	if i <= prng.maxMask {
+		return i ^ prng.mask
+	} else {
+		return i
+	}
+}
+
+func (prng *QuadraticResiduePRNG) permuteQPR(i uint64) uint64 {
+	residue := (i * i) % prng.prime
+	if i <= (prng.prime / 2) {
+		return residue
+	} else {
+		return prng.prime - residue
+	}
+}
+
+// QuadraticResiduePRNG implements math/rand.Source
+
+func (prng *QuadraticResiduePRNG) Int63() int64 {
+	return int64(prng.Uint64() & INT63_MASK)
+}
+
+// Seed changes the seed of the PRNG instance and resets the internal state of the generator.
+func (prng *QuadraticResiduePRNG) Seed(seed int64) {
+	if seed >= 0 {
+		prng.intermediateOffset = uint64(seed)
+	} else {
+		prng.intermediateOffset = math.MaxUint64 - uint64(-1*seed)
+	}
+	prng.idx = 0
+}
+
+// QuadraticResiduePRNG implements math/rand.Source64
+
+func (prng *QuadraticResiduePRNG) Uint64() uint64 {
+	n, _ := prng.Index(prng.idx)
+	prng.idx++
+	return n
+}
+
+// ===== Private functions =====
+
+func validate(prime, offset, intermediateOffset uint64) error {
+	if prime%4 != 3 {
+		return fmt.Errorf("invalid prime %d: must be 3 mod 4", prime)
+	} else if intermediateOffset >= prime {
+		return fmt.Errorf("invalid intermediate offset %d: must be less than chosen prime", intermediateOffset)
+	} else if p := bigIntFromUint64(prime); !p.ProbablyPrime(0) {
+		return fmt.Errorf("invalid prime %d: number is not prime", prime)
+	}
+	return nil
+}
+
+func calculateMaxMask(prime uint64) uint64 {
+	primeBits := bits.Len64(prime - 1)
+	return (1 << (primeBits - 1)) - 1
+}
+
+func calculateMask(prime, intermediateOffset, maxMask uint64) uint64 {
+	min := uint64(1 << (bits.Len64(maxMask) - 1))
+	return min + ((prime + intermediateOffset) % (maxMask - min))
+}
+
+func bigIntFromUint64(n uint64) *big.Int {
+	var result *big.Int
+	if n <= math.MaxInt64 {
+		result = big.NewInt(int64(n))
+	} else {
+		result = big.NewInt(math.MaxInt64)
+		result.Add(result, big.NewInt(int64(n-math.MaxInt64)))
+	}
+	return result
+}