RE: [Caml-list] Efficient and canonical set representation?

RE: [Caml-list] Efficient and canonical set representation?

[Fred Smith]

Ocaml has a structural comparison function (<=) : 'a -> 'a -> bool.  

But the bigger problem with this approach (as Jean-Baptiste Rouquier) kindly pointed out to me is that insertion takes O(n), not O(log(n)) time. Duh.

-Fred


> -----Original Message-----
> From: owner-caml-list@pauillac.inria.fr 
> [mailto:owner-caml-list@pauillac.inria.fr] On Behalf Of Samuel Lacas
> Sent: Friday, November 07, 2003 10:44 AM
> To: caml-list@inria.fr
> Subject: Re: [Caml-list] Efficient and canonical set representation?
> 
> 
> Fred Smith a écrit 2.2K le Fri, Nov 07, 2003 at 10:27:25AM -0500: # 
> # I guess what you're looking for are sorted arrays:
> #   1) O(log n) lookup and insertion via binary search
> #   2) O(n) union and intersection are simple
> #   3) Equal sets are represented by structurally equivalent objects.
> # 
> # -Fred
> 
> Hmm, except that, if I'm not wrong, it was required the 
> structure to hold any kind of object. Sorted arrays require 
> the elements to be sortable. Using the hash of the objects 
> may be an answer ?
> 
> sL
> 
> -------------------
> To unsubscribe, mail caml-list-request@inria.fr Archives: 
http://caml.inria.fr Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/ Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

RE: [Caml-list] Efficient and canonical set representation?

[Diego Olivier Fernandez Pons]

    Bonjour,

> After your remarks and Brian's, I'm starting to wonder if it is
> possible at all to do what I want. Maybe I should be looking for an
> impossibility proof instead...

Patricia sets seem to be what you are looking for.
 (1). Efficient usual operations (lookup, insertion, union)
 (2). Structural equality

Their only problem is that they cannot handle polymorphic orderable
types but only integers...

Hash the data, use this key to insert it in a patricia map and solve
the collisions by chaining in an ordered list (with the polymorphic
[compare] function). (1) and (2) still hold under usual hypothesis on
the rate of collisions.

A few changes to JCF's implementation should be enough.

        Diego Olivier



-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

[Caml-list] Rounding mode

[Christophe Raffalli]

I saw a previous discussion about rounding mode for OCaml (in 2000).

As anyone implemented the functions to change the rounding mode for 
floating point from OCaml (possibly as a patch to add primitives and 
save a C function call) ?

This is necessary to implement interval arithmetic ...

Note: I found the assembly code for 6 architectures on the web:
http://www-anp.lip6.fr/cadna/Rounding_mode_Dir/Accueil.php

 From this, it would be easy to add it in OCaml ?


-- 
Christophe Raffalli
Université de Savoie
Batiment Le Chablais, bureau 21
73376 Le Bourget-du-Lac Cedex

tél: (33) 4 79 75 81 03
fax: (33) 4 79 75 87 42
mail: Christophe.Raffalli@univ-savoie.fr
www: http://www.lama.univ-savoie.fr/~RAFFALLI
---------------------------------------------
IMPORTANT: this mail is signed using PGP/MIME
At least Enigmail/Mozilla, mutt or evolution
can check this signature
---------------------------------------------

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Eric Dahlman]

Christophe Raffalli wrote:

>
> I saw a previous discussion about rounding mode for OCaml (in 2000).
>
> As anyone implemented the functions to change the rounding mode for 
> floating point from OCaml (possibly as a patch to add primitives and 
> save a C function call) ?
>
> This is necessary to implement interval arithmetic ...

Somewhat off topic but why is this necessary from a numerical math type 
of perspective.  I am honestly curious as I don't see how this would 
interact with the calculation in a meaningful way.

-Eric

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Brian Hurt]

On Mon, 10 Nov 2003, Eric Dahlman wrote:

> Christophe Raffalli wrote:
> 
> >
> > I saw a previous discussion about rounding mode for OCaml (in 2000).
> >
> > As anyone implemented the functions to change the rounding mode for 
> > floating point from OCaml (possibly as a patch to add primitives and 
> > save a C function call) ?
> >
> > This is necessary to implement interval arithmetic ...
> 
> Somewhat off topic but why is this necessary from a numerical math type 
> of perspective.  I am honestly curious as I don't see how this would 
> interact with the calculation in a meaningful way.
>

