Skip to content

Instantly share code, notes, and snippets.

@jburse

jburse/Matrice.java

Last active Mar 26, 2017
Embed
What would you like to do?
Matrice Operations via Prolog Object Orientation
/**
* Vector/Matrice Operatios Prolog OO Prototype for Jekejeke Prolog.
*
* Warranty & Liability
* To the extent permitted by applicable law and unless explicitly
* otherwise agreed upon, XLOG Technologies GmbH makes no warranties
* regarding the provided information. XLOG Technologies GmbH assumes
* no liability that any problems might be solved with the information
* provided by XLOG Technologies GmbH.
*
* Rights & License
* All industrial property rights regarding the information - copyright
* and patent rights in particular - are the sole property of XLOG
* Technologies GmbH. If the company was not the originator of some
* excerpts, XLOG Technologies GmbH has at least obtained the right to
* reproduce, change and translate the information.
*
* Reproduction is restricted to the whole unaltered document. Reproduction
* of the information is only allowed for non-commercial uses. Selling,
* giving away or letting of the execution of the library is prohibited.
* The library can be distributed as part of your applications and libraries
* for execution provided this comment remains unchanged.
*
* Trademarks
* Jekejeke is a registered trademark of XLOG Technologies GmbH.
*/
% ?- sys_add_path('file:/Projects/Jekejeke/Prototyping/experiment/other/').
:- package(library(diffuse/cas/matrice)).
:- module(generic, []).
/* for the residue hooks */
:- sys_auto_load(matrice).
:- sys_auto_load(number).
:- sys_auto_load(vector).
/*********************************************************************/
/* Arithmetic */
/*********************************************************************/
/**
* X is E:
* The predicate succeeds in evaluating E by using polymorphism.
*/
:- override is/2.
:- public is/2.
X is E :- number(E), !,
X = E.
X is [] :- !,
X = vector.
X is {F | G} :- !,
findall(Y, (G, Y is F), L),
X is L.
X is E :- atom(E), !,
user:(X is E).
X is F :- functor(F, vector, _), !,
X = F.
X is F :- functor(F, matrice, _), !,
X = F.
X is E :- compound(E), !,
E =.. [F|L],
sys_eval_list(L, X, [Y|T]),
sys_poly_send(Y, F, T).
/**
* sys_eval_list(L, X, R):
* The predicate succeeds in R for evaluating the elements of L
* and adding the place holder X.
*/
:- private sys_eval_list/3.
sys_eval_list([X|Y], H, [Z|T]) :-
Z is X,
sys_eval_list(Y, H, T).
sys_eval_list([], H, [H]).
/*********************************************************************/
/* Comparison */
/*********************************************************************/
/**
* E =:= F:
* The predicate succeeds when evaluating E and F by using
* polymorphism gives the same result.
*/
:- override (=:=)/2.
:- public (=:=)/2.
E =:= F :-
X is E, Y is F,
sys_poly_send(X, gen_eq, [Y]).
/**
* E =\= F:
* The predicate succeeds when evaluating E and F by using
* polymorphism dont give the same result.
*/
:- override (=\=)/2.
:- public (=\=)/2.
E =\= F :-
X is E, Y is F,
\+ sys_poly_send(X, gen_eq, [Y]).
/**
* E < F:
* The preicate succeeds when evaluating E by using polymorphism
* is less than evaluating F by using polymorphism.
*/
:- override (<)/2.
:- public (<)/2.
E < F :-
X is E, Y is F,
sys_poly_send(X, gen_ls, [Y]).
/**
* E =< F:
* The preicate succeeds when evaluating E by using polymorphism
* is less or equal than evaluating F by using polymorphism.
*/
:- override (=<)/2.
:- public (=<)/2.
E =< F :-
X is E, Y is F,
\+ sys_poly_send(Y, gen_ls, [X]).
/**
* E > F:
* The preicate succeeds when evaluating E by using polymorphism
* is greater than evaluating F by using polymorphism.
*/
:- override (>)/2.
:- public (>)/2.
E > F :-
X is E, Y is F,
sys_poly_send(Y, gen_ls, [X]).
/**
* E >= F:
* The preicate succeeds when evaluating E by using polymorphism
* is greater or equal than evaluating F by using polymorphism.
*/
:- override (>=)/2.
:- public (>=)/2.
E >= F :-
X is E, Y is F,
\+ sys_poly_send(X, gen_ls, [Y]).
/**
* sys_poly_send(Y, F, T):
* The predicate succeeds for sending a message with funtor F
* and arguments T to the object Y.
*/
:- private sys_poly_send/3.
sys_poly_send(Y, F, T) :- number(Y), !,
M =.. [F,Y|T],
number:M.
sys_poly_send(Y, F, T) :-
M =.. [F|T],
Y::M.
package omonia.util.math;
/**
* <p>This class provides a double matrice. The matrice has not
* to be quadratic. You can use the constructor to construct
* a matrice of a given size with nulls in it. The methodes
* add(), mul(), neg() and inv() perform the usual arithmentic
* operations on matrices. The matrice supports elements which
* are double. A matrice with complex values is currently not
* yet supported.
* </p>
* <p>For the inversion the sweep algorithm is used as it is
* for example documented in the following booklet:
* </p>
* <p> Hoidn, H.P., Kirchgraber, U. and Marti, J.: Linear Algebra,
* <br> 3. Auflage, Verlag der Fachvereine Zürich, 1983
* </p>
* <p>This algorithm is cubic in its complexity. For special
* forms of matrices, such as positive definite ones or
* band matrices, better algorithms, such as cholesky
* decomposition respective scarce respresentations, are
* available. Nevertheless the given algorithm performs quite
* well in Java, for example a 100x100 matrice is inverted in
* about 50ms.</p>
*
* @version Omonia 1.3 (a simulation environment)
* @author Copyright 1997-2007 XLOG Technologies GmbH, Zürich Switzerland
*/
public class Matrice {
/* the number of columns of the matrice*/
public int cols;
/* the number of rows of the matrice */
public int rows;
/* the values of the matrice */
public double[][] data;
/**
* Construct a singleton matrice.
*/
public Matrice() {
this(1, 1);
}
/**
* Construct a matrice with the given number of columns and rows.
*
* @param m The number of columns.
* @param n The number of rows.
*/
public Matrice(int m, int n) {
cols = m;
rows = n;
data = new double[m][n];
}
/**
* Compute the sum of this matrice and the given matrice.
* This matrice and the given matrice must have the same
* dimension, otherwise an exception is thrown.
*
* @param arg The given matrice.
* @return The sum.
*/
public Matrice add(Matrice arg) {
if (cols != arg.cols || rows != arg.rows) {
throw new IllegalArgumentException("Matrices not of same size.");
}
Matrice temp = new Matrice(cols, rows);
for (int i = 0; i < cols; i++) {
for (int j = 0; j < rows; j++) {
temp.data[i][j] = data[i][j] + arg.data[i][j];
}
}
return temp;
}
/**
* Compute the product of this matrice and the given matrice.
* This matrice number of columns must correspond with the
* number of rows of the given matrice, otherwise an exception
* is thrown.
*
* @param arg The given matrice.
* @return The product.
*/
public Matrice mul(Matrice arg) {
if (cols != arg.rows) {
throw new IllegalArgumentException("Matrices cannot be aligned.");
}
Matrice temp = new Matrice(arg.cols, rows);
for (int i = 0; i < arg.cols; i++) {
for (int j = 0; j < rows; j++) {
double s = 0;
for (int k = 0; k < cols; k++) {
s += data[k][j] * arg.data[i][k];
}
temp.data[i][j] = s;
}
}
return temp;
}
/**
* Computes the negation of this matrice.
*
* @return The negation.
*/
public Matrice neg() {
Matrice temp = new Matrice(cols, rows);
for (int i = 0; i < cols; i++) {
for (int j = 0; j < rows; j++) {
temp.data[i][j] = -data[i][j];
}
}
return temp;
}
/**
* Make a copy of this matrice.
*
* @return A copy.
*/
public Matrice copy() {
Matrice temp = new Matrice(cols, rows);
for (int i = 0; i < cols; i++) {
//noinspection ManualArrayCopy
for (int j = 0; j < rows; j++) {
temp.data[i][j] = data[i][j];
}
}
return temp;
}
/**
* Computes the inversion of this matrice. If the matrice is not quadratic
* an exception is thrown. If during the inversion a sweep step cannot be
* performed any more an exception is thrown which indicates that the
* matrice is singular.
*
* @return The inversion.
*/
public Matrice inv() {
if (cols != rows)
throw new IllegalArgumentException("Matrice not quadratic.");
Matrice temp = copy();
int[] pivot = new int[cols];
for (int i = 0; i < cols; i++)
pivot[i] = -1;
for (int i = 0; i < cols; i++) {
int current;
for (current = 0; current < cols &&
(pivot[current] != -1 || temp.data[current][i] == 0); current++) ;
if (!(current < cols))
throw new ArithmeticException("Singular matrice " + i + "/" + cols);
pivot[current] = i;
temp.data[current][i] = 1 / temp.data[current][i];
for (int j = 0; j < cols; j++) {
if (j != current)
temp.data[j][i] = -temp.data[j][i] * temp.data[current][i];
}
for (int j = 0; j < cols; j++) {
for (int k = 0; k < cols; k++) {
if (j != current && k != i)
temp.data[j][k] += temp.data[j][i] * temp.data[current][k];
}
}
for (int j = 0; j < cols; j++) {
if (j != i)
temp.data[current][j] = temp.data[current][j] * temp.data[current][i];
}
}
for (int i = 0; i < cols; i++) {
int current = pivot[i];
if (current != i) {
for (int j = 0; j < cols; j++) {
double h = temp.data[j][i];
temp.data[j][i] = temp.data[j][current];
temp.data[j][current] = h;
}
int find;
for (find = i + 1; pivot[find] != i; find++) ;
pivot[find] = current;
double[] h = temp.data[i];
temp.data[i] = temp.data[current];
temp.data[current] = h;
}
}
return temp;
}
/**
* @param c Column
* @param r Row
* @param d Data
*/
public void setElement(int c, int r, double d) {
data[c][r] = d;
}
/**
* @param c Column
* @param r Row
* @return Data
*/
public double getElement(int c, int r) {
return data[c][r];
}
/**
*/
public String toString() {
String res = "[\n";
for (int r = 0; r < rows; r++) {
res += "[";
for (int c = 0; c < cols; c++) {
res += data[c][r] + " ";
}
res += "]\n";
}
return res + "]";
}
}
/**
* Vector/Matrice Operatios Prolog OO Prototype for Jekejeke Prolog.
*
* Warranty & Liability
* To the extent permitted by applicable law and unless explicitly
* otherwise agreed upon, XLOG Technologies GmbH makes no warranties
* regarding the provided information. XLOG Technologies GmbH assumes
* no liability that any problems might be solved with the information
* provided by XLOG Technologies GmbH.
*
* Rights & License
* All industrial property rights regarding the information - copyright
* and patent rights in particular - are the sole property of XLOG
* Technologies GmbH. If the company was not the originator of some
* excerpts, XLOG Technologies GmbH has at least obtained the right to
* reproduce, change and translate the information.
*
* Reproduction is restricted to the whole unaltered document. Reproduction
* of the information is only allowed for non-commercial uses. Selling,
* giving away or letting of the execution of the library is prohibited.
* The library can be distributed as part of your applications and libraries
* for execution provided this comment remains unchanged.
*
* Trademarks
* Jekejeke is a registered trademark of XLOG Technologies GmbH.
*/
:- package(library(diffuse/cas/matrice)).
:- use_package(library(jekpro/frequent/misc)).
:- module(matrice, []).
:- use_module(generic).
:- use_module(library(advanced/arith)).
/***********************************************************/
/* Array Builder & Access */
/***********************************************************/
/**
* X[Y, Z]:
* The predicate unifies Z with the Y-the vector
* of the matrice X.
*/
:- override []/3.
:- public []/3.
[](X, Y, Z) :- integer(Y), !,
arg(Y, X, Z).
/**
* X[Y, Z, T]:
* The predicate unifies T with the Z-the element
* of the Y-the vector of the matrice X.
*/
:- override []/4.
:- public []/4.
[](X, Y, Z, T) :- integer(Y), integer(Z), !,
arg(Y, X, H),
arg(Z, H, T).
/**
* len(X, Y):
* The predicate unifies Y with the number of vectors
* in the matrice X.
*/
:- public len/2.
len(X, Y) :-
functor(X, _, Y).
/***********************************************************/
/* Basic Arithmetic */
/***********************************************************/
/**
* -(X, Y):
* The predicate unifies Y with the sign changed matrice X.
*/
:- override (-)/2.
:- public (-)/2.
-(X, Y) :-
L is len(X), Y is {-(X[I]) | between(1,L,I) }.
/**
* +(X, Y, Z):
* The predicate unifies Z with the sum of the matrice X and
* the matrice Y.
*/
:- override (+)/3.
:- public (+)/3.
+(X, Y, Z) :- functor(Y, matrice, _), L is len(X), L =:= len(Y), !,
Z is { X[I]+Y[I] | between(1,L,I) }.
/**
* -(X, Y, Z):
* The predicate unifizes Z with the the matrice X subtracted
* by the matrice Y.
*/
:- override (-)/3.
:- public (-)/3.
-(X, Y, Z) :- functor(Y, matrice, _), L is len(X), L =:= len(Y), !,
Z is { X[I]-Y[I] | between(1,L,I) }.
/**
* *(X, Y, Z):
* The predicate unifies Z with the product of the matrice X followed
* by tbe matrice Y.
*/
:- override (*)/3.
:- public (*)/3.
*(X, Y, Z) :- functor(Y, matrice, _), L is len(X[1]), L =:= len(Y), !,
M is len(X), N is len(Y[1]),
Z is { { sum({ X[J,K] * Y[K,I] | between(1,L,K) })
| between(1,N,I) }
| between(1,M,J) }.
/**
* ^(X, Y, Z):
* The predicate unifies Z with the Y-the power of the matrice X.
*/
:- override (^)/3.
:- public (^)/3.
^(X, Y, R) :- integer(Y), user: (Y =:= 0), L is len(X), L =:= len(X[1]), !,
R is { { V | between(1,L,I), (user: (I =:= J) -> V = 1; V = 0) }
| between(1,L,J) }.
^(X, Y, R) :- integer(Y), user: (Y =:= 1), !,
R = X.
^(X, Y, R) :- integer(Y), user: (Y > 1), user: (Y rem 2 =:= 1), !,
user: -(Y, 1, Z),
R is X*X^Z.
^(X, Y, R) :- integer(Y), user: (Y > 1),
user: div(Y, 2, Z),
H is X^Z,
R is H*H.
/***********************************************************/
/* Basic Comparison */
/***********************************************************/
:- public gen_eq/2.
gen_eq(X, Y) :- functor(Y, matrice, _), L is len(X), L =:= len(Y), !,
\+ (between(1,L,I), X[I]=\=Y[I]).
/***********************************************************/
/* CAS Display Hook */
/***********************************************************/
:- public residue:sys_portray_eq/2.
:- multifile residue:sys_portray_eq/2.
:- meta_predicate residue:sys_portray_eq(0, 0).
residue:sys_portray_eq(_ = X, _) :- var(X), !, fail.
residue:sys_portray_eq(X = F, X is R) :- functor(F, matrice, _), !,
F =.. [_|L],
sys_portray_list(L, R).
:- private sys_portray_list/2.
sys_portray_list([F|L], [G|R]) :-
F =.. [_|G],
sys_portray_list(L, R).
sys_portray_list([], []).
/**
* Vector/Matrice Operatios Prolog OO Prototype for Jekejeke Prolog.
*
* Warranty & Liability
* To the extent permitted by applicable law and unless explicitly
* otherwise agreed upon, XLOG Technologies GmbH makes no warranties
* regarding the provided information. XLOG Technologies GmbH assumes
* no liability that any problems might be solved with the information
* provided by XLOG Technologies GmbH.
*
* Rights & License
* All industrial property rights regarding the information - copyright
* and patent rights in particular - are the sole property of XLOG
* Technologies GmbH. If the company was not the originator of some
* excerpts, XLOG Technologies GmbH has at least obtained the right to
* reproduce, change and translate the information.
*
* Reproduction is restricted to the whole unaltered document. Reproduction
* of the information is only allowed for non-commercial uses. Selling,
* giving away or letting of the execution of the library is prohibited.
* The library can be distributed as part of your applications and libraries
* for execution provided this comment remains unchanged.
*
* Trademarks
* Jekejeke is a registered trademark of XLOG Technologies GmbH.
*/
:- package(library(diffuse/cas/matrice)).
:- use_package(library(jekpro/frequent/misc)).
:- module(number, []).
:- use_module(generic).
:- use_module(library(advanced/arith)).
/***********************************************************/
/* Array Builder & Access */
/***********************************************************/
/**
* .(X, Y, Z):
* The predicate unifies in Z the consing of X and Y.
*/
:- public '.'/3.
'.'(X, Y, Z) :- functor(Y, vector, _), !,
Y =.. [_|L],
Z =.. [vector,X|L].
/***********************************************************/
/* Basic Arithmetic */
/***********************************************************/
:- public bitlength/2.
bitlength(X, Y) :-
misc/bits: bitlength(X, Y).
/***********************************************************/
/* Basic Comparison */
/***********************************************************/
:- public gen_eq/2.
gen_eq(X, Y) :- number(Y), !,
user: =:=(X, Y).
:- public gen_ls/2.
gen_ls(X, Y) :- number(Y), !,
user: <(X, Y).
/***********************************************************/
/* CAS Display Hook */
/***********************************************************/
:- public residue:sys_portray_eq/2.
:- multifile residue:sys_portray_eq/2.
:- meta_predicate residue:sys_portray_eq(0, 0).
residue:sys_portray_eq(_ = X, _) :- number(X), !, fail.
/**
* Vector/Matrice Operatios Prolog OO Prototype for Jekejeke Prolog.
*
* Warranty & Liability
* To the extent permitted by applicable law and unless explicitly
* otherwise agreed upon, XLOG Technologies GmbH makes no warranties
* regarding the provided information. XLOG Technologies GmbH assumes
* no liability that any problems might be solved with the information
* provided by XLOG Technologies GmbH.
*
* Rights & License
* All industrial property rights regarding the information - copyright
* and patent rights in particular - are the sole property of XLOG
* Technologies GmbH. If the company was not the originator of some
* excerpts, XLOG Technologies GmbH has at least obtained the right to
* reproduce, change and translate the information.
*
* Reproduction is restricted to the whole unaltered document. Reproduction
* of the information is only allowed for non-commercial uses. Selling,
* giving away or letting of the execution of the library is prohibited.
* The library can be distributed as part of your applications and libraries
* for execution provided this comment remains unchanged.
*
* Trademarks
* Jekejeke is a registered trademark of XLOG Technologies GmbH.
*/
:- use_module(generic).
/**
* pell(N, R, S):
* Compute the N-the solution to p*p+1=2*q*q via matrices and return
* p/q in decimal format in R. An upper bound 1/(2*p*p) for the relative
* error is returned as logarithm base 10 in S.
*/
pell(N, R, S) :-
Z is [[3,4],[2,3]]^N*[[1],[1]],
S is ceiling((-2*bitlength(Z[1,1])-1)*log(2)/log(10)),
D is 10^(-S),
H is Z[1,1]*D//Z[2,1],
J is 2*D*D,
nudge(H, J, R).
/**
* nudge(H, J, R):
* Do a last correction an return deciomal format in R. H is the possibly
* too low result and J is the value 2 in decimal format.
*/
nudge(H, J, R) :-
K is H+1, J > K*K, !,
nudge(K, J, R).
nudge(H, _, H).
% ?- use_module(library(advanced/arith)).
% ?- between(1, 23, N), pell(N, R, S), write('R='), write(R), write(', S='), write(S), nl, fail; true.
% R=141, S=-2
% R=1414, S=-3
% R=141421, S=-5
% R=1414213, S=-6
% R=141421356, S=-8
% R=1414213562, S=-9
% R=141421356237, S=-11
% R=1414213562373, S=-12
% R=141421356237309, S=-14
% R=1414213562373095, S=-15
% R=141421356237309504, S=-17
% R=1414213562373095048, S=-18
% R=141421356237309504880, S=-20
% R=1414213562373095048801, S=-21
% R=141421356237309504880168, S=-23
% R=1414213562373095048801688, S=-24
% R=141421356237309504880168872, S=-26
% R=14142135623730950488016887242, S=-28
% R=141421356237309504880168872420, S=-29
% R=14142135623730950488016887242096, S=-31
% R=141421356237309504880168872420969, S=-32
% R=14142135623730950488016887242096980, S=-34
% R=141421356237309504880168872420969807, S=-35
% ?- time(pell(653124, _, S)).
% % Uptime 1,646 ms, GC Time 30 ms, Thread Cpu Time 1,640 ms
% S = -1000000
/**
* Vector/Matrice Operatios Prolog OO Prototype for Jekejeke Prolog.
*
* Warranty & Liability
* To the extent permitted by applicable law and unless explicitly
* otherwise agreed upon, XLOG Technologies GmbH makes no warranties
* regarding the provided information. XLOG Technologies GmbH assumes
* no liability that any problems might be solved with the information
* provided by XLOG Technologies GmbH.
*
* Rights & License
* All industrial property rights regarding the information - copyright
* and patent rights in particular - are the sole property of XLOG
* Technologies GmbH. If the company was not the originator of some
* excerpts, XLOG Technologies GmbH has at least obtained the right to
* reproduce, change and translate the information.
*
* Reproduction is restricted to the whole unaltered document. Reproduction
* of the information is only allowed for non-commercial uses. Selling,
* giving away or letting of the execution of the library is prohibited.
* The library can be distributed as part of your applications and libraries
* for execution provided this comment remains unchanged.
*
* Trademarks
* Jekejeke is a registered trademark of XLOG Technologies GmbH.
*/
:- package(library(diffuse/cas/matrice)).
:- use_package(library(jekpro/frequent/misc)).
:- module(vector, []).
:- use_module(generic).
:- use_module(library(advanced/arith)).
/***********************************************************/
/* Array Builder & Access */
/***********************************************************/
/**
* .(X, Y, Z):
* The predicate unifies in Z the consing of X and Y. If
* Y is a vector it is promoted to a matrice.
*/
:- public '.'/3.
'.'(X, Y, Z) :- functor(Y, matrice, _), !,
Y =.. [_|L],
Z =.. [matrice,X|L].
'.'(X, Y, Z) :- functor(Y, vector, _),
Y =.. [_|L],
sys_promote_vector(L, R),
Z =.. [matrice,X|R].
:- private sys_promote_vector/2.
sys_promote_vector([X|L], [vector(X)|R]) :-
sys_promote_vector(L, R).
sys_promote_vector([], []).
/**
* X[Y, Z]:
* The predicate unifies Z with the Y-the cell
* of the vector X.
*/
:- override []/3.
:- public []/3.
[](X, Y, Z) :- integer(Y),
arg(Y, X, Z).
/**
* len(X, Y):
* The predicate unifies Y with the number of cells
* in the vector X.
*/
:- public len/2.
len(X, Y) :-
functor(X, _, Y).
/**
* sum(X, Y):
* The predicate unifies Y with the sum of cells
* in the vector X.
*/
:- public sum/2.
sum(X, Y) :-
X =.. [_|L], sys_sum_list(L, Y).
:- private sys_sum_list/2.
sys_sum_list([X|L], R) :-
sys_sum_list(L, H), R is X+H.
sys_sum_list([], 0).
/***********************************************************/
/* Basic Arithmetic */
/***********************************************************/
/**
* -(X, Y):
* The predicate unifies Y with the sign changed vector X.
*/
:- override (-)/2.
:- public (-)/2.
-(X, Y) :-
L is len(X), Y is { -(X[I]) | between(1,L,I) }.
/**
* +(X, Y, Z):
* The predicate unifies Z with the sum of the vector X and
* the vector Y.
*/
:- override (+)/3.
:- public (+)/3.
+(X, Y, Z) :- functor(Y, vector, _), L is len(X), L =:= len(Y), !,
Z is { X[I]+Y[I] | between(1,L,I) }.
/**
* -(X, Y, Z):
* The predicate unifizes Z with the the vector X subtracted
* by the vector Y.
*/
:- override (-)/3.
:- public (-)/3.
-(X, Y, Z) :- functor(Y, vector, _), L is len(X), L =:= len(Y), !,
Z is { X[I]-Y[I] | between(1,L,I) }.
/***********************************************************/
/* Basic Comparison */
/***********************************************************/
:- public gen_eq/2.
gen_eq(X, Y) :- functor(Y, vector, _), L is len(X), L =:= len(Y), !,
\+ (between(1,L,I), X[I] =\= Y[I]).
/***********************************************************/
/* CAS Display Hook */
/***********************************************************/
:- public residue:sys_portray_eq/2.
:- multifile residue:sys_portray_eq/2.
:- meta_predicate residue:sys_portray_eq(0, 0).
residue:sys_portray_eq(_ = X, _) :- var(X), !, fail.
residue:sys_portray_eq(X = F, X is G) :- functor(F, vector, _), !,
F =.. [_|G].
@jburse

