From 106226c443f2ee6871b2fb6c5f81ccda9c799708 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20Lindstr=C3=B8m?= Date: Mon, 8 Jan 2024 09:59:25 +0100 Subject: [PATCH 1/2] Add sampling of qf's and fix fft --- .../ruffini/common/abstractions/Set.java | 6 +- .../ruffini/common/algorithms/BinaryGCD.java | 25 ++++ .../ruffini/common/util/MathUtils.java | 4 + .../ruffini/common/util/SamplingUtils.java | 40 ++++++ common/src/test/java/AlgorithmsTests.java | 48 +++++++ .../elliptic/structures/bn254/BN254.java | 132 ++++++++++++++++++ .../quadraticform/ClassGroup.java | 51 +++++++ .../quadraticform/QuadraticForm.java | 9 +- .../src/test/java/QuadraticFormTests.java | 14 ++ .../integers/structures/limbs/BigElement.java | 16 +++ .../structures/limbs/BigElements.java | 70 ++++++++++ .../polynomials/algorithms/Modulus.java | 56 ++++++++ .../structures/PolynomialRingFFT.java | 129 ++++++++++++++--- .../src/test/java/PolynomialTests.java | 51 +++++-- 14 files changed, 616 insertions(+), 35 deletions(-) create mode 100644 common/src/main/java/dk/jonaslindstrom/ruffini/common/algorithms/BinaryGCD.java create mode 100644 elliptic/src/main/java/dk/jonaslindstrom/ruffini/elliptic/structures/bn254/BN254.java create mode 100644 integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElement.java create mode 100644 integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElements.java create mode 100644 polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/algorithms/Modulus.java diff --git a/common/src/main/java/dk/jonaslindstrom/ruffini/common/abstractions/Set.java b/common/src/main/java/dk/jonaslindstrom/ruffini/common/abstractions/Set.java index 58a0a243..885a7252 100644 --- a/common/src/main/java/dk/jonaslindstrom/ruffini/common/abstractions/Set.java +++ b/common/src/main/java/dk/jonaslindstrom/ruffini/common/abstractions/Set.java @@ -8,9 +8,11 @@ public interface Set { /** - * Returns a human readable string representation of an element in this set. + * Returns a human-readable string representation of an element in this set. */ - String toString(E a); + default String toString(E a) { + return a.toString(); + } /** * Returns true if and only if a = b as elements of this set. diff --git a/common/src/main/java/dk/jonaslindstrom/ruffini/common/algorithms/BinaryGCD.java b/common/src/main/java/dk/jonaslindstrom/ruffini/common/algorithms/BinaryGCD.java new file mode 100644 index 00000000..4419c7c4 --- /dev/null +++ b/common/src/main/java/dk/jonaslindstrom/ruffini/common/algorithms/BinaryGCD.java @@ -0,0 +1,25 @@ +package dk.jonaslindstrom.ruffini.common.algorithms; + +import java.math.BigInteger; + +public class BinaryGCD { + + public static BigInteger apply(BigInteger x, BigInteger y) { + int beta = Math.min(twoValence(x), twoValence(y)); + x = x.shiftRight(beta); + y = y.shiftRight(beta); + + while (!y.equals(x)) { + BigInteger z = x.min(y); + BigInteger d = x.subtract(y).abs(); + y = d.shiftRight(twoValence(d)); + x = z; + } + return x.shiftLeft(beta); + } + + private static int twoValence(BigInteger x) { + return x.getLowestSetBit(); + } + +} diff --git a/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/MathUtils.java b/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/MathUtils.java index 635f8933..da8941df 100644 --- a/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/MathUtils.java +++ b/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/MathUtils.java @@ -32,6 +32,10 @@ public static int floorLog2(int n) { return 31 - Integer.numberOfLeadingZeros(n); } + public static int bigLog2(BigInteger n) { + return n.bitLength() - 1; + } + public static int ceilLog2(int n) { return 32 - Integer.numberOfLeadingZeros(n - 1); } diff --git a/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/SamplingUtils.java b/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/SamplingUtils.java index ae64c105..68ed3a2f 100644 --- a/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/SamplingUtils.java +++ b/common/src/main/java/dk/jonaslindstrom/ruffini/common/util/SamplingUtils.java @@ -1,7 +1,15 @@ package dk.jonaslindstrom.ruffini.common.util; +import org.checkerframework.checker.units.qual.A; + import java.math.BigInteger; +import java.util.ArrayList; +import java.util.List; +import java.util.Optional; import java.util.Random; +import java.util.stream.Collectors; + +import static dk.jonaslindstrom.ruffini.common.util.MathUtils.bigLog2; public class SamplingUtils { @@ -17,4 +25,36 @@ public static BigInteger sample(BigInteger n, Random random) { return r; } + public static Optional>> sampleFactoredNumber(BigInteger upperBound, Random random) { + List s = generateDecreasingSequence(upperBound, random); + List r = s.stream().filter(si -> si.isProbablePrime(30)).toList(); + BigInteger result = r.stream().reduce(BigInteger.ONE, BigInteger::multiply); + + boolean biasCorrected = random.nextDouble() < result.doubleValue() / upperBound.doubleValue(); + if (result.compareTo(upperBound) > 0 || !biasCorrected) { + return Optional.empty(); + } else { + return Optional.of(Pair.of(result, r)); + } + } + + /** + * Sample a random BigInteger in the range [0,n] + */ + public static BigInteger sampleInclusive(BigInteger n, Random random) { + return sample(n.add(BigInteger.ONE), random); + } + + /** + * Generate a decreasing list of numbers s_0, s_1, ..., s_k such that s_0 \leq n and s_k = 1. + */ + public static List generateDecreasingSequence(BigInteger n, Random random) { + List result = new ArrayList<>(); + while (n.compareTo(BigInteger.ONE) > 0) { + n = sampleInclusive(n, random); + result.add(n); + } + return result; + } + } diff --git a/common/src/test/java/AlgorithmsTests.java b/common/src/test/java/AlgorithmsTests.java index 43a3cee8..fa8f66b9 100644 --- a/common/src/test/java/AlgorithmsTests.java +++ b/common/src/test/java/AlgorithmsTests.java @@ -2,6 +2,7 @@ import dk.jonaslindstrom.ruffini.common.abstractions.EuclideanDomain; import dk.jonaslindstrom.ruffini.common.abstractions.Ring; import dk.jonaslindstrom.ruffini.common.algorithms.*; +import dk.jonaslindstrom.ruffini.common.util.Pair; import dk.jonaslindstrom.ruffini.common.util.SamplingUtils; import dk.jonaslindstrom.ruffini.common.util.TestUtils; import dk.jonaslindstrom.ruffini.common.vector.Vector; @@ -11,8 +12,12 @@ import java.math.BigInteger; import java.util.ArrayList; import java.util.List; +import java.util.Optional; import java.util.Random; +import static dk.jonaslindstrom.ruffini.common.util.SamplingUtils.generateDecreasingSequence; +import static dk.jonaslindstrom.ruffini.common.util.SamplingUtils.sampleFactoredNumber; + public class AlgorithmsTests { @Test @@ -142,4 +147,47 @@ public void testEuclideanAlgorithm() { } } + @Test + public void testBinaryGCD() { + Random random = new Random(1234); + + int tests = 100; + for (int i = 0; i < tests; i++) { + BigInteger a = new BigInteger(16, random); + BigInteger b = new BigInteger(16, random); + BigInteger gcd = BinaryGCD.apply(a,b); + Assert.assertEquals(a.gcd(b), gcd); + } + } + + + @Test + public void countQuadraticResidues() { + int n = 1000001; + // 1000001 = 101 * 9901 + int c = 0; + for (int i = 1; i < n; i++) { + int j = new JacobiSymbol().applyAsInt(i, n); + } + } + + @Test + public void testDecreasingSequence() { + Random random = new Random(1234); + BigInteger upperBound = BigInteger.valueOf(100); + List result = generateDecreasingSequence(upperBound, random); + System.out.println(result); + } + + @Test + public void testSampleFactored() { + Random random = new Random(1234); + BigInteger upperBound = BigInteger.valueOf(100000); + + for (int i = 0; i < 1000; i++) { + Optional>> result = sampleFactoredNumber(upperBound, random); + result.ifPresent(System.out::println); + } + + } } diff --git a/elliptic/src/main/java/dk/jonaslindstrom/ruffini/elliptic/structures/bn254/BN254.java b/elliptic/src/main/java/dk/jonaslindstrom/ruffini/elliptic/structures/bn254/BN254.java new file mode 100644 index 00000000..482a4994 --- /dev/null +++ b/elliptic/src/main/java/dk/jonaslindstrom/ruffini/elliptic/structures/bn254/BN254.java @@ -0,0 +1,132 @@ +package dk.jonaslindstrom.ruffini.elliptic.structures.bn254; + +import dk.jonaslindstrom.ruffini.common.abstractions.Group; +import dk.jonaslindstrom.ruffini.common.util.SamePair; +import dk.jonaslindstrom.ruffini.elliptic.algorithms.OptimalAtePairing; +import dk.jonaslindstrom.ruffini.elliptic.elements.AffinePoint; +import dk.jonaslindstrom.ruffini.elliptic.structures.ShortWeierstrassCurveAffine; +import dk.jonaslindstrom.ruffini.finitefields.AlgebraicFieldExtension; +import dk.jonaslindstrom.ruffini.finitefields.BigPrimeField; +import dk.jonaslindstrom.ruffini.polynomials.elements.Polynomial; + +import java.math.BigInteger; +import java.util.List; +import java.util.function.BiFunction; +import java.util.stream.Collectors; + +import static dk.jonaslindstrom.ruffini.common.util.MathUtils.binaryExpansion; + +/** + * Implementation of the BN254 (aka BN128) pairing-friendly elliptic curve construction. + */ +public class BN254 { + + /** + * Modulus of the base field. + */ + public static BigInteger p = new BigInteger("21888242871839275222246405745257275088696311157297823662689037894645226208583", 10); + + /** + * The base field FP = Fp. + */ + public static BigPrimeField FP = new BigPrimeField(p); + + /** + * FP2 = FP(u) / (u2 + 1) is a quadratic field extension of base field FP. + */ + public static AlgebraicFieldExtension FP2 = + new AlgebraicFieldExtension<>(FP, "u", Polynomial.of( + FP.identity(), + FP.zero(), + FP.identity())); + + /** + * FP6 = FP2(v) / (v3 - (u + 9)) is a cubic field extension of FP2. + */ + public static AlgebraicFieldExtension, AlgebraicFieldExtension> FP6 = + new AlgebraicFieldExtension<>(FP2, "v", Polynomial.of( + FP2.negate(Polynomial.of(FP.integer(9), FP.identity())), + FP2.zero(), + FP2.zero(), + FP2.identity())); + + /** + * FP12 = FP6(w) / (w2 - v)) is a quadratic field extension of FP6. + */ + public static AlgebraicFieldExtension>, + AlgebraicFieldExtension, AlgebraicFieldExtension>> FP12 = + new AlgebraicFieldExtension<>(FP6, "w", Polynomial.of( + FP6.negate(Polynomial.of(FP2.zero(), FP2.identity())), + FP6.zero(), + FP6.identity())); + public static Group>>> GT = FP12; + + /** + * Curve over FP2 containing the G2 subgroup. + */ + public static ShortWeierstrassCurveAffine, ?> G2 = + new ShortWeierstrassCurveAffine<>(FP2, FP2.zero(), + Polynomial.of(new BigInteger("19485874751759354771024239261021720505790618469301721065564631296452457478373"), + new BigInteger("266929791119991161246907387137283842545076965332900288569378510910307636690"))); + + /** + * Curve over FP containing the G1 subgroup. + */ + public static ShortWeierstrassCurveAffine G1 = + new ShortWeierstrassCurveAffine<>(FP, BigInteger.ZERO, BigInteger.valueOf(2)); + public static AffinePoint> G2_GENERATOR = new AffinePoint<>( + Polynomial.of( + new BigInteger("61A10BB519EB62FEB8D8C7E8C61EDB6A4648BBB4898BF0D91EE4224C803FB2B", 16), + new BigInteger("516AAF9BA737833310AA78C5982AA5B1F4D746BAE3784B70D8C34C1E7D54CF3", 16)), + Polynomial.of( + new BigInteger("21897A06BAF93439A90E096698C822329BD0AE6BDBE09BD19F0E07891CD2B9A", 16), + new BigInteger("EBB2B0E7C8B15268F6D4456F5F38D37B09006FFD739C9578A2D1AEC6B3ACE9B", 16)) + ); + /** + * Order of subgroups of G1, G2 and GT + */ + public static BigInteger q = new BigInteger("21888242871839275222246405745257275088548364400416034343698204186575808495617", 10); + + /** + * Prime field of order q. + */ + public static BigPrimeField FQ = new BigPrimeField(q); + + /** + * Generator for the G1 subgroup of order q. + */ + public static AffinePoint G1_GENERATOR = new AffinePoint<>( + p.subtract(BigInteger.ONE), + BigInteger.ONE); + + /** + * The parameter t + */ + private static BigInteger t = new BigInteger("-4080000000000001", 16); + + /** + * Binary expansion of t + */ + private static List ci = binaryExpansion(t.negate()).stream().map(x -> -x).collect(Collectors.toList()); + /** + * The optimal Ate pairing which is a bilinear function e: G1 x G2 → GT. + */ + public static BiFunction, AffinePoint>, Polynomial>>> PAIRING = + (a, b) -> new OptimalAtePairing<>( + g1 -> FP2.embed(g1), + G2, + g2 -> FP12.embed(FP6.embed(g2)), + FP12, + BN254::twist, + p, q, 12).pairing(a, b, ci); + + public static SamePair>>> twist(AffinePoint p) { + return new SamePair<>( + Polynomial.of( + Polynomial.of(FP2.zero(), FP2.zero(), Polynomial.of(p.x())), + Polynomial.constant(FP2.zero())), + Polynomial.of( + Polynomial.constant(FP2.zero()), + Polynomial.constant(Polynomial.of(p.y())))); + } +} diff --git a/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/ClassGroup.java b/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/ClassGroup.java index c0992a0f..200c160e 100644 --- a/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/ClassGroup.java +++ b/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/ClassGroup.java @@ -1,9 +1,15 @@ package dk.jonaslindstrom.ruffini.finitefields.quadraticform; import dk.jonaslindstrom.ruffini.common.abstractions.Group; +import dk.jonaslindstrom.ruffini.common.algorithms.BigLegendreSymbol; +import dk.jonaslindstrom.ruffini.common.algorithms.ChineseRemainderTheorem; +import dk.jonaslindstrom.ruffini.common.util.SamplingUtils; +import dk.jonaslindstrom.ruffini.common.vector.Vector; +import dk.jonaslindstrom.ruffini.integers.algorithms.ModularSquareRoot; import dk.jonaslindstrom.ruffini.integers.structures.BigIntegers; import java.math.BigInteger; +import java.util.Random; /** * This implements the ideal class group for a negative discriminant. @@ -17,6 +23,7 @@ public ClassGroup(BigInteger discriminant) { if (discriminant.compareTo(BigInteger.ZERO) >= 0) { throw new IllegalArgumentException("Discriminant must be negative"); } + this.discriminant = discriminant; // Compute identity @@ -61,4 +68,48 @@ public String toString(QuadraticForm a) { public boolean equals(QuadraticForm a, QuadraticForm b) { return a.getA().equals(b.getA()) && a.getB().equals(b.getB()) && a.getC().equals(b.getC()); } + + public QuadraticForm sample(Random random) { + if (discriminant.mod(BigInteger.valueOf(4)).intValue() != 1 || !discriminant.abs().isProbablePrime(100)) { + throw new IllegalStateException("Discriminant must be prime and 1 mod 4"); + } + + // Note: The number of candidates can be computed with the prime number theorem. + + // If a is smaller than this bound and |b| < a, the form is guaranteed to be reduced. + BigInteger bound = discriminant.abs().sqrt().shiftRight(1); + + // Sample a prime a such that a = 3 mod 4 and the discriminant is a quadratic residue mod a + // TODO: Sample a with factorization and use this to compute the legendre symbol of the discriminant (https://pages.cs.wisc.edu/~cs812-1/pfrn.pdf or https://inst.eecs.berkeley.edu/~cs276/fa20/notes/Kalai_generating_factored.pdf) + BigInteger a; + int iterations = 0; + do { + a = SamplingUtils.sample(bound.shiftRight(2), random).shiftLeft(2).add(BigInteger.valueOf(3)); + // TODO: Once a square root algorithm is implemented, we may allow a to be 1 mod 4 also + iterations++; + // Note: We could also check legendre(a, discriminant.negate()) != 1 due to the law of quadratic reciprocity + } while (legendre(discriminant, a) != 1 || !a.isProbablePrime(100)); + System.out.println("Iterations: " + iterations); + + // Compute square root of the discriminant mod a + BigInteger bPrime = new ModularSquareRoot(a).apply(discriminant); + BigInteger b = bPrime.multiply(a.add(BigInteger.ONE)).subtract(a).mod(a.shiftLeft(2)); + + // The "normalization operator" replaces b with b-2a until b <= a. + // Here, b is < 4a, so we only need to do this at most twice. + while (b.compareTo(a) > 0) { + b = b.subtract(a.shiftLeft(1)); + } + + BigInteger c = b.multiply(b).subtract(discriminant).divide(a.shiftLeft(2)); + QuadraticForm result = new QuadraticForm<>(BigIntegers.getInstance(), a, b, c); + System.out.println("Reduced: " + result.isReduced()); + System.out.println("Discriminant: " + result.discriminant().equals(discriminant)); + return result.reduce(); + } + + private static int legendre(BigInteger a, BigInteger p) { + BigInteger l = a.modPow(p.subtract(BigInteger.ONE).divide(BigInteger.TWO), p); + return l.equals(BigInteger.ONE) ? 1 : -1; + } } diff --git a/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/QuadraticForm.java b/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/QuadraticForm.java index 33bfae2e..312668cf 100644 --- a/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/QuadraticForm.java +++ b/finite-fields/src/main/java/dk/jonaslindstrom/ruffini/finitefields/quadraticform/QuadraticForm.java @@ -55,17 +55,18 @@ public boolean isPositiveDefinite() { } public boolean isReduced() { - E absB = ring.abs(b, ordering); - if (ring.greaterThan(absB, a)) { + if (!isNormal()) { return false; } + if (ring.greaterThan(a, c)) { return false; } - if (ring.equals(a, absB) || ring.equals(a, c)) { + + if (ring.equals(a, c)) { return ring.greaterThanOrEqual(b, ring.zero()); } - return false; + return true; } private boolean isNormal() { diff --git a/finite-fields/src/test/java/QuadraticFormTests.java b/finite-fields/src/test/java/QuadraticFormTests.java index c57d948c..352c6b49 100644 --- a/finite-fields/src/test/java/QuadraticFormTests.java +++ b/finite-fields/src/test/java/QuadraticFormTests.java @@ -6,6 +6,7 @@ import org.junit.Test; import java.math.BigInteger; +import java.util.Random; public class QuadraticFormTests { @@ -85,4 +86,17 @@ public void testLargeGroup() { } Assert.assertEquals(group.identity(), power.apply(P, 7)); } + + @Test + public void testSampling() { + BigInteger d = BigInteger.valueOf(-49513019); + ClassGroup group = new ClassGroup(d); + + Random random = new Random(1234); + for (int i = 0; i < 100; i++) { + System.out.println(); + QuadraticForm x = group.sample(random); + System.out.println(x); + } + } } diff --git a/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElement.java b/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElement.java new file mode 100644 index 00000000..e9dc613d --- /dev/null +++ b/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElement.java @@ -0,0 +1,16 @@ +package dk.jonaslindstrom.ruffini.integers.structures.limbs; + +import java.util.List; + +public class BigElement { + + final List limbs; + + public BigElement(List limbs) { + this.limbs = limbs; + } + + public String toString() { + return limbs.toString(); + } +} diff --git a/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElements.java b/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElements.java new file mode 100644 index 00000000..faf92c89 --- /dev/null +++ b/integers/src/main/java/dk/jonaslindstrom/ruffini/integers/structures/limbs/BigElements.java @@ -0,0 +1,70 @@ +package dk.jonaslindstrom.ruffini.integers.structures.limbs; + +import dk.jonaslindstrom.ruffini.common.abstractions.EuclideanDomain; +import dk.jonaslindstrom.ruffini.common.abstractions.Ring; +import dk.jonaslindstrom.ruffini.common.util.Pair; + +import java.util.ArrayList; +import java.util.List; + +public class BigElements implements Ring> { + + private final EuclideanDomain ring; + private final E modulus; + + public BigElements(EuclideanDomain ring, E modulus) { + this.ring = ring; + this.modulus = modulus; + } + + @Override + public BigElement negate(BigElement a) { + return new BigElement<>(a.limbs.stream().map(ai -> ring.subtract(modulus, ai)).toList()); + } + + @Override + public BigElement add(BigElement a, BigElement b) { + E carry = ring.zero(); + List limbs = new ArrayList<>(); + for (int i = 0; i < Math.max(a.limbs.size(), b.limbs.size()); i++) { + E ai = i < a.limbs.size() ? a.limbs.get(i) : ring.zero(); + E bi = i < b.limbs.size() ? b.limbs.get(i) : ring.zero(); + E di = ring.add(ai, bi, carry); + Pair qr = ring.divide(di, modulus); + carry = qr.first; + limbs.add(qr.second); + } + if (!ring.isZero(ring.zero())) { + limbs.add(carry); + } + return new BigElement<>(limbs); + } + + @Override + public BigElement zero() { + return new BigElement<>(List.of(ring.zero())); + } + + @Override + public BigElement identity() { + return new BigElement<>(List.of(ring.identity())); + } + + @Override + public BigElement multiply(BigElement a, BigElement b) { + return null; + } + + @Override + public boolean equals(BigElement a, BigElement b) { + if (a.limbs.size() != b.limbs.size()) { + return false; + } + for (int i = 0; i < a.limbs.size(); i++) { + if (!ring.equals(a.limbs.get(i), b.limbs.get(i))) { + return false; + } + } + return true; + } +} diff --git a/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/algorithms/Modulus.java b/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/algorithms/Modulus.java new file mode 100644 index 00000000..7b691b7f --- /dev/null +++ b/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/algorithms/Modulus.java @@ -0,0 +1,56 @@ +//package dk.jonaslindstrom.ruffini.polynomials.algorithms; +// +//import dk.jonaslindstrom.ruffini.common.util.Pair; +//import dk.jonaslindstrom.ruffini.common.vector.Vector; +//import dk.jonaslindstrom.ruffini.polynomials.elements.Polynomial; +//import dk.jonaslindstrom.ruffini.polynomials.structures.PolynomialRingFFT; +//import dk.jonaslindstrom.ruffini.polynomials.structures.PolynomialRingOverRing; +// +//import java.util.function.BiFunction; +//import java.util.function.UnaryOperator; +// +///** +// * Algorthm 9.5 from Modern Computer Algebra. This is asymptotically faster than the usual division algorithm, and in +// * particular if the degree of the divisor is almost as large as the degree of the dividend. +// */ +//public class Modulus implements UnaryOperator> { +// +// private final PolynomialRingFFT ring; +// private final PolynomialRingFFT.TransformedPolynomial b; +// private final Vector.TransformedPolynomial> inverses; +// +// public Modulus(PolynomialRingFFT ring, PolynomialRingFFT.TransformedPolynomial b, int maxDegree) { +// this.ring = ring; +// this.b = b; +// // Compute the inverse of b(x) mod x^{m+1} +// this.inverses = Vector.of(maxDegree, i -> ring.fromPolynomial(new Inversion<>(ring.getField()).apply(ring.toPolynomial(b).reverse(), i+1))); +// } +// +// @Override +// public Polynomial apply(Polynomial a) { +// if (a.degree() < b.degree()) { +// return a; +// } +// +// int m = a.degree() - b.degree(); +// +// Polynomial bReverseInverse = inverses.get(m + b.degree()); +// +// Polynomial.Builder builder = new Polynomial.Builder<>(ring.getRing()); +// a.reverse().forEach((i, ai) -> { +// bReverseInverse.forEach((j, bj) -> { +// if (i+j < m+1) { +// builder.addTo(i + j, ring.getRing().multiply(ai, bj)); +// } +// }); +// }); +// Polynomial q = builder.build().reverse(); +// +// return ring.subtract(a, ring.multiply(b, q)); +// } +// +// @Override +// public Polynomial apply(Polynomial polynomial) { +// return null; +// } +//} diff --git a/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/structures/PolynomialRingFFT.java b/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/structures/PolynomialRingFFT.java index 51fcc4c6..39e836a1 100644 --- a/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/structures/PolynomialRingFFT.java +++ b/polynomials/src/main/java/dk/jonaslindstrom/ruffini/polynomials/structures/PolynomialRingFFT.java @@ -1,37 +1,132 @@ package dk.jonaslindstrom.ruffini.polynomials.structures; import dk.jonaslindstrom.ruffini.common.abstractions.Field; +import dk.jonaslindstrom.ruffini.common.abstractions.Ring; import dk.jonaslindstrom.ruffini.common.algorithms.DiscreteFourierTransform; import dk.jonaslindstrom.ruffini.common.algorithms.InverseDiscreteFourierTransform; import dk.jonaslindstrom.ruffini.common.vector.Vector; import dk.jonaslindstrom.ruffini.polynomials.elements.Polynomial; -import java.util.SortedMap; +import java.util.ArrayList; +import java.util.List; /** * This class implements the ring of polynomials K[x] over a field K. */ -public class PolynomialRingFFT extends PolynomialRing { +public class PolynomialRingFFT implements Ring.TransformedPolynomial> { - private final SortedMap rootsOfUnity; + private final int n; + private final DiscreteFourierTransform fft; + private final InverseDiscreteFourierTransform ifft; + private final Field field; + private final E root; - public PolynomialRingFFT(Field field, SortedMap rootsOfUnity) { - super(field); - this.rootsOfUnity = rootsOfUnity; + public PolynomialRingFFT(Field field, E rootOfUnity, int n) { + this.field = field; + this.n = n; + this.root = rootOfUnity; + this.fft = new DiscreteFourierTransform<>(field, rootOfUnity, n); + this.ifft = new InverseDiscreteFourierTransform<>(field, rootOfUnity, n); + } + + public Field getField() { + return field; + } + + /** Transformation of x^m */ + public TransformedPolynomial monomial(int m) { + E am = field.power(root, m); + List coefficients = new ArrayList<>(); + coefficients.add(field.identity()); + coefficients.add(am); + E current = am; + for (int i = 2; i < n; i++) { + current = field.multiply(current, am); + coefficients.add(current); + } + return new TransformedPolynomial(Vector.ofList(coefficients)); + } + + public TransformedPolynomial fromConstant(E a) { + return new TransformedPolynomial(Vector.of(n, i -> a)); + } + + public TransformedPolynomial scale(TransformedPolynomial a, E b) { + return new TransformedPolynomial(a.coefficients.map(c -> field.multiply(c, b))); + } + + public TransformedPolynomial fromPolynomial(Polynomial a) { + return new TransformedPolynomial(a); + } + + public TransformedPolynomial fromFourierCoefficients(Vector a) { + return new TransformedPolynomial(a); + } + + public Polynomial toPolynomial(TransformedPolynomial a) { + return a.toPolynomial(); + } + + @Override + public TransformedPolynomial negate(TransformedPolynomial a) { + return new TransformedPolynomial(a.coefficients.map(field::negate)); + } + + @Override + public TransformedPolynomial add(TransformedPolynomial a, TransformedPolynomial b) { + return new TransformedPolynomial(a.coefficients.coordinateWise(b.coefficients, field::add)); + } + + @Override + public TransformedPolynomial zero() { + return new TransformedPolynomial(Vector.of(n, i -> field.zero())); } @Override - public Polynomial multiply(Polynomial a, Polynomial b) { - int n = a.degree() + b.degree() + 1; - if (rootsOfUnity.containsKey(n)) { - DiscreteFourierTransform DFT = new DiscreteFourierTransform<>(field, rootsOfUnity.get(n), n); - Vector â = DFT.apply(a.vector(this.field.zero())); - Vector b̂ = DFT.apply(b.vector(this.field.zero())); - Vector ĉ = â.coordinateWise(b̂, this.field::multiply); - InverseDiscreteFourierTransform inverseDFT = new InverseDiscreteFourierTransform<>(field, rootsOfUnity.get(n), n); - return new Polynomial<>(inverseDFT.apply(ĉ), this.field); - } else { - return super.multiply(a, b); + public TransformedPolynomial identity() { + return new TransformedPolynomial(Vector.of(n, i -> field.identity())); + } + + @Override + public TransformedPolynomial multiply(TransformedPolynomial a, TransformedPolynomial b) { + return new TransformedPolynomial(a.coefficients.coordinateWise(b.coefficients, field::multiply)); + } + + @Override + public String toString(TransformedPolynomial a) { + return a.toString(); + } + + @Override + public boolean equals(TransformedPolynomial a, TransformedPolynomial b) { + return a.coefficients.equals(b.coefficients); + } + + + public class TransformedPolynomial { + + private final Vector coefficients; + + public TransformedPolynomial(Polynomial from) { + this(fft.apply(from.vector(field.zero()))); + } + + private TransformedPolynomial(Vector coefficients) { + this.coefficients = coefficients; + } + + public Vector getCoefficients() { + return coefficients; } + + public Polynomial toPolynomial() { + return new Polynomial<>(ifft.apply(this.coefficients), field); + } + + public String toString() { + return "TransformedPolynomial(" + coefficients + ")"; + } + } + } diff --git a/polynomials/src/test/java/PolynomialTests.java b/polynomials/src/test/java/PolynomialTests.java index 7457aeb7..59616190 100644 --- a/polynomials/src/test/java/PolynomialTests.java +++ b/polynomials/src/test/java/PolynomialTests.java @@ -2,6 +2,7 @@ import dk.jonaslindstrom.ruffini.common.helpers.PerformanceLoggingField; import dk.jonaslindstrom.ruffini.common.util.Pair; import dk.jonaslindstrom.ruffini.common.util.TestUtils; +import dk.jonaslindstrom.ruffini.common.vector.Vector; import dk.jonaslindstrom.ruffini.polynomials.algorithms.*; import dk.jonaslindstrom.ruffini.polynomials.elements.Polynomial; import dk.jonaslindstrom.ruffini.polynomials.structures.PolynomialRing; @@ -63,22 +64,20 @@ public void testInterpolation() { @Test public void testFFT() { - PerformanceLoggingField field = new PerformanceLoggingField<>(new TestUtils.TestField(13)); - - Polynomial p = Polynomial.of(2, 4, 5); - Polynomial q = Polynomial.of(1, 2, 3, 4); + PerformanceLoggingField field = new PerformanceLoggingField<>(new TestUtils.TestField(17)); - PolynomialRing polynomialRingFft = new PolynomialRingFFT<>(field, ImmutableSortedMap.of(3, 3, 2, 12, 6, 4, 4, 5)); - Polynomial actual = polynomialRingFft.multiply(p, q); - System.out.println("With FFT"); - System.out.println(field); - field.reset(); + Polynomial p = Polynomial.of(1, 2); + Polynomial q = Polynomial.of(4, 1, 2); PolynomialRing polynomialRing = new PolynomialRing<>(field); Polynomial expected = polynomialRing.multiply(p, q); - System.out.println(); - System.out.println("Without FFT"); - System.out.println(field); + + PolynomialRingFFT polynomialRingFft = new PolynomialRingFFT<>(field, 3, 16); + PolynomialRingFFT.TransformedPolynomial pHat = polynomialRingFft.fromPolynomial(p); + PolynomialRingFFT.TransformedPolynomial qHat = polynomialRingFft.fromPolynomial(q); + PolynomialRingFFT.TransformedPolynomial actualHat = polynomialRingFft.multiply(pHat, qHat); + Polynomial actual = polynomialRingFft.toPolynomial(actualHat); + Assert.assertEquals(expected, actual); } @@ -137,4 +136,32 @@ public void testFastDivision() { } + @Test + public void testFftDivision() { + PerformanceLoggingField field = new PerformanceLoggingField<>(new TestUtils.TestField(17)); + + Polynomial p = Polynomial.of(1, 2, 3, 4, 5); + Polynomial q = Polynomial.of(2, 3, 4, 5, 1); + + PolynomialRing polynomialRing = new PolynomialRing<>(field); + Pair, Polynomial> expected = polynomialRing.divide(p, q); + + PolynomialRingFFT polynomialRingFft = new PolynomialRingFFT<>(field, 3, 16); + + PolynomialRingFFT.TransformedPolynomial pHat = polynomialRingFft.fromPolynomial(p); + PolynomialRingFFT.TransformedPolynomial qHat = polynomialRingFft.fromPolynomial(q); + + System.out.println(qHat); + + Vector actualCoefficients = pHat.getCoefficients().coordinateWise(qHat.getCoefficients(), field::divide); + PolynomialRingFFT.TransformedPolynomial actualHat = polynomialRingFft.fromFourierCoefficients(actualCoefficients); + Polynomial actual = polynomialRingFft.toPolynomial(actualHat); + + System.out.println(polynomialRingFft.multiply(qHat, actualHat)); + + System.out.println(expected); + System.out.println(actual); + + } + } From 9b6b566ef0ca99fd65225824494f97b73d8a0f30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20Lindstr=C3=B8m?= Date: Mon, 8 Jan 2024 10:03:48 +0100 Subject: [PATCH 2/2] Add citation.cff --- CITATION.cff | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 00000000..ee310111 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,12 @@ +# This CITATION.cff file was generated with cffinit. +# Visit https://bit.ly/cffinit to generate yours today! + +cff-version: 1.2.0 +title: Ruffini +message: Computations over algebraic structures in Java made easy +type: software +authors: + - given-names: Jonas + family-names: Lindstrøm + email: mail@jonaslindstrom.dk + orcid: 'https://orcid.org/0000-0002-1989-3019'