[Maxima-discuss] ordering in a product

[Barton Willis]

What is the term ordering in a product? Consider

(%i88) (2*x+1)^2 * (x^2-1)$

(%i89) args(%);
(%o89) [(2*x+1)^2,x^2-1]

(%i90) sort(%, orderlessp);
(%o90) [x^2-1,(2*x+1)^2]

(%i91) a*b*c;
(%o91) a*b*c

(%i92) args(%);
(%o92) [a,b,c]

(%i93) sort(args(%), orderlessp);
(%o93) [a,b,c]

For a*b*c, the term order is the same as orderlessp, for (2*x+1)^2 * (x^2-1) it isn't? Or am I being deceived?

As I understand the simptimes code, the main routine is timesin. The function timesin does its own sorting without calling sort---that makes it challenging to know how terms are sorted.

Arguably understanding the term ordering in a product isn't particularly important. A robust algorithm should be somewhat immune to changes in the way expressions are simplified--generally I would say that very few algorithms should rely on a particular term ordering for sum or products. But if timesin isn't sorting consistently, that could lead to bug--I haven't found any.


--Barton

[Robert Dodier]

On 2017-07-21, Barton Willis <***@unk.edu> wrote:

> (%i88) (2*x+1)^2 * (x^2-1)$
>
> (%i89) args(%);
> (%o89) [(2*x+1)^2,x^2-1]
>
> (%i90) sort(%, orderlessp);
> (%o90) [x^2-1,(2*x+1)^2]
>
> (%i91) a*b*c;
> (%o91) a*b*c
>
> (%i92) args(%);
> (%o92) [a,b,c]
>
> (%i93) sort(args(%), orderlessp);
> (%o93) [a,b,c]
>
> For a*b*c, the term order is the same as orderlessp, for (2*x+1)^2 * (x^2-1) it isn't? Or am I being deceived?

TIMESIN calls GREAT to compare terms, the same as 'sort'. But in this
case, TIMESIN compares only 2*x + 1, not (2*x + 1)^2.

(%i18) (2*x - 1)*(x + 1)^2;
1 Enter great [x+1,2*x-1]
2 Enter great [x,2*x]
3 Enter great [1,2]
3 Exit great false
2 Exit great false
1 Exit great false
(%o18) (x+1)^2*(2*x-1)

But calling sort compares (2*x + 1)^2.

(%i19) sort(args(%));
1 Enter great [(x+1)^2,2*x-1]
2 Enter great [(x+1)^2,2*x]
3 Enter great [(x+1)^2,x]
4 Enter great [x+1,x]
5 Enter great [0,1]
5 Exit great false
4 Exit great true
3 Exit great true
2 Exit great true
1 Exit great true
(%o19) [2*x-1,(x+1)^2]

Now I'm aware that GREAT looks inside of nonatomic expressions when it's
figuring out how to sort them -- e.g. a(x)*b(w)*c(z) --> b(w)*a(x)*c(z).
But in this case, we have GREAT(f(e1), g(e2)) but not GREAT(e1, e2). I
don't know how common such cases are. We could try to discover other
cases by e.g. fishing "*" (or "+" for good measure) expressions out of
the test suite and determining whether args(e) = sort(args(e)).

I see that COMMUTATIVE1 in src/asum.lisp (invoked after declare(foo,
commutative)) calls GREAT without peeking into expressions.

Perhaps it would be interesting to experiment to modify TIMESIN to call
GREAT on whole terms, with no peeking? Just wondering out loud here.

Hope this helps,

Robert Dodier

[Stavros Macrakis (Σταῦρος Μακράκης)]

Presumably the reason that timesin sorts the base of a^b rather than the
whole expression is to make a^b and a^c adjacent, so that it can later
simplify that to a^(b+c).

Without experimenting or looking at the code, I predict that modifying
timesin to call great on whole terms will break that very important
simplification.

On Sat, Jul 22, 2017 at 1:58 AM, Robert Dodier <***@gmail.com>
wrote:

