someone wrote:
>
> > > 2.) Which fix do we use for the wrong values of "Shi" and "Chi"?
> > > I think you wanted to try another solution that what I suggested
> > > first. What is the state of this?
> >
> > AFAICS the "wrong" examples posted on the list are actually OK,
> > that is difference between true and computed result is within
> > error bounds (see my last message in the thread). We should
> > increase accuracy of those functions because for some other
> > arguments error is slightly larger than promised.
>
>
> Maybe one should interface to other libraries for numerics.
> AFAIK Axiom can use GMP and related libraries for arbitrary
> precision computations.
>
> There is a really nice new library for ball arithmetics:
>
> http://fredrikj.net/arb/
> http://fredrikj.net/blog/
>
> The most interesting point (beside being one of the fastest
> libraries of its kind) is the rigorous error control.
>
Well, at first glance reusing work looks attractive. However,
when I looked at existing codes I have found number of implemented
functions is limited (for example I do not know of any free code
outside of FriCAS implementing arbitrary precision Weierstrass
elliptic functions). Some codes are numerically unstable
(elliptic functions in Reduce and at least older Maxima
elliptic functions). There are restrictions on arguments:
I have found _no_ free implementation of incoplete Gamma function
for complex arguments. Sometimes code only works for very
limited range of arguments. So even _selecting_ functions to
link requires nontrivial work. Then there is problem of
actual interfacing which typically requires data convertions,
and consequently limits performance at low precisions.
Once we have external parts complexity of the whole system
grows quite a lot: there are multiple languages and different
build systems. Typically there is serious code bloat.
As a somewhat extreme example let me mention that few days
ago I built FriCAS on Gentoo. I wanted to do this in
"Gentoo" way so I installed Gentoo 'noweb' package.
'noweb' is tiny, but it pulled more than 500 Mb of other
software. So I prefer to be selective in what we link
to. Basically we want establised standalone packages.
Currently we use GMP. I plan to use BLAS and Lapack,
however ATM main blocker is limited availablity of
BLAS and Lapack and potential build problems when
they are not available as standard component.
Next thing probably would be fplll.
Concering rigorous precision, for Ci in Wikipedia notation
we have
Ci(x) = -Cin(x) + gamma + log(x)
we compute all term (including Cin(x)) with rigorous error
bounds. But then there is cancellation between Cin and the
other terms. So absolute precision is as expected, but
relative precision is worse. ATM I do not know how much
effort spend on getting better precision. In principle
we could detect cancellation and repeat computation at
higher precision. But cancellation is a problem only
when Ci(x) is small (basically close to zeros). If
we add such small number to someting else the precision
will be essentially wasted.
BTW, book "Precise Numerical Methods" by Aberth indicates
how to compute large class of functions with rigorous
error bound. It would be not hard to implement such
method. However, I have doubts about performance. For
simple cases (and Ci almost surely is such a case)
performance at moderate precision will be OK. But
I expect that in "interesting" case performance will
be so bad that they will be impossible to do in
reasonable time. In other words getting rigorous
error bounds tends to be expensive, either in
computer time or in programmer time. Interval
style arithmetic typically has significant overhead
maintaining error bounds and obtained bounds tends
to be pessimistic. Conseqently, interval computation
usually are much slower than non-interval multiple
precision computations (which in turn are much slower
than machine floats). AFAICS one way to speed up
interval computations is to have bigger blocks
with rigorous error estimates (inside block
computations would be done using non-interval
arithmetic). Of course this requires effort.
FriCAS functions are step in such direction, for
functions that I implemented I did error analysis
and most functions should deliver results accurate
up to designed precision. Of course there are bugs
(for example in earlier version of Gamma my error
estimates were too optimistic). And there are
deliberate shotcuts (for example Beta and Ci).
>
> > So ATM, what you do looks at attempt to optimize
> > cases when we need sort + parity of permutation.
> > More precisely, it looks like premature optimization
>
> I only need sort + parity, right. This is the best
> and easiest way to compute that. At least I have
> no better idea to come up with.
As I wrote parity is linear. Good sorting algoritms
perform of order n*log(n) operations. So it looks
better to compute parity separately then have
per operation overhead in sort.
More precisely, if you have permutation of integers
between 1 and n then there is easy linear algorithm
to compute parity: follow cycles and multiply parities
of cycles. If you have general sequence sort
pairs (x, i) where i increases from 1 to n. Then
extract second components and compute their parity.
>
> For the moment, I'll factor out the new code
> into an own package. This way VectorAlgebra can
> depend on it but it need not to be part of Fricas
> and we could review this point later.
>
> For VectorAlgebra any way to compute sort + parity
> easily would be enough. (Maybe I should review the
> proposals on wrapping around SortPackage.)
>
>
> > Concerning merge sort: the default sort is merge
> > sort, so we already have one implementation.
>
> Sure, but not a counting one.
>
>
> > BTW: AFAICS SortPackage is not used by any other part
> > of algebra. It is included because it is
> > used as an example in chapter 11.7 of the Axiom
> > book.
>
> Well, this might be true but I think that Axiom
> definitely should include generic sorting algorithms.
> So SortPackage is a valid component in my view.
Merge sort _is_ generic. If you want something which
works on arrays in place I have on my disk heapsort
(non-generic but could be easily genericized by changing
type declarations). Both perform less comparisons than
bubble sort. So, there are really no reasons to
use bubble sort.
Note that in actual use lists are quite frequent. Using
in place sort on lists is going to have very bad performance.
OTOH merge sort on arrays leads to modarate overhead due to
copies and memory allocation. So from performance point of
view merge sort is the only really generic sort.
Let me reiterate my view about performance. Most of the
time we want simple code that works. We want to avoid
obvious performance traps, in particular if there is simple
method with better asymptotic complexity then we use it.
Having running code we can measure which parts are
heavily used and determine bottlenecks. Then we try
to optimize parts that matter. But even in critical
part we have to balance performance and code
complexity.
For me generic code is largely to avoid duplicationg code
(or putting unnatural restrictions on users). Once we
have _one_ generic routine doing something there is no
need for another (except if we find better replacement).
To have better chance of reuse operations are frequently
split into smaller parts. Multiple versions of routine
typically are needed for optimization (and freqently
type specific version is a performance win). But before
we start optimizing we first want to know that operation
matters enough to justify optimization.
--
Waldek Hebisch
[email protected]
--
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 [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/fricas-devel?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