Simplistic discussion of why changing the rounding modes is usefull for 
numeric programming.  If you know numerical analysis, stop reading now 
:-).

FP numbers are only precise to so many digits.  So every time you do an 
operation, a small amount of error is introduced- rather than being the 
exact answer, it's the closest number to the exact answer that FP numbers 
can represent.  This is why the following code does not return the answer 
you expect:

let unexpected () =
    let loop i n =
        if (i < 10) then
            loop (i + 1) (n +. 0.1)
        else
            n
    in
    loop 0 0.
;;

1/10 in binary is infinitely repeating the way 1/3 is in decimal.  So the 
expression n +. 0.1 isn't adding *exactly* 0.1, but instead 0.1 - e for 
some small e.  Where e is ~1E-16.  But this is large enough that 
unexpected() returns  0.999999999999999889 and not 1.0.

Note also that this is nothing to do with Ocaml- C, Java, Perl, Python, 
FORTRAN, and every other language out there that does floating point has 
this problem- it's a problem with floating point.

Now, for a single operation this error is negligable- one the order of
1e-15 (assuming the answer is 1.0).  But over large numbers of
calculations, this error can accumulate, and become large enough to affect
the results.

Consider, for example, multiplying two numbers- x and y.  But you don't 
have x and y exactly, you have (x + e) and (y + e) for some error amount 
e.  You multiply them, and you get (x + e)*(y + e) = x*y + x*e + y*e + 
e*e.  Now, e*e is vanishingly small, we can ignore it.  And, assuming x 
and y are both about 1.0, our answer is x*y + ~2*e- in other words, we've 
about doubled the amount of error.  Even starting at 1e-15, doubling it 
every operation and you eat up all your precision surprisingly fast.

How can we tell how much error we have?  Well, you can either do a whole
bunch of graduate-level mathematics, or you can run the same algorithm
with different rounding modes.  If you always round towards positive
infinity, you'll always be adding an error term- i.e. it'll be (x + e) and
(y + e) with a result of x*y + 2e.  If you always round to negative
infinity, you'll be subtracting the error term- (x - e) and (y - e) with a 
result x*Y - 2e.  Comparing the results gives you both a good idea of 
where the answer really is, and how big your errors are.

Note, all is not necessarily lost.  If instead we have (x + e) and (y - e)  
a little bit of quick algebra shows that the result is actually *more*
accurate than either of the original inputs.  The errors tend to cancel
instead of accumulate.  Numerical analysts call algorithms whose errors
tend to cancel "stable", while algorithsm whose errors tend to accumulate
are "unstable".  A good textbook on numerical methods will generally only 
discuss stable algorithms- use 'em.

The effect can also be sort of acheived by using a random rounding method.  
For each operation, the CPU picks one of it's rounding methods at random 
to use (well, pseudo-random anyways).  Now, while this will in general 
increase the accuracy of floating point programs, it also means that every 
time you run a numerical program, even over the same data, you will get a 
slightly different answer.  Newbies wondering why their program keeps 
changing it's answer will probably rival newbies wondering why their stack 
overflows pretty quickly.  Which is why I'm hesitant to advocate doing 
this to all programs.  I still would like to be able to select the 
rounding mode- either at compile time or at run time, which ever is 
easier.  Actually, run time is probably easier to implement and more 
usefull (I can switch rounding modes for different parts of the program).
 
-- 
"Usenet is like a herd of performing elephants with diarrhea -- massive,
difficult to redirect, awe-inspiring, entertaining, and a source of
mind-boggling amounts of excrement when you least expect it."
                                - Gene Spafford 
Brian


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Eric Dahlman]

Brian Hurt wrote:

