Kernels in AlgebraicNumber

[Ralf Hemmecke]

Dear Waldek,

Basic question: How can I see what kernels are cached?

Explanation follows.

In some program I recursively traverse the operator/kernel/ratfun tree
of an algebraic number and rebuild another algebraic number. I seemingly
have a bug there, but cannot easily find it.

The lines

dbgPrint("csfX z", z)
dbgPrint("csfXxx", x)
t := z*x
dbgPrint("CSFX t", t)


gave me the following output

[:> , csfX z, sqrt(-1)*nthRoot(-1, 3)*sqrt(3)]
[:> , csfXxx, nthRoot(-sqrt(-1)*sqrt(3)+1, 3)]
[:> , CSFX t, 3*sqrt(-1)*nthRoot(-1, 3)]

The value of t is completely wrong.

In another place I see two algebraic numbers that print in my program as

-- aa:=(sqrt(-1)*nthRoot(-1, 3)*sqrt(3)+nthRoot(-1,
3))*nthRoot(-sqrt(-1)*sqrt(3)+1, 3)-2*nthRoot(2, 3)
-- ab:=sqrt(-1)*nthRoot(-1, 3)*sqrt(3)*nthRoot(-sqrt(-1)*sqrt(3)+1,
3)+nthRoot(-1, 3)*nthRoot(-sqrt(-1)*sqrt(3)+1, 3)-2*nthRoot(2, 3)

with respective towers

-- ta:=[sqrt(-1),_
nthRoot(-1, 3),_
nthRoot(2, 3),_
sqrt(3),_
nthRoot(-sqrt(-1)*sqrt(3)+1, 3)]
-- tb:=[sqrt(-1),_
nthRoot(-1, 3),_
sqrt(3),_
nthRoot(2, 3),_
nthRoot(-sqrt(-1)*sqrt(3)+1, 3),_
nthRoot(-sqrt(-1)*sqrt(3)+1, 3)]

The surprising fact for me is in this double entry in tb.
Entering ab and computing "tower ab" in the same session gives:

%%% (6) -> ab:=sqrt(-1)*nthRoot(-1,
3)*sqrt(3)*nthRoot(-sqrt(-1)*sqrt(3)+1, 3)+nthRoot(-1,
3)*nthRoot(-sqrt(-1)*sqrt(3)+1, 3)-2*nthRoot(2, 3)

+----------------+
+---+3+---+ +-+ 3+---+ 3| +---+ +-+ 3+-+
(6) (\|- 1 \|- 1 \|3 + \|- 1 )\|- \|- 1 \|3 + 1 - 2 \|2
Type: AlgebraicNumber
%%% (7) -> tower ab

+----------------+
+---+ 3+---+ 3+-+ +-+ 3| +---+ +-+
(7) [\|- 1 , \|- 1 , \|2 , \|3 , \|- \|- 1 \|3 + 1 ]
Type: List(Kernel(AlgebraicNumber))

That is, no duplication of the last kernel.

Long story short... it seems that my program creates a kernel and the
cache does not recognize that it is the same. In order to find exactly
where this happens, I would like to know how I can look at the kernel
cache from SPAD.

Thank you.
Ralf

[Waldek Hebisch]

On Sun, Mar 16, 2025 at 07:46:43AM +0100, 'Ralf Hemmecke' via FriCAS - computer algebra system wrote:
> Dear Waldek,
>
> Basic question: How can I see what kernels are cached?

When I needed to debug kernel related functions I added bunch
of printouts to SortedCache and Kernel. They were not suitable
for normal use, so I did not commit them. If you want to
debug this you may add your own printouts.

> Explanation follows.
>
> In some program I recursively traverse the operator/kernel/ratfun tree of an
> algebraic number and rebuild another algebraic number. I seemingly have a
> bug there, but cannot easily find it.

I would suggest to avoid creating new kernels in AlgebraicNumber.
At implementation level AlgebraicNumber is essentially a slightly
changed view of Expression(Integer) and IIRC corresponding
kernels live in cache for Kernel(Expression(Integer)). If you
manage to put kernel into cache for Kernel(AlgebraicNumber),
this is likely to lead to confusion, possibly with effects
like below.

It is better to keep kernel transformations in Expression(Integer).
I think that currently all operations on kernels within
AlgebraicNumber are delegated to Expression(Integer).

BTW: It probably would be better to change our declarations so
that AlgebraicNumber does not export FunctionSpace, but doing
that without removing useful operations looks tricky.

