rootSplit

13 views
Skip to first unread message

Ralf Hemmecke

unread,
Jul 2, 2024, 6:25:25 AMJul 2
to fricas-devel
I wonder whether rootSplit in manip.spad is implemented correctly.

rootkernels l == select!((z1 : K) : Boolean +->
is?(operator z1, 'nthRoot), l)
rootSplit x ==
lk := rootkernels tower x
eval(x, lk, [rootExpand k for k in lk])

rootExpand k ==
x := first argument k
n := second argument k
op := operator k
op(numer(x)::F, n) / op(denom(x)::F, n)

AFAIU, it evaluates the rational function that represents x
by the respective "rootExpand"ed stuff.

With that implementation, it is no surprise that
for expr := sqrt(y+sqrt((z+1)/(z-1))))

rootSplit expr

gives back the same expression.

I wonder whether this indeed was the intention.

%%% (29) -> rootSplit((y + sqrt((z+1)/(z-1))))

+-----+ +-----+
\|z + 1 + y\|z - 1
(29) --------------------
+-----+
\|z - 1

%%% (30) -> rootSplit(sqrt(y + sqrt((z+1)/(z-1))))

+------------+
| +-----+
| |z + 1
(30) | |----- + y
\|\|z - 1

Otherwise I wonder why tower(x) is used and not just kernels(x).

Any opinions?

Ralf

Ralf Hemmecke

unread,
Jul 3, 2024, 8:34:49 AMJul 3
to fricas-devel
It looks like smpeval

smpeval(p, lk, lv) ==
map(x +-> match(lk, lv, x,
notfound((z : K) : % +-> map(s +-> eval(s, lk, lv), z),
lk))$ListToMap(K, %),
y +-> y::%,
p)$PolynomialCategoryLifting(IndexedExponents K, K, R, MP, %)

is supposed to apply "eval" recursively.

However, when one analyses what match(lk,lv,x,f) actually does, then it
means, an extension of the mapping k+->v for any value k in lk by
applying f if k is not in lk.

Now the f here is this "notfound" function.

notfound(fn, lk) ==
k +->
empty? setIntersection(tower(f := k::%), lk) => f
fn k

and that says: apply fn whenever in the tower of k is an element from lk.

In particular this function will never be applied to any value that is
in lk, because these are handled directly by "match(lk,lv,x)".

Clearly, "rootSplit" starts with "lk := rootkernels tower x".
and in "eval(x, lk, [rootExpand k for k in lk])" the values
are (unfortunately) NOT recursivly computed in case lk is something like

[sqrt(6), sqrt(sqrt(6))]

Thus, in rootFactor(sqrt(sqrt(6)), smpeval is called with

lk = [sqrt(6), sqrt(sqrt(6))]]
lv = [sqrt(2)*sqrt(3), sqrt(sqrt(6))]]

and it is clear what will happen. The toplevel rational function
representing sqrt(sqrt(6)) is just k/1 for the kernel sqrt(sqrt(6)).
So eventually sqrt(sqrt(6)) is replaced by itself and not a recursively
evaluated value.

I think that this is not what was intended by Bronstein, otherwise, I do
not (yet) see, why he would have put the "eval" inside the "notfound" call.

That case certainly appears here:

%%% (1) -> rootFactor(log(sqrt 6)*sqrt(sqrt(6)))
+----+
+-+ +-+ | +-+
(3) log(\|2 \|3 )\|\|6

I find it strange that substitution inside log does happen, but not
inside sqrt.

Opinions? Should I try to come up with a proper fix?
I somehow feel such a change will have quite some consequences in
different part of FriCAS.

Ralf

Waldek Hebisch

unread,
Jul 4, 2024, 10:41:51 AM (13 days ago) Jul 4
to fricas...@googlegroups.com
Impact should be limited. Functions in AlgebraicManipulations are
mostly "user level" functions. It seems that 'ratDenom' and
'rootSimp' are called from rest of algebra. IMO we should get
rid of calls to 'rootSimp' (instead extend normalization to handle
algebraic kernels). As you noticed 'ratDenom' may be unreliable,
so probably we should replace is by something else in most of algebra.

Concerning 'rootSplit', it make sense to create better implementation.
My impression is that original AlgebraicManipulations just used
very simple tricks to get functionality similar to what other CAS-es
offerd. But it was very limited. There are now some extentions
and enhancements, but it is still limited.

--
Waldek Hebisch

Ralf Hemmecke

unread,
Jul 4, 2024, 12:28:11 PM (13 days ago) Jul 4
to fricas...@googlegroups.com
On 7/4/24 16:41, Waldek Hebisch wrote:
>> Opinions? Should I try to come up with a proper fix? I somehow
>> feel such a change will have quite some consequences in different
>> part of FriCAS.
>
> Impact should be limited. Functions in AlgebraicManipulations are
> mostly "user level" functions. It seems that 'ratDenom' and
> 'rootSimp' are called from rest of algebra.

