mersenne: add Sqr

goos: linux goarch: amd64 pkg: modernc.org/mathutil/mersenne cpu: AMD Ryzen 9 3900X 12-Core Processor BenchmarkSqr BenchmarkSqr/M1-mersenne-sqr BenchmarkSqr/M1-mersenne-sqr-24 13628048 79.38 ns/op 40 B/op 2 allocs/op BenchmarkSqr/M1-math-big-sqr BenchmarkSqr/M1-math-big-sqr-24 8802544 131.7 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M2-mersenne-sqr BenchmarkSqr/M2-mersenne-sqr-24 16625383 81.08 ns/op 40 B/op 2 allocs/op BenchmarkSqr/M2-math-big-sqr BenchmarkSqr/M2-math-big-sqr-24 9768525 131.8 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M4-mersenne-sqr BenchmarkSqr/M4-mersenne-sqr-24 14647800 85.37 ns/op 40 B/op 2 allocs/op BenchmarkSqr/M4-math-big-sqr BenchmarkSqr/M4-math-big-sqr-24 9056818 130.9 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M8-mersenne-sqr BenchmarkSqr/M8-mersenne-sqr-24 11580205 98.52 ns/op 40 B/op 2 allocs/op BenchmarkSqr/M8-math-big-sqr BenchmarkSqr/M8-math-big-sqr-24 9426870 128.2 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M16-mersenne-sqr BenchmarkSqr/M16-mersenne-sqr-24 10476908 116.7 ns/op 40 B/op 2 allocs/op BenchmarkSqr/M16-math-big-sqr BenchmarkSqr/M16-math-big-sqr-24 8949414 134.0 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M32-mersenne-sqr BenchmarkSqr/M32-mersenne-sqr-24 6495727 173.0 ns/op 48 B/op 2 allocs/op BenchmarkSqr/M32-math-big-sqr BenchmarkSqr/M32-math-big-sqr-24 9296271 132.0 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M64-mersenne-sqr BenchmarkSqr/M64-mersenne-sqr-24 4802661 241.9 ns/op 56 B/op 2 allocs/op BenchmarkSqr/M64-math-big-sqr BenchmarkSqr/M64-math-big-sqr-24 10145700 133.9 ns/op 88 B/op 3 allocs/op BenchmarkSqr/M128-mersenne-sqr BenchmarkSqr/M128-mersenne-sqr-24 3003795 413.0 ns/op 80 B/op 2 allocs/op BenchmarkSqr/M128-math-big-sqr BenchmarkSqr/M128-math-big-sqr-24 5549552 209.3 ns/op 144 B/op 3 allocs/op BenchmarkSqr/M256-mersenne-sqr BenchmarkSqr/M256-mersenne-sqr-24 1659831 722.0 ns/op 112 B/op 2 allocs/op BenchmarkSqr/M256-math-big-sqr BenchmarkSqr/M256-math-big-sqr-24 4880001 251.4 ns/op 192 B/op 3 allocs/op BenchmarkSqr/M512-mersenne-sqr BenchmarkSqr/M512-mersenne-sqr-24 858310 1349 ns/op 176 B/op 2 allocs/op BenchmarkSqr/M512-math-big-sqr BenchmarkSqr/M512-math-big-sqr-24 3462825 329.5 ns/op 288 B/op 3 allocs/op BenchmarkSqr/M1024-mersenne-sqr BenchmarkSqr/M1024-mersenne-sqr-24 420891 2565 ns/op 320 B/op 2 allocs/op BenchmarkSqr/M1024-math-big-sqr BenchmarkSqr/M1024-math-big-sqr-24 1964814 615.2 ns/op 480 B/op 3 allocs/op BenchmarkSqr/M2048-mersenne-sqr BenchmarkSqr/M2048-mersenne-sqr-24 206782 5177 ns/op 608 B/op 2 allocs/op BenchmarkSqr/M2048-math-big-sqr BenchmarkSqr/M2048-math-big-sqr-24 1000000 1255 ns/op 896 B/op 3 allocs/op BenchmarkSqr/M4096-mersenne-sqr BenchmarkSqr/M4096-mersenne-sqr-24 108934 10175 ns/op 1184 B/op 2 allocs/op BenchmarkSqr/M4096-math-big-sqr BenchmarkSqr/M4096-math-big-sqr-24 377995 3055 ns/op 1762 B/op 3 allocs/op BenchmarkSqr/M8192-mersenne-sqr BenchmarkSqr/M8192-mersenne-sqr-24 59828 20511 ns/op 2336 B/op 2 allocs/op BenchmarkSqr/M8192-math-big-sqr BenchmarkSqr/M8192-math-big-sqr-24 148688 8677 ns/op 3492 B/op 3 allocs/op BenchmarkSqr/M16384-mersenne-sqr BenchmarkSqr/M16384-mersenne-sqr-24 30609 38871 ns/op 4896 B/op 2 allocs/op BenchmarkSqr/M16384-math-big-sqr BenchmarkSqr/M16384-math-big-sqr-24 48075 25592 ns/op 7211 B/op 3 allocs/op BenchmarkSqr/M32768-mersenne-sqr BenchmarkSqr/M32768-mersenne-sqr-24 15621 82111 ns/op 9504 B/op 2 allocs/op BenchmarkSqr/M32768-math-big-sqr BenchmarkSqr/M32768-math-big-sqr-24 14744 87142 ns/op 32205 B/op 3 allocs/op BenchmarkSqr/M65536-mersenne-sqr BenchmarkSqr/M65536-mersenne-sqr-24 7050 154968 ns/op 18464 B/op 2 allocs/op BenchmarkSqr/M65536-math-big-sqr BenchmarkSqr/M65536-math-big-sqr-24 6040 237927 ns/op 66936 B/op 3 allocs/op BenchmarkSqr/M131072-mersenne-sqr BenchmarkSqr/M131072-mersenne-sqr-24 3741 324530 ns/op 40992 B/op 2 allocs/op BenchmarkSqr/M131072-math-big-sqr BenchmarkSqr/M131072-math-big-sqr-24 1988 615355 ns/op 125121 B/op 3 allocs/op BenchmarkSqr/M262144-mersenne-sqr BenchmarkSqr/M262144-mersenne-sqr-24 1950 620119 ns/op 73760 B/op 2 allocs/op BenchmarkSqr/M262144-math-big-sqr BenchmarkSqr/M262144-math-big-sqr-24 688 1755040 ns/op 246108 B/op 3 allocs/op BenchmarkSqr/M524288-mersenne-sqr BenchmarkSqr/M524288-mersenne-sqr-24 956 1290968 ns/op 139296 B/op 2 allocs/op BenchmarkSqr/M524288-math-big-sqr BenchmarkSqr/M524288-math-big-sqr-24 246 4932569 ns/op 476234 B/op 3 allocs/op BenchmarkSqr/M1048576-mersenne-sqr BenchmarkSqr/M1048576-mersenne-sqr-24 470 2560887 ns/op 270368 B/op 2 allocs/op BenchmarkSqr/M1048576-math-big-sqr BenchmarkSqr/M1048576-math-big-sqr-24 84 14315411 ns/op 935650 B/op 4 allocs/op BenchmarkSqr/M2097152-mersenne-sqr BenchmarkSqr/M2097152-mersenne-sqr-24 225 5156848 ns/op 532512 B/op 2 allocs/op BenchmarkSqr/M2097152-math-big-sqr BenchmarkSqr/M2097152-math-big-sqr-24 26 42251646 ns/op 1854422 B/op 4 allocs/op BenchmarkSqr/M4194304-mersenne-sqr BenchmarkSqr/M4194304-mersenne-sqr-24 117 10486764 ns/op 1056800 B/op 2 allocs/op BenchmarkSqr/M4194304-math-big-sqr BenchmarkSqr/M4194304-math-big-sqr-24 9 118619660 ns/op 3689229 B/op 4 allocs/op BenchmarkSqr/M8388608-mersenne-sqr BenchmarkSqr/M8388608-mersenne-sqr-24 61 19271563 ns/op 2105377 B/op 2 allocs/op BenchmarkSqr/M8388608-math-big-sqr BenchmarkSqr/M8388608-math-big-sqr-24 3 347243627 ns/op 7360130 B/op 5 allocs/op BenchmarkSqr/M16777216-mersenne-sqr BenchmarkSqr/M16777216-mersenne-sqr-24 27 37213931 ns/op 4202528 B/op 2 allocs/op BenchmarkSqr/M16777216-math-big-sqr BenchmarkSqr/M16777216-math-big-sqr-24 1 1044838615 ns/op 14704448 B/op 7 allocs/op BenchmarkSqr/M33554432-mersenne-sqr BenchmarkSqr/M33554432-mersenne-sqr-24 15 73673501 ns/op 8396838 B/op 2 allocs/op BenchmarkSqr/M33554432-math-big-sqr BenchmarkSqr/M33554432-math-big-sqr-24 1 3144805101 ns/op 29384512 B/op 7 allocs/op BenchmarkSqr/M67108864-mersenne-sqr BenchmarkSqr/M67108864-mersenne-sqr-24 8 145513519 ns/op 16785440 B/op 2 allocs/op BenchmarkSqr/M67108864-math-big-sqr BenchmarkSqr/M67108864-math-big-sqr-24 1 9167346999 ns/op 58744640 B/op 7 allocs/op BenchmarkSqr/M134217728-mersenne-sqr BenchmarkSqr/M134217728-mersenne-sqr-24 4 282906478 ns/op 33562656 B/op 2 allocs/op BenchmarkSqr/M134217728-math-big-sqr BenchmarkSqr/M134217728-math-big-sqr-24 1 27725921736 ns/op 117464896 B/op 7 allocs/op BenchmarkSqr/M268435456-mersenne-sqr BenchmarkSqr/M268435456-mersenne-sqr-24 2 584291840 ns/op 67117088 B/op 2 allocs/op BenchmarkSqr/M268435456-math-big-sqr BenchmarkSqr/M268435456-math-big-sqr-24 1 83496304443 ns/op 234905408 B/op 7 allocs/op PASS ok modernc.org/mathutil/mersenne 203.388s