This comment has been minimized.

Copy link
Owner Author

@jburse jburse commented Dec 5, 2016

There are a lot of things to do:

  • Range checks.
  • Type checks.
  • Matrix inversion.
  • Integrate with rationals, symbolic expressions or complex.
  • Micro engines.
  • Parallelization.
@jburse

This comment has been minimized.

Copy link
Owner Author

@jburse jburse commented Dec 5, 2016

We just added the positive integer power function over matrices and did a little fun
experiment in calculating a million digits of sqrt(2) in less than a minute:

?- time(pell(10000000,_)).
% Uptime 49,641 ms, GC Time 1,329 ms, Thread Cpu Time 48,203 ms
Yes 

Here are some few of the first digits. Compare the below with
http://apod.nasa.gov/htmltest/gifcity/sqrt2.10mil:

?- time(pell(10000000,R)).
% Uptime 48,469 ms, GC Time 1,233 ms, Thread Cpu Time 47,110 ms
R = 14142135623730950488016887242096980785696718753769480731766797379907324784621
07038850387534327641572735013846230912297024924836055850737212644121497099935831
41322266592750559275579995050115278206057147010955997160597027453459686201472851
74186408891986095523292304843087143214508397626036279952514079896872533965463318
08829640620615258352395054745750287759961729835575220337531857011354374603408498
84716038689997069900481503054402779031645424782306849293691862158057846311159666
87130130156185689872372352885092648612494977154218334204285686060146824720771435
85487415565706967765372022648544701585880162075847492265722600208558446652145839
88939443709265918003113882464681570826301005948587040031864803421948972782906410
45072636881313739855256117322040245091227700226941127573627280495738108967504018
36986836845072579936472906076299694138047565482372899718032680247442062926912485
90521810044598421505911202494413417285314781058036033710773091828693147101711116
83916581726889419758716582152128229518488472089694633862891562882765952635140542
26765323969461751129160240871551013515045538128756005263146801712740265396947024
03005174953188629256313851881634780015693691768818523786840522878376293892143006
55869568685964595155501644724509836896036887323114389415576651040883914292338113
20605243362948531704991577175622854974143899918802176243096520656421182731672625
75395947172559346372386322614827426222086711558395999265211762526989175409881593
48640083457085181472231814204070426509056532333398436457865796796519267292399875
36661721598257886026336361782749599421940377775368142621773879919455139723127406
68983299898953867288228563786977496625199665835257761989393228453447356947949629
52168891485492538904755828834526096524096542889394538646625744927556381964410316
97983306185201937938494005715633372054806854057586799967012137223947582142630658
@jburse

