Attribute Variables and Computer Algebra? (Gröbner Basis)

[burs...@gmail.com]

What do have attribute variables and computer algebra
in common? First it seems not so much. But then observe
a nice property of attribute variables.

They usually buoy inside a term, at least many implementations
do that. They prefer a non-attribute variable to be instantiated
with a attribute variable and not vice versa.

Since instantiating an attribute variable fires the attribute
variable hooks. So when we can delay this the better. But here
is another application for computer algebra.

So we first started a little polynom manipulation library,
which allows doing polynom arithmetic +, -, * as follows,
just using Prolog variables as the symbolic variables:

?- S is (X-1)*(X+1).
S is X^2-1

Under the hood something more is going on. We wrap attribute
variables into object references so that we can write Prolog
rules that see these variables as atoms.

Interesting applications of this approach are for example
use cases where variables need to be freshly generated. We
can do what every Prolog programmer always does, just mention
a fresh variable, thats all:

sys_poly_lcm(A, B, C) :-
S is A*Z,
T is B*(1-Z),
sys_poly_groeb([S,T],[(0,1)],[C|_]).

In the above code Z is a fresh variable. We compute a polynomial
LCM via a Gröbner Basis trick mentioned in a Python documentation.
As a result we can normalize polynomial fractions, even
multivariate ones:

?- A is (X^2*Y-X)^2, B is (X*Y^2-Y)^2, R is A/B.
A is X^4*Y^2-2*X^3*Y+X^2,
B is X^2*Y^4-2*X*Y^3+Y^2,
R is X^2/Y^2

For more information:

Preview: Gröbner Basis GCD via OO Prolog. (Jekejeke)
https://plus.google.com/+JekejekeCh/posts/aiGNsA8QY9J

Gröbner bases and their applications - Mateusz, P. (2010)
https://mattpap.github.io/masters-thesis/html/src/groebner.html

[burs...@gmail.com]

Dear All,

It turned out that our realization of the polynomial GCD
had a little glitch. The better code for the LCM that is
used in the computation of the GCD is as follows:

sys_poly_lcm(A, B, C) :-
S is A*Z,
T is B*(1-Z),

sys_poly_groeb([S,T],L),
sys_poly_min(L, C).

The fix is necessary since sys_poly_groeb/2 doesn't return
the Gröbner basis as a sorted list, sorted by the lexical
ordering that puts Z at last. With this little fix some
border cases such as this one now work correctly:

?- S is 1, T is 2/3*A+4/5*B, R is S/T.
S = 1,
T is 4/5*B+2/3*A,
R is 15/(12*B+10*A)

The algorithm itself works with rational monomial coefficients,
but there is also finishing step to make the numerator and
denumerator of the polynomial fraction with integer monomial
coefficients. This is all seen in the open source gist.

Now that we have polynomial GCD, what can we do with it?
Well we just made a little experiment in providing symbolic
matrices in Prolog. We made such an experiment already in
december 2016, but we added matrice inversion now.

We implemented an inline matrice inversion algorithm. Its
name is "Austauschschritt Verfahren" or AT-Algorithm. It
is in fact realized as Prolog-ish tail recursive, but our
implementation is derived from an originally inline approach.

Hilbert matrices are numerically difficult,
but we can do them now exact:

?- hilbert_matrice(4, X), Y is X^(-1).
X is [[1,1/2,1/3,1/4],[1/2,1/3,1/4,1/5],[1/3,1/4,1/5,1/6],[1/4,1/5,1/6,1/7]],
Y is [[16,-120,240,-140],[-120,1200,-2700,1680],[240,-2700,6480,-4200],[-140,1680,-4200,2800]]

The problem matrice can also have
fully symbolic entries:

?- X is [[1,B],[B^2,B]], Y is X^(-1).
X is [[1,B],[B^2,B]],
Y is [[-1/(B^2-1),1/(B^2-1)],[B/(B^2-1),-1/(B^3-B)]]

The AT-Algorithm uses less steps than for example the
Cramer method where we would compute a couple of
determinants. At the moment we only implemented a version
without pivoting.

Upon code inspection one will see that we used the same
programming approach for matrice inversion as in our
december 2016 matrice multiplication prototype, i.e
findall/3 and arg/3 that can be used inside is/2.