> On 2017-07-21, Barton Willis <***@unk.edu> wrote:
>
> > (%i88) (2*x+1)^2 * (x^2-1)$
> >
> > (%i89) args(%);
> > (%o89) [(2*x+1)^2,x^2-1]
> >
> > (%i90) sort(%, orderlessp);
> > (%o90) [x^2-1,(2*x+1)^2]
> >
> > (%i91) a*b*c;
> > (%o91) a*b*c
> >
> > (%i92) args(%);
> > (%o92) [a,b,c]
> >
> > (%i93) sort(args(%), orderlessp);
> > (%o93) [a,b,c]
> >
> > For a*b*c, the term order is the same as orderlessp, for (2*x+1)^2 *
> (x^2-1) it isn't? Or am I being deceived?
>
> TIMESIN calls GREAT to compare terms, the same as 'sort'. But in this
> case, TIMESIN compares only 2*x + 1, not (2*x + 1)^2.
>
> (%i18) (2*x - 1)*(x + 1)^2;
> 1 Enter great [x+1,2*x-1]
> 2 Enter great [x,2*x]
> 3 Enter great [1,2]
> 3 Exit great false
> 2 Exit great false
> 1 Exit great false
> (%o18) (x+1)^2*(2*x-1)
>
> But calling sort compares (2*x + 1)^2.
>
> (%i19) sort(args(%));
> 1 Enter great [(x+1)^2,2*x-1]
> 2 Enter great [(x+1)^2,2*x]
> 3 Enter great [(x+1)^2,x]
> 4 Enter great [x+1,x]
> 5 Enter great [0,1]
> 5 Exit great false
> 4 Exit great true
> 3 Exit great true
> 2 Exit great true
> 1 Exit great true
> (%o19) [2*x-1,(x+1)^2]
>
> Now I'm aware that GREAT looks inside of nonatomic expressions when it's
> figuring out how to sort them -- e.g. a(x)*b(w)*c(z) --> b(w)*a(x)*c(z).
> But in this case, we have GREAT(f(e1), g(e2)) but not GREAT(e1, e2). I
> don't know how common such cases are. We could try to discover other
> cases by e.g. fishing "*" (or "+" for good measure) expressions out of
> the test suite and determining whether args(e) = sort(args(e)).
>
> I see that COMMUTATIVE1 in src/asum.lisp (invoked after declare(foo,
> commutative)) calls GREAT without peeking into expressions.
>
> Perhaps it would be interesting to experiment to modify TIMESIN to call
> GREAT on whole terms, with no peeking? Just wondering out loud here.
>
> Hope this helps,
>
> Robert Dodier
>
>
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Maxima-discuss mailing list
> Maxima-***@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
>

[Robert Dodier]

On 2017-07-22, Stavros Macrakis <***@alum.mit.edu> wrote:

> Presumably the reason that timesin sorts the base of a^b rather than the
> whole expression is to make a^b and a^c adjacent, so that it can later
> simplify that to a^(b+c).
>
> Without experimenting or looking at the code, I predict that modifying
> timesin to call great on whole terms will break that very important
> simplification.

Yes, good point. I think that's right. Incidentally I trawled through
the test suite looking for "*" expressions for which sort(args(e)) #
args(e). These are not too common (or perhaps there's a bug in the way I
was searchng for them). Anyway for the record I found these:

z^((-a)-1)*((a^2+a)*''gamma_greek(a,z)*z+((-a^2)-a)*''gamma_greek(a+1,z))$
a^p*gamma_incomplete(-p,a)$
a^(n-1)*gamma_incomplete(1-n,-a*z)*(-1)^n$
a^(d*z)*h^(c*sqrt(z))$
a^(d*z+e)*h^(c*sqrt(z))$
a^(d*z)*h^(c*sqrt(z)+g)$
a^(d*z+e)*h^(c*sqrt(z)+g)$
a^v*gamma_incomplete(-v,a)$

I see these all involve exponents (and not other identifiable terms).

What about modifying GREAT to implement the same policy as TIMESIN? The
ordering of expressions by GREAT mostly matters because of
simplification, right? Perhaps it makes sense to have GREAT match the de
facto order imposed by simplification. But SIMPLUS also wants to order
terms, and probably differently than TIMESIN. So I don't know if we can
make both of them happy by redefining GREAT.

How about implementing the ordering predicate as a separate function, so
that we can at least say, "The ordering imposed on '*' expressions is
orderlessp_times" and "... '+' expressions is orderlessp_plus". Then at
least one could order terms the same way, if it matters. Maybe it
doesn't. 1/x and -x will throw a monkey wrench into this too.

Just thinking out loud here. I don't have any intention of changing
GREAT or TIMESIN or anything at this time.

best,

Robert Dodier

[Richard Fateman]

There was a burst of publications on simplification/ multiplication without
ordering, a few decades ago.

It is actually not required that terms be ordered in a canonical fashion.
Maple decides on the order based on which symbol names are encountered
first.
What is required is that terms that "match" and can be combined, are
suitably identified.
A hash table with an appropriate hashing function can do the matching
in (probably) linear time, rather than sorting in (probably) n^2 or n^3 or
worse, as done by Maxima's simplifier.

There are implementations of such schemes, but in reality, most expressions
in Maxima are rather short.

RJF


On 7/22/2017 2:54 PM, Robert Dodier wrote:
> On 2017-07-22, Stavros Macrakis <***@alum.mit.edu> wrote:
>
>> Presumably the reason that timesin sorts the base of a^b rather than the
>> whole expression is to make a^b and a^c adjacent, so that it can later
>> simplify that to a^(b+c).
>>
>> Without experimenting or looking at the code, I predict that modifying
>> timesin to call great on whole terms will break that very important
>> simplification.
> Yes, good point. I think that's right. Incidentally I trawled through
> the test suite looking for "*" expressions for which sort(args(e)) #
> args(e). These are not too common (or perhaps there's a bug in the way I
> was searchng for them). Anyway for the record I found these:
>
> z^((-a)-1)*((a^2+a)*''gamma_greek(a,z)*z+((-a^2)-a)*''gamma_greek(a+1,z))$
> a^p*gamma_incomplete(-p,a)$
> a^(n-1)*gamma_incomplete(1-n,-a*z)*(-1)^n$
> a^(d*z)*h^(c*sqrt(z))$
> a^(d*z+e)*h^(c*sqrt(z))$
> a^(d*z)*h^(c*sqrt(z)+g)$
> a^(d*z+e)*h^(c*sqrt(z)+g)$
> a^v*gamma_incomplete(-v,a)$
>
> I see these all involve exponents (and not other identifiable terms).
>
> What about modifying GREAT to implement the same policy as TIMESIN? The
> ordering of expressions by GREAT mostly matters because of
> simplification, right? Perhaps it makes sense to have GREAT match the de
> facto order imposed by simplification. But SIMPLUS also wants to order
> terms, and probably differently than TIMESIN. So I don't know if we can
> make both of them happy by redefining GREAT.
>
> How about implementing the ordering predicate as a separate function, so
> that we can at least say, "The ordering imposed on '*' expressions is
> orderlessp_times" and "... '+' expressions is orderlessp_plus". Then at
> least one could order terms the same way, if it matters. Maybe it
> doesn't. 1/x and -x will throw a monkey wrench into this too.
>
> Just thinking out loud here. I don't have any intention of changing
> GREAT or TIMESIN or anything at this time.
>
> best,
>
> Robert Dodier
>
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Maxima-discuss mailing list
> Maxima-***@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss

[Stavros Macrakis (Σταῦρος Μακράκης)]

I'm not sure I understand what the point of exposing orderlessp_times and
_plus would be.

After all, the orderings are not normally seen by the user, since some
kinds of adjacent terms *do* get combined.

And anyway, I don't think the order of additive or multiplicative terms is
part of the "contract" with the user. The contract is simply that there is
a single consistent order. Users — and other parts of the system —
shouldn't care what that order is.

As for orderlessp, it is just an implementation decision to use the
simplifier's ordering in orderlessp. I don't think we promise anywhere that
additive terms are in orderlessp order.

As RJF says, we could just as well decide to hash instead of sort someday,
and the ordering that induces is pretty opaque.

-s
-s


On Sat, Jul 22, 2017 at 5:54 PM, Robert Dodier <***@gmail.com>
wrote:

> On 2017-07-22, Stavros Macrakis <***@alum.mit.edu> wrote:
>
> > Presumably the reason that timesin sorts the base of a^b rather than the
> > whole expression is to make a^b and a^c adjacent, so that it can later
> > simplify that to a^(b+c).
> >
> > Without experimenting or looking at the code, I predict that modifying
> > timesin to call great on whole terms will break that very important
> > simplification.
>
> Yes, good point. I think that's right. Incidentally I trawled through
> the test suite looking for "*" expressions for which sort(args(e)) #
> args(e). These are not too common (or perhaps there's a bug in the way I
> was searchng for them). Anyway for the record I found these:
>
> z^((-a)-1)*((a^2+a)*''gamma_greek(a,z)*z+((-a^2)-a)*''gamma_greek(a+1,z))$
> a^p*gamma_incomplete(-p,a)$
> a^(n-1)*gamma_incomplete(1-n,-a*z)*(-1)^n$
> a^(d*z)*h^(c*sqrt(z))$
> a^(d*z+e)*h^(c*sqrt(z))$
> a^(d*z)*h^(c*sqrt(z)+g)$
> a^(d*z+e)*h^(c*sqrt(z)+g)$
> a^v*gamma_incomplete(-v,a)$
>
> I see these all involve exponents (and not other identifiable terms).
>
> What about modifying GREAT to implement the same policy as TIMESIN? The
> ordering of expressions by GREAT mostly matters because of
> simplification, right? Perhaps it makes sense to have GREAT match the de
> facto order imposed by simplification. But SIMPLUS also wants to order
> terms, and probably differently than TIMESIN. So I don't know if we can
> make both of them happy by redefining GREAT.
>
> How about implementing the ordering predicate as a separate function, so
> that we can at least say, "The ordering imposed on '*' expressions is
> orderlessp_times" and "... '+' expressions is orderlessp_plus". Then at
> least one could order terms the same way, if it matters. Maybe it
> doesn't. 1/x and -x will throw a monkey wrench into this too.
>
> Just thinking out loud here. I don't have any intention of changing
> GREAT or TIMESIN or anything at this time.
>
> best,
>
> Robert Dodier
>
>
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Maxima-discuss mailing list
> Maxima-***@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
>