>On Mon, 10 Nov 2003, Eric Dahlman wrote:
>
>  
>
>>Christophe Raffalli wrote:
>>
>>    
>>
>>>I saw a previous discussion about rounding mode for OCaml (in 2000).
>>>
>>>As anyone implemented the functions to change the rounding mode for 
>>>floating point from OCaml (possibly as a patch to add primitives and 
>>>save a C function call) ?
>>>
>>>This is necessary to implement interval arithmetic ...
>>>      
>>>
>>Somewhat off topic but why is this necessary from a numerical math type 
>>of perspective.  I am honestly curious as I don't see how this would 
>>interact with the calculation in a meaningful way.
>>
>>    
>>
>
>Simplistic discussion of why changing the rounding modes is usefull for 
>numeric programming.  If you know numerical analysis, stop reading now 
>:-).
>  
>
I read it anyway ;-)

[Snip a description of how floating point works.]

>Consider, for example, multiplying two numbers- x and y.  But you don't 
>have x and y exactly, you have (x + e) and (y + e) for some error amount 
>e.  You multiply them, and you get (x + e)*(y + e) = x*y + x*e + y*e + 
>e*e.  Now, e*e is vanishingly small, we can ignore it.  And, assuming x 
>and y are both about 1.0, our answer is x*y + ~2*e- in other words, we've 
>about doubled the amount of error.  Even starting at 1e-15, doubling it 
>every operation and you eat up all your precision surprisingly fast.
>  
>
>How can we tell how much error we have?  Well, you can either do a whole
>bunch of graduate-level mathematics, or you can run the same algorithm
>with different rounding modes.  If you always round towards positive
>infinity, you'll always be adding an error term- i.e. it'll be (x + e) and
>(y + e) with a result of x*y + 2e.  If you always round to negative
>infinity, you'll be subtracting the error term- (x - e) and (y - e) with a 
>result x*Y - 2e.  Comparing the results gives you both a good idea of 
>where the answer really is, and how big your errors are.
>  
>
Now I don't know that it is this simple.  I understand that by always 
rounding in one direction you are in effect  calculating a bound on the 
error in that direction but that will only work in the most simplistic 
of cases.  For instance if half of your values had their sign changed 
before being summed the sign of half  the associated error terms would 
also change.  The result would not give you any information about the 
effective lower/upper bound you thought you were calculating.  It seems 
to me that this approach is only applicable to problems which would be 
trivially solvable with graduate level mathematics and of absolutely no 
value in any interesting cases.

>Note, all is not necessarily lost.  If instead we have (x + e) and (y - e)  
>a little bit of quick algebra shows that the result is actually *more*
>accurate than either of the original inputs.  The errors tend to cancel
>instead of accumulate.  Numerical analysts call algorithms whose errors
>tend to cancel "stable", while algorithsm whose errors tend to accumulate
>are "unstable".  A good textbook on numerical methods will generally only 
>discuss stable algorithms- use 'em.
>  
>
This is incorrect because you are thinking about real numbers while 
working with floating point numbers, this is especially dangerous when 
you are operating near or below the epsilon for the particular floating 
point representation like we are now.  If you were to look at the set of 
numbers which are precisely represented by a given floating point 
representation you would see that it is concentrated near 0 and that it 
tends to thin out as you move away from there.   The distance  between 
any  two adjacent floating point numbers is highly variable and not 
uniform or "nice" by almost any measure.  I really am having a hard time 
figuring out what you can say about the tiny e terms in this case as we 
move between regions with differing "floating point resolutions".  

Just to illustrate this the different "resolutions" the smallest 
positive double float which can be represented is  2.225 x 10^(-308) 
while the epsilon for doubles in general is 1.11 x 10^(-16) and can be 
thought of at the coarsest resolution.  That is a big difference and is 
because floats are packed in tightly around 0.

As for a numerically stable algorithm being able to produce a result 
which is more accurate than the original inputs that is total hogwash.  
You cannot ever have more accuracy than was present in your original 
input, ever. Nihil ex nihilo.

>The effect can also be sort of acheived by using a random rounding method.  
>For each operation, the CPU picks one of it's rounding methods at random 
>to use (well, pseudo-random anyways).  Now, while this will in general 
>increase the accuracy of floating point programs, it also means that every 
>time you run a numerical program, even over the same data, you will get a 
>slightly different answer.  Newbies wondering why their program keeps 
>changing it's answer will probably rival newbies wondering why their stack 
>overflows pretty quickly.  Which is why I'm hesitant to advocate doing 
>this to all programs.  I still would like to be able to select the 
>rounding mode- either at compile time or at run time, which ever is 
>easier.  Actually, run time is probably easier to implement and more 
>usefull (I can switch rounding modes for different parts of the program). 
>  
>
As an oldby I would be disconcerted by this as well.  All you are doing 
is trading one sort of  numerical noise for another which you guess 
should be better behaved but in the absence of a neat proof to the 
contrary I don't believe is.  There are just *SO* many sources of huger 
error terms in numerical calculations that the tiny rounding errors are 
really lost in the noise.  I even tried to come up with a case where it 
would matter but that got so farcical that I gave up on it.

What I was really wondering about when I first asked this question was 
if the other person had exactly that kind of example where they were 
using some clever trick to usefully misapply floating point calculations 
in their problem domain.

-Eric

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Brian Hurt]

On Mon, 10 Nov 2003, Eric Dahlman wrote:

> Brian Hurt wrote:
> 
> >
> >Simplistic discussion of why changing the rounding modes is usefull for 
> >numeric programming.  If you know numerical analysis, stop reading now 
> >:-).
> >  
> >
> I read it anyway ;-)

The reason I asked you to stop reading was that I was going to be 
massively simplifying a subject that is in effect it's own branch of PhD 
level mathematics.  The result will be somewhat like those books on 
relativity aimed at grade school kids.  Physicists who read them are 
regularly horrified that important subtlities are glossed over, hand waved 
past, or simply blatantly ignored.  The reason is, of course, I was 
writting an email not a book :-).

Also, I'm far from an expert in this field myself (although I know a 
couple).  I know just enough to know how much I don't know.

[Snip happens]

> Now I don't know that it is this simple.  

You're right, it's not at all that simple.  For any given two rounding 
modes, it's possible to construct algorithms for which both rounding modes 
will error in the same direction- i.e. the two results will not bracket 
the "real" answer, and will also misdirect you as to how large your error 
bounds really are.  Also, doing this trick only tells you about the data 
you are using.  Algorithms can have unexpected failure modes where they 
lose a lot more presision.  A classic example of this is subtraction.  You 
may know that x = 1.00000000000002 to 15 digits of precision, and y = 
1.00000000000001 to 15 digits of precision, but you only know x-y = 
0.00000000000001 to 1 digit of precision.  Most of your precision just got 
canceled.  Opps.

Interval arithmetic has similiar problems.  There are no mechanistic 
replacements for numerically analyzing your algorithms.  Both are usefull 
as experimental evidence that your analysis is correct.  But numerically 
analyzing your algorithms is a subtle and tricky- which is why my advice 
is to leave to the numerical analysts, and only use algorithms they have 
blessed (if possible).

Hmm.  Actually, I wonder if type systems could be extended to detect 
numerical instabilities?

> This is incorrect because you are thinking about real numbers while 
> working with floating point numbers, this is especially dangerous when 
> you are operating near or below the epsilon for the particular floating 
> point representation like we are now.  If you were to look at the set of 
> numbers which are precisely represented by a given floating point 
> representation you would see that it is concentrated near 0 and that it 
> tends to thin out as you move away from there.   The distance  between 
> any  two adjacent floating point numbers is highly variable and not 
> uniform or "nice" by almost any measure.  I really am having a hard time 
> figuring out what you can say about the tiny e terms in this case as we 
> move between regions with differing "floating point resolutions".  

Most of the numerical proofs I have seen assume error is a scale factor,
not an additive factor.  I.e., it's x * (1 + e), not x + e.  Then x*y with
error terms becomes x*(1+e_x)*y*(1+e_y) = x*y*(1 + e_x + e_y + e_x*e_y).  
This way I sidestep any assumptions about the "sizes" of x and y.  You
still have to deal with the corner cases where x and y become too small,
and you start hitting denormalized numbers.  This is more accurate, but
more complex math.  Thus the simplification to (x + e), while making sure
the conclusions are still more or less correct.

> As for a numerically stable algorithm being able to produce a result 
> which is more accurate than the original inputs that is total hogwash.  
> You cannot ever have more accuracy than was present in your original 
> input, ever. Nihil ex nihilo.

Yep.  But you can have less, even none what so ever.  Stable algorithms
are stable precisely because the errors they introduce tend to cancel each
other out instead of accumulating.  Which was the point I was trying to
make- I'll admit I was speaking inaccurately when I implied that stable
algorithms somehow magically added precision that wasn't there already.

You can make both type I and type II errors (warping statistic
terminology), and both are incorrect.  Type I errors would here be
assuming that because the computer printed out 15 digits of output, all 15
digits are to be beleived (despite the fact that the inputs are only
accurate to 3 digits).  Type II errors are assuming that errors only
accumulate.  This leads to demands for, for example, 2048-bit floating
point types.  The hope here is that if you can make the errors small
enough, you may still have some precision when you're done.

In my experience, type II errors are most often made by programmers who 
have recently "discovered" that floating point introduces errors into 
calculations.  Which is why I was addressing it early- the solution isn't 
to throw more precision at the problem, the solution is to make the errors 
introduced cancel each other.  And in so doing, unwittingly played to 
making type I errors.

[ snip: randomized rounding as the default ]

> As an oldby I would be disconcerted by this as well.  All you are doing
> is trading one sort of numerical noise for another which you guess
> should be better behaved but in the absence of a neat proof to the
> contrary I don't believe is.  There are just *SO* many sources of huger
> error terms in numerical calculations that the tiny rounding errors are
> really lost in the noise.  I even tried to come up with a case where it
> would matter but that got so farcical that I gave up on it.

All stable algorithms I know of remain stable in the face of randomized 
rounding.  There are some algorithms which are instable with any fixed 
rounding scheme, but which are stable with randomized rounding.

Consider the following case: adding a list of three numbers together-
1.0, 0.5e-15, and 0.5e-15.  Adding them in that order is unstable no 
matter what fixed rounding you use.  The correct answer is, obviously, 
1.000000000000001.  If you always round up, you get 1.000000000000002, 
while if you always round down, you get 1.000000000000000.  With 
randomized rounding, you have a 50% shot of getting the right answer- more 
likely than any other answer.  Now, replace that three-element list with 
1.0 followed by 0.5e-15 repeated a million times.

The stable algorithm here is to sort your long list of integers to add 
first, and always add from smallest to largest.  If you do this with the 
three element list, you'll add 0.5e-15 to itself, getting 1e-15, a number 
large enough not get totally rounded off when added to 1.0.  Adding 
0.5e-15 to itself a million times gets you 0.5e-9, which when added to 1.0 
always gets the right answer, with very little dependence on what rounding 
happens.  Replacing an instable algorithm with a stable algorithm is 
always a good idea.  Of course, in this case, you're replacing an O(n) 
algorithm with an O(n log n) algorithm.

There are reasons why people who know what they're doing may want to play 
with rounding modes.  And while randomized rounding may in some sense be 
"better", I agree that making it the default it a bad idea.  I'm not even 
sure if all architectures offer it.

I've got to admit to a certain fondness for the writtings of Prof. Kahan:
http://www.cs.berkeley.edu/~wkahan/

In fact, I'll go so far as saying the following should be required 
reading for the Ocaml developers:
http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf

(Well, except for operator overloading, I'll disagree with Kahan on that. 
 Although for what Kahan wants, I wonder if Ocaml's ability to define new
infix operators isn't enough, or even better.)

> 
> What I was really wondering about when I first asked this question was 
> if the other person had exactly that kind of example where they were 
> using some clever trick to usefully misapply floating point calculations 
> in their problem domain.

I agree- this is an incredibly bad idea.

-- 
"Usenet is like a herd of performing elephants with diarrhea -- massive,
difficult to redirect, awe-inspiring, entertaining, and a source of
mind-boggling amounts of excrement when you least expect it."
                                - Gene Spafford 
Brian


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Xavier Leroy]

> Hmm.  Actually, I wonder if type systems could be extended to detect 
> numerical instabilities?

It's an extremely tough problem, but you might be interested in the
following papers:

"Static Analyses of the Precision of Floating-Point Operations",
Eric Goubault, Static Analysis Symposium 2001.

"Propagation of Roundoff Errors in Finite Precision Computations: a
Semantics Approach", Matthieu Martel, European Symposium on
Programming 2002.

- Xavier Leroy

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Diego Olivier Fernandez Pons]

    Bonjour,

> Somewhat off topic but why is this necessary from a numerical math
> type of perspective. I am honestly curious as I don't see how this
> would interact with the calculation in a meaningful way.

You are right when you say that there are many sources of errors in
numerical computations and that rounding errors are usually
insignificant with respect to them.

The point is that stochastic arithmetic (and its deterministic variant
interval arithmetic) are useful to find where the accurancy of your
computation is falling drastically (e.g. cancellations)

I really haven't the place to explain extensively how CESTAC works but
there are a few explanations in the ANP website

   http://anp.lip6.fr/cadna/Accueil.php

(CADNA for C/C++ source codes, user's guide. Chapter 4. Survey of the
CESTAC method. Many examples also on the homepages).

The main idea is that in a first order approximation, the number of
significant digits of a result can be estimated with respects to the
dispersion of the different values it can take using several rounding
modes.

Then, you can avoid doing unstable computations like dividing by a
small number (epsilon) very noised which makes you believe it is a
good 'pivot' in a gaussian resolution, etc. The whole computation will
then give a more accurate value.

The website gives an example where usual gauss method finds

   x1 = 60 x2 = - 8.9 x3 = 0.0 and x4 = 1.0

when you estimate the errors, you find

   x1 = 1.0 x2 = 1.0 x3 = 0.1 e-07 and x4 = 1.0

exact values are

   x1 = 1 x2 = 1 x3 = 0.1 e-07 x4 = 1

The difference is only due to a 'bad' pivot succesfully detected and
therefore avoided.


        Diego Olivier


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Eric Dahlman]

Hi Diego,

I looked at that page thanks a bunch!  That was just what I was hoping 
for and it answered many of my questions.

-Eric


Diego Olivier Fernandez Pons wrote:

>     Bonjour,
> 
> 
>>Somewhat off topic but why is this necessary from a numerical math
>>type of perspective. I am honestly curious as I don't see how this
>>would interact with the calculation in a meaningful way.
> 
> 
> You are right when you say that there are many sources of errors in
> numerical computations and that rounding errors are usually
> insignificant with respect to them.
> 
> The point is that stochastic arithmetic (and its deterministic variant
> interval arithmetic) are useful to find where the accurancy of your
> computation is falling drastically (e.g. cancellations)
> 
> I really haven't the place to explain extensively how CESTAC works but
> there are a few explanations in the ANP website
> 
>    http://anp.lip6.fr/cadna/Accueil.php
> 
> (CADNA for C/C++ source codes, user's guide. Chapter 4. Survey of the
> CESTAC method. Many examples also on the homepages).
> 
> The main idea is that in a first order approximation, the number of
> significant digits of a result can be estimated with respects to the
> dispersion of the different values it can take using several rounding
> modes.
> 
> Then, you can avoid doing unstable computations like dividing by a
> small number (epsilon) very noised which makes you believe it is a
> good 'pivot' in a gaussian resolution, etc. The whole computation will
> then give a more accurate value.
> 
> The website gives an example where usual gauss method finds
> 
>    x1 = 60 x2 = - 8.9 x3 = 0.0 and x4 = 1.0
> 
> when you estimate the errors, you find
> 
>    x1 = 1.0 x2 = 1.0 x3 = 0.1 e-07 and x4 = 1.0
> 
> exact values are
> 
>    x1 = 1 x2 = 1 x3 = 0.1 e-07 x4 = 1
> 
> The difference is only due to a 'bad' pivot succesfully detected and
> therefore avoided.
> 
> 
>         Diego Olivier
> 
> 
> -------------------
> To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
> Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
> Beginner's list: http://groups.yahoo.com/group/ocaml_beginners
> 


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Diego Olivier Fernandez Pons]

    Bonjour,

> Actually, I wonder if type systems could be extended to detect
> numerical instabilities ?

I forgot to answer to that question...

There was a discusion on the types forum

http://www.cis.upenn.edu/~bcpierce/types/archives/current/msg01419.html

The main problem could be stated as "what properties do you want to
guarantee for your arithmetics" : algebraic properties ? topological
properties ? subtyping properties ?

Several examples violated algebraic properties were given by Tim
Sweeney in his post (Leibniz equalities, uniqueness of the zero,
forall x. x = x, etc.)

Several examples of topological properties were given by Joe Darcy :

