Developement plans

[Waldek Hebisch]

This is to let you know about my developement plans. With
little details, as details would make it much longer.

1) Numeric special functions. I wrote few month ago about
Bessel functions, there is some progress here, few pieces
are commited, few other wait as they depend on to be
written parts. Eventually this will replace 'sfsfun.boot'.
1a) There is also old work I did on numeric elliptic functions.
While what is commited is quite good, there are some problems
to fix and a few functions are missing.
2) I would like to rework our polynomial system, more preciely
factorization and GCD. I have found scheme which I consider
quite promising. It starts from some old ideas of von zur
Gathen and Zippel with some newer additions. The exact
combination is probably new. I hope that this will give
quite good perfomance. Drawback is that it is substantial
work to implement and speed strongly depends on speed of
low level routines.
2a) There is new univariate factorizer over finite fields. Recently
commited version in principle can handle all finite
fields. There is some benchmarking to do, but there
is good chance that 'moddfact.spad' is reduntant now
(new factorizer can factorizations done by 'moddfact.spad',
but benchmaring is need to verify check for possible
slowdowns). Also, for general finite fields new
factorizer should be more or less as good as 'ddfact.spad'.
But in case of algebraic extentions there is faster method
which I want to add. And there is some work needed on
low level routines for important special cases.
Actually, some low level routines should be usable
in general multivariate code, so I will probaly finish
this before doing serious work on multivariate case.
3) I would like to reorganize our expression machinery.
Unfortunately my original idea is blocked by bugs
in Spad compiler. So so I am looking at alternative
plans.
3a) Change in expression handling is needed for integrator,
it should improve speed and make solution to some
problems easier. In particular with changes expression
structure it would be easier to solve problems caused
by 'real'
3b) Our limit code needs improvement. Currenly mrv_limit
uses Gruntz method, which has some performance problems
and does not handle all functions. I think that I can
improve speed and handle more functions (in particular
all correctly handled by old limit). But this needs
reorganization which would be easier if changes to
expressions machinery are done first.
3c) I would like to improve our ODE solver and transcendental
equation solver. For nonlinear ODE-s I have proof
of concept code which handles well some popular classes
of equations. But it needs work to finish. In case
of transcendental equations we probably could easily
improve equations solver quite a lot, but still,
there is some work to do.
3d) As you probably noticed I spent a lot of effort on
integrator. I plan several changes, but I also
want to improve other parts of FriCAS.
4) Spad compiler and interpreter. I have long term plans for
both, but major improvement to Spad compiler can not be
done quickly and bigger changes to interpreter should be
coordinated with complier. In short term I limit work
on compiler and interpreter to bugs fixes, instead spending
most time on points above.

BTW: In last severel years various things caused me to spent
less time on FriCAS that I would like to. I hope for best,
but just when I hoped that things will go better we got
Covid...

--
Waldek Hebisch

[Qian Yun]

I'm very interested in 3), expression handling is core to
the design of CAS. I wonder what's your original idea
and blocked by what compiler bugs?

I'd like to see improvements on expression, but I have
no ideas on where to start...

- Qian

[Waldek Hebisch]

On Sun, Aug 28, 2022 at 11:42:10AM +0800, Qian Yun wrote:
> I'm very interested in 3), expression handling is core to
> the design of CAS. I wonder what's your original idea
> and blocked by what compiler bugs?
>
> I'd like to see improvements on expression, but I have
> no ideas on where to start...

One thing to to have "per expression" kernels. More
precisely, currently Expression serves as a kind of
universal differential field inside which we embed
other fields. This is good as user interface, because
user are spared burden of creating each needed
differential field. But integration, limits and
solvers really should work with differenatial
fields. And they need to build new fields.
Differential fields need more information than
provided in current kernels. For this reason
and for efficiency I would like to avoid current
kernel and use a new one. Basically:
- I want more info in kernels to avoid some searches
and simplify coding
- when creating kernels I want to deal only with
given differential field (and not with all
kernels that happen to live in kernel cache).