[Barton Willis]

>There are implementations of such schemes, but in reality, most expressions in Maxima are rather short.


Yes--running the testsuite of 8 195 344 calls to simptimes, the prodand has 1 term 412 842 times and

two terms 7 040 552 times. In only 61 calls to simptimes, the number of terms in the product is ten or more. It is, I think, hard to win by building a better O algorithm. The distribution of terms in a product for running the testsuite is


[[1, 412842], [2, 7040552], [3, 568802], [4, 97374], [5, 57450], [6, 13551], [7, 3525], [8, 995], [9, 192], [10, 36], [11, 21], [14, 1], [16, 1], [43, 2]]

Multiplication of long lists of complex floats is pokey slow--this was reported to this list earlier this summer:

(%i4) L : makelist(random(1.00) + %i*random(1.00), k,1,10^3)$

(%i6) expand(xreduce("*",L));
Evaluation took 1.2630 seconds (1.2650 elapsed) using 182.193 MB.

(%o6) 4.883033506584155e-166 %i + 1.033163490174675e-166



--Barton



--Barton
________________________________
From: Stavros Macrakis (ΣταῊρος Μακράκης) <***@alum.mit.edu>
Sent: Saturday, July 22, 2017 6:10:04 PM
To: Robert Dodier
Cc: <maxima-***@lists.sourceforge.net>
Subject: Re: [Maxima-discuss] ordering in a product

I'm not sure I understand what the point of exposing orderlessp_times and _plus would be.

After all, the orderings are not normally seen by the user, since some kinds of adjacent terms do get combined.

And anyway, I don't think the order of additive or multiplicative terms is part of the "contract" with the user. The contract is simply that there is a single consistent order. Users — and other parts of the system — shouldn't care what that order is.

As for orderlessp, it is just an implementation decision to use the simplifier's ordering in orderlessp. I don't think we promise anywhere that additive terms are in orderlessp order.

As RJF says, we could just as well decide to hash instead of sort someday, and the ordering that induces is pretty opaque.

-s
-s


On Sat, Jul 22, 2017 at 5:54 PM, Robert Dodier <***@gmail.com<mailto:***@gmail.com>> wrote:
On 2017-07-22, Stavros Macrakis <***@alum.mit.edu<mailto:***@alum.mit.edu>> wrote:

> Presumably the reason that timesin sorts the base of a^b rather than the
> whole expression is to make a^b and a^c adjacent, so that it can later
> simplify that to a^(b+c).
>
> Without experimenting or looking at the code, I predict that modifying
> timesin to call great on whole terms will break that very important
> simplification.

Yes, good point. I think that's right. Incidentally I trawled through
the test suite looking for "*" expressions for which sort(args(e)) #
args(e). These are not too common (or perhaps there's a bug in the way I
was searchng for them). Anyway for the record I found these:

z^((-a)-1)*((a^2+a)*''gamma_greek(a,z)*z+((-a^2)-a)*''gamma_greek(a+1,z))$
a^p*gamma_incomplete(-p,a)$
a^(n-1)*gamma_incomplete(1-n,-a*z)*(-1)^n$
a^(d*z)*h^(c*sqrt(z))$
a^(d*z+e)*h^(c*sqrt(z))$
a^(d*z)*h^(c*sqrt(z)+g)$
a^(d*z+e)*h^(c*sqrt(z)+g)$
a^v*gamma_incomplete(-v,a)$

I see these all involve exponents (and not other identifiable terms).

What about modifying GREAT to implement the same policy as TIMESIN? The
ordering of expressions by GREAT mostly matters because of
simplification, right? Perhaps it makes sense to have GREAT match the de
facto order imposed by simplification. But SIMPLUS also wants to order
terms, and probably differently than TIMESIN. So I don't know if we can
make both of them happy by redefining GREAT.

