New life

README.md

@@ -1,21 +1,17 @@
 # modular
-A modular arithmetic library with an emphasis on speed
-[![Status](https://github.com/stewi1014/modular/actions/workflows/go.yml/badge.svg)](https://github.com/stewi1014/modular/actions/workflows/go.yml)
-
-Float64
-[![GoDoc](https://godoc.org/github.com/stewi1014/modular/modular64?status.svg)](https://godoc.org/github.com/stewi1014/modular/modular64)
+A math library with an emphasis on speed.
 
-Float32
-[![GoDoc](https://godoc.org/github.com/stewi1014/modular/modular32?status.svg)](https://godoc.org/github.com/stewi1014/modular/modular32)
+[![Status](https://github.com/stewi1014/modular/actions/workflows/go.yml/badge.svg)](https://github.com/stewi1014/modular/actions/workflows/go.yml)
+[![Go Reference](https://pkg.go.dev/badge/github.com/stewi1014/modular.svg)](https://pkg.go.dev/github.com/stewi1014/modular)
 
 
-Modular tries to leverage pre-computation as much as possible to allow direct computation in Congruent() and Index(), using [fastdiv] and pre-computed lookup tables. I can't test it on all hardware, but in principle should perform better than traditional modulo functions on all but the strangest of hardware.
+Modular tries to leverage pre-computation as much as possible to allow direct computation in Mod() and Index(), using [fastdiv] and pre-computed lookup tables. I can't test it on all hardware, but in principle should perform better than traditional modulo functions on all but the strangest of hardware.
 
 For example, on my machine with a modulo of 1e-25;
 
 **float64**
 
-| Number | Math.Mod | Modulus.Congruent | Indexer.Index |
+| Number | Math.Mod | Modulus.Mod | Indexer.Index |
 | ------ | ------ | ------ | ------ |
 | 0 | 14.1 ns/op | 6.07 ns/op | 12.4 ns/op |
 | 2.5e-25 | 29.2 ns/op | 14.4 ns/op | 15.7 ns/op |
@@ -25,7 +21,7 @@ For example, on my machine with a modulo of 1e-25;
 
 **float32**
 
-| Number | Math.Mod | Modulus.Congruent | Indexer.Index |
+| Number | Math.Mod | Modulus.Mod | Indexer.Index |
 | ------ | ------ | ------ | ------ |
 | 0 | 15.6 ns/op | 5.49 ns/op | 11.3 ns/op |
 | 2.5e-25 | 33.5 ns/op | 12.1 ns/op | 11.3 ns/op |
@@ -34,7 +30,7 @@ For example, on my machine with a modulo of 1e-25;
 
 ***
 
-I've use the name 'Congruent' as it's a more explicit definition of the function, and helps avoid confusion with other functions. It is the same as a euclidian 'modulo' function; that is, it finds the number satisfying '0 <= n < modulus' that is representative of the given number's congruency class, hence the name Congruent.
+Modular implements the euclidian modulo function; that is, it finds the number satisfying '0 <= n < modulus' that is representative of the given number's congruency class.
 
 
 [fastdiv]: <https://github.com/bmkessler/fastdiv>

go.mod

@@ -1,9 +1,7 @@
 module github.com/stewi1014/modular
 
-go 1.12
+go 1.21
 
 require (
 	github.com/bmkessler/fastdiv v0.0.0-20190227075523-41d5178f2044
-	github.com/chewxy/math32 v1.0.0
-	github.com/go-gl/mathgl v0.0.0-20190713194549-592312d8590a
 )

go.sum

@@ -1,9 +1,2 @@
 github.com/bmkessler/fastdiv v0.0.0-20190227075523-41d5178f2044 h1:8Rz0TcIbkvU+x53bDQgezQ3tbjrQSpZRr6h9JnR9lZU=
 github.com/bmkessler/fastdiv v0.0.0-20190227075523-41d5178f2044/go.mod h1:OI0uaNyGvxANSxteY6/mFRZs9EcQGqK30Bd1wqQj9zQ=
-github.com/chewxy/math32 v1.0.0 h1:RTt2SACA7BTzvbsAKVQJLZpV6zY2MZw4bW9L2HEKkHg=
-github.com/chewxy/math32 v1.0.0/go.mod h1:Miac6hA1ohdDUTagnvJy/q+aNnEk16qWUdb8ZVhvCN0=
-github.com/go-gl/mathgl v0.0.0-20190713194549-592312d8590a h1:yoAEv7yeWqfL/l9A/J5QOndXIJCldv+uuQB1DSNQbS0=
-github.com/go-gl/mathgl v0.0.0-20190713194549-592312d8590a/go.mod h1:yhpkQzEiH9yPyxDUGzkmgScbaBVlhC06qodikEM0ZwQ=
-golang.org/x/image v0.0.0-20190321063152-3fc05d484e9f h1:FO4MZ3N56GnxbqxGKqh+YTzUWQ2sDwtFQEZgLOxh9Jc=
-golang.org/x/image v0.0.0-20190321063152-3fc05d484e9f/go.mod h1:kZ7UVZpmo3dzQBMxlp+ypCbDeSB+sBbTgSJuh5dn5js=
-golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ=

integer/sqrt.go

@@ -0,0 +1,34 @@
+package integer
+
+// Sqrt returns the square root of n, rounding down and returning the remainder.
+func Sqrt(n int) (sqrt int, remainder int) {
+	// TODO; optimise
+
+	if n < 0 {
+		panic("square root of negative number")
+	}
+
+	if n < 2 {
+		return n, 0
+	}
+
+	shift := 2
+	for (n >> shift) != 0 {
+		shift += 2
+	}
+
+	result := 0
+	for shift >= 0 {
+		result = result << 1
+
+		large_cand := result + 1
+
+		if large_cand*large_cand <= n>>shift {
+			result = large_cand
+		}
+
+		shift -= 2
+	}
+
+	return result, n - (result * result)
+}

integer/sqrt_test.go

@@ -0,0 +1,43 @@
+package integer
+
+import "testing"
+
+// TODO; add more tests
+func TestSqrt(t *testing.T) {
+	tests := []struct {
+		name          string
+		arg           int
+		wantSqrt      int
+		wantRemainder int
+	}{
+		{
+			name:          "sqrt(4) = 2,0",
+			arg:           4,
+			wantSqrt:      2,
+			wantRemainder: 0,
+		},
+		{
+			name:          "sqrt(9) = 3,0",
+			arg:           9,
+			wantSqrt:      3,
+			wantRemainder: 0,
+		},
+		{
+			name:          "sqrt(10) = 3,1",
+			arg:           10,
+			wantSqrt:      3,
+			wantRemainder: 1,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			gotSqrt, gotRemainder := Sqrt(tt.arg)
+			if gotSqrt != tt.wantSqrt {
+				t.Errorf("Sqrt() gotSqrt = %v, want %v", gotSqrt, tt.wantSqrt)
+			}
+			if gotRemainder != tt.wantRemainder {
+				t.Errorf("Sqrt() gotRemainder = %v, want %v", gotRemainder, tt.wantRemainder)
+			}
+		})
+	}
+}

integeru64/sqrt.go

@@ -0,0 +1,30 @@
+package integeru64
+
+// Sqrt returns the square root of n, rounding down and returning the remainder.
+func Sqrt(n uint64) (sqrt uint64, remainder uint64) {
+	// TODO; optimise
+
+	if n < 2 {
+		return n, 0
+	}
+
+	shift := 2
+	for (n >> shift) != 0 {
+		shift += 2
+	}
+
+	result := uint64(0)
+	for shift >= 0 {
+		result = result << 1
+
+		large_cand := result + 1
+
+		if large_cand*large_cand <= n>>shift {
+			result = large_cand
+		}
+
+		shift -= 2
+	}
+
+	return result, n - (result * result)
+}

integeru64/sqrt_test.go

@@ -0,0 +1,43 @@
+package integeru64
+
+import "testing"
+
+// TODO; add more tests
+func TestSqrt(t *testing.T) {
+	tests := []struct {
+		name          string
+		arg           uint64
+		wantSqrt      uint64
+		wantRemainder uint64
+	}{
+		{
+			name:          "sqrt(4) = 2,0",
+			arg:           4,
+			wantSqrt:      2,
+			wantRemainder: 0,
+		},
+		{
+			name:          "sqrt(9) = 3,0",
+			arg:           9,
+			wantSqrt:      3,
+			wantRemainder: 0,
+		},
+		{
+			name:          "sqrt(10) = 3,1",
+			arg:           10,
+			wantSqrt:      3,
+			wantRemainder: 1,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			gotSqrt, gotRemainder := Sqrt(tt.arg)
+			if gotSqrt != tt.wantSqrt {
+				t.Errorf("Sqrt() gotSqrt = %v, want %v", gotSqrt, tt.wantSqrt)
+			}
+			if gotRemainder != tt.wantRemainder {
+				t.Errorf("Sqrt() gotRemainder = %v, want %v", gotRemainder, tt.wantRemainder)
+			}
+		})
+	}
+}

modular32/float32.go

@@ -2,11 +2,16 @@ package modular32
 
 import (
 	"math/bits"
-
-	math "github.com/chewxy/math32"
+	"unsafe"
 )
 
 const (
+	nan    = 0x7FE00000
+	posinf = 0x7F800000
+	neginf = 0xFF800000
+
+	MaxFloat = 0x1p127 * (1 + (1 - 0x1p-23))
+
 	fExponentBits = 8
 	fFractionBits = 23
 	fTotalBits    = 32
@@ -27,22 +32,22 @@ func shiftSub(up, down uint, n uint32) uint32 {
 	return n >> (down - up)
 }
 
-// frexp splits a float into it's exponent and fraction component. Sign bit is discarded.
+// Frexp splits a float into it's exponent and fraction component. Sign bit is discarded.
 // The 24th implied bit is placed in the fraction if appropriate
-func frexp(f float32) (uint32, uint) {
-	fbits := math.Float32bits(f)
+func Frexp(f float32) (uint32, uint) {
+	fbits := ToBits(f)
 	exp := uint((fbits & fExponentMask) >> fFractionBits)
 	if exp == 0 {
 		return fbits & fFractionMask, 0
 	}
 	return (fbits & fFractionMask) | (1 << fFractionBits), exp
 }
 
-// ldexp assembles a float from an exponent and fraction component. Sign is ignored.
+// Ldexp assembles a float from an exponent and fraction component. Sign is ignored.
 // Expects the 24th implied bit to be set if appropriate.
-func ldexp(fr uint32, exp uint) float32 {
+func Ldexp(fr uint32, exp uint) float32 {
 	if exp == 0 || fr == 0 {
-		return math.Float32frombits(fr & fFractionMask)
+		return FromBits(fr & fFractionMask)
 	}
 	shift := uint(bits.LeadingZeros32(fr) - fExponentBits)
 	if shift >= exp {
@@ -52,5 +57,35 @@ func ldexp(fr uint32, exp uint) float32 {
 		exp -= uint(shift)
 	}
 	fr = fr << shift
-	return math.Float32frombits((uint32(exp) << fFractionBits) | (fr & fFractionMask))
+	return FromBits((uint32(exp) << fFractionBits) | (fr & fFractionMask))
+}
+
+func ToBits(n float32) uint32 {
+	return *(*uint32)(unsafe.Pointer(&n))
+}
+
+func FromBits(bits uint32) float32 {
+	return *(*float32)(unsafe.Pointer(&bits))
+}
+
+func Abs(n float32) float32 {
+	if n <= 0 {
+		return -n
+	}
+	return n
+}
+
+func NaN() float32 {
+	return FromBits(nan)
+}
+
+func IsInf(n float32) bool {
+	return n > MaxFloat || n < -MaxFloat
+}
+
+func Inf(sign int) float32 {
+	if sign < 0 {
+		return FromBits(neginf)
+	}
+	return FromBits(posinf)
 }

modular32/indexer.go

@@ -4,7 +4,6 @@ import (
 	"errors"
 
 	"github.com/bmkessler/fastdiv"
-	math "github.com/chewxy/math32"
 )
 
 // Error types
@@ -18,27 +17,28 @@ var (
 // index must not be larger than 2**16, and modulus must be a normalised float
 //
 // Special cases:
-// 		NewIndexer(m, 0) = panic(integer divide by zero)
-// 		NewIndexer(m, i > 2**16) = ErrBadIndex
-// 		NewIndexer(0, i) = ErrBadModulo
-// 		NewIndexer(±Inf, i) = ErrBadModulo
-// 		NewIndexer(NaN, i) = ErrBadModulo
-// 		NewIndexer(m, i) = ErrBadModulo for |m| < 2**-126
+//
+//	NewIndexer(m, 0) = panic(integer divide by zero)
+//	NewIndexer(m, i > 2**16) = ErrBadIndex
+//	NewIndexer(0, i) = ErrBadModulo
+//	NewIndexer(±Inf, i) = ErrBadModulo
+//	NewIndexer(NaN, i) = ErrBadModulo
+//	NewIndexer(m, i) = ErrBadModulo for |m| < 2**-126
 func NewIndexer(modulus float32, index int) (Indexer, error) {
 	mod := NewModulus(modulus)
 	return mod.NewIndexer(index)
 }
 
 // NewIndexer creates a new indexer from the Modulus.
 func (m Modulus) NewIndexer(index int) (Indexer, error) {
-	if math.IsInf(m.mod, 0) || math.IsNaN(m.mod) || m.exp == 0 {
+	if IsInf(m.mod) || m.mod != m.mod || m.exp == 0 {
 		return Indexer{}, ErrBadModulo
 	}
 	if index > (1<<16) || index < 1 {
 		return Indexer{}, ErrBadIndex
 	}
 
-	modfr, _ := frexp(m.mod)
+	modfr, _ := Frexp(m.mod)
 	r := modfr << fExponentBits //r - range; is shifted fExponentBits to get a little more
 	rDivisor := r / uint32(index)
 	return Indexer{
@@ -63,14 +63,15 @@ type Indexer struct {
 // Otherwise, it always satisfies 0 <= num < index
 //
 // Special cases:
-// 		Index(NaN) = index
-// 		Index(±Inf) = index
+//
+//	Index(NaN) = index
+//	Index(±Inf) = index
 func (i Indexer) Index(n float32) int {
-	if math.IsNaN(n) || math.IsInf(n, 0) || i.i == 0 {
+	if n != n || IsInf(n) || i.i == 0 {
 		return i.i
 	}
 
-	nfr, nexp := frexp(n)
+	nfr, nexp := Frexp(n)
 	var nr uint32
 	switch {
 	case n > i.mod: