JasPerformance.java

import cc.r2.core.poly.Integers;
import cc.r2.core.poly.IntegersModulo;
import cc.r2.core.poly.multivar.DegreeVector;
import cc.r2.core.poly.multivar.MultivariateGCD;
import cc.r2.core.poly.multivar.MultivariatePolynomial;
import cc.r2.core.poly.univar.Factorization;
import cc.r2.core.poly.univar.UnivariatePolynomial;
import com.wolfram.jlink.KernelLink;
import com.wolfram.jlink.MathLinkFactory;
import edu.jas.arith.BigInteger;
import edu.jas.arith.ModInteger;
import edu.jas.arith.ModIntegerRing;
import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.ufd.FactorAbstract;
import edu.jas.ufd.FactorFactory;
import edu.jas.ufd.GCDFactory;
import edu.jas.ufd.GreatestCommonDivisorAbstract;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well44497b;

import java.util.SortedMap;

public class Main {

    static KernelLink ml;

    private static boolean suppressOutput = false;

    static {
        try {
            ml = MathLinkFactory.createKernelLink("-linkmode launch -linkname '\"/Applications/Mathematica.app/Contents/MacOS/MathKernel\" -mathlink'");
            ml.discardAnswer();
        } catch (Exception e) {}
    }

    private static void print() {
        if (!suppressOutput)
            System.out.println();
    }

    private static void print(Object o) {
        if (!suppressOutput)
            System.out.println(o);
    }

    public static void main(String[] args) throws Exception {
        String[] vars;

        //warm up JVM
        runJAS = false; // otherwise it will took infinite time
        suppressOutput = true;
        vars = new String[]{"x"};
        testFactorizationInZ(vars, 5, 5, 3, 1000);
        vars = new String[]{"x", "y", "z"};
        testGCDInZ(vars, 5, 5, 1000);

        // run benchmark
        runJAS = true;
        suppressOutput = false;
        print("Warmed up");

        vars = new String[]{"x"};
        testFactorizationInZ(vars, 70, 70, 3, 100);

        vars = new String[]{"x"};
        testFactorizationInZp(vars, 17, 250, 250, 3, 100);

        vars = new String[]{"x", "y"};
        testFactorizationInZ(vars, 4, 10, 3, 100);

        vars = new String[]{"x", "y", "z"};
        testGCDInZ(vars, 25, 125, 100);
    }

    private static boolean runJAS = true;

    static void testFactorizationInZ(String[] vars,
                                     int factorsSize, int factorsDegree, int nFactors, int nIterations) {
        GenPolynomialRing<BigInteger> factory = new GenPolynomialRing<BigInteger>(BigInteger.ONE, vars);

        Well44497b rnd = new Well44497b(123);//fix seed
        FactorAbstract<BigInteger> jasFactor = FactorFactory.getImplementation(BigInteger.ONE);
        for (int i = 0; i < nIterations; i++) {
            print("Iteration: " + i);
            GenPolynomial<BigInteger> poly = factory.getONE();
            MultivariatePolynomial<cc.r2.core.number.BigInteger> rPoly = MultivariatePolynomial.one(vars.length, Integers.Integers, DegreeVector.LEX);
            for (int j = 0; j < nFactors; j++) {
                String fString = randomPoly(rnd, factorsSize, factorsDegree, vars);
                GenPolynomial<BigInteger> factor = factory.parse(fString);
                if (factor.isConstant()) {
                    --j;
                    continue;
                }
                MultivariatePolynomial<cc.r2.core.number.BigInteger> rFactor = MultivariatePolynomial.parse(fString, vars);
                poly = poly.multiply(factor);
                rPoly = rPoly.multiply(rFactor);
            }
            print("Input degree: " + poly.degree());

            // run Mathematica
            long[] mma = factorWithMathematica(poly.toString(), -1);
            print("Mathematica time: " + mma[0] + "ms");

            if (vars.length == 1) {
                // run r20x
                long[] r20x = factorWithR20x(rPoly);
                print("r20x time: " + r20x[0] + "ms");
            }

            if (runJAS) {
                // run JAS
                print("Now running JAS...");
                print("Input: "  + poly);
                long start = System.currentTimeMillis();
                SortedMap<GenPolynomial<BigInteger>, Long> map = jasFactor.factors(poly);
                print("JAS #factors: " + map.size());
                print("JAS time: " + (System.currentTimeMillis() - start) + "ms");
                assert map.size() >= nFactors;
            }
            print();
        }
    }