> The lines
>
> dbgPrint("csfX z", z)
> dbgPrint("csfXxx", x)
> t := z*x
> dbgPrint("CSFX t", t)
>
>
> gave me the following output
>
> [:> , csfX z, sqrt(-1)*nthRoot(-1, 3)*sqrt(3)]
> [:> , csfXxx, nthRoot(-sqrt(-1)*sqrt(3)+1, 3)]
> [:> , CSFX t, 3*sqrt(-1)*nthRoot(-1, 3)]
>
> The value of t is completely wrong.

Look at position fields of kernels. Kernel equality assumes
that kernels with the same nonzero positions are equal:

k1 = k2 ==
p1 := position(k1)
p2 := position(k2)
p1 ~= 0 and p2 ~= 0 => p1 = p2
if p1 = 0 then k1 := kernelEnterInCache(k1)
if p2 = 0 then k2 := kernelEnterInCache(k2)
position(k1) = position(k2)

If you manage to mix kernels from different caches, then this
leads to unexpected results. Similarly, when kernel is removed
from the cache, but is position filed in not set to 0 (I hope
that problems with removal are fixed now, but it happende in
the past).

--
Waldek Hebisch

[Ralf Hemmecke]

On 3/16/25 12:21, Waldek Hebisch wrote:
> On Sun, Mar 16, 2025 at 07:46:43AM +0100, 'Ralf Hemmecke' via FriCAS - computer algebra system wrote:

> I would suggest to avoid creating new kernels in AlgebraicNumber.
> At implementation level AlgebraicNumber is essentially a slightly
> changed view of Expression(Integer) and IIRC corresponding
> kernels live in cache for Kernel(Expression(Integer)).

Aaaahh.... I was just about to test and wondered why the cache for
SortedCache(Kernel(AlgebraicNumber)) printed as empty though I entered
sqrt(2) into the session. Additionally, I get

%%% (1) -> a2:=sqrt(2)

+-+
(1) \|2
Type: AlgebraicNumber
%%% (2) -> k2 := kernels(a2).1

+-+
(2) \|2
Type: Kernel(AlgebraicNumber)
%%% (3) -> position(k2)

(3) 2048
Type: PositiveInteger
%%% (4) -> position(k2)$Kernel(AlgebraicNumber)

(4) 2048
Type: NonNegativeInteger

I added
getCache(): PrimitiveArray S == cache
to SortedCache.

%%% (5) -> getCache()$SortedCache(Kernel(Expression(Integer)))

+-+
(5) [%%var, \|2 , %%var, %%var, %%var, %%var, %%var, %%var, %%var,
%%var]
Type: PrimitiveArray(Kernel(Expression(Integer)))
%%% (6) -> getCache()$SortedCache(Kernel(AlgebraicNumber))

(6) []
Type: PrimitiveArray(Kernel(AlgebraicNumber))

So certainly I must create kernels in Kernel(Expression Integer). That's
a good hint.

However, why position gives such a big number whereas the true position
in the primitive array should actually be 2, is a mystery for me.

Anyhow, if I could ensure all my kernels are truely in Kernel(AN) (which
I probably cannot), then all would be fine?

> BTW: It probably would be better to change our declarations so
> that AlgebraicNumber does not export FunctionSpace, but doing
> that without removing useful operations looks tricky.

Erm, according to https://fricas.github.io/api/AlgebraicNumber.html, it
doesn´t. You probably meant ExpressionSpace.
My first thought was to implement explicitly all those exports in AN via
coercion to E==>Expression(INT), calling the respective function from E
and then converting back the result. But that cannot work without
confusion, because if a result contains Kernel(AN), then an entry in
SortedCache(Kernel(AN)) would be created, and two kernels from AN might
have different positions in SortedCache(Kernel(AN)) and
SortedCache(Kernel(E)). So when k1 < k2 in Kernel(AN) may give a
different result then the one in Kernel(E).
Wait, something must be wrong with my thinking. How would otherwise
getCache()$SortedCache(Kernel(AlgebraicNumber)) still return the empty
list although I clearly have created an element of Kernel(AN) in (2) above?

Thank you
Ralf

[Waldek Hebisch]

Positions have gaps, that is deliberate to limit renumbering.

> Anyhow, if I could ensure all my kernels are truely in Kernel(AN) (which I
> probably cannot), then all would be fine?
>
> > BTW: It probably would be better to change our declarations so
> > that AlgebraicNumber does not export FunctionSpace, but doing
> > that without removing useful operations looks tricky.
>
> Erm, according to https://fricas.github.io/api/AlgebraicNumber.html, it
> doesn´t. You probably meant ExpressionSpace.

Yes, ExpressionSpace.