mersenne: add Sqr
8ae7a21d · cznic · 19b10d78 · 8ae7a21d · 8ae7a21d
Commit 8ae7a21d authored 3 years ago by cznic
--- a/mersenne/all_test.go
+++ b/mersenne/all_test.go
@@ -5,6 +5,7 @@
 package mersenne // import "modernc.org/mathutil/mersenne"
 import (
+	"fmt"
 	"math"
 	"math/big"
 	"math/rand"
@@ -941,3 +942,39 @@ func TestModPow2(t *testing.T) {
 		f(uint32(rg.Next()), uint32(rg.Next()))
 	}
 }
+func TestSqr(t *testing.T) {
+	const N = 25000
+	for i := 0; i <= N; i++ {
+		g := Sqr(uint32(i))
+		e := New(uint32(i))
+		e.Mul(e, e)
+		if g.Cmp(e) != 0 {
+			t.Errorf("Sqr of M%d incorrect\ngot %s %#x\nexp %s %#x", i, g, g.Bits(), e, e.Bits())
+		}
+	}
+}
+var bigResult *big.Int
+func BenchmarkSqr(b *testing.B) {
+	for i := 0; i <= 28; i++ {
+		n := uint32(1) << i
+		b.Run(fmt.Sprintf("M%d-mersenne-sqr", n), func(b *testing.B) {
+			b.ReportAllocs()
+			for i := 0; i < b.N; i++ {
+				bigResult = Sqr(n)
+			}
+		})
+		b.Run(fmt.Sprintf("M%d-math-big-sqr", n), func(b *testing.B) {
+			a := New(n)
+			b.ReportAllocs()
+			b.ResetTimer()
+			for i := 0; i < b.N; i++ {
+				b := big.NewInt(0).Set(a)
+				b.Mul(b, b)
+				bigResult = b
+			}
+		})
+	}
+}
--- a/mersenne/mersenne.go
+++ b/mersenne/mersenne.go
@@ -24,8 +24,8 @@ import (
 	"math"
 	"math/big"
-	"modernc.org/mathutil"
 	"github.com/remyoudompheng/bigfft"
+	"modernc.org/mathutil"
 )
 var (
@@ -34,7 +34,7 @@ var (
 	_2 = big.NewInt(2)
 )
-// Knowns list the exponent of currently (March 2012) known Mersenne primes
+// Knowns list the exponent of currently (May 2021) known Mersenne primes
 // exponents in order.  See also: http://oeis.org/A000043 for a partial list.
 var Knowns = []uint32{
 	2,  // #1
@@ -91,10 +91,12 @@ var Knowns = []uint32{
 	57885161, // #48
 	74207281, // #49
 	77232917, // #50
+	82589933, // #51
 }
 // Known maps the exponent of known Mersenne primes its ordinal number/rank.
-// Ranks > 45 are currently provisional.
+// Ranks > 47 are currently provisional.
 var Known map[uint32]int
 func init() {
@@ -296,3 +298,63 @@ func ProbablyPrime(n, a uint32) bool {
 	x := ModPow(a, n-1, n)
 	return x.Cmp(_1) == 0 || x.Cmp(nMinus1) == 0
 }
+// Sqr returns the square of Mn.
+func Sqr(n uint32) *big.Int {
+	if n == 0 {
+		// (2^(0)-1)^2 = 0
+		return big.NewInt(0)
+	}
+	bitSize := 2 * int(n)
+	wordSize := bitSize/mathutil.IntBits + 1
+	bits := make([]big.Word, 0, wordSize)
+	var d, c, w big.Word
+	m := big.Word(1)
+	for i := uint32(0); i < n; i++ {
+		d++
+		c += d
+		if c&1 != 0 {
+			w |= m
+		}
+		m <<= 1
+		c >>= 1
+		if m == 0 {
+			bits = append(bits, w)
+			w = 0
+			m = 1
+		}
+	}
+	for i := uint32(1); i < n; i++ {
+		d--
+		c += d
+		if c&1 != 0 {
+			w |= m
+		}
+		m <<= 1
+		c >>= 1
+		if m == 0 {
+			bits = append(bits, w)
+			w = 0
+			m = 1
+		}
+	}
+	for c != 0 {
+		if c&1 != 0 {
+			w |= m
+		}
+		m <<= 1
+		c >>= 1
+		if m == 0 {
+			bits = append(bits, w)
+			w = 0
+			m = 1
+		}
+	}
+	if w != 0 {
+		bits = append(bits, w)
+	}
+	var r big.Int
+	r.SetBits(bits)
+	return &r
+}