Re[3]: floats

[Bruno Haible <haible@ilog.fr>]

Sam Steingold <sshteingold@cctrading.com> writes:
> 
>      1. I do set *read-default-float-format* = double-float. The numbers 
>      are still read as singles.

Just a guess why this happens: Maybe you have floating point constants
without precision suffix in your program. When you load a compiled file
containing

       (setq *read-default-float-format* 'double-float)
       (defvar x 1.0)

then x will be a single-float nevertheless. To avoid this, wrap the
setq form inside  (eval-when (compile load eval) ....).

>      To use your impnotes example: "(- (+ 1.7 pi) pi)  should not return  
>      1.700000726342836417234L0,"  This 1.700000whatever is still 1.7, for 
>      purposes of practical computations.

But when a program prints out a 21-digit number, how could you (as a user)
guess that this number has only 7 correct places?

>      But when I am getting *negative* 
>      standard deviations, this is not correct in any possible sense!!!!

If you are getting negative radicands when computing standard deviations,
you either

   need to take a formula with better numerical stability [i.e. replace
   n*sum(i,x_i^2) - sum(i,x_i)^2  by  n*sum(i,(x_i-sum(j,x_j)/n)^2) ],

or

   need to work with larger precision. In this case, if somehow a
   single-float enters the computations, you would still see negative
   radicands and wonder where they come from. With CLISP you'll see
   a single-float negative radicand, and that will point you to the
   problem.

>      You think that a float is always imprecise. You are *wrong*. 123456.7 
>      is *precise*, and I should not have to use ratios to tell clisp that!

If you want 123456.7 to be precise, you either need to use ratios (or,
equivalently, a fixed-point representation), or a floating-point
representation with base 10 (e.g. GNU bc uses this). CLISP chose, for
efficiency and for compliance with the IEEE standards, a floating-point
representation with base 2.

>      Anyway, can you name another computer system where 
>      single+double=single?!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Yes. PARI 2.0.alpha for example:

$ gp
                      GP/PARI CALCULATOR Version 2.0.alpha
                      (i586 running linux 32-bit version)
                  (readline disabled, extended help available)

                             Copyright 1989-1997 by
          C. Batut, K. Belabas, D. Bernardi, H. Cohen and M. Olivier.

Type ? for help.

   realprecision = 28 significant digits
   seriesprecision = 16 significant terms
   format = g0.28

parisize = 4000000, primelimit = 500000, buffersize = 30000
? x = sqrt(2)
%1 = 1.414213562373095048801688724
? precision(x)
%2 = 28
? default(realprecision,50)
   realprecision = 50 significant digits
? precision(x)
%3 = 28
? y = sqrt(3)
%4 = 1.7320508075688772935274463415058723669428052538103
? precision(y)
%5 = 57
? precision(x+y)
%6 = 28
? precision(y+x)
%7 = 28
? precision(x*y)
%8 = 28

                            Bruno