For more information:

Preview: Symbolic Matrice Inversion in OO Prolog. (Jekejeke)
https://plus.google.com/+JekejekeCh/posts/8jRmzruDMgs

Code Golf: Symbolic Matrix Multiplication in Prolog
http://codegolf.stackexchange.com/questions/106570/symbolic-matrix-multiplication/113875#113875

[kint...@gmail.com]

> ?- S is (X-1)*(X+1).
> S is X^2-1

~~~
That code does not work for me .
Are you using some kind of equation rewriter ?

First I tried :

?- S is (X-1)*(X+1) .
ERROR: is/2: Arguments are not sufficiently instantiated
?-

Then I tried , giving X a value , still no success :

?- S is (X-1)*(X+1) , X = 7.
?-

Then the little birdies started chattering in my mind again and I am off my meds so I tried this :

?- use_module(library(clpfd)) .
true.

?-
(
MACRO_INSTALL_EXPR = system:asserta(goal_expansion(OLD_EXPR_0,NEW_EXPR_9)) ,
OLD_EXPR_0=(S is WHATEVER) ,
NEW_EXPR_9=(S #= WHATEVER) ,
system:call(MACRO_INSTALL_EXPR)
)
.
true.

With that in place , the test providing X is successful :

?- S is (X-1)*(X+1) , X=7 .
S = 48,
X = 7.

The test not providing X now "succeeds" (.i.e.|it does not return an (implicit) false|.i.e.),
but now there seem to be a bunch of N0NS3N5E for the answer :

?- S is (X-1)*(X+1) .
_4072*_4074#=S,
_4072+1#=X,
X+1#=_4074.

?-

I tried a test where I provided an S instead of an X .
I expected the answer to be '7' , a simple reversal of the X=7 test above .
I have noticed in the past that in some circumstances prolog can do this . You can make a variable on the left hand-side of your equation have a real value , and then leave a variable open on the right-side of your equation , and prolog rewrites the equation and puts the open variable on the left-hand side so it can receive the value .

Instead of an answer I got this weird junk :

? S is (X-1)*(X+1) , S=48 .
S = 48,
X in -7.. -2\/0\/2..7,
X+1#=_6992,
_7020+1#=X,
_6992 in -6.. -1\/1\/3..8,
_7020*_6992#=48,
_7020 in -8.. -3\/ -1\/1..6.

I started to doubt myself about that equation rewriting thing ,
so I fired up a new instance of swipl and tried this again :

S is (X-1)*(X+1) , S=48 .

ERROR: is/2: Arguments are not sufficiently instantiated

I guess I was wrong about equation rewritting .,
Then I realized something obvious .
Now I feel like an idiot .
Of course there is no sensible answer the computer can provide for this question !

When the computer is asked :

?- S is (X-1)*(X+1) .

- There are 2 unknowns , S and X .
- The computer can't calculate S unless an X value is known .
- The computer can't calculate X unless an S value is known .
- **Therefore** either an S or an X must be provided or the computer cannot calculate the answer .

The little birdies always tell me "use clpfd , use clpfd" . They say lots of weird stuff like that . At this point I am starting to wonder if they are deliberately trying to make me insane .

Sorry if this is off-topic from your post but it seems like it mnight somehow be related .

~~kintalken~~
~~~
~~~

[kint...@gmail.com]

call
====
---

> You can make a variable on the left hand-side of your equation have a real
> value , and then leave a variable open on the right-side of your equation ,
> and prolog rewrites the equation and puts the open variable on the left-hand
> side so it can receive the value .

---
---

response
========
---

It seems that you are an imperativ programmer, which is not surprising because 99.99999% of programmers are .
This is betrayed by the way your answer contains action-oriented commands , like ">>>put<< the variable" and ">>>receive<< the value .

Prolog is a meta-meta-logical programming language .

It is a declarativ approach to programming that describes relationships ; in constrast to the imperativ approach that describes operations to perform .

In fact Prolog completely eschews the operaational approach , and as you will learn if you try to learn Prolog it is a completely baffling task to get Prolog to >>>do<<< anything at all . This is because >>do<< is an action word , and ideally there is no action , everything is a noun .

You might have more success posting your question elsewhere .
Nobody around here is interested in communicating to you because we believe that your approach is wrong .

---
astimony
---
---

[burs...@gmail.com]

Its a preview of, there is still some testing and improvement done:
https://gist.github.com/jburse/3a14fbae8d7e14b238407556d7835e7b

I will probably bundle it with the next release 1.2.1 of Jekejeke Prolog.
Its not some CLP(*) at the moment, not sure whether it will move in this direction.

[burs...@gmail.com]

Hint: Open the Gist, scroll down till you find the comment
section, there are more links, since the code is distributed
over two Gists (one Gist for matrices, one Gist for polynomials).

[burs...@gmail.com]

As a next step maybe will do something in this direction:

Symbolic QR decomposition via Isabelle/HOL - Aransay & Divasónn 2016
https://t.co/24uNNDS2Et

[burs...@gmail.com]

There was an interesting experiment by Samer Abdallah
on SWI-Prolog google groups with CHR and derive.

Which is not only degrading attribute variables to atoms,
as I do here, as a convenience for sympy symbol('x'),

but using attribute variables with instantiation for
evaluation. But I have my doubts that this works

satisfactory, a use case being taylor expansion. But
one never knows, maybe it works? See here:
https://groups.google.com/d/msg/swi-prolog/PONN0X1HD_w/StspvbMkBAAJ

Am Samstag, 18. März 2017 21:55:27 UTC+1 schrieb burs...@gmail.com:

[burs...@gmail.com]

I shouldn't have deleted the below citation. Just
stepped over a new use case, for generating Prolog
variables on the fly in some CAS application.

Basically I couldn't resist, implemented a new expression
deriv(P, X) where X is a binder, and a new expression
subst(P, X, R) where again X is a binder. There
was already a prototype for that in december 2016.

That one of the arguments is a binder doesn't give much
trouble for CAS, just maybe check that the argument has
the type of a variable, in our case the same type as
the here discussed wrapped attribute variables.

Putting deriv and subst together, we could
taylor expand:

?- X is 1/(1+A), Y is taylor(X, A, 5).
X is 1/(1+A),
Y is 1-A+A^2-A^3+A^4-A^5

The surprise came when we tackled laurent series. The
mathematical presentations are usually quite contrived
and don't cut through the thicket. The code reads
as simple as, the fresh variable is Y:

:- public element:laurent/4.
element:laurent(P, X, N, R) :-
H is subst(P, X, 1/Y),
J is taylor(H, Y, N),
R is subst(J, Y, 1/X).

So developing a laurent series for the same
example as befor yields:

?- X is 1/(1+A), Y is laurent(X, A, 5).
X is 1/(1+A),
Y is (1-A+A^2-A^3+A^4)/A^5

The prototype also does not yet support special functions
and it does not yet evaluate the given function with the
help of limes. Finally a multi-variable version of the
expansion is not yet supported.

See also: Symbolic Series Development in OO Prolog.
https://plus.google.com/+JekejekeCh/posts/Y5QMvvW5emP

Jan Burse wrote:
> Interesting applications of this approach are for example
> use cases where variables need to be freshly generated. We
> can do what every Prolog programmer always does, just mention

> a fresh variable, thats all, the fresh variable is Z:

[burs...@gmail.com]

This link has now a proof of concept that CHR
can be also used for taylor expansion:
https://groups.google.com/d/msg/swi-prolog/PONN0X1HD_w/zYvJGrksAAAJ

All done by Samer Abdallah with a neat little example,
taylor expansion of sqrt(1+X):

?- add(1,X,X1), pow(-1 rdiv 2,X1,Y), time(taylor(50,X,Y,Cs)), X=0.
% 971,174 inferences, 0.196 CPU in 0.198 seconds (99% CPU, 4966118 Lips)
X = 0,
X1 = Y, Y = 1,
Cs = [1, -1 rdiv 2, 3 rdiv 8, -5 rdiv 16, 35 rdiv 128, -63 rdiv 256, 231 rdiv 1024, -429 rdiv 2048, ...rdiv...|…].

Am Samstag, 15. April 2017 16:14:56 UTC+2 schrieb burs...@gmail.com:

[burs...@gmail.com]

Corr.
The taylor expansion is for f(x)=1/sqrt(1+x).

[burs...@gmail.com]

Dear All,

We made a little addition to our Gröbner Basis
algorithm. We now also support a number datatype that
has the form:

a = p + s1*sqrt(p1) + .. + sn*sqrt(pn)

where p, p1, .., pn are rational numbers and where
s1, .., sn are signs in the set {-1,1}. We managed to
provide fully symbolic basic arithmetic, including
reciprocal:

?- use_module(library(groebner/generic)).
% 21 consults and 0 unloads in 719 ms.
Yes

?- X is (1+sqrt(2))^3.
X is 7+sqrt(50)

?- X is (1+sqrt(2))^(-1).
X is - 1+sqrt(2)

Interestingly this was all that was needed to plug
in the new datatype into the Gröbner basis algorithm,
so that now the following things work:

?- X is 1-10*A^2+A^4, Y is -5-sqrt(24)-(4+sqrt(24))*A^2+A^4, Z is X/Y.
X is 1-10*A^2+A^4,
Y is - 5-sqrt(24)-(4+sqrt(24))*A^2+A^4,
Z is 1-(6-sqrt(24))/(1+A^2)

More recently we also considered exact symbolic ordering
the radical datatype and computing an integer floor for
the radical datatype. Both with success:

?- X is sqrt(11)-sqrt(5)-sqrt(2),
(X < 0 -> write('Jau'), nl; write('Nö'), nl).
Jau
X is -sqrt(2)-sqrt(5)+sqrt(11)

?- X is floor((sqrt(77/3)+sqrt(3/77))*10^100).
X is 526361355967815190707142656486221473163024293045
96489408819263671700931193265582373375139798275704672

The basic arithmetic makes use of an isqrt function to
determine linear dependence of sqrt summands. The
reciprocal and the ordering make use of Swinnerton-Dyer
Polynomials.

Best Regards

See also:

Preview: Floor for Symbolic Square Root Sums. (Jekejeke)
https://plus.google.com/+JekejekeCh/posts/ghQhMmDR9fA

Efficient Computation of Swinnerton-Dyer Polynomials
https://math.stackexchange.com/questions/1121954/efficient-computation-of-swinnerton-dyer-polynomials

Integer Square Root
https://en.wikipedia.org/wiki/Integer_square_root

Am Samstag, 18. März 2017 21:55:27 UTC+1 schrieb burs...@gmail.com:

[burs...@gmail.com]

Here is a nice case where CLP(FD) cannot help,
but CAS does indeed help here. Especially when
it has a radical datatype:

Without CAS, wrong output:
==========================

?- eval_keys([(sqrt(98428513)+sqrt(101596577))-foo,sqrt(400025090)-bar],L).
L = [20000.62724016424-foo, 20000.627240164245-bar].

5 ?- min_key([20000.62724016424-foo, 20000.627240164245-bar], X).
X = foo.

With CAS, correct output:
=========================

?- eval_keys([(sqrt(98428513)+sqrt(101596577))-foo,sqrt(400025090)-bar],L).
L = [radical(0,[98428513-1,101596577-1])-foo,
radical(0,[400025090-1])-bar]

?- min_key([radical(0,[98428513-1,101596577-1])-foo,
radical(0,[400025090-1])-bar],X).
X = bar

Neither the standard predicate is/2 nor a CLP(FD)
predicate (#=)/2 can do mathematics here. So that in
the end, for certain applications such as exact geometry,
they might not be suitable.

On the other hand many CAS systems support radicals.
Our Prototype even supports comparison of them so that
we can obtain the exact result. That bar is the exact
result can be seen for example by using a multi-precision
calculator.

?- use_module(library(decimal/multi)).
% 7 consults and 0 unloads in 319 ms.
Yes

?- X is mp(sqrt(98428513)+sqrt(101596577), 32).
X = 0d20000.627240164244658958331341095

?- X is mp(sqrt(400025090), 32).
X = 0d20000.627240164244408966171597261

But a CAS need not to proceed this way. For example
our Prolog implementation uses a Swinnerton-Dyer
polynomial inspired method to compare radical
expressions, which works purely symbolic.