    static void testFactorizationInZp(String[] vars,
                                      long modulus,
                                      int factorsSize, int factorsDegree, int nFactors, int nIterations) {
        ModIntegerRing ring = new ModIntegerRing(modulus);
        GenPolynomialRing<ModInteger> factory = new GenPolynomialRing<ModInteger>(ring, vars);

        Well44497b rnd = new Well44497b(123);//fix seed
        FactorAbstract<ModInteger> jasFactor = FactorFactory.getImplementation(ring);
        for (int i = 0; i < nIterations; i++) {
            print("Iteration: " + i);
            GenPolynomial<ModInteger> poly = factory.getONE();
            IntegersModulo domain = new IntegersModulo(modulus);
            MultivariatePolynomial<cc.r2.core.number.BigInteger> rPoly = MultivariatePolynomial.one(vars.length, domain, DegreeVector.LEX);
            for (int j = 0; j < nFactors; j++) {
                String fString = randomPoly(rnd, factorsSize, factorsDegree, vars);
                GenPolynomial<ModInteger> factor = factory.parse(fString);
                if (factor.isConstant()) {
                    --j;
                    continue;
                }
                MultivariatePolynomial<cc.r2.core.number.BigInteger> rFactor = MultivariatePolynomial.parse(fString, domain, DegreeVector.LEX, vars);
                poly = poly.multiply(factor);
                rPoly = rPoly.multiply(rFactor);
            }
            print("Input degree: " + poly.degree());

            // run Mathematica
            long[] mma = factorWithMathematica(poly.toString(), modulus);
            print("Mathematica time: " + mma[0] + "ms");

            if (vars.length == 1) {
                // run r20x
                long[] r20x = factorWithR20x(rPoly);
                print("r20x time: " + r20x[0] + "ms");
            }

            if (runJAS) {
                // run JAS
                print("Now running JAS...");
                print("Input: "  + poly);
                long start = System.currentTimeMillis();
                SortedMap<GenPolynomial<ModInteger>, Long> map = jasFactor.factors(poly);
                print("JAS #factors: " + map.size());
                print("Time: " + (System.currentTimeMillis() - start) + "ms");
                assert map.size() >= nFactors;
            }
            print();
        }
    }

    static void testGCDInZ(String[] vars,
                           int factorsSize, int factorsDegree, int nIterations) {
        GenPolynomialRing<BigInteger> factory = new GenPolynomialRing<BigInteger>(BigInteger.ONE, vars);

        Well44497b rnd = new Well44497b(123);//fix seed
        for (int i = 0; i < nIterations; i++) {
            print("Iteration: " + i);
            String sPoly = randomPoly(rnd, factorsSize, factorsDegree, vars);
            GenPolynomial<BigInteger> a = factory.parse(sPoly);
            MultivariatePolynomial<cc.r2.core.number.BigInteger> ra = MultivariatePolynomial.parse(sPoly, vars);

            sPoly = randomPoly(rnd, factorsSize, factorsDegree, vars);
            GenPolynomial<BigInteger> b = factory.parse(sPoly);
            MultivariatePolynomial<cc.r2.core.number.BigInteger> rb = MultivariatePolynomial.parse(sPoly, vars);

            sPoly = randomPoly(rnd, factorsSize, factorsDegree, vars);
            GenPolynomial<BigInteger> gcd = factory.parse(sPoly);
            MultivariatePolynomial<cc.r2.core.number.BigInteger> rgcd = MultivariatePolynomial.parse(sPoly, vars);

            a = a.multiply(gcd);
            b = b.multiply(gcd);

            ra = ra.multiply(rgcd);
            rb = rb.multiply(rgcd);

            GreatestCommonDivisorAbstract<BigInteger> gcdFactory = GCDFactory.getImplementation(BigInteger.ONE);


            print("Input a: " + a);
            print("Input b: " + b);
            print("Input gcd: " + gcd);

            // run Mathematica
            long[] mma = gcdWithMathematica(a.toString(), b.toString(), -1);
            print("Mathematica time: " + mma[0] + "ms");

            long[] r20x = gcdWithR20x(ra, rb);
            print("r20x time: " + r20x[0] + "ms");
            assert r20x[1] >= rgcd.size();

            // run JAS
            if (runJAS) {
                print("Now running JAS...");
                print("Input a: "  + a);
                print("Input b: "  + b);
                long start = System.currentTimeMillis();
                GenPolynomial<BigInteger> actual = gcdFactory.gcd(a, b);
                print("JAS gcd: " + actual);
                print("Time: " + (System.currentTimeMillis() - start) + "ms");
                assert actual.degree() >= gcd.degree();
            }
            print();

        }
    }

    static long[] gcdWithMathematica(String a, String b, long modulus) {
        String eval =
                modulus == -1 ? String.format("Timing[PolynomialGCD[%s, %s];]", a, b)
                        : String.format("Timing[PolynomialGCD[%s, %s, Modulus->%s]]", a, b, modulus);
        String r = ml.evaluateToInputForm(eval, -1);
        String[] split = r.replace("{", "").replace("}", "").replace(" ", "").split(",");
        long milliseconds = (long) (1000 * Double.parseDouble(split[0]));
        return new long[]{milliseconds, 0};
    }

    static long[] factorWithMathematica(String expr, long modulus) {
        String eval =
                modulus == -1 ? String.format("Timing[Length[Factor[%s]]]", expr)
                        : String.format("Timing[Length[Factor[%s, Modulus->%s]]]", expr, modulus);
        String r = ml.evaluateToInputForm(eval, -1);
        String[] split = r.replace("{", "").replace("}", "").replace(" ", "").split(",");
        long milliseconds = (long) (1000 * Double.parseDouble(split[0]));
        int nFactors = -1;//Integer.parseInt(split[1]);

        return new long[]{milliseconds, nFactors};
    }

    static long[] factorWithR20x(MultivariatePolynomial<cc.r2.core.number.BigInteger> poly) {
        if (poly.nUsedVariables() > 1)
            throw new IllegalArgumentException("multivariate factorization is not supported: " + poly);
        try {
            long start = System.currentTimeMillis();
            UnivariatePolynomial<cc.r2.core.number.BigInteger> uPoly = poly.asUnivariate();
            int nFactors;
            if (!poly.isOverField())
                nFactors = Factorization.factor(uPoly).size();
            else
                nFactors = Factorization.factor(UnivariatePolynomial.asLongPolyZp(
                        uPoly)).size();
            long milliseconds = System.currentTimeMillis() - start;
            return new long[]{milliseconds, nFactors};
        } catch (Exception e) {
            System.out.println(poly);
            throw new RuntimeException(e);
        }
    }

    static long[] gcdWithR20x(MultivariatePolynomial<cc.r2.core.number.BigInteger> a, MultivariatePolynomial<cc.r2.core.number.BigInteger> b) {
        long start = System.currentTimeMillis();
        int gcdSize = MultivariateGCD.PolynomialGCD(a, b).size();
        long milliseconds = System.currentTimeMillis() - start;
        return new long[]{milliseconds, gcdSize};
    }

    static String randomPoly(RandomGenerator rnd, int factorSize, int factorMaxDegree, String... vars) {
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < factorSize; i++) {
            if (i != 0)
                sb.append((rnd.nextBoolean() && rnd.nextBoolean()) ? "-" : "+");
            int n = rnd.nextInt(100);
            sb.append(n);
            for (String var : vars)
                sb.append("*").append(var).append("^").append(rnd.nextInt(factorMaxDegree));
        }
        return sb.toString();
    }
}