- 'compactification' "IEEE 754 floating-point arithmetics is a
complete system; that is for all inputs, every arithmetic operation in
the standard has a defined IEEE 754 value as a result [...] There are
two ways to add infinite values to the field of real numbers : 1. a
single infinity 2. two affine infinities ..." Actually, there are much
more ways of compactifying a locally compact space and
'compactification' hasn't the same mathematical meaning than
'completion' or 'inner operations' but that is not the problem here.

- funny theorems on limits and series like 'if the general term of a
serie converges, then the whole serie converges'

The subtyping question seems to be related to commutative diagrams :
we want the surjection (e.g. 64 -> 32 bits) to commute with arithmetic
operations.

"Though one can losslessly convert from 32-bit to 64-bit and back,
performing a 32-bit arithmetic operation on such 32-bit numbers can
produce a different result than performing a 64-bit arithmetic
operation on their 64-bit extensions, and then converting back."


We surely can answer to some of these question : with respect to group
properties there is a well known theorem 'any subgroup of R either is
dense or has the form aZ'. In both cases they are infinite. Therefore
you cannot expect any floting points arithmetic to keep all group
properties.

For the typing question, we can give part of an answer :

(I use stochastic arithmetic because interval arithmetic is usually
too pessimistic. I am not a statistician but they all seem to agree
with the fact that in many cases being 'almost sure' is much more easy
than being 'absolutely sure')

Cestac allows you to compute the 'meaningful part' of a computation
(in a statistical sense). Then you can see a result like a couple
float * int where the integers gives the number of exact digits.

If you see a coercion like a surjection (f, n) -> (f, m) with m < n
(which just means that you loose accuracy) then :

(i) We can construct a 'clean' family of subtyping operators that is
such as there is at least one surjection that commutes with your
computations (which one of course depends on the stability of the
given computations).

(ii) We can construct a 'type system' which ensures dynamically that
every result is 'well typed'. Actually, this is exactly what CADNA
does but as a library.


        Diego Olivier


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode

[Christophe Raffalli]