That was motivation for KernelCategory, instead
of Kernel one should be able to use any domain
that offers kernel interface. But currently
it breaks when I try to generilize ExpressionSpace
and FunctionSpace, my current try leads to compile
failure later. This is problematic, as I would like
to reuse several packages that currently require
FunctionSpace (so de facto ATM are limited to
Expression).

Another things is "factored expressions". Currently
expressions are rational functions of kernels, normalised
so that algebraic kernels are reduced by their defining
equations. Rational functions are fully expanded.
I would like variant where rational functions are
partially factored. If you look at Rubi testsuite
many epressions are already in factored form.
When input as Expression they get expanded.
That in many cases leads to much bigger expressions
later. At the input in many cases we could restore
factored form and take advantage of known factors
during computation. Of course, to recognize zero
we need to expand. But I feel that there will be
net gain.

I wrote above about differential fields. EFSTRUC
currently helps dealing with fields of elementary
functions (with some extentions to Liovillian
functions). I would like similar functionality
but doing more and covering more functions.
Basically to produce all data needed to convert
given set of expressions into differential
field generated by them and back (plus some
other info to allow more efficient queries
about structure of the field).

--
Waldek Hebisch

[Tim Daly]

Despite past issues, I have to say I am impressed with your

ability to improve the system.

Keep up the excellent work.

Tim

[Ralf Hemmecke]

>> I'd like to see improvements on expression, but I have
>> no ideas on where to start...
>
> One thing to to have "per expression" kernels. More
> precisely, currently Expression serves as a kind of
> universal differential field inside which we embed
> other fields. This is good as user interface, because
> user are spared burden of creating each needed
> differential field. But integration, limits and
> solvers really should work with differenatial
> fields.

Wouldn't it make sense to use Expression just as user interface and when
it comes to "integrate", take the input expression and construct a
proper differential field out of if. As far as I understand the
transcendental part of the Risch algorithm, one works with differential
field s of the form K(t1)(t2)...(tn) where tk is an indeterminate such
that D(tk) = rational function in t{k-1}, ..., t1. So at each time
there is a concrete differential field and it is clear which tk are
involve and what these tk model. Isn´t it done this way in FriCAS so
that the tk are represented by the kernels? Or is internally always
Expression(...) employed and the actual differential field must be read
from the expression that is involved?

Ralf

[Waldek Hebisch]

On Sun, Aug 28, 2022 at 09:04:26PM +0200, Ralf Hemmecke wrote:
> >>I'd like to see improvements on expression, but I have
> >>no ideas on where to start...
> >
> >One thing to to have "per expression" kernels. More
> >precisely, currently Expression serves as a kind of
> >universal differential field inside which we embed
> >other fields. This is good as user interface, because
> >user are spared burden of creating each needed
> >differential field. But integration, limits and
> >solvers really should work with differenatial
> >fields.
>
> Wouldn't it make sense to use Expression just as user interface and when it
> comes to "integrate", take the input expression and construct a proper
> differential field out of if.

My plans goes in this direction. But I am not sure if we agree what
"proper" means.

> As far as I understand the transcendental part
> of the Risch algorithm, one works with differential field s of the form
> K(t1)(t2)...(tn) where tk is an indeterminate such that D(tk) = rational
> function in t{k-1}, ..., t1.

Well, that is for elementary functions. In the recent example
we have Q(x, log(x), log(x+1)) and role of log(x) and log(x+1)
is really symmetric. When dealing with polylogs having to
put one of kernels at the top is inconvenient. And even
in case of elementary functions it is useful to re-arrage
tower in cases when this is possible.

> So at each time
> there is a concrete differential field and it is clear which tk are involve
> and what these tk model. Isn´t it done this way in FriCAS so that the tk are
> represented by the kernels? Or is internally always Expression(...) employed
> and the actual differential field must be read from the expression that is
> involved?

We have mostly expression (algebraic integrator has special type for
curves). In algebraic code there is field of coefficients, main
transcendental and algebraic generator. In legacy code (coming from
Axiom) at each stage structure is read from expression, more
precisely we take field generated by kernels containd in the
expression. This has led to a bug: when recursively calling
integrator we need to find integral in smaller field. But
function to integrate missed some kernels, so Axiom looked
for integral in too small field and falsely declared whole
function as unintegrable. intpar.spad uses different technique:
differential fields are represented by list of kernels which
generate them, so is free from such problems.

