fyi, Fricas 1.3.10 still give division by zero on some integrals.

[Nasser M. Abbasi]

FYI,

When first used, Fricas gives division by zero on some integrals. When called again, the error goes away. This still happens in V 1.3.10. So all these count as failed since not possible to try again same integral as each call uses new instance of Fricas process each time.

>fricas
FRICAS="/usr/local/lib/fricas/target/x86_64-linux-gnu"
spad-lib="/usr/local/lib/fricas/target/x86_64-linux-gnu/lib/libspad.so"
                       FriCAS Computer Algebra System
                Version: FriCAS 1.3.10 built with sbcl 2.3.11
                 Timestamp: Wed Jan 10 09:37:52 PM CST 2024

(1) -> setSimplifyDenomsFlag(true)
(2) -> integrate(1/(2^(2/3)-x)/(x^3-1)^(1/2),x) 
   >> Error detected within library code:
   catdef: division by zero

(2) -> integrate(1/(2^(2/3)-x)/(x^3-1)^(1/2),x)
   (2)
          +-+
         \|3
etc...

Same for  (with minus sign)

1) -> setSimplifyDenomsFlag(true)
   (1)  false
(2) ->  integrate(1/(2^(2/3)-x)/(-x^3-1)^(1/2),x)
   >> Error detected within library code:
   catdef: division by zero

I noticed some of these go away when not using  setSimplifyDenomsFlag(true). From new session:

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

           x
         ++              1
   (1)   |   ------------------------ d%B
        ++                +---------+
              3+-+2       |    3
             (\|2   - %B)\|- %B  - 1

No division by zero.

--Nasser

[Qian Yun]

Hi Waldek,

Since the method used by "possibleOrder" is random in nature,
so it doesn't harm to try a few times right?

With following patch it gives another try when encountering "division
by zero". And now the failure rate drops to lower than 10%. Of course
we can try more times.

- Qian

diff --git a/src/algebra/pfo.spad b/src/algebra/pfo.spad
index 22deba96..ba961c56 100644
--- a/src/algebra/pfo.spad
+++ b/src/algebra/pfo.spad
@@ -376,7 +376,12 @@ PointsOfFiniteOrder(R0, F, UP, UPUP, R) : Exports
== Implementation where
kbad3Num(h, m) == lcm [kbadBadNum(c, m) for c in coefficients h]

torsionIfCan d ==
- zero?(n := possibleOrder(d := reduce d)) => "failed"
+ d := reduce d
+ nu : Union(N, "failed") := trappedSpadEval(possibleOrder(d))$Lisp
+ if nu case "failed" then nu := trappedSpadEval(possibleOrder(d))$Lisp
+ if nu case "failed" then error "tried twice"
+ n := nu::N
+ zero?(n) => "failed"
(g := generator reduce(n::Z * d)) case "failed" => "failed"
[n, g@R]



On 1/24/24 22:05, 'Nasser M. Abbasi' via FriCAS - computer algebra

[Waldek Hebisch]

On Thu, Jan 25, 2024 at 06:36:12PM +0800, Qian Yun wrote:
> Hi Waldek,
>
> Since the method used by "possibleOrder" is random in nature,
> so it doesn't harm to try a few times right?
>
> With following patch it gives another try when encountering "division
> by zero". And now the failure rate drops to lower than 10%. Of course
> we can try more times.

Well, there is a harm: 'possibleOrder' is potentially quite
expensive and it may be particularly expensive in failing
cases.
A little theory: 'possibleOrder' reduces divisor, first
plugging in integers for parameters and then reducing
coefficients modulo a prime. With current code this part
may cause division by zero even on "good" divisors.

The result of this part is a divisor over finite field.
Over finite field we try to compute order of divisor
essentially by brute force: we compute all powers of
divisor until we get a divisor of a function.

Theory says that it is enough to compute order over finite
field once: we either get actual order (if divisor is of
finite order) or divisor is of infinite order. So
after result of 'possibleOrder' we need to check to
exclude cases of infinite order (this is done by the
line with 'generator'). This check again may be quite
expensive. To reduce chance of failure in the final
check in many cases we run main part of 'possibleOrder'
twice. If second run produces division by zero, then
the results of first run is effectively lost.

So, IMO 'torsionIfCan' is not the right place to catch
errors. Rather we should catch them at lower level, that
is check if reduction to finite field worked. That way
we could reuse computations.

Also, there is another aspect: in our code randomness comes
from integers which are substituted for paramemeters. If
there are no parameters, then there is no randomness.
Theoretically prime could be rendomized, but cost of
computation grows (potentially quite fast) with size of
the prime, so currently we try smallest suitable prime.
But if changing parameters does not help we could try
bigger prime.

BTW: I wrote "with current code", because at least theoretically
there is a way to avoid division by zero when reducing
divisors to obtain divisor over finite field. But
this requires checking at various stages and more important,
a different representaion of divisors: our current code can only
represet divisors which are zero "at infinity", but
division by zero means that we will get part at infinity.

BTW2: Keeping reduction to finite field fixed at least in
principle allows various efficiency improvements. And different
representation of divisors should give us better efficiency.
Also, different representation is needed to complete Risch
algorithm.

--
Waldek Hebisch

[Ralf Hemmecke]

Hi Qian,

just a side remark...

On 1/25/24 11:36, Qian Yun wrote:

> +      nu : Union(N, "failed") := trappedSpadEval(possibleOrder(d))$Lisp

Oh, trappedSpadEval $ Lisp looks like a useful function.
Sounds like a "try [CODE] catch ..." construction.

What exactly can it do?

Hmmmm... unfortunately it is another "foo $ Lisp" appearing in code.
Maybe that function can be put into a package in order to reduce the
"Lisp" occurrences in ordinary SPAD code.

Ralf

[Qian Yun]

On 1/25/24 22:46, Ralf Hemmecke wrote:
> Hi Qian,
>
> just a side remark...
>
> On 1/25/24 11:36, Qian Yun wrote:
>
>> +      nu : Union(N, "failed") := trappedSpadEval(possibleOrder(d))$Lisp
>
> Oh, trappedSpadEval $ Lisp looks like a useful function.
> Sounds like a "try [CODE] catch ..." construction.
>
> What exactly can it do?

If "exprA" gives result of typeA, then "trappedSpadEval(exprA)$Lisp"
gives type Union(typeA, "failed"). Its value takes typeA branch if
"exprA" successes without error, takes "failed" branch if "exprA"
fails with error.

> Hmmmm... unfortunately it is another "foo $ Lisp" appearing in code.
> Maybe that function can be put into a package in order to reduce the
> "Lisp" occurrences in ordinary SPAD code.
>
> Ralf
>

Maybe we can put it into a package named LispMacrosPackage.

BTW, trappedSpadEval is already used in spad files, you can search for
its usage.

- Qian

[Waldek Hebisch]

'trappedSpadEval' is a macro, this makes a difference in use.
Wrapping it in a function would defeat the purpose: it is
suppossed to catch errors in evalutation of its argument.
Functions get evaluated arguments, so any possible error
in evalutation happens before function can do anything about it.

Currently we can not import/export Spad macros and using Spad macro
as a wrapper would not help with main drawback, that is lack
of typechecking.

One could do version which takes a function as an argument, but
such version is less convenient to use.

--
Waldek Hebisch

[Qian Yun]

So in general you agree with this approach (before more
complete method is implemented)?

The right place to catch errors and try again is
"algcurve(d, selectIntegers first la, first la)"
"ratcurve(d, selIntegers())", or we need to refactor
and catch errors at a deeper level?

The following integrals take the "algcurve" branch.
integrate((1-2^(1/3)*x)/(2^(2/3)+x)/(x^3+1)^(1/2),x)
integrate(1/(2^(2/3)-x)/(x^3-1)^(1/2),x)
integrate(1/(2^(2/3)+x)/(-x^3-1)^(1/2),x)
integrate((2+3*x)/(2^(2/3)-x)/(x^3-1)^(1/2),x)
integrate((1+x+3^(1/2))/(d*x+c)/(-x^3-1)^(1/2),x)
integrate((1+x-3^(1/2))/(d*x+c)/(-x^3-1)^(1/2),x)

This one takes the "ratcurve" branch, but still fails
after repeated tries:
integrate((x^6+x^2+1)^(1/2)*(2*x^6-1)/(x^6+1)/(2*x^6-x^2+2),x)

- Qian

[Waldek Hebisch]

On Sat, Jan 27, 2024 at 10:39:33AM +0800, Qian Yun wrote:
> So in general you agree with this approach (before more
> complete method is implemented)?

Yes.

> The right place to catch errors and try again is
> "algcurve(d, selectIntegers first la, first la)"
> "ratcurve(d, selIntegers())", or we need to refactor
> and catch errors at a deeper level?

I think that some refactor is needed. Look at
FunctionSpaceReduce. Reductions go as three nested loops,
and we need to catch failure at outer loop. There is
a buch of similarly named functions, so finding good place
needs some effort.

Extra thing: there is first step which substitutes integers,
it can fail but failures are quite rare. We should catch
errors here and try different values. There is second
stage which reduces modulo a prime, most failures happen
there. As I wrote, there is a little subtlety that in
case of no parameters second try with the same prime will
do the same caldulations, so retry needs some way to
detect case of no parameters.

> The following integrals take the "algcurve" branch.
> integrate((1-2^(1/3)*x)/(2^(2/3)+x)/(x^3+1)^(1/2),x)
> integrate(1/(2^(2/3)-x)/(x^3-1)^(1/2),x)
> integrate(1/(2^(2/3)+x)/(-x^3-1)^(1/2),x)
> integrate((2+3*x)/(2^(2/3)-x)/(x^3-1)^(1/2),x)
> integrate((1+x+3^(1/2))/(d*x+c)/(-x^3-1)^(1/2),x)
> integrate((1+x-3^(1/2))/(d*x+c)/(-x^3-1)^(1/2),x)
>
> This one takes the "ratcurve" branch, but still fails
> after repeated tries:
> integrate((x^6+x^2+1)^(1/2)*(2*x^6-1)/(x^6+1)/(2*x^6-x^2+2),x)

No parameters there, so each try is doing the same thing.

--
Waldek Hebisch