From 2a834f5d2839e2213fce4e0a8dc20c2bcd26614f Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 20 Sep 2024 16:58:49 -0400 Subject: [PATCH 01/11] feat: implement Eisenstein integers arithmetic --- ecc/eisenstein.go | 103 ++++++++++++++++++ ecc/eisenstein_test.go | 232 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 335 insertions(+) create mode 100644 ecc/eisenstein.go create mode 100644 ecc/eisenstein_test.go diff --git a/ecc/eisenstein.go b/ecc/eisenstein.go new file mode 100644 index 000000000..7d6d223b0 --- /dev/null +++ b/ecc/eisenstein.go @@ -0,0 +1,103 @@ +// The Eisenstein integers form a commutative ring of algebraic integers in the +// algebraic number field Q(ω) – the third cyclotomic field. These are of the +// form z = a + bω, where a and b are integers and ω is a primitive third root +// of unity i.e. ω²+ω+1 = 0. +package ecc + +import ( + "math/big" +) + +// A ComplexNumber represents an arbitrary-precision Eisenstein integer. +type ComplexNumber struct { + A0, A1 big.Int +} + +// Equal returns true if z equals x, false otherwise +func (z *ComplexNumber) Equal(x *ComplexNumber) bool { + return z.A0.Cmp(&x.A0) == 0 && z.A1.Cmp(&x.A1) == 0 +} + +// Set sets z equal to x, and returns z. +func (z *ComplexNumber) Set(x *ComplexNumber) *ComplexNumber { + z.A0.Set(&x.A0) + z.A1.Set(&x.A1) + return z +} + +// Neg sets z to the negative of x, and returns z. +func (z *ComplexNumber) Neg(x *ComplexNumber) *ComplexNumber { + z.A0.Neg(&x.A0) + z.A1.Neg(&x.A1) + return z +} + +// Conjugate sets z to the conjugate of x, and returns z. +func (z *ComplexNumber) Conjugate(x *ComplexNumber) *ComplexNumber { + z.A0.Sub(&x.A0, &x.A1) + z.A1.Neg(&x.A1) + return z +} + +// Add sets z to the sum of x and y, and returns z. +func (z *ComplexNumber) Add(x, y *ComplexNumber) *ComplexNumber { + z.A0.Add(&x.A0, &y.A0) + z.A1.Add(&x.A1, &y.A1) + return z +} + +// Sub sets z to the difference of x and y, and returns z. +func (z *ComplexNumber) Sub(x, y *ComplexNumber) *ComplexNumber { + z.A0.Sub(&x.A0, &y.A0) + z.A1.Sub(&x.A1, &y.A1) + return z +} + +// Mul sets z to the product of x and y, and returns z. +// +// Given that ω²+ω+1=0, the explicit formula is: +// +// (x0+x1ω)(y0+y1ω) = (x0y0-x1y1) + (x0y1+x1y0-x1y1)ω +func (z *ComplexNumber) Mul(x, y *ComplexNumber) *ComplexNumber { + var t [3]big.Int + var z0, z1 big.Int + t[0].Mul(&x.A0, &y.A0) + t[1].Mul(&x.A1, &y.A1) + z0.Sub(&t[0], &t[1]) + t[0].Mul(&x.A0, &y.A1) + t[2].Mul(&x.A1, &y.A0) + t[0].Add(&t[0], &t[2]) + z1.Sub(&t[0], &t[1]) + z.A0.Set(&z0) + z.A1.Set(&z1) + return z +} + +// Norm returns the norm of z. +// +// The explicit formula is: +// +// N(x0+x1ω) = x0² + x1² - x0*x1 +func (z *ComplexNumber) Norm() *big.Int { + norm := new(big.Int) + temp := new(big.Int) + norm.Add( + norm.Mul(&z.A0, &z.A0), + temp.Mul(&z.A1, &z.A1), + ) + norm.Sub( + norm, + temp.Mul(&z.A0, &z.A1), + ) + return norm +} + +// Quo sets z to the quotient of x and y, and returns z. +func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { + norm := y.Norm() + z.Conjugate(y) + z.Mul(x, z) + z.A0.Quo(&z.A0, norm) + z.A1.Quo(&z.A1, norm) + return z +} diff --git a/ecc/eisenstein_test.go b/ecc/eisenstein_test.go new file mode 100644 index 000000000..bfe8381ed --- /dev/null +++ b/ecc/eisenstein_test.go @@ -0,0 +1,232 @@ +package ecc + +import ( + "math/big" + "math/rand" + "testing" + + "github.com/leanovate/gopter" + "github.com/leanovate/gopter/prop" +) + +const ( + nbFuzzShort = 10 + nbFuzz = 50 +) + +func TestEisensteinReceiverIsOperand(t *testing.T) { + + t.Parallel() + parameters := gopter.DefaultTestParameters() + if testing.Short() { + parameters.MinSuccessfulTests = nbFuzzShort + } else { + parameters.MinSuccessfulTests = nbFuzz + } + + properties := gopter.NewProperties(parameters) + + genE := GenComplexNumber() + + properties.Property("Having the receiver as operand (addition) should output the same result", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c, d ComplexNumber + d.Set(a) + c.Add(a, b) + a.Add(a, b) + b.Add(&d, b) + return a.Equal(b) && a.Equal(&c) && b.Equal(&c) + }, + genE, + genE, + )) + + properties.Property("Having the receiver as operand (sub) should output the same result", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c, d ComplexNumber + d.Set(a) + c.Sub(a, b) + a.Sub(a, b) + b.Sub(&d, b) + return a.Equal(b) && a.Equal(&c) && b.Equal(&c) + }, + genE, + genE, + )) + + properties.Property("Having the receiver as operand (mul) should output the same result", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c, d ComplexNumber + d.Set(a) + c.Mul(a, b) + a.Mul(a, b) + b.Mul(&d, b) + return a.Equal(b) && a.Equal(&c) && b.Equal(&c) + }, + genE, + genE, + )) + + properties.Property("Having the receiver as operand (neg) should output the same result", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + b.Neg(a) + a.Neg(a) + return a.Equal(&b) + }, + genE, + )) + + properties.Property("Having the receiver as operand (conjugate) should output the same result", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + b.Conjugate(a) + a.Conjugate(a) + return a.Equal(&b) + }, + genE, + )) + + properties.TestingRun(t, gopter.ConsoleReporter(false)) +} + +func TestEisensteinArithmetic(t *testing.T) { + + t.Parallel() + parameters := gopter.DefaultTestParameters() + if testing.Short() { + parameters.MinSuccessfulTests = nbFuzzShort + } else { + parameters.MinSuccessfulTests = nbFuzz + } + + properties := gopter.NewProperties(parameters) + + genE := GenComplexNumber() + + properties.Property("sub & add should leave an element invariant", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c ComplexNumber + c.Set(a) + c.Add(&c, b).Sub(&c, b) + return c.Equal(a) + }, + genE, + genE, + )) + + properties.Property("neg twice should leave an element invariant", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + b.Neg(a).Neg(&b) + return a.Equal(&b) + }, + genE, + )) + + properties.Property("conj twice should leave an element invariant", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + b.Conjugate(a).Conjugate(&b) + return a.Equal(&b) + }, + genE, + )) + + properties.Property("add zero should leave element invariant", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + zero := new(ComplexNumber) + b.Add(a, zero) + return a.Equal(&b) + }, + genE, + )) + + properties.Property("mul by one should leave element invariant", prop.ForAll( + func(a *ComplexNumber) bool { + var b ComplexNumber + one := &ComplexNumber{ + *big.NewInt(1), + *big.NewInt(0), + } + b.Mul(a, one) + return a.Equal(&b) + }, + genE, + )) + + properties.Property("add should be commutative", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c, d ComplexNumber + c.Add(a, b) + d.Add(b, a) + return c.Equal(&d) + }, + genE, + genE, + )) + + properties.Property("add should be assiocative", prop.ForAll( + func(a, b, c *ComplexNumber) bool { + var d, e ComplexNumber + d.Add(a, b).Add(&d, c) + e.Add(c, b).Add(&e, a) + return e.Equal(&d) + }, + genE, + genE, + genE, + )) + + properties.Property("mul should be commutative", prop.ForAll( + func(a, b *ComplexNumber) bool { + var c, d ComplexNumber + c.Mul(a, b) + d.Mul(b, a) + return c.Equal(&d) + }, + genE, + genE, + )) + + properties.Property("mul should be assiocative", prop.ForAll( + func(a, b, c *ComplexNumber) bool { + var d, e ComplexNumber + d.Mul(a, b).Mul(&d, c) + e.Mul(c, b).Mul(&e, a) + return e.Equal(&d) + }, + genE, + genE, + genE, + )) + + properties.Property("norm should always be positive", prop.ForAll( + func(a *ComplexNumber) bool { + return a.Norm().Sign() > 0 + }, + genE, + )) + + properties.TestingRun(t, gopter.ConsoleReporter(false)) +} + +// GenNumber generates a random integer +func GenNumber() gopter.Gen { + return func(genParams *gopter.GenParameters) *gopter.GenResult { + elmt := *big.NewInt(rand.Int63()) + genResult := gopter.NewGenResult(elmt, gopter.NoShrinker) + return genResult + } +} + +// GenComplexNumber generates a random integer +func GenComplexNumber() gopter.Gen { + return gopter.CombineGens( + GenNumber(), + GenNumber(), + ).Map(func(values []interface{}) *ComplexNumber { + return &ComplexNumber{A0: values[0].(big.Int), A1: values[1].(big.Int)} + }) +} From dda7a27023eb77d3e00e2a0ab54ea194617115f6 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Sat, 21 Sep 2024 00:44:57 -0400 Subject: [PATCH 02/11] feat: half-GCD for Eisenstein integers --- ecc/eisenstein.go | 58 +++++++++++++++++++++++++++++++++++++++--- ecc/eisenstein_test.go | 30 ++++++++++++++++++++++ 2 files changed, 85 insertions(+), 3 deletions(-) diff --git a/ecc/eisenstein.go b/ecc/eisenstein.go index 7d6d223b0..a950a51fc 100644 --- a/ecc/eisenstein.go +++ b/ecc/eisenstein.go @@ -13,18 +13,37 @@ type ComplexNumber struct { A0, A1 big.Int } +// String implements Stringer interface for fancy printing +func (z *ComplexNumber) String() string { + return z.A0.String() + "+(" + z.A1.String() + "*ω)" +} + // Equal returns true if z equals x, false otherwise func (z *ComplexNumber) Equal(x *ComplexNumber) bool { return z.A0.Cmp(&x.A0) == 0 && z.A1.Cmp(&x.A1) == 0 } -// Set sets z equal to x, and returns z. +// Set sets z to x, and returns z. func (z *ComplexNumber) Set(x *ComplexNumber) *ComplexNumber { z.A0.Set(&x.A0) z.A1.Set(&x.A1) return z } +// Set sets z to 0, and returns z. +func (z *ComplexNumber) SetZero() *ComplexNumber { + z.A0 = *big.NewInt(0) + z.A1 = *big.NewInt(0) + return z +} + +// Set sets z to 1, and returns z. +func (z *ComplexNumber) SetOne() *ComplexNumber { + z.A0 = *big.NewInt(1) + z.A1 = *big.NewInt(0) + return z +} + // Neg sets z to the negative of x, and returns z. func (z *ComplexNumber) Neg(x *ComplexNumber) *ComplexNumber { z.A0.Neg(&x.A0) @@ -97,7 +116,40 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { norm := y.Norm() z.Conjugate(y) z.Mul(x, z) - z.A0.Quo(&z.A0, norm) - z.A1.Quo(&z.A1, norm) + z.A0.Div(&z.A0, norm) + z.A1.Div(&z.A1, norm) return z } + +// HalfGCD returns the rational reconstruction of l mod r. +// This outputs u, v s.t. l = u/v mod r +func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { + + var aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber + + aRun.Set(a) + bRun.Set(b) + u.SetOne() + v.SetZero() + u_.SetZero() + v_.SetOne() + nbits := a.Norm().BitLen() / 2 + + for aRun.Norm().BitLen() > nbits { + quotient.Quo(&aRun, &bRun) + t.Mul(&bRun, "ient) + remainder.Sub(&aRun, &t) + t.Mul(&u_, "ient) + t1.Sub(&u, &t) + t.Mul(&v_, "ient) + t2.Sub(&v, &t) + aRun.Set(&bRun) + u.Set(&u_) + v.Set(&v_) + bRun.Set(&remainder) + u_.Set(&t1) + v_.Set(&t2) + } + + return [3]*ComplexNumber{&aRun, &v, &u} +} diff --git a/ecc/eisenstein_test.go b/ecc/eisenstein_test.go index bfe8381ed..794b2caeb 100644 --- a/ecc/eisenstein_test.go +++ b/ecc/eisenstein_test.go @@ -212,6 +212,36 @@ func TestEisensteinArithmetic(t *testing.T) { properties.TestingRun(t, gopter.ConsoleReporter(false)) } +func TestEisensteinHalfGCD(t *testing.T) { + + t.Parallel() + parameters := gopter.DefaultTestParameters() + if testing.Short() { + parameters.MinSuccessfulTests = nbFuzzShort + } else { + parameters.MinSuccessfulTests = nbFuzz + } + + properties := gopter.NewProperties(parameters) + + genE := GenComplexNumber() + + properties.Property("half-GCD", prop.ForAll( + func(a, b *ComplexNumber) bool { + res := HalfGCD(a, b) + var c, d ComplexNumber + c.Mul(b, res[1]) + d.Mul(a, res[2]) + d.Add(&c, &d) + return d.Equal(res[0]) + }, + genE, + genE, + )) + + properties.TestingRun(t, gopter.ConsoleReporter(false)) +} + // GenNumber generates a random integer func GenNumber() gopter.Gen { return func(genParams *gopter.GenParameters) *gopter.GenResult { From 6f69c7e12ed97fdeebd625f86410cee1ede88abf Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Sat, 21 Sep 2024 00:56:16 -0400 Subject: [PATCH 03/11] test: half-GCD test with bigger integers --- ecc/eisenstein_test.go | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/ecc/eisenstein_test.go b/ecc/eisenstein_test.go index 794b2caeb..40374cdbf 100644 --- a/ecc/eisenstein_test.go +++ b/ecc/eisenstein_test.go @@ -1,8 +1,8 @@ package ecc import ( + "crypto/rand" "math/big" - "math/rand" "testing" "github.com/leanovate/gopter" @@ -245,8 +245,9 @@ func TestEisensteinHalfGCD(t *testing.T) { // GenNumber generates a random integer func GenNumber() gopter.Gen { return func(genParams *gopter.GenParameters) *gopter.GenResult { - elmt := *big.NewInt(rand.Int63()) - genResult := gopter.NewGenResult(elmt, gopter.NoShrinker) + var prime, _ = new(big.Int).SetString("7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed", 16) // 2^255 - 19 + elmt, _ := rand.Int(rand.Reader, prime) + genResult := gopter.NewGenResult(*elmt, gopter.NoShrinker) return genResult } } From 9b3186918b44a860918374a2c91665df3cb924d7 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Sat, 21 Sep 2024 01:13:54 -0400 Subject: [PATCH 04/11] style: clean comments --- ecc/eisenstein.go | 8 ++++---- ecc/eisenstein_test.go | 15 +++++++++++++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/ecc/eisenstein.go b/ecc/eisenstein.go index a950a51fc..e24dbc165 100644 --- a/ecc/eisenstein.go +++ b/ecc/eisenstein.go @@ -30,14 +30,14 @@ func (z *ComplexNumber) Set(x *ComplexNumber) *ComplexNumber { return z } -// Set sets z to 0, and returns z. +// SetZero sets z to 0, and returns z. func (z *ComplexNumber) SetZero() *ComplexNumber { z.A0 = *big.NewInt(0) z.A1 = *big.NewInt(0) return z } -// Set sets z to 1, and returns z. +// SetOne sets z to 1, and returns z. func (z *ComplexNumber) SetOne() *ComplexNumber { z.A0 = *big.NewInt(1) z.A1 = *big.NewInt(0) @@ -121,8 +121,8 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { return z } -// HalfGCD returns the rational reconstruction of l mod r. -// This outputs u, v s.t. l = u/v mod r +// HalfGCD returns the rational reconstruction of a, b. +// This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { var aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber diff --git a/ecc/eisenstein_test.go b/ecc/eisenstein_test.go index 40374cdbf..fdd9f2543 100644 --- a/ecc/eisenstein_test.go +++ b/ecc/eisenstein_test.go @@ -261,3 +261,18 @@ func GenComplexNumber() gopter.Gen { return &ComplexNumber{A0: values[0].(big.Int), A1: values[1].(big.Int)} }) } + +// bench +func BenchmarkHalfGCD(b *testing.B) { + var prime, _ = new(big.Int).SetString("7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed", 16) // 2^255 - 19 + a0, _ := rand.Int(rand.Reader, prime) + a1, _ := rand.Int(rand.Reader, prime) + c0, _ := rand.Int(rand.Reader, prime) + c1, _ := rand.Int(rand.Reader, prime) + a := ComplexNumber{A0: *a0, A1: *a1} + c := ComplexNumber{A0: *c0, A1: *c1} + b.ResetTimer() + for i := 0; i < b.N; i++ { + HalfGCD(&a, &c) + } +} From 705a31ade028608759d7198becd54ed943c2f565 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Tue, 24 Sep 2024 15:11:22 -0400 Subject: [PATCH 05/11] refactor: move eisenstein under field/ --- {ecc => field/eisenstein}/eisenstein.go | 2 +- {ecc => field/eisenstein}/eisenstein_test.go | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) rename {ecc => field/eisenstein}/eisenstein.go (99%) rename {ecc => field/eisenstein}/eisenstein_test.go (99%) diff --git a/ecc/eisenstein.go b/field/eisenstein/eisenstein.go similarity index 99% rename from ecc/eisenstein.go rename to field/eisenstein/eisenstein.go index e24dbc165..8276a416b 100644 --- a/ecc/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -2,7 +2,7 @@ // algebraic number field Q(ω) – the third cyclotomic field. These are of the // form z = a + bω, where a and b are integers and ω is a primitive third root // of unity i.e. ω²+ω+1 = 0. -package ecc +package eisenstein import ( "math/big" diff --git a/ecc/eisenstein_test.go b/field/eisenstein/eisenstein_test.go similarity index 99% rename from ecc/eisenstein_test.go rename to field/eisenstein/eisenstein_test.go index fdd9f2543..197036f1e 100644 --- a/ecc/eisenstein_test.go +++ b/field/eisenstein/eisenstein_test.go @@ -1,4 +1,4 @@ -package ecc +package eisenstein import ( "crypto/rand" From 7febb567664cdbb5d4214d3cb3b7e47120be10e4 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Wed, 25 Sep 2024 14:24:14 -0400 Subject: [PATCH 06/11] fix: apply review suggestions --- field/eisenstein/doc.go | 7 +++++ field/eisenstein/eisenstein.go | 44 ++++++++++++++++++++++------- field/eisenstein/eisenstein_test.go | 39 +++++++++++++------------ 3 files changed, 60 insertions(+), 30 deletions(-) create mode 100644 field/eisenstein/doc.go diff --git a/field/eisenstein/doc.go b/field/eisenstein/doc.go new file mode 100644 index 000000000..3d353e2ab --- /dev/null +++ b/field/eisenstein/doc.go @@ -0,0 +1,7 @@ +// Package Eisenstein provides Eisenstein integer arithmetic. +// +// The Eisenstein integers form a commutative ring of algebraic integers in the +// algebraic number field Q(ω) – the third cyclotomic field. These are of the +// form z = a + bω, where a and b are integers and ω is a primitive third root +// of unity i.e. ω²+ω+1 = 0. +package eisenstein diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index 8276a416b..eddb54994 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -1,7 +1,3 @@ -// The Eisenstein integers form a commutative ring of algebraic integers in the -// algebraic number field Q(ω) – the third cyclotomic field. These are of the -// form z = a + bω, where a and b are integers and ω is a primitive third root -// of unity i.e. ω²+ω+1 = 0. package eisenstein import ( @@ -114,6 +110,9 @@ func (z *ComplexNumber) Norm() *big.Int { // Quo sets z to the quotient of x and y, and returns z. func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { norm := y.Norm() + if norm.Cmp(big.NewInt(0)) == 0 { + panic("division by zero") + } z.Conjugate(y) z.Mul(x, z) z.A0.Div(&z.A0, norm) @@ -121,11 +120,37 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { return z } +// QuoRem sets z to the quotient of x and y, r to the remaind, and returns z and r. +func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *ComplexNumber) { + norm := y.Norm() + if norm.Cmp(big.NewInt(0)) == 0 { + panic("division by zero") + } + z.Conjugate(y) + z.Mul(x, z) + z.A0.Div(&z.A0, norm) + z.A1.Div(&z.A1, norm) + var t ComplexNumber + r.Sub(x, t.Mul(y, z)) + return z, r +} + +// Min returns the minimum of z... with respect to the norm. +func Min(z ...*ComplexNumber) *ComplexNumber { + min := z[0] + for _, v := range z { + if v.Norm().Cmp(min.Norm()) == -1 { + min = v + } + } + return min +} + // HalfGCD returns the rational reconstruction of a, b. // This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { - var aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber + var check, aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber aRun.Set(a) bRun.Set(b) @@ -133,12 +158,11 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { v.SetZero() u_.SetZero() v_.SetOne() - nbits := a.Norm().BitLen() / 2 - for aRun.Norm().BitLen() > nbits { - quotient.Quo(&aRun, &bRun) - t.Mul(&bRun, "ient) - remainder.Sub(&aRun, &t) + // Eisenstein integers form an Euclidean domain for the norm + aNorm := a.Norm() + for check.Mul(&aRun, &aRun).Norm().Cmp(aNorm) == 1 { + quotient.QuoRem(&aRun, &bRun, &remainder) t.Mul(&u_, "ient) t1.Sub(&u, &t) t.Mul(&v_, "ient) diff --git a/field/eisenstein/eisenstein_test.go b/field/eisenstein/eisenstein_test.go index 197036f1e..224497201 100644 --- a/field/eisenstein/eisenstein_test.go +++ b/field/eisenstein/eisenstein_test.go @@ -12,6 +12,7 @@ import ( const ( nbFuzzShort = 10 nbFuzz = 50 + boundSize = 256 ) func TestEisensteinReceiverIsOperand(t *testing.T) { @@ -26,7 +27,7 @@ func TestEisensteinReceiverIsOperand(t *testing.T) { properties := gopter.NewProperties(parameters) - genE := GenComplexNumber() + genE := GenComplexNumber(boundSize) properties.Property("Having the receiver as operand (addition) should output the same result", prop.ForAll( func(a, b *ComplexNumber) bool { @@ -102,7 +103,7 @@ func TestEisensteinArithmetic(t *testing.T) { properties := gopter.NewProperties(parameters) - genE := GenComplexNumber() + genE := GenComplexNumber(boundSize) properties.Property("sub & add should leave an element invariant", prop.ForAll( func(a, b *ComplexNumber) bool { @@ -135,9 +136,9 @@ func TestEisensteinArithmetic(t *testing.T) { properties.Property("add zero should leave element invariant", prop.ForAll( func(a *ComplexNumber) bool { - var b ComplexNumber - zero := new(ComplexNumber) - b.Add(a, zero) + var b, zero ComplexNumber + zero.SetZero() + b.Add(a, &zero) return a.Equal(&b) }, genE, @@ -145,12 +146,9 @@ func TestEisensteinArithmetic(t *testing.T) { properties.Property("mul by one should leave element invariant", prop.ForAll( func(a *ComplexNumber) bool { - var b ComplexNumber - one := &ComplexNumber{ - *big.NewInt(1), - *big.NewInt(0), - } - b.Mul(a, one) + var b, one ComplexNumber + one.SetOne() + b.Mul(a, &one) return a.Equal(&b) }, genE, @@ -204,7 +202,7 @@ func TestEisensteinArithmetic(t *testing.T) { properties.Property("norm should always be positive", prop.ForAll( func(a *ComplexNumber) bool { - return a.Norm().Sign() > 0 + return a.Norm().Sign() >= 0 }, genE, )) @@ -224,7 +222,7 @@ func TestEisensteinHalfGCD(t *testing.T) { properties := gopter.NewProperties(parameters) - genE := GenComplexNumber() + genE := GenComplexNumber(boundSize) properties.Property("half-GCD", prop.ForAll( func(a, b *ComplexNumber) bool { @@ -243,20 +241,21 @@ func TestEisensteinHalfGCD(t *testing.T) { } // GenNumber generates a random integer -func GenNumber() gopter.Gen { +func GenNumber(boundSize int64) gopter.Gen { return func(genParams *gopter.GenParameters) *gopter.GenResult { - var prime, _ = new(big.Int).SetString("7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed", 16) // 2^255 - 19 - elmt, _ := rand.Int(rand.Reader, prime) + var bound big.Int + bound.Exp(big.NewInt(2), big.NewInt(boundSize), nil) + elmt, _ := rand.Int(genParams.Rng, &bound) genResult := gopter.NewGenResult(*elmt, gopter.NoShrinker) return genResult } } // GenComplexNumber generates a random integer -func GenComplexNumber() gopter.Gen { +func GenComplexNumber(boundSize int64) gopter.Gen { return gopter.CombineGens( - GenNumber(), - GenNumber(), + GenNumber(boundSize), + GenNumber(boundSize), ).Map(func(values []interface{}) *ComplexNumber { return &ComplexNumber{A0: values[0].(big.Int), A1: values[1].(big.Int)} }) @@ -273,6 +272,6 @@ func BenchmarkHalfGCD(b *testing.B) { c := ComplexNumber{A0: *c0, A1: *c1} b.ResetTimer() for i := 0; i < b.N; i++ { - HalfGCD(&a, &c) + _ = HalfGCD(&a, &c) } } From 12d53bfaad6f1706d8f940137f2512fb27a29d25 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Wed, 25 Sep 2024 14:30:30 -0400 Subject: [PATCH 07/11] fix: makes linter happy --- field/eisenstein/eisenstein.go | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index eddb54994..54bfde8af 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -120,7 +120,7 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { return z } -// QuoRem sets z to the quotient of x and y, r to the remaind, and returns z and r. +// QuoRem sets z to the quotient of x and y, r to the remainder, and returns z and r. func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *ComplexNumber) { norm := y.Norm() if norm.Cmp(big.NewInt(0)) == 0 { From fa1b905322a14c23c71718bf4a361ea53490ab75 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 26 Sep 2024 16:01:43 -0400 Subject: [PATCH 08/11] fix: consider all possible remainders --- field/eisenstein/eisenstein.go | 35 ++++++++++++++++++++++------- field/eisenstein/eisenstein_test.go | 2 +- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index 54bfde8af..b77aa9a39 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -120,6 +120,23 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { return z } +func (z *ComplexNumber) Rem(x, y, q *ComplexNumber) *ComplexNumber { + var t ComplexNumber + z.Sub(x, t.Mul(y, q)) + var r1, r2, r3, r4, r5, r6, w, wy, w2y ComplexNumber + w.A1.SetUint64(1) + wy.Mul(&w, y) + w2y.Mul(&w, &wy) + r1.Sub(z, y) + r2.Add(z, y) + r3.Sub(z, &wy) + r4.Add(z, &wy) + r5.Sub(z, &w2y) + r6.Add(z, &w2y) + z = Min(z, &r1, &r2, &r3, &r4, &r5, &r6) + return z +} + // QuoRem sets z to the quotient of x and y, r to the remainder, and returns z and r. func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *ComplexNumber) { norm := y.Norm() @@ -130,8 +147,8 @@ func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *Complex z.Mul(x, z) z.A0.Div(&z.A0, norm) z.A1.Div(&z.A1, norm) - var t ComplexNumber - r.Sub(x, t.Mul(y, z)) + r.Rem(x, y, z) + return z, r } @@ -150,7 +167,7 @@ func Min(z ...*ComplexNumber) *ComplexNumber { // This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { - var check, aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber + var check, aRun, bRun, u, v, u_, v_, quotient, remainder, t1, t2 ComplexNumber aRun.Set(a) bRun.Set(b) @@ -161,12 +178,10 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { // Eisenstein integers form an Euclidean domain for the norm aNorm := a.Norm() - for check.Mul(&aRun, &aRun).Norm().Cmp(aNorm) == 1 { + for check.Mul(&aRun, &aRun).Norm().Cmp(aNorm) >= 0 { quotient.QuoRem(&aRun, &bRun, &remainder) - t.Mul(&u_, "ient) - t1.Sub(&u, &t) - t.Mul(&v_, "ient) - t2.Sub(&v, &t) + t1.Rem(&u, &u_, "ient) + t2.Rem(&v, &v_, "ient) aRun.Set(&bRun) u.Set(&u_) v.Set(&v_) @@ -175,5 +190,9 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { v_.Set(&t2) } + if v_.Norm().Cmp(v.Norm()) == -1 { + return [3]*ComplexNumber{&bRun, &v_, &u_} + } + return [3]*ComplexNumber{&aRun, &v, &u} } diff --git a/field/eisenstein/eisenstein_test.go b/field/eisenstein/eisenstein_test.go index 224497201..89c7dc5d4 100644 --- a/field/eisenstein/eisenstein_test.go +++ b/field/eisenstein/eisenstein_test.go @@ -12,7 +12,7 @@ import ( const ( nbFuzzShort = 10 nbFuzz = 50 - boundSize = 256 + boundSize = 128 ) func TestEisensteinReceiverIsOperand(t *testing.T) { From 13fa612255da2bcf57b77a24510c69a17b5476e6 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 27 Sep 2024 17:21:45 -0400 Subject: [PATCH 09/11] fix: use sqrt in eisenstein halfgcd condition --- field/eisenstein/eisenstein.go | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index b77aa9a39..9fe12ad83 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -167,7 +167,8 @@ func Min(z ...*ComplexNumber) *ComplexNumber { // This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { - var check, aRun, bRun, u, v, u_, v_, quotient, remainder, t1, t2 ComplexNumber + var aRun, bRun, u, v, u_, v_, quotient, remainder, t1, t2 ComplexNumber + var sqrt big.Int aRun.Set(a) bRun.Set(b) @@ -177,8 +178,8 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { v_.SetOne() // Eisenstein integers form an Euclidean domain for the norm - aNorm := a.Norm() - for check.Mul(&aRun, &aRun).Norm().Cmp(aNorm) >= 0 { + sqrt.Sqrt(a.Norm()) + for bRun.Norm().Cmp(&sqrt) == 1 { quotient.QuoRem(&aRun, &bRun, &remainder) t1.Rem(&u, &u_, "ient) t2.Rem(&v, &v_, "ient) @@ -190,9 +191,5 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { v_.Set(&t2) } - if v_.Norm().Cmp(v.Norm()) == -1 { - return [3]*ComplexNumber{&bRun, &v_, &u_} - } - - return [3]*ComplexNumber{&aRun, &v, &u} + return [3]*ComplexNumber{&bRun, &v_, &u_} } From 507de623a0cd11ebe1fd3e858e12612090e38670 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 3 Oct 2024 15:30:25 -0400 Subject: [PATCH 10/11] perf(eisentein/half-GCD): only 1 remainder option --- field/eisenstein/eisenstein.go | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index 9fe12ad83..1b1f2a066 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -120,6 +120,7 @@ func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { return z } +// Rem sets z to the minimum of all possible remainders of x by y, and returns z. func (z *ComplexNumber) Rem(x, y, q *ComplexNumber) *ComplexNumber { var t ComplexNumber z.Sub(x, t.Mul(y, q)) @@ -147,7 +148,8 @@ func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *Complex z.Mul(x, z) z.A0.Div(&z.A0, norm) z.A1.Div(&z.A1, norm) - r.Rem(x, y, z) + r.Mul(y, z) + r.Sub(x, r) return z, r } @@ -167,7 +169,7 @@ func Min(z ...*ComplexNumber) *ComplexNumber { // This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { - var aRun, bRun, u, v, u_, v_, quotient, remainder, t1, t2 ComplexNumber + var aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber var sqrt big.Int aRun.Set(a) @@ -179,10 +181,12 @@ func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { // Eisenstein integers form an Euclidean domain for the norm sqrt.Sqrt(a.Norm()) - for bRun.Norm().Cmp(&sqrt) == 1 { + for bRun.Norm().Cmp(&sqrt) >= 0 { quotient.QuoRem(&aRun, &bRun, &remainder) - t1.Rem(&u, &u_, "ient) - t2.Rem(&v, &v_, "ient) + t.Mul(&u_, "ient) + t1.Sub(&u, &t) + t.Mul(&v_, "ient) + t2.Sub(&v, &t) aRun.Set(&bRun) u.Set(&u_) v.Set(&v_) From 2e01a9f52a522c3630fec533ec5db066751ccd08 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 3 Oct 2024 15:53:16 -0400 Subject: [PATCH 11/11] refactor: apply review suggestions --- field/eisenstein/eisenstein.go | 109 +++++++++++----------------- field/eisenstein/eisenstein_test.go | 22 +++--- 2 files changed, 54 insertions(+), 77 deletions(-) diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go index 1b1f2a066..e80915413 100644 --- a/field/eisenstein/eisenstein.go +++ b/field/eisenstein/eisenstein.go @@ -6,7 +6,18 @@ import ( // A ComplexNumber represents an arbitrary-precision Eisenstein integer. type ComplexNumber struct { - A0, A1 big.Int + A0, A1 *big.Int +} + +func (z *ComplexNumber) init() { + if z.A0 == nil { + z.A0 = new(big.Int) + + } + if z.A1 == nil { + z.A1 = new(big.Int) + + } } // String implements Stringer interface for fancy printing @@ -16,55 +27,60 @@ func (z *ComplexNumber) String() string { // Equal returns true if z equals x, false otherwise func (z *ComplexNumber) Equal(x *ComplexNumber) bool { - return z.A0.Cmp(&x.A0) == 0 && z.A1.Cmp(&x.A1) == 0 + return z.A0.Cmp(x.A0) == 0 && z.A1.Cmp(x.A1) == 0 } // Set sets z to x, and returns z. func (z *ComplexNumber) Set(x *ComplexNumber) *ComplexNumber { - z.A0.Set(&x.A0) - z.A1.Set(&x.A1) + z.init() + z.A0.Set(x.A0) + z.A1.Set(x.A1) return z } // SetZero sets z to 0, and returns z. func (z *ComplexNumber) SetZero() *ComplexNumber { - z.A0 = *big.NewInt(0) - z.A1 = *big.NewInt(0) + z.A0 = big.NewInt(0) + z.A1 = big.NewInt(0) return z } // SetOne sets z to 1, and returns z. func (z *ComplexNumber) SetOne() *ComplexNumber { - z.A0 = *big.NewInt(1) - z.A1 = *big.NewInt(0) + z.A0 = big.NewInt(1) + z.A1 = big.NewInt(0) return z } // Neg sets z to the negative of x, and returns z. func (z *ComplexNumber) Neg(x *ComplexNumber) *ComplexNumber { - z.A0.Neg(&x.A0) - z.A1.Neg(&x.A1) + z.init() + z.A0.Neg(x.A0) + z.A1.Neg(x.A1) return z } // Conjugate sets z to the conjugate of x, and returns z. func (z *ComplexNumber) Conjugate(x *ComplexNumber) *ComplexNumber { - z.A0.Sub(&x.A0, &x.A1) - z.A1.Neg(&x.A1) + z.init() + z.A0.Sub(x.A0, x.A1) + z.A1.Neg(x.A1) return z } // Add sets z to the sum of x and y, and returns z. func (z *ComplexNumber) Add(x, y *ComplexNumber) *ComplexNumber { - z.A0.Add(&x.A0, &y.A0) - z.A1.Add(&x.A1, &y.A1) + z.init() + z.A0.Add(x.A0, y.A0) + z.A1.Add(x.A1, y.A1) return z } // Sub sets z to the difference of x and y, and returns z. func (z *ComplexNumber) Sub(x, y *ComplexNumber) *ComplexNumber { - z.A0.Sub(&x.A0, &y.A0) - z.A1.Sub(&x.A1, &y.A1) + z.init() + z.A0.Sub(x.A0, y.A0) + z.A1.Sub(x.A1, y.A1) return z } @@ -74,13 +90,14 @@ func (z *ComplexNumber) Sub(x, y *ComplexNumber) *ComplexNumber { // // (x0+x1ω)(y0+y1ω) = (x0y0-x1y1) + (x0y1+x1y0-x1y1)ω func (z *ComplexNumber) Mul(x, y *ComplexNumber) *ComplexNumber { + z.init() var t [3]big.Int var z0, z1 big.Int - t[0].Mul(&x.A0, &y.A0) - t[1].Mul(&x.A1, &y.A1) + t[0].Mul(x.A0, y.A0) + t[1].Mul(x.A1, y.A1) z0.Sub(&t[0], &t[1]) - t[0].Mul(&x.A0, &y.A1) - t[2].Mul(&x.A1, &y.A0) + t[0].Mul(x.A0, y.A1) + t[2].Mul(x.A1, y.A0) t[0].Add(&t[0], &t[2]) z1.Sub(&t[0], &t[1]) z.A0.Set(&z0) @@ -97,47 +114,16 @@ func (z *ComplexNumber) Norm() *big.Int { norm := new(big.Int) temp := new(big.Int) norm.Add( - norm.Mul(&z.A0, &z.A0), - temp.Mul(&z.A1, &z.A1), + norm.Mul(z.A0, z.A0), + temp.Mul(z.A1, z.A1), ) norm.Sub( norm, - temp.Mul(&z.A0, &z.A1), + temp.Mul(z.A0, z.A1), ) return norm } -// Quo sets z to the quotient of x and y, and returns z. -func (z *ComplexNumber) Quo(x, y *ComplexNumber) *ComplexNumber { - norm := y.Norm() - if norm.Cmp(big.NewInt(0)) == 0 { - panic("division by zero") - } - z.Conjugate(y) - z.Mul(x, z) - z.A0.Div(&z.A0, norm) - z.A1.Div(&z.A1, norm) - return z -} - -// Rem sets z to the minimum of all possible remainders of x by y, and returns z. -func (z *ComplexNumber) Rem(x, y, q *ComplexNumber) *ComplexNumber { - var t ComplexNumber - z.Sub(x, t.Mul(y, q)) - var r1, r2, r3, r4, r5, r6, w, wy, w2y ComplexNumber - w.A1.SetUint64(1) - wy.Mul(&w, y) - w2y.Mul(&w, &wy) - r1.Sub(z, y) - r2.Add(z, y) - r3.Sub(z, &wy) - r4.Add(z, &wy) - r5.Sub(z, &w2y) - r6.Add(z, &w2y) - z = Min(z, &r1, &r2, &r3, &r4, &r5, &r6) - return z -} - // QuoRem sets z to the quotient of x and y, r to the remainder, and returns z and r. func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *ComplexNumber) { norm := y.Norm() @@ -146,25 +132,14 @@ func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *Complex } z.Conjugate(y) z.Mul(x, z) - z.A0.Div(&z.A0, norm) - z.A1.Div(&z.A1, norm) + z.A0.Div(z.A0, norm) + z.A1.Div(z.A1, norm) r.Mul(y, z) r.Sub(x, r) return z, r } -// Min returns the minimum of z... with respect to the norm. -func Min(z ...*ComplexNumber) *ComplexNumber { - min := z[0] - for _, v := range z { - if v.Norm().Cmp(min.Norm()) == -1 { - min = v - } - } - return min -} - // HalfGCD returns the rational reconstruction of a, b. // This outputs w, v, u s.t. w = a*u + b*v. func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber { diff --git a/field/eisenstein/eisenstein_test.go b/field/eisenstein/eisenstein_test.go index 89c7dc5d4..6aff795f9 100644 --- a/field/eisenstein/eisenstein_test.go +++ b/field/eisenstein/eisenstein_test.go @@ -246,7 +246,7 @@ func GenNumber(boundSize int64) gopter.Gen { var bound big.Int bound.Exp(big.NewInt(2), big.NewInt(boundSize), nil) elmt, _ := rand.Int(genParams.Rng, &bound) - genResult := gopter.NewGenResult(*elmt, gopter.NoShrinker) + genResult := gopter.NewGenResult(elmt, gopter.NoShrinker) return genResult } } @@ -257,21 +257,23 @@ func GenComplexNumber(boundSize int64) gopter.Gen { GenNumber(boundSize), GenNumber(boundSize), ).Map(func(values []interface{}) *ComplexNumber { - return &ComplexNumber{A0: values[0].(big.Int), A1: values[1].(big.Int)} + return &ComplexNumber{A0: values[0].(*big.Int), A1: values[1].(*big.Int)} }) } // bench +var benchRes [3]*ComplexNumber + func BenchmarkHalfGCD(b *testing.B) { - var prime, _ = new(big.Int).SetString("7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed", 16) // 2^255 - 19 - a0, _ := rand.Int(rand.Reader, prime) - a1, _ := rand.Int(rand.Reader, prime) - c0, _ := rand.Int(rand.Reader, prime) - c1, _ := rand.Int(rand.Reader, prime) - a := ComplexNumber{A0: *a0, A1: *a1} - c := ComplexNumber{A0: *c0, A1: *c1} + var n, _ = new(big.Int).SetString("100000000000000000000000000000000", 16) // 2^128 + a0, _ := rand.Int(rand.Reader, n) + a1, _ := rand.Int(rand.Reader, n) + c0, _ := rand.Int(rand.Reader, n) + c1, _ := rand.Int(rand.Reader, n) + a := ComplexNumber{A0: a0, A1: a1} + c := ComplexNumber{A0: c0, A1: c1} b.ResetTimer() for i := 0; i < b.N; i++ { - _ = HalfGCD(&a, &c) + benchRes = HalfGCD(&a, &c) } }