[-- Attachment #1: Type: text/plain, Size: 1338 bytes --]

Eric Dahlman wrote:
> Christophe Raffalli wrote:
> 
>>
>> I saw a previous discussion about rounding mode for OCaml (in 2000).
>>
>> As anyone implemented the functions to change the rounding mode for 
>> floating point from OCaml (possibly as a patch to add primitives and 
>> save a C function call) ?
>>
>> This is necessary to implement interval arithmetic ...
> 
> 
> Somewhat off topic but why is this necessary from a numerical math type 
> of perspective.  I am honestly curious as I don't see how this would 
> interact with the calculation in a meaningful way.
> 
> -Eric
> 

interval arithmetic proceeds by computing not one appriximation of a 
real number x but an interval containing x. To add two intervals, you 
need to add both ends ... but one end with near to plus infinity round 
and the other end with near to minus infinity (if you care about correct 
program :-).

-- 
Christophe Raffalli
Université de Savoie
Batiment Le Chablais, bureau 21
73376 Le Bourget-du-Lac Cedex

tél: (33) 4 79 75 81 03
fax: (33) 4 79 75 87 42
mail: Christophe.Raffalli@univ-savoie.fr
www: http://www.lama.univ-savoie.fr/~RAFFALLI
---------------------------------------------
IMPORTANT: this mail is signed using PGP/MIME
At least Enigmail/Mozilla, mutt or evolution
can check this signature
---------------------------------------------

[-- Attachment #2: Type: application/pgp-signature, Size: 252 bytes --]

Re: [Caml-list] Rounding mode + extended

[Christophe Raffalli]

At this address, you will find a very small library for changing the
rounding mode from OCaml (works on 6 processors, but only tested on i486)

ftp://www.lama.univ-savoie.fr/pub/users/RAFFALLI/OCaml-round.tgz

...

Is it possible to use "extended" instead of "double" in OCaml for float
(I am ready to recompile OCaml for this)

What was the reason to have only one type of "float" ? Three types of
float (float, double and extended) do not explode the complexity of the
optimization for float array or float in record (I think).

Ok, you should have three additions +. +.. and +... :-(

Is this the only reason ?

-- 
Christophe Raffalli
Université de Savoie
Batiment Le Chablais, bureau 21
73376 Le Bourget-du-Lac Cedex

tél: (33) 4 79 75 81 03
fax: (33) 4 79 75 87 42
mail: Christophe.Raffalli@univ-savoie.fr
www: http://www.lama.univ-savoie.fr/~RAFFALLI
---------------------------------------------
IMPORTANT: this mail is signed using PGP/MIME
At least Enigmail/Mozilla, mutt or evolution
can check this signature
---------------------------------------------


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

Re: [Caml-list] Rounding mode + extended

[Xavier Leroy]

> At this address, you will find a very small library for changing the
> rounding mode from OCaml (works on 6 processors, but only tested on i486)

A general comment on the rounding mode issue and your question about
extended-precision floats.

While I agree it would be nice to have full IEEE 754 support in OCaml,
one has to keep in mind that the OCaml standard libraries are severely
constrained by the general requirement that the bytecoded part should
work on any "reasonable" platform, where "reasonable" is currently
defined as "ANSI C plus Unix 98".

Unfortunately, this combination of standards has very little to offer
in terms of advanced features of IEEE float arithmetic.  In
particular, there's no standardized API for controlling rounding
modes.

ISO C9X does provide functions (fegetround, fesetround) to control
rounding modes, and it might be a good idea to use them in your
library when available.  The problem is that very few platforms
implement these functions, let alone full ISO C9X.  For instance, glibc
is the only C library that I know of that provides fe[gs]etround.

So, between standards that few implement and processor-specific asm
statements like you did, there is very little hope of achieving portability.

> Is it possible to use "extended" instead of "double" in OCaml for float
> (I am ready to recompile OCaml for this)

You mean "long double"? :-)  (That's the ISO C9X name for
extended-precision floats.)  

The core bytecode system might keep working if you do this change,
except output_value/input_value on floats, which require 8-byte floats.
As for the native-code compiler, that would require non-trivial
changes throughout.

> What was the reason to have only one type of "float" ? Three types of
> float (float, double and extended) do not explode the complexity of the
> optimization for float array or float in record (I think).

If you want all three array types (float array, double array, extended array)
to be unboxed, it will certainly explode that complexity, turning (in
the polymorphic case) a 2-way branch into a 4-way branch.

> Ok, you should have three additions +. +.. and +... :-(

And much more, notably conversions between the various float types,
conversions to/from integers, to/from strings, printf formats, and so on.
Just look at the Int32, Int64 and Nativeint modules for an idea of the
combinatorics involved.

> Is this the only reason ?

No, there are two others.  The first is that single-precision floats
(32 bits) are nearly useless, since modern processors compute just as
fast over double-precision floats.  The only case where you'd want to
use single-precision floats is as a storage type, to halve the size of
large arrays of numbers, but that's an option supported by bigarrays.

The second is that extended-precision floats are not fully specified
(these can use 80, 96 or 128 bits if I remember correctly), and few
processors implement them efficiently in hardware: IA32, perhaps IA64,
maybe some IBM 360 series processors.  Even on IA32, the trend is to
use SSE2 for floating-point computation, which don't support 80-bit
floats.  So, again, if you program for the common denominator, you'd
use 64-bit floats and nothing else.

- Xavier Leroy

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners

RE: [Caml-list] Efficient and canonical set representation?

[Julien Signoles]

On Mon, 10 Nov 2003, Diego Olivier Fernandez Pons wrote:

> Patricia sets seem to be what you are looking for.
>  (1). Efficient usual operations (lookup, insertion, union)
>  (2). Structural equality
>
> Their only problem is that they cannot handle polymorphic orderable
> types but only integers...
>
> Hash the data, use this key to insert it in a patricia map and solve
> the collisions by chaining in an ordered list (with the polymorphic
> [compare] function). (1) and (2) still hold under usual hypothesis on
> the rate of collisions.
>
> A few changes to JCF's implementation should be enough.

I think JCF's Hmap module is what you want.
A hmap is a map over hash-consed values implemented as Patricia Trees.
See http://www.lri.fr/~filliatr/software.en.html for more details.

Julien Signoles
-- 
mailto:Julien.Signoles@lri.fr ; http://www.lri.fr/~signoles
"In theory, practice and theory are the same,
but in practice they are different" (Larry McVoy)

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners