[BUG] r1976: Remove argumentless 'random'; pfo.spad

[oldk1331]

I discovered this bug while I am digging another one:

Martin's orignal integral: integrate((a + b*x)/((3 + x^2)*(1 - x^2)^(1/3)), x)
fails. However integrate(1/((3 + x^2)*(1 - x^2)^(1/3)), x) succeeds
while integrate(a/((3 + x^2)*(1 - x^2)^(1/3)), x) fails:

(1) -> integrate(a/((3 + x^2)*(1 - x^2)^(1/3)), x)

>> Error detected within library code:
catdef: division by zero

Note that this is not an "integration error", it comes from catdef.

More interestingly, the third time you invok this expression,
it stucks instead of giving this error.

")set break break" and the backtrace leads me to bringDown
from FunctionSpaceReduce from pfo.spad.

After some "print", I find that there's:
nn : Integer := 1
......
K2Z k == ... setelt!(redmap, k, random(nn)$Z)::F
......
eval(f, lk := kernels f, [K2Z k for k in lk])

The nn is supposed to be the initial seed, however
random(nn) will always return 0 before nn get modified
elsewhere. So when f includes a kernel in denominator,
the divesion by zero error will occur.

This part of code is modified by r1976, which seems to
have bugs:

1) The default value for nn is incorrect, and nn is not
updated after being used in K2Z.

Actually, I think argumentless 'random' is a good idea.

2) Seems the type of random is Integer -> NNI, in this
case it should require a PI. But I'm not very clear with
FunctionSpaceReduce

3) Before r1976, seems this integration still doesn't return
in a long time, maybe there's a bug somewhere else.

[oldk1331]

> 1) The default value for nn is incorrect, and nn is not
> updated after being used in K2Z.
>
> Actually, I think argumentless 'random' is a good idea.

I misunderstood a little, my point is:
1) The default value of nn should not be 1.
2) We should use "1+random(nn)" in this case to avoid 0.

[Waldek Hebisch]

oldk1331 wrote:
>
> I discovered this bug while I am digging another one:
>
> Martin's orignal integral: integrate((a + b*x)/((3 + x^2)*(1 - x^2)^(1/3)), x)
> fails. However integrate(1/((3 + x^2)*(1 - x^2)^(1/3)), x) succeeds
> while integrate(a/((3 + x^2)*(1 - x^2)^(1/3)), x) fails:
>
> (1) -> integrate(a/((3 + x^2)*(1 - x^2)^(1/3)), x)
>
> >> Error detected within library code:
> catdef: division by zero
>
> Note that this is not an "integration error", it comes from catdef.
>
> More interestingly, the third time you invok this expression,
> it stucks instead of giving this error.
>
> ")set break break" and the backtrace leads me to bringDown
> from FunctionSpaceReduce from pfo.spad.
>
> After some "print", I find that there's:
> nn : Integer := 1
> ......
> K2Z k == ... setelt!(redmap, k, random(nn)$Z)::F
> ......
> eval(f, lk := kernels f, [K2Z k for k in lk])
>
> The nn is supposed to be the initial seed, however
> random(nn) will always return 0 before nn get modified
> elsewhere. So when f includes a kernel in denominator,
> the divesion by zero error will occur.

nn is set by 'newReduc'. 'newReduc' should be called before
'bringDown' -- if it is called too late this is a bug.

> This part of code is modified by r1976, which seems to
> have bugs:
>
> 1) The default value for nn is incorrect, and nn is not
> updated after being used in K2Z.

The FSRED code is somewhat fishy. First we estimate
how much randomness we need and call 'newReduc'.
Then there are several calls to 'bringDown'. If we
have the same kernel in two calls the kernel should
be replaced by the same integer. In other words
we should get consistent reduction in several calls.
After we are done with one divisor we call 'newReduc'
and start new series of calls to bringDown'.

AFAICS code using FSRED has design problem: due to random
choice we may reduce denominator to 0. I believe that
in such case we should restart computations using different
reduction. But current code just signals error (which is
definitely wrong).

> 3) Before r1976, seems this integration still doesn't return
> in a long time, maybe there's a bug somewhere else.

Computations with divisors may be quite heavy and take very
long time (days, months and possibly years). It may also
need substantial memory, so after long computation you
may funally run out of memory. Speeding up computations
with divisors is essentially a research problem: there
are some papers giving theoretical algorithms, but complexity
estimates are high enough that it is not clear if such
algorithms will work well in practice. And I am not
aware of descriptions of practical implementations.


--
Waldek Hebisch

[Waldek Hebisch]

In general we may have something like '1/(x + k)' where 'k'
is arbitrary integer. So you can not avoid zero denominator
by adding some bias. In other words: you may reduce probablity
of zero denominator but you can not elimintate it. To
deal with zeros we need to handle failing reduction (restart
computation).

--
Waldek Hebisch

[oldk1331]

> AFAICS code using FSRED has design problem: due to random
> choice we may reduce denominator to 0. I believe that
> in such case we should restart computations using different
> reduction. But current code just signals error (which is
> definitely wrong).

Yes, that should be the correct solution.

> Computations with divisors may be quite heavy and take very
> long time (days, months and possibly years). It may also
> need substantial memory, so after long computation you
> may funally run out of memory.

The problem is that if a=1, this integral can be computed in seconds,
I though replace the constant numerator with a symbol should
not make it takes much longer time.

[oldk1331]

Another integral with exact same problem:

(1) -> integrate((k*x^2 - 1)/((a*k*x + b)*(b*x + a)*sqrt((1 - x^2)*(1
- k^2*x^2))), x)

>> Error detected within library code:
catdef: division by zero

Under FriCAS 1.2.7, this integral returns after 1 minute.

[oldk1331]

Another integral affected by this change:

(1) -> integrate(1/sqrt(1 + b*sec(x))^3,x)

>> Error detected within library code:

RadicalFunctionField: curve is not irreducible

It worked in 1.2.7.

A simple change can make most of these problems go away,
but as discussed before, there's a better solution.

diff --git a/src/algebra/pfo.spad b/src/algebra/pfo.spad
index 816fe22..5904fee 100644
--- a/src/algebra/pfo.spad
+++ b/src/algebra/pfo.spad
@@ -235,7 +235,7 @@ FunctionSpaceReduce(R, F) : Exports == Implementation where
K2Z : K -> F

redmap := table()$AssociationList(K, Z)
- nn : Integer := 1
+ nn : Integer := 10^6

newReduc(n) ==
nn := n