How about implementing the ordering predicate as a separate function, so
that we can at least say, "The ordering imposed on '*' expressions is
orderlessp_times" and "... '+' expressions is orderlessp_plus". Then at
least one could order terms the same way, if it matters. Maybe it
doesn't. 1/x and -x will throw a monkey wrench into this too.

Just thinking out loud here. I don't have any intention of changing
GREAT or TIMESIN or anything at this time.

best,

Robert Dodier


------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot<https://urldefense.proofpoint.com/v2/url?u=http-3A__sdm.link_slashdot&d=DwMFaQ&c=9ZQuuHhOefNvAzlN-3viIA&r=Y-7rxY5GkJ0PrHulpV2qSA&m=epOosRD9eGApq5l1C1mLmjZliI1QwJWfAYP3M0rM9HM&s=6BSoeZyOY3rSrNXxEa5GPJuorznOTdZSVMIIpL-wEeI&e=>
_______________________________________________
Maxima-discuss mailing list
Maxima-***@lists.sourceforge.net<mailto:Maxima-***@lists.sourceforge.net>
https://lists.sourceforge.net/lists/listinfo/maxima-discuss<https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.sourceforge.net_lists_listinfo_maxima-2Ddiscuss&d=DwMFaQ&c=9ZQuuHhOefNvAzlN-3viIA&r=Y-7rxY5GkJ0PrHulpV2qSA&m=epOosRD9eGApq5l1C1mLmjZliI1QwJWfAYP3M0rM9HM&s=SAvKKAsVZTUQddyCdv7yNr4fceKfAoPGpefFLaAuA6Q&e=>

[Stavros Macrakis (Σταῦρος Μακράκης)]

Re

On Sun, Jul 23, 2017 at 2:07 PM, Barton Willis <***@unk.edu> wrote:

> Multiplication of long lists of complex floats is pokey slow--this was
> reported to this list earlier this summer:
> (%i4) L : makelist(random(1.00) + %i*random(1.00), k,1,10^3)$
> (%i6) expand(xreduce("*",L));
> Evaluation took 1.2630 seconds (1.2650 elapsed) using 182.193 MB.
>

​Actually, you can multiply a long list of floats much more efficiently
than that, first of all by not simplifying it when you construct it (40% of
the time), and secondly by using the appropriate function (rectform rather
than expand).

load(simplifying)$
L: makelist(random(1.00) + %i*random(1.00), k,1,10^3)$

0.0080 sec
M: xreduce("*",L)$
0.1370 sec <<< simplification
N: simpfunmake(verbify("*"),L)$
0.0000 sec <<< no simplification
rectform(M)$ == rectform(N)
0.0080 sec <<< fast

expand(M)$ == expand(N)

0.2150 sec <<< expand is slow

​Because of the way rat simplifies things like %i^2 with algebraic=true,
unfortunately it is also slow, though in other cases it is much faster than
expand:

rat(M),keepfl
​​
oat:true,algebraic
​$​

​
0.2200 sec
​


If you can avoid generating the multiplicative expression in the first
place, you also do pretty well:

(P:1.0, thru 1000 do P:rectform(P*(random(1.0)+%i*random(1.0))))$
0.0170 sec

[Robert Dodier]

On 2017-07-22, Stavros Macrakis <***@alum.mit.edu> wrote:

> I'm not sure I understand what the point of exposing orderlessp_times and
> _plus would be.

Well, the only point would be to be able to state explicitly what's the
ordering of terms in expressions. While I think that's desirable, it
seems a little clumsy as stated.

> And anyway, I don't think the order of additive or multiplicative terms is
> part of the "contract" with the user. The contract is simply that there is
> a single consistent order. Users — and other parts of the system —
> shouldn't care what that order is.

I dunno. As it stands, the ordering imposed by TIMESIN is almost but not
quite the same as orderlessp. One easily slips into believing that they
understand the ordering, but they'd be wrong. The actual ordering is
currently unstated and possibly inexplicable. I don't think this is a
desirable state of affairs.

> As for orderlessp, it is just an implementation decision to use the
> simplifier's ordering in orderlessp. I don't think we promise anywhere that
> additive terms are in orderlessp order.

There's no pressing reason for consistency, except that it makes the
whole system more nearly comprehensible. That seems like a good idea.

> As RJF says, we could just as well decide to hash instead of sort someday,
> and the ordering that induces is pretty opaque.

That would be great. TIMESIN or its replacement could do whatever it
wants, and then the terms could be ordered by orderlessp.

best,

Robert Dodier