> My first thought was to implement explicitly all those exports in AN via
> coercion to E==>Expression(INT), calling the respective function from E and
> then converting back the result. But that cannot work without confusion,
> because if a result contains Kernel(AN), then an entry in
> SortedCache(Kernel(AN)) would be created, and two kernels from AN might have
> different positions in SortedCache(Kernel(AN)) and SortedCache(Kernel(E)).
> So when k1 < k2 in Kernel(AN) may give a different result then the one in
> Kernel(E).
> Wait, something must be wrong with my thinking. How would otherwise
> getCache()$SortedCache(Kernel(AlgebraicNumber)) still return the empty list
> although I clearly have created an element of Kernel(AN) in (2) above?

AFAICS this is really Kernel(Expression(Integer)). They share
representation, so in many case it does not matter if it
is Kernel(AlgebraicNumber) or Kernel(Expression(Integer)). But
it matters for cache. And once you mix kernels from two caches
you may get weird results.

ATM I do not have good solution to this. Correct solution would
probably remove all traces of Kernel(AlgebraicNumber) from
algebra code, but currently it is used in several important
operations.

BTW: Are you recompiling some domains? More precisely, do you
get wrong results when you )LIB files that you recompiled
separately and avoid any operations that could clear kernel
cache?

--
Waldek Hebisch

[Ralf Hemmecke]

On 3/16/25 16:42, Waldek Hebisch wrote:
>> %%% (3) -> position(k2)
>>
>> (3) 2048

> Positions have gaps, that is deliberate to limit renumbering.

OK. But then I must have missed the place where these gaps are
introduced. In fact, I do not quite get the ordering of the kernels.
I saw the triage function that is used for comparison of kernels.

triage(k1, k2) ==
height(k1) ~= height(k2) => B2Z(height(k1) < height(k2))
operator(k1) ~= operator(k2) => B2Z(operator(k1) < operator(k2))
(n1 := #(argument k1)) ~= (n2 := #(argument k2)) => B2Z(n1 < n2)
-- Handled by linear search earlier
-- ((func := property(operator k1, SPECIALEQUAL)) case None) and
-- (((func::None) pretend ((%, %) -> Boolean)) (k1, k2)) => 0
for x1 in argument(k1) for x2 in argument(k2) repeat
x1 ~= x2 => return B2Z(smaller?(x1, x2))
0

That seems to suggest that if the same kernels are introduced in two
different sessions they would inserted into the cache array in the same
order. Maybe not in the same position, but comparing the kernels with <
should be the same in both sessions. (OK, except when "smaller?" behaves
differently.)

> BTW: Are you recompiling some domains? More precisely, do you
> get wrong results when you )LIB files that you recompiled
> separately and avoid any operations that could clear kernel
> cache?

The behaviour I experienced was not done with )lib. But in my original
setup I started a fresh fricas, compiled my code then ran it. It used
"trueEqual $ AN" with my seemingly wrongly constructed algebraic numbers
and trueEqual crashed, because there appeared a sqrt(3) in the
polynomial computed by "norm" which couldn't be retracted to
Fraction(INT). Interestingly, when I afterwards ran the same code for
the second time, it went through, giving even the right result.
I meanwhile know that compiling is dangerous for the cache.

Anyway, than you for your hint to Kernel(Expression INT). I wouldn't
have seen this so easily. With that I can probably continue.

Ralf

[Waldek Hebisch]

On Sun, Mar 16, 2025 at 05:05:02PM +0100, 'Ralf Hemmecke' via FriCAS - computer algebra system wrote:
> On 3/16/25 16:42, Waldek Hebisch wrote:
> > > %%% (3) -> position(k2)
> > >
> > > (3) 2048
>
> > Positions have gaps, that is deliberate to limit renumbering.
>
> OK. But then I must have missed the place where these gaps are introduced.

See 'enterInCache' in 'kl.spad (second signature) and constant DIFF.

> In fact, I do not quite get the ordering of the kernels.
> I saw the triage function that is used for comparison of kernels.
>
> triage(k1, k2) ==
> height(k1) ~= height(k2) => B2Z(height(k1) < height(k2))
> operator(k1) ~= operator(k2) => B2Z(operator(k1) < operator(k2))
> (n1 := #(argument k1)) ~= (n2 := #(argument k2)) => B2Z(n1 < n2)
> -- Handled by linear search earlier
> -- ((func := property(operator k1, SPECIALEQUAL)) case None) and
> -- (((func::None) pretend ((%, %) -> Boolean)) (k1, k2)) => 0
> for x1 in argument(k1) for x2 in argument(k2) repeat
> x1 ~= x2 => return B2Z(smaller?(x1, x2))
> 0
>
> That seems to suggest that if the same kernels are introduced in two
> different sessions they would inserted into the cache array in the same
> order. Maybe not in the same position, but comparing the kernels with <
> should be the same in both sessions. (OK, except when "smaller?" behaves
> differently.)

Roughly yes. There are two irregularites. Some kernels have special
equality and 'smaller?' is not an order. So if you insert kernel
in a different order also order on kernels may be different.
And extra kernels may affect result (in real ordered set order
on a subset is independent on extra elements in the set, but
we have only "almost order").

--
Waldek Hebisch

[Ralf Hemmecke]

> See 'enterInCache' in 'kl.spad (second signature) and constant DIFF.

Aha. Thanks.

Since we are close to SortedCache... The following can do without the x
parameter.

linearSearch : (S, S -> Boolean) -> Union(S, "failed")
++ linearSearch(x, f) searches x in the cache, calling \spad{f(y)}
++ to determine whether x is equal to y. It returns y from cache
++ if f(y) is true or "failed" if no such y exists.


linearSearch(x : S, equal? : S -> Boolean) ==
k : Integer := 0
-- Can not use for loop because equal? can insert new elements
-- and change cache_use
while k < cache_use repeat
vscan := cache
y := vscan(k)
equal?(y) => return y
vscan := cache
-- skip over elements possibly inserted by equal?
while not(EQ(y, vscan(k))$Lisp) repeat k := k + 1
k := k + 1
return "failed"

Bur maybe it's not worth to change it.

Ralf

[Waldek Hebisch]

I have now removed most uses of Kernel(AlgebraicNumber). But
remaining ones (in AlgebraicNumber) are more tricky: if I remove
several exported functions and no longer export ExpressionSpace,
then AlgebraicNumber compiles, but lacks important
functionality.

--
Waldek Hebisch

[Ralf Hemmecke]

> I have now removed most uses of Kernel(AlgebraicNumber). But
> remaining ones (in AlgebraicNumber) are more tricky: if I remove
> several exported functions and no longer export ExpressionSpace,
> then AlgebraicNumber compiles, but lacks important
> functionality.

Thanks. I have also rewriten my code and now construct kernels as
Kernel(Expression Integer). That solved all my issues as far as I can see.

Obviously, if in AlgebraicNumber functions are called that are inherited
from the Expression Integer representation, I think there is anyway no
issue, because of the way inheritance works. (OK, that means, as far as
I believe it works.) So if in those functions kernels are created they
are automatically Kernel(Expression(INT)) and the inheritance mechanism
then just claims/pretends that thes kernel are in Kernel(AN).

Problem comes if a user create a Kernel(AN) element that is not actually
an element of Kernel(Expression(INT)) that is pretended to be Kernel(AN).

So for the particular case of AN, maybe Kernal(X) checks whether X is AN
and den does the kernel creation via Expression(INT).
Maybe this would work, but I would find this very, very bad design if
the Kernel constructor has to know about AN and Expression(INT)
explicitly and thus, do certainly NOT support such an approach.

In fact, I find it a bit strange that in the design of AN it was decided
to internally build on Kernel(Expression(INT)). Well, maybe it was not
decided, but just happened (for most of the cases), because many
functions from ExpressionSpace are just inherited from Rep:=Expression(INT).

A tricky issue, indeed.

Ralf

[Waldek Hebisch]

On Tue, Mar 18, 2025 at 01:59:47PM +0100, 'Ralf Hemmecke' via FriCAS - computer algebra system wrote:
>
> In fact, I find it a bit strange that in the design of AN it was decided to
> internally build on Kernel(Expression(INT)). Well, maybe it was not decided,
> but just happened (for most of the cases), because many functions from
> ExpressionSpace are just inherited from Rep:=Expression(INT).

It was natural to inherit representation and operations from
Expression(Integer). But operations in Expression(Integer)
produce kernels from Kernel(Expression(Integer)). Also,
AlgebraicNumber has operations like norm and subst which
take kernels as arguments. Those operations are defined
in ExpressionSpace, so we make AlgebraicNumber into
ExpressionSpace. Oops, now the kernels live in two caches :(

There seem to be two possible ways to solve this. We have
now ExpressionSpace2 which can express fact that kernels
are not from Kernel(AlgebraicNumber), so we could declare
that AlgebraicNumber is ExpressionSpace2(SomeNewKernel)
and implement SomeNewKernel as Kernel(Expression(Integer)).

Alternatively, we can give some new signatures to 'norm'
and 'subst'. Actually, 'norm' could be declared in
some Category and implemented in Expression. But trouble
is that currently there in Kernel in signature of 'norm'.

--
Waldek Hebisch