I think that using list of kernels is resonable, but I want
add more info and this would be simpler with something
different than Kernel.


--
Waldek Hebisch

[Qian Yun]

On 8/28/22 21:51, Waldek Hebisch wrote:
> Another things is "factored expressions". Currently
> expressions are rational functions of kernels, normalised
> so that algebraic kernels are reduced by their defining
> equations. Rational functions are fully expanded.
> I would like variant where rational functions are
> partially factored. If you look at Rubi testsuite
> many epressions are already in factored form.
> When input as Expression they get expanded.
> That in many cases leads to much bigger expressions
> later. At the input in many cases we could restore
> factored form and take advantage of known factors
> during computation. Of course, to recognize zero
> we need to expand. But I feel that there will be
> net gain.
>

I have thought about "factored expressions" as well.
My rough idea is using "paren/box" to prevent expand,
and adding checks for "paren/box" to expand when needed.
Do you think this is a viable solution?

Similar to "factored expressions", I think pattern matching,
integration by part, heuristic algorithm also relies on
"do not auto expand expression beyond recognition."

- Qian

[Waldek Hebisch]

On Mon, Aug 29, 2022 at 09:45:20AM +0800, Qian Yun wrote:
>
> On 8/28/22 21:51, Waldek Hebisch wrote:
> >Another things is "factored expressions". Currently
> >expressions are rational functions of kernels, normalised
> >so that algebraic kernels are reduced by their defining
> >equations. Rational functions are fully expanded.
> >I would like variant where rational functions are
> >partially factored. If you look at Rubi testsuite
> >many epressions are already in factored form.
> >When input as Expression they get expanded.
> >That in many cases leads to much bigger expressions
> >later. At the input in many cases we could restore
> >factored form and take advantage of known factors
> >during computation. Of course, to recognize zero
> >we need to expand. But I feel that there will be
> >net gain.
> >
>
> I have thought about "factored expressions" as well.
> My rough idea is using "paren/box" to prevent expand,
> and adding checks for "paren/box" to expand when needed.
> Do you think this is a viable solution?

Well, peren and box allow keeping factors/terms when
we use Expression as representation. I am not sure
if it would work well for more extensive computations.
One trouble is that our zero test depends on expansion

> Similar to "factored expressions", I think pattern matching,
> integration by part, heuristic algorithm also relies on
> "do not auto expand expression beyond recognition."

Shortest description of Risch algorithm is "integrate by
part when this makes integral simpler, otherwise use
special methods". It depends on having expression
in specific form and keeping it in this form during
computation. So "smart" automatic transformations
are harmful for wide range of computation.

OTOH Risch is insensitve to expansion (except for
cost of computations) while pattern matching
is very sensitive to syntactic form. I must
admit that after experience in FriCAS I am
very sceptical about pattern matching. Namely,
when one wants reliabilty then other methods
usually were simpler.

--
Waldek Hebisch

[Qian Yun]

I wonder if there are simpler issues that can be separated from
this grand plan, so that me and others can help with that.

- Qian

On 8/28/22 09:46, Waldek Hebisch wrote:

[Waldek Hebisch]

On Tue, Aug 30, 2022 at 06:01:08PM +0800, Qian Yun wrote:
> I wonder if there are simpler issues that can be separated from
> this grand plan, so that me and others can help with that.

Let me first explain what takes most time. Tim Daly says
that Axiom codebase represents more than 140 man years of
work. In original source there is about 350 thousends of
wc lines in code files (I do no write LOC, because LOC
count should be done removing comments). Resonable coder
should be able to write about 100 lines daily, really
fast one maybe 2-3 times more. Assuming 200 working
days in a year we get 17.5 man years of coding. So
why so much time to do this. This 100 lines already
includes debugging. One thing that take time is documentation.
But main factor is research. Developers must invent
new methods. To do this they must do experiments and
measurements. Partially written code may be discarded
because there is better alternative.

