Merge pull request #9 from jonas-lj/fft

Fft

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: [email protected]
+    orcid: 'https://orcid.org/0000-0002-1989-3019'

common/src/main/java/dk/jonaslindstrom/ruffini/common/abstractions/Set.java

@@ -8,9 +8,11 @@
 public interface Set<E> {
 
     /**
-     * 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 <code>true</code> if and only if <i>a = b</i> as elements of this set.

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();
+    }
+
+}

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);
     }

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<Pair<BigInteger, List<BigInteger>>> sampleFactoredNumber(BigInteger upperBound, Random random) {
+        List<BigInteger> s = generateDecreasingSequence(upperBound, random);
+        List<BigInteger> 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<BigInteger> generateDecreasingSequence(BigInteger n, Random random) {
+        List<BigInteger> result = new ArrayList<>();
+        while (n.compareTo(BigInteger.ONE) > 0) {
+            n = sampleInclusive(n, random);
+            result.add(n);
+        }
+        return result;
+    }
+
 }

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<BigInteger> 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<Pair<BigInteger, List<BigInteger>>> result = sampleFactoredNumber(upperBound, random);
+            result.ifPresent(System.out::println);
+        }
+
+    }
 }

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 <code>FP</code><i> = F<sub>p</sub></i>.
+     */
+    public static BigPrimeField FP = new BigPrimeField(p);
+
+    /**
+     * <code>FP2</code><i> = FP(u) / (u<sup>2</sup> + 1)</i> is a quadratic field extension of base field <code>FP</code>.
+     */
+    public static AlgebraicFieldExtension<BigInteger, BigPrimeField> FP2 =
+            new AlgebraicFieldExtension<>(FP, "u", Polynomial.of(
+                    FP.identity(),
+                    FP.zero(),
+                    FP.identity()));
+
+    /**
+     * <code>FP6</code><i> = FP2(v) / (v<sup>3</sup> - (u + 9))</i> is a cubic field extension of <code>FP2</code>.
+     */
+    public static AlgebraicFieldExtension<Polynomial<BigInteger>, AlgebraicFieldExtension<BigInteger, BigPrimeField>> FP6 =
+            new AlgebraicFieldExtension<>(FP2, "v", Polynomial.of(
+                    FP2.negate(Polynomial.of(FP.integer(9), FP.identity())),
+                    FP2.zero(),
+                    FP2.zero(),
+                    FP2.identity()));
+
+    /**
+     * <code>FP12</code><i> = FP6(w) / (w<sup>2</sup> - v))</i> is a quadratic field extension of <code>FP6</code>.
+     */
+    public static AlgebraicFieldExtension<Polynomial<Polynomial<BigInteger>>,
+            AlgebraicFieldExtension<Polynomial<BigInteger>, AlgebraicFieldExtension<BigInteger, BigPrimeField>>> FP12 =
+            new AlgebraicFieldExtension<>(FP6, "w", Polynomial.of(
+                    FP6.negate(Polynomial.of(FP2.zero(), FP2.identity())),
+                    FP6.zero(),
+                    FP6.identity()));
+    public static Group<Polynomial<Polynomial<Polynomial<BigInteger>>>> GT = FP12;
+
+    /**
+     * Curve over <code>FP2</code> containing the G2 subgroup.
+     */
+    public static ShortWeierstrassCurveAffine<Polynomial<BigInteger>, ?> G2 =
+            new ShortWeierstrassCurveAffine<>(FP2, FP2.zero(),
+                    Polynomial.of(new BigInteger("19485874751759354771024239261021720505790618469301721065564631296452457478373"),
+                            new BigInteger("266929791119991161246907387137283842545076965332900288569378510910307636690")));
+
+    /**
+     * Curve over <code>FP</code> containing the G1 subgroup.
+     */
+    public static ShortWeierstrassCurveAffine<BigInteger, ?> G1 =
+            new ShortWeierstrassCurveAffine<>(FP, BigInteger.ZERO, BigInteger.valueOf(2));
+    public static AffinePoint<Polynomial<BigInteger>> 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 <i>q</i>.
+     */
+    public static BigPrimeField FQ = new BigPrimeField(q);
+
+    /**
+     * Generator for the G1 subgroup of order <i>q</i>.
+     */
+    public static AffinePoint<BigInteger> 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<Integer> ci = binaryExpansion(t.negate()).stream().map(x -> -x).collect(Collectors.toList());
+    /**
+     * The optimal Ate pairing which is a bilinear function <i>e: G1 x G2 &rarr; GT</i>.
+     */
+    public static BiFunction<AffinePoint<BigInteger>, AffinePoint<Polynomial<BigInteger>>, Polynomial<Polynomial<Polynomial<BigInteger>>>> 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<Polynomial<Polynomial<Polynomial<BigInteger>>>> twist(AffinePoint<BigInteger> 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()))));
+    }
+}

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<BigInteger, BigIntegers> a) {
     public boolean equals(QuadraticForm<BigInteger, BigIntegers> a, QuadraticForm<BigInteger, BigIntegers> b) {
         return a.getA().equals(b.getA()) && a.getB().equals(b.getB()) && a.getC().equals(b.getC());
     }
+
+    public QuadraticForm<BigInteger, BigIntegers> 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<BigInteger, BigIntegers> 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;
+    }
 }