diff --git a/src/main/java/palaiologos/kamilalisp/runtime/FunctionRegistry.java b/src/main/java/palaiologos/kamilalisp/runtime/FunctionRegistry.java index fbb55a05..655c60ef 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/FunctionRegistry.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/FunctionRegistry.java @@ -23,15 +23,12 @@ import palaiologos.kamilalisp.runtime.math.bit.*; import palaiologos.kamilalisp.runtime.math.cmp.*; import palaiologos.kamilalisp.runtime.math.hyperbolic.*; -import palaiologos.kamilalisp.runtime.math.num.Det; -import palaiologos.kamilalisp.runtime.math.num.PLUDecomposition; +import palaiologos.kamilalisp.runtime.math.num.*; import palaiologos.kamilalisp.runtime.math.prime.IsPrime; import palaiologos.kamilalisp.runtime.math.prime.NextPrime; import palaiologos.kamilalisp.runtime.math.prime.PrimeFactors; import palaiologos.kamilalisp.runtime.math.prime.PrimeNo; import palaiologos.kamilalisp.runtime.math.trig.*; -import palaiologos.kamilalisp.runtime.math.num.LUDecomposition; -import palaiologos.kamilalisp.runtime.math.num.Trace; import palaiologos.kamilalisp.runtime.cas.MatrixTrace; import palaiologos.kamilalisp.runtime.meta.*; import palaiologos.kamilalisp.runtime.net.*; @@ -256,6 +253,8 @@ public static void registerDefault(Environment env) { env.setPrimitive("num:PLU", "⎕⍉↙↗", new Atom(new PLUDecomposition())); env.setPrimitive("num:trace", "⎕∑", new Atom(new Trace())); env.setPrimitive("num:det", "⎕∆", new Atom(new Det())); + env.setPrimitive("num:I", "⎕I", new Atom(new I())); + env.setPrimitive("num:invert", "⎕¯¹", new Atom(new Inv())); env.setPrimitive("cas:matrix:trace", new Atom(new MatrixTrace())); diff --git a/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Base.java b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Base.java index aba76051..a557947c 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Base.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Base.java @@ -338,6 +338,7 @@ public void registerFlt64(Environment env) { env.setPrimitive("flt64:LU", new Atom(new Flt64LU())); env.setPrimitive("flt64:PLU", new Atom(new Flt64PLU())); env.setPrimitive("flt64:det", new Atom(new Flt64Det())); + env.setPrimitive("flt64:invert", new Atom(new Flt64Inv())); env.setPrimitive("flt64:trace", new Atom(new Flt64Trace())); env.setPrimitive("flt64:e", toAtom(Math.E)); env.setPrimitive("flt64:pi", toAtom(Math.PI)); diff --git a/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Det.java b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Det.java index f49226b3..4f3d523c 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Det.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Det.java @@ -9,41 +9,28 @@ import java.util.List; public class Flt64Det extends PrimitiveFunction implements Lambda { - @Override - public Atom apply(Environment env, List args) { - Atom a1 = args.get(0); - if (Rank.computeRank(a1) != 2) { - throw new RuntimeException("Expected a matrix of rank 2."); - } - - List> l1 = a1.getList().stream().map(Atom::getList).toList(); - if (l1.stream().anyMatch(x -> x.size() != l1.get(0).size())) { - throw new RuntimeException("Expected a square matrix."); - } - - double[][] A = l1.stream().map(x -> x.stream().mapToDouble(Flt64Base::toFlt64).toArray()).toArray(double[][]::new); - + public static double det(double[][] A) { if(A.length <= 1) throw new RuntimeException("Expected at least a 2x2 matrix."); else if(A.length == 2) { // 2x2 determinant. - return Flt64Base.toAtom(A[0][0] * A[1][1] - A[0][1] * A[1][0]); + return A[0][0] * A[1][1] - A[0][1] * A[1][0]; } else if(A.length == 3) { // 3x3 determinant. - return Flt64Base.toAtom( + return A[0][0] * A[1][1] * A[2][2] + - A[0][1] * A[1][2] * A[2][0] + - A[0][2] * A[1][0] * A[2][1] - - A[0][2] * A[1][1] * A[2][0] - - A[0][1] * A[1][0] * A[2][2] - - A[0][0] * A[1][2] * A[2][1] - ); + A[0][1] * A[1][2] * A[2][0] + + A[0][2] * A[1][0] * A[2][1] - + A[0][2] * A[1][1] * A[2][0] - + A[0][1] * A[1][0] * A[2][2] - + A[0][0] * A[1][2] * A[2][1] + ; } else { double[][][] lup; try { lup = Flt64PLU.lu(A); } catch(ArithmeticException e) { - return Atom.FALSE; + return 0; } double sumDiagP = 0; @@ -56,12 +43,29 @@ else if(A.length == 2) { prodU *= lup[1][i][i]; double detp = (lup[2].length - sumDiagP - 1) % 2 == 0 ? 1 : -1; - return Flt64Base.toAtom(detp * prodL * prodU); + return detp * prodL * prodU; + } + } + + @Override + public Atom apply(Environment env, List args) { + Atom a1 = args.get(0); + if (Rank.computeRank(a1) != 2) { + throw new RuntimeException("Expected a matrix of rank 2."); } + + List> l1 = a1.getList().stream().map(Atom::getList).toList(); + if (l1.stream().anyMatch(x -> x.size() != l1.get(0).size())) { + throw new RuntimeException("Expected a square matrix."); + } + + double[][] A = l1.stream().map(x -> x.stream().mapToDouble(Flt64Base::toFlt64).toArray()).toArray(double[][]::new); + + return Flt64Base.toAtom(det(A)); } @Override protected String name() { - return null; + return "flt64:det"; } } diff --git a/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Inv.java b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Inv.java new file mode 100644 index 00000000..b7047f1c --- /dev/null +++ b/src/main/java/palaiologos/kamilalisp/runtime/flt64/Flt64Inv.java @@ -0,0 +1,41 @@ +package palaiologos.kamilalisp.runtime.flt64; + +import palaiologos.kamilalisp.atom.Atom; +import palaiologos.kamilalisp.atom.Environment; +import palaiologos.kamilalisp.atom.Lambda; +import palaiologos.kamilalisp.atom.PrimitiveFunction; +import palaiologos.kamilalisp.runtime.array.Rank; + +import java.util.Arrays; +import java.util.List; + +public class Flt64Inv extends PrimitiveFunction implements Lambda { + @Override + public Atom apply(Environment env, List args) { + Atom a1 = args.get(0); + if (Rank.computeRank(a1) != 2) { + throw new RuntimeException("Expected a matrix of rank 2."); + } + + List> l1 = a1.getList().stream().map(Atom::getList).toList(); + if (l1.stream().anyMatch(x -> x.size() != l1.get(0).size())) { + throw new RuntimeException("Expected a square matrix."); + } + + double[][] A = l1.stream().map(x -> x.stream().mapToDouble(Flt64Base::toFlt64).toArray()).toArray(double[][]::new); + + double invdet = 1 / Flt64Det.det(A); + for(int i = 0; i < A.length; i++) { + for(int j = 0; j < A.length; j++) { + A[i][j] *= invdet; + } + } + + return new Atom(Arrays.stream(A).map(x -> new Atom(Arrays.stream(x).mapToObj(Flt64Base::toAtom).toList())).toList()); + } + + @Override + protected String name() { + return "flt64:invert"; + } +} diff --git a/src/main/java/palaiologos/kamilalisp/runtime/math/Slash.java b/src/main/java/palaiologos/kamilalisp/runtime/math/Slash.java index 41576391..d33c3d3a 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/math/Slash.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/math/Slash.java @@ -13,7 +13,7 @@ import java.util.stream.Collectors; public class Slash extends PrimitiveFunction implements Lambda { - private static Atom quot2(Environment e, Atom a, Atom b) { + public static Atom quot2(Environment e, Atom a, Atom b) { a.assertTypes(Type.INTEGER, Type.REAL, Type.COMPLEX, Type.STRING, Type.LIST); b.assertTypes(Type.INTEGER, Type.REAL, Type.COMPLEX, Type.STRING, Type.LIST); if (a.getType() == Type.COMPLEX && b.getType() == Type.COMPLEX) { @@ -41,7 +41,7 @@ private static Atom quot2(Environment e, Atom a, Atom b) { } } - private static Atom quot1(Environment e, Atom a) { + public static Atom quot1(Environment e, Atom a) { a.assertTypes(Type.INTEGER, Type.REAL, Type.COMPLEX, Type.LIST); if (a.getType() == Type.COMPLEX) { return new Atom(a.getComplex().reciprocal(e.getMathContext())); diff --git a/src/main/java/palaiologos/kamilalisp/runtime/math/Star.java b/src/main/java/palaiologos/kamilalisp/runtime/math/Star.java index a3decf6c..7abf2314 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/math/Star.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/math/Star.java @@ -11,7 +11,7 @@ import java.util.stream.Collectors; public class Star extends PrimitiveFunction implements Lambda { - private static Atom multiply2(Atom a, Atom b) { + public static Atom multiply2(Atom a, Atom b) { a.assertTypes(Type.INTEGER, Type.REAL, Type.COMPLEX, Type.STRING, Type.LIST); b.assertTypes(Type.INTEGER, Type.REAL, Type.COMPLEX, Type.STRING, Type.LIST); if (a.getType() == Type.COMPLEX && b.getType() == Type.COMPLEX) { diff --git a/src/main/java/palaiologos/kamilalisp/runtime/math/num/Det.java b/src/main/java/palaiologos/kamilalisp/runtime/math/num/Det.java index e4e71c76..05d2a95d 100644 --- a/src/main/java/palaiologos/kamilalisp/runtime/math/num/Det.java +++ b/src/main/java/palaiologos/kamilalisp/runtime/math/num/Det.java @@ -6,17 +6,11 @@ import java.math.BigDecimal; import java.math.BigInteger; +import java.math.MathContext; import java.util.List; public class Det extends PrimitiveFunction implements Lambda { - @Override - public Atom apply(Environment env, List args) { - Atom a1 = args.get(0); - if (Rank.computeRank(a1) != 2) { - throw new RuntimeException("Expected a matrix of rank 2."); - } - - List> l1 = a1.getList().stream().map(Atom::getList).toList(); + public static Atom det(MathContext mc, List> l1) { if (l1.stream().anyMatch(x -> x.size() != l1.get(0).size())) { throw new RuntimeException("Expected a square matrix."); } @@ -48,7 +42,7 @@ public Atom apply(Environment env, List args) { BigDecimal[][] A = l1.stream().map(x -> x.stream().map(Atom::getReal).toArray(BigDecimal[]::new)).toArray(BigDecimal[][]::new); BigDecimal[][][] lup; try { - lup = PLUDecomposition.lu(env.getMathContext(), A); + lup = PLUDecomposition.lu(mc, A); } catch (ArithmeticException e) { return Atom.FALSE; } @@ -63,7 +57,7 @@ public Atom apply(Environment env, List args) { BigDecimal detp = BigInteger.valueOf(lup[2].length - 1).subtract(sumDiagP).remainder(BigInteger.TWO).compareTo(BigInteger.ZERO) == 0 ? BigDecimal.ONE : BigDecimal.ONE.negate(); BigDecimal result = detp.multiply(prodL).multiply(prodU); - result = result.setScale(env.getMathContext().getPrecision(), env.getMathContext().getRoundingMode()); + result = result.setScale(mc.getPrecision(), mc.getRoundingMode()); return new Atom(result); } } else { @@ -93,7 +87,7 @@ public Atom apply(Environment env, List args) { BigComplex[][] A = l1.stream().map(x -> x.stream().map(Atom::getComplex).toArray(BigComplex[]::new)).toArray(BigComplex[][]::new); BigComplex[][][] lup; try { - lup = PLUDecomposition.lu(env.getMathContext(), A); + lup = PLUDecomposition.lu(mc, A); } catch (ArithmeticException e) { return Atom.FALSE; } @@ -108,13 +102,24 @@ public Atom apply(Environment env, List args) { BigComplex detp = BigInteger.valueOf(lup[2].length - 1).subtract(sumDiagP).remainder(BigInteger.TWO).compareTo(BigInteger.ZERO) == 0 ? BigComplex.ONE : BigComplex.ONE.negate(); BigComplex result = detp.multiply(prodL).multiply(prodU); - BigDecimal re = result.re.setScale(env.getMathContext().getPrecision(), env.getMathContext().getRoundingMode()); - BigDecimal im = result.im.setScale(env.getMathContext().getPrecision(), env.getMathContext().getRoundingMode()); + BigDecimal re = result.re.setScale(mc.getPrecision(), mc.getRoundingMode()); + BigDecimal im = result.im.setScale(mc.getPrecision(), mc.getRoundingMode()); return new Atom(BigComplex.valueOf(re, im)); } } } + @Override + public Atom apply(Environment env, List args) { + Atom a1 = args.get(0); + if (Rank.computeRank(a1) != 2) { + throw new RuntimeException("Expected a matrix of rank 2."); + } + + List> l1 = a1.getList().stream().map(Atom::getList).toList(); + return det(env.getMathContext(), l1); + } + @Override protected String name() { return "num:det"; diff --git a/src/main/java/palaiologos/kamilalisp/runtime/math/num/I.java b/src/main/java/palaiologos/kamilalisp/runtime/math/num/I.java new file mode 100644 index 00000000..047e562e --- /dev/null +++ b/src/main/java/palaiologos/kamilalisp/runtime/math/num/I.java @@ -0,0 +1,31 @@ +package palaiologos.kamilalisp.runtime.math.num; + +import palaiologos.kamilalisp.atom.Atom; +import palaiologos.kamilalisp.atom.Environment; +import palaiologos.kamilalisp.atom.Lambda; +import palaiologos.kamilalisp.atom.PrimitiveFunction; + +import java.math.BigInteger; +import java.util.ArrayList; +import java.util.List; + +public class I extends PrimitiveFunction implements Lambda { + @Override + public Atom apply(Environment env, List args) { + assertArity(args, 1); + int size = args.get(0).getInteger().intValueExact(); + List result = new ArrayList<>(); + for(int i = 0; i < size; i++) { + List row = new ArrayList<>(); + for(int j = 0; j < size; j++) + row.add(new Atom(BigInteger.valueOf(i == j ? 1 : 0))); + result.add(new Atom(row)); + } + return new Atom(result); + } + + @Override + protected String name() { + return "num:I"; + } +} diff --git a/src/main/java/palaiologos/kamilalisp/runtime/math/num/Inv.java b/src/main/java/palaiologos/kamilalisp/runtime/math/num/Inv.java new file mode 100644 index 00000000..816071c2 --- /dev/null +++ b/src/main/java/palaiologos/kamilalisp/runtime/math/num/Inv.java @@ -0,0 +1,27 @@ +package palaiologos.kamilalisp.runtime.math.num; + +import palaiologos.kamilalisp.atom.*; +import palaiologos.kamilalisp.runtime.array.Rank; +import palaiologos.kamilalisp.runtime.math.Slash; +import palaiologos.kamilalisp.runtime.math.Star; + +import java.util.List; + +public class Inv extends PrimitiveFunction implements Lambda { + @Override + public Atom apply(Environment env, List args) { + Atom a1 = args.get(0); + if (Rank.computeRank(a1) != 2) { + throw new RuntimeException("Expected a matrix of rank 2."); + } + + List> l1 = a1.getList().stream().map(Atom::getList).toList(); + Atom result = Det.det(env.getMathContext(), l1); + return Star.multiply2(Slash.quot1(env, result), a1); + } + + @Override + protected String name() { + return "num:invert"; + } +}