To illustrate this let me mention few examples:
- in case of ffact.spad IIRC corrctly I coded it better than
100 lines per day. Before I started coding theory was
mostly worked out so I know exactly how to go.
- in case cyclo.spad also coding was fast thank to theory
in Arnold and Monagan paper (in fact googling for hints
and reading literature took more time than coding).
- in case of intpar.spad coding also was resonably fast,
but before coding was rather long period of analysis
and planning.
- in case of rdeefx.spad there was quite long period when
I just looked at problem from theoretical point of view.
Coding was slowed down because at first I got too
ambitious and tried too complicated methods for some
substeps (now rdeefx.spad handles main cases but still
need work to cover some important special cases).

OK, what needs to be done:
1. we need packages with interface like ModularAlgebraicGcdTools2
and ModularAlgebraicGcdTools3 but faster. Expected use case
for both is dense, so polynomials should be kept in arrays.
Already a single routine, like polynomial multiplications
would be good step forward. Once we figure out how to do
fast multiplications other routines are likely to be similar
or can use multiplication routine.
2. we need package like ModularFactorizationTools1, but
for algebraic extention fields. Main operation, that is
multiplication is the same as in ModularAlgebraicGcdTools2.
In fact, the two tasks are closely related.
3. in case of multiple algebraic extentions it is not clear
if doing arithmetic directly can be really fast, all
schemes seem give extra power factor in worst case. As
an alternative one could search for primitive element.
In fact, finding primitve elements should be easy, but
there is cost of changing representation.
4. For factoring and GCD we need fast Hensel lifting. Bernardin
in his PhD thesis describe scheme in two variables which
should run resonably well. So probably we should write
routine using his scheme (but some variation may give
improvement).
5. Magma literature says that for moderately large extentions
of small finite field they use Zech logaritms. We could
try similar thing, that is write version of routines
from previous points based on Zech logaritms.
6. ffact.spad currently uses sligtly better algorithm than
ddfact.spad and ffact.spad. Version over Z_p for small p
has significant low level advantage. But generic version
has some disadvantage compared to ddfact.spad and should
be similar to ffact.spad. It would be simpler to just
work with ffact.spad (which is single generic code that
can be specialized to fast verisions) instead of several
packages. But we should do bencharking, to make sure
that there is no unexpected slowdown (and possibly
improve ffact.spad if there is). So interesting cases
are Z_p for large p (so that we can not use machine
sized integers) comparing to ddfact.spad and factorization
over extention fields comparing to ffact.spad.
7. Part of work on numerics is checking accuracy of
routines that we have. As I mentioned, AiryAi had problem
for positive argument, so it is replaced by new version
which still have problems, but should be better than old
one. One issue now it to find out if/where our complex
polygamma works. I have code to compute Bessel functions
and polygamma to resonably high accuracy. This code
is quite slow so not appropriate as general routine.
But if somebody wants to look at accuracy of existing code,
then I can share it and it can serve os source of good
values.

My example with cyclo.spad indicate a class of possiblities:
find promising and well described algorithm in literature
that is not implemented in FriCAS and add it.

On less algorithmic ground: AFAICS current method of computing
Laplace transform and inverse Laplace transform basically
try to match to known cases. Here general hypergeometric
functions help, making part of it more systematic, but
still, this is mostly business of matching to special cases.
It would help if somebody implemented more cases for Laplace
and inverse Laplace transform. And possibly also for
Fourier and Z transform (there is old .input file by
Alasdair McAndrew but I would need convertion to Spad
constructs).

Let me add that usual form of contribution is fixing bugs.
My policy is to fix easily fixable bugs as soon as posible
(known bugs can mask other problems and if left unfixed
can lead to lowering quality), which means that fixes
for known bugs are probably not so easy. OTOH you can
look for unkown bugs, some of them may be easily fixable.
When searching for bugs one possiblity is to repeat
what Vladimir Bondarenko did. Start from large collection
of expression, possibly randomly generated, possibly
random variation of known formula and use them as input
to the system. For indefinite integration things are
easy: derivative of an expression is supposed to be
integrable and up to additive constant should agree with
original expression. For factoring one could multiply
known factors. Getting good scheme requires some
invention, for example my first generator of random
expressions generated too deep expressions, variant
inspired by Facebook folks is much better. But one
can try other ways to bias generated expression towards
some interesting properties.

Let me add that first obstacle in random testing is
to have some criterion for validity of answers. When
testing routines using pseudo-random number simple
approch is to run the same input several times. In
many cases, when answer differs this indicates a bug
(there are cases when output in naturally non-unique
like in case of primitieve elements).

Another thing is fuzzing: one deliberately introduces changes
(probably bugs). Tests that give different result after
change (detect change) are good tests and candidates to
add to testsuite.

--
Waldek Hebisch

[Tim Daly]

If anyone is looking for test suites, this might be useful:

http://axiom-developer.org/axiom-website/CATS/index.html

In the integration tests I tried to use derivatives and various

normalizations to compare the derivatives of each integration

matched the original expression within a constant. Some of

them did not succeed. Some of those failures are likely due

to bugs.

Tim

[Grégory Vanuxem]

Hello Tim,

(trying to find your mail)

Do you happen to know if "FortranFunctionCategory" was working ?

In other words the Fortran interface?

Cheers,

__

Greg

[Tim Daly]

My email is axiomcas at gmail.com

Domains handling Fortran from NAG generally are found using

)apropos Asp

Conversions to Fortran work, e.g

outputAsFortran(3*x) ==> R8=3*x

but the code that invokes NAG routines such as f02aff from NagEigenPackage

compile to a Lisp call 'invokeNagman' which wanders down to a lisp 'nagCall'.

That invokes 'sockSendString' to talk to a NAG co-process.

I don't have a NAG license so I don't have access to NAG libraries to test this.

It would be possible to write such a process and have it attach to a socket

that handled the messages.

Some of the functionality you need may be in BLAS and LAPACK.

https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms

https://en.wikipedia.org/wiki/LAPACK

Axiom has converted the BLAS and LAPACK source code to Common Lisp.

Portions of this have been added to Axiom directly (more to come) so there

is no need for a process/socket.

For example, from BlasLevelOne:

(9) -> t2:Complex DoubleFloat := complex(1.0,1.0)

   (9)  1. + %i
                                                   Type: Complex(DoubleFloat)

( https://netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga53914429b84cb315da483e71e27ed5c1.html )

(10) -> dcabs1(t2)

   (10)  2.
                                                            Type: DoubleFloat

Tim

[Grégory Vanuxem]

Am astothised .  Will respond later

--
You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/f7ae2284-ac19-46a2-80a1-76b845751c27n%40googlegroups.com.

[Tim Daly]

The embedded version of BLAS and LAPACK is delayed as the code

should use a GPU if it is available on your hardware (quite common these days).

So, for instance, SAXPY (a*x[i]+y[i]) should use CUDA if available or default to

the CPU if not (Torch will do this, for example). I just recently got a Google TPU

which will be even faster though much less common. The fun never ends.

[Grégory Vanuxem]

Take into account that the graph in the right pane is animated:
https://docs.juliaplots.org/latest/

Le mer. 28 déc. 2022 à 10:22, Tim Daly <axio...@gmail.com> a écrit :
>
> The embedded version of BLAS and LAPACK is delayed as the code
> should use a GPU if it is available on your hardware (quite common these days).

I broke two 1200$ laptop. With GPU ... One was a present  from my mother. NVIDIA.

No more than 6 months alive. I find myself very silly sometimes. I preferred lesse.

This is in the pipeline for me.

> So, for instance, SAXPY (a*x[i]+y[i]) should use CUDA if available or default to
> the CPU if not (Torch will do this, for example). I just recently got a Google TPU
> which will be even faster though much less common.

"Oh, c'est du lourd". A literally french affective expression that means it's a lot heavy ;)

Later...More...

  > fun never ends..

I _love _ to read that, many thanks