This comment has been minimized.

Copy link
Owner Author

@jburse jburse commented Dec 20, 2016

We improved the code in that we estimate an upper bound for a relative
error, and further compute stable digits by nudging the result.

This is the result when using only the relative error bound:

?- between(1, 23, N), pell(N, R, S), write('R='), write(R), write(', S='), write(S), nl,  fail; true.
R=140, S=-2
R=1413, S=-3
R=141420, S=-5
R=1414213, S=-6
R=141421355, S=-8
R=1414213562, S=-9
R=141421356236, S=-11
R=1414213562372, S=-12
R=141421356237308, S=-14
R=1414213562373094, S=-15
R=141421356237309504, S=-17
R=1414213562373095048, S=-18
R=141421356237309504879, S=-20
R=1414213562373095048801, S=-21
R=141421356237309504880168, S=-23
R=1414213562373095048801688, S=-24
R=141421356237309504880168871, S=-26
R=14142135623730950488016887240, S=-28
R=141421356237309504880168872420, S=-29
R=14142135623730950488016887242095, S=-31
R=141421356237309504880168872420969, S=-32
R=14142135623730950488016887242096979, S=-34
R=141421356237309504880168872420969807, S=-35

This is the result when nudging is also done:

?- between(1, 23, N), pell(N, R, S), write('R='), write(R), write(', S='), write(S), nl,  fail; true.
R=141, S=-2
R=1414, S=-3
R=141421, S=-5
R=1414213, S=-6
R=141421356, S=-8
R=1414213562, S=-9
R=141421356237, S=-11
R=1414213562373, S=-12
R=141421356237309, S=-14
R=1414213562373095, S=-15
R=141421356237309504, S=-17
R=1414213562373095048, S=-18
R=141421356237309504880, S=-20
R=1414213562373095048801, S=-21
R=141421356237309504880168, S=-23
R=1414213562373095048801688, S=-24
R=141421356237309504880168872, S=-26
R=14142135623730950488016887242, S=-28
R=141421356237309504880168872420, S=-29
R=14142135623730950488016887242096, S=-31
R=141421356237309504880168872420969, S=-32
R=14142135623730950488016887242096980, S=-34
R=141421356237309504880168872420969807, S=-35

And we can bring down the time to compute a million digits to less than 2 secs:

% Jekejeke Runtime 1.1.7
?- time(pell(653124, _, S)).
% Uptime 1,646 ms, GC Time 30 ms, Thread Cpu Time 1,640 ms
S = -1000000
@jburse

This comment has been minimized.

Copy link
Owner Author

@jburse jburse commented Mar 26, 2017

Meanwhile we have forked our matrice development. There is now a second branch based on our new polynomial GCD which also supports symbolic matrice inversion. You find the sources of this new development here:

Symbolic Matrices with Inversion
https://gist.github.com/jburse/3a14fbae8d7e14b238407556d7835e7b

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.