OK, then I propose the attached patch. Admittedly, I have not looked
deep into the Expression related stuff, so I hope that someone more
familiar with Expression looks over it.

The basic idea of the change I have already analysed in an earlier
message. There is the variable lv of kernels of the input expression the
tower of kernels with nthRoot as operator.

One call to root_factor_k basically only changes the leaves of the
kernel tree. So if there are nested roots the unchanged value would be
substituted back during eval.

This patch changes it such that first step by step every element of lv
get the evaluations with the elements from before.
This, of course, assumes that "rootkernels tower x" returns the kernels
from the leaves to the root.
Glancing at the implementation of SortedCache, tells me that this is
might be always true. Is that correct?

> IMO we should get rid of calls to 'rootSimp' (instead extend
> normalization to handle algebraic kernels).

Yes. The more we have properly specified code, the better. rootSimp
sounds like an unspecifiable function.

> As you noticed 'ratDenom' may be unreliable, so probably we should
> replace is by something else in most of algebra.

That would be good. The reason why I am currently looking into it is
that I want some kind of normal form for AlgebraicNumber. In my
computations there appear two of them in different form and I would like
to be able to verify that they are equal (and, of course, it looks
better, if the output is visually small). I don't need to solve every
case but nested sqrt and nthRoot(..., 3) would already be a good progress.

> Concerning 'rootSplit', it make sense to create better
> implementation.

It seems that rootFactor actually also covers rootSplit. But probably
rootSplit was supposed to do things without factorization.
Anyway, if you agree with the way I "solve" the problem, I can look at
the other cases.

> My impression is that original AlgebraicManipulations just used very
> simple tricks to get functionality similar to what other CAS-es
> offerd. But it was very limited.

Yes. But I think, if we agree that AlgebraicManipulations should stay at
top-level, i.e., convenience functions for users and not otherwise used
in the algebra library, then we should make those functions as
convenient as possible to assist the user in transforming expressions in
the way they want them.

> There are now some extentions and enhancements, but it is still
> limited.

What do you mean here?

Ralf
0001-fix-rootFactor-function-to-apply-recursively.patch

Waldek Hebisch

unread,
Jul 4, 2024, 3:30:44 PM (13 days ago) Jul 4
to fricas...@googlegroups.com
On Thu, Jul 04, 2024 at 06:28:01PM +0200, Ralf Hemmecke wrote:
> On 7/4/24 16:41, Waldek Hebisch wrote:
> > > Opinions? Should I try to come up with a proper fix? I somehow
> > > feel such a change will have quite some consequences in different
> > > part of FriCAS.
> >
> > Impact should be limited. Functions in AlgebraicManipulations are
> > mostly "user level" functions. It seems that 'ratDenom' and 'rootSimp'
> > are called from rest of algebra.
>
> OK, then I propose the attached patch. Admittedly, I have not looked
> deep into the Expression related stuff, so I hope that someone more
> familiar with Expression looks over it.
<snip>
> This, of course, assumes that "rootkernels tower x" returns the kernels
> from the leaves to the root.

Actually, I am not sure if this is true in corner cases. But
basically all expresion code expects such behaviour, so if needed
we should fix tower so that this holds.

> > As you noticed 'ratDenom' may be unreliable, so probably we should
> > replace is by something else in most of algebra.
>
> That would be good. The reason why I am currently looking into it is
> that I want some kind of normal form for AlgebraicNumber. In my
> computations there appear two of them in different form and I would like
> to be able to verify that they are equal (and, of course, it looks
> better, if the output is visually small). I don't need to solve every
> case but nested sqrt and nthRoot(..., 3) would already be a good progress.

So, basically you want 'normalize' to work with algebraic expression.
That is high priority item on my toto list.

> > Concerning 'rootSplit', it make sense to create better implementation.
>
> It seems that rootFactor actually also covers rootSplit. But probably
> rootSplit was supposed to do things without factorization.

Yes.

> Anyway, if you agree with the way I "solve" the problem, I can look at
> the other cases.

I would do this somewhat differently. Untested code below:

rootExpand(lk, lv, k) ==
x := first argument k
n := second argument k
op := operator k
x := eval(x, lk, lv) --> eval here
op(numer(x)::F, n) / op(denom(x)::F, n)

rootSplit x ==
lk0 := rootkernels tower x
lk : List(K) := []
lv : List(F) := []
for k in lk0 repeat
v := rootExpand(lk, lv, k)
lk := cons(k, lk)
lv := cons(v, lv)
eval(x, lk, lv)

What is the difference? Main is that eval is done _before_ solitting,
that allows cancellation between numerator and denominator. And
little differences like using 'cons' to _prepend_ to a list. Compared
to concat this avoids quadratic behaviour. Not very important here,
but IMO it is good to keep style consistent and this is normal
style in FriCAS code.

> > My impression is that original AlgebraicManipulations just used very
> > simple tricks to get functionality similar to what other CAS-es offerd.
> > But it was very limited.
>
> Yes. But I think, if we agree that AlgebraicManipulations should stay at
> top-level, i.e., convenience functions for users and not otherwise used
> in the algebra library, then we should make those functions as
> convenient as possible to assist the user in transforming expressions in
> the way they want them.

Yes it makes sense to improve them.

> > There are now some extentions and enhancements, but it is still
> > limited.
>
> What do you mean here?

I added rootFactor, there were various improvements, especially
by Qian. But the approach looks limited, I do not see how to
extend it to do things which and be done using different approach.
And to that matter, unstated assumption is that user can do
desired transformation by sequence of simpler transformations.
Currently I think that _single_ say 'transform' function with
many options would be better. The point is that multiple options
are simpler both for user and implementation than multiple
functions. And interactions of options can produce results
that would be hard to do as sequence of transformations.

--
Waldek Hebisch

Ralf Hemmecke

unread,
Jul 5, 2024, 4:20:26 AM (13 days ago) Jul 5
to fricas...@googlegroups.com
On 7/4/24 21:30, Waldek Hebisch wrote:

> So, basically you want 'normalize' to work with algebraic
> expression.

Yes, that would be nice, but more urgent for me is to have that just for
algebraic numbers, i.e. no variables involved. So there might be even
some properties that can be used. For example, I've realized that
sometimes there appears a sqrt(-1) factor while applying rootFactor.
(The following is with the current trunk code.)

%%% (3) -> rootFactor((sqrt(sqrt(2)-3*sqrt(3))))

+-------------+
+---+ | +-+ +-+
(3) \|- 1 \|3 \|3 - \|2

%%% (4) -> rootFactor((sqrt(5*sqrt(2)-sqrt(3))))

+-------------+
+---+ | +-+ +-+
(4) \|- 1 \|\|3 - 5 \|2


That looks OK for the first case, but definitely not for the second.
Also

%%% (57) -> a := sqrt((3*sqrt(2)-sqrt(3)::E)*(7*sqrt(5)-sqrt(7)::E))

+----------------------------------------------+
| +-+ +-+ +-+ +-+ +-+ +-+
(57) \|(\|3 - 3 \|2 )\|7 + (- 7 \|3 + 21 \|2 )\|5

%%% (58) -> [a::Complex(Float), rootFactor(a)::Complex(Float)]

(58) [5.7144160659_439266652, - 5.7144160659_439266652]

does not seem to be desirable even though it matches the current
specification of rootFactor.

At least for AlgebraicNumber, I would like to see an implementation that
does not change the value of the algebraic number if the value under the
root is positive. Yes, that is a bit more work, but if rootFactor
changes the sign somewhere in a nested root expression (that is in fact
positive at all places) that can change the floating point value of that
number in a pretty unpredictable way and thus becomes useless.

I somehow think of checking positiveness of the expression and their
factors with

algnum::Complex(DoubleFloat)::Float.

Similarly, one could have a rootCombine that doesn't chang the value.
I guess all that should live in AlgebraicNumber directly and fall under
the name "normalize" (or is even silently applied)/

ad ((your suggested code))

> What is the difference? Main is that eval is done _before_
> splitting, that allows cancellation between numerator and
> denominator.

OK, that makes sense.

> And little differences like using 'cons' to _prepend_ to a list.
> Compared to concat this avoids quadratic behaviour.

I actually knew about that quadratic issue and wanted to do it with
cons, but wasn't completely sure whether eval wants its arguments in a
specific sorting order.

> I added rootFactor, there were various improvements, especially by
> Qian. But the approach looks limited, I do not see how to extend it
> to do things which and be done using different approach. And to that
> matter, unstated assumption is that user can do desired
> transformation by sequence of simpler transformations. Currently I
> think that _single_ say 'transform' function with many options would
> be better. The point is that multiple options are simpler both for
> user and implementation than multiple functions. And interactions
> of options can produce results that would be hard to do as sequence
> of transformations.

I stayed away from expression transformation since long. I feel, it's a
hard problem and I began to love AXIOM, because it worked with canonical
representations in most places.

Introducing a 'transform' function with options is probably a good idea.
Honestly, I would have liked that such a function already existed in
FriCAS. Programming it myself, is currently not my main interest. I go
for simple improvements that help me with what I currently need.

Ralf
Reply all
Reply to author
Forward
0 new messages