emacs-devel
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Code for converting between Elisp and Calc floats


From: Stefan Monnier
Subject: Re: Code for converting between Elisp and Calc floats
Date: Thu, 22 Oct 2009 20:50:53 -0400
User-agent: Gnus/5.13 (Gnus v5.13) Emacs/23.1.50 (gnu/linux)

> Vincent Belaïche has been working on making Calc more useful as an
> Emacs library.  Among the things he's working on is efficient
> conversions between Calc floats and Elisp floats.  This is best done
> using some lisp functions that are written in C.  Since it's C code, I
> wanted an okay before I committed it.  I have a copy of it at the end
> of this message.

The coding style needs to be adjusted: be careful with indentation,
spacing around parentheses, placement of function names (at beginning
of line, with the return type on the previous line), ...

As for the feature itself, I'm not sure whether I like it.
Could someone explain exactly what's the rationale behind this
(e.g. a typical use case), what alternatives were considered, etc...?
E.g. an "obvious" alternative would be to add builtin "big floats" via
libgmp or some such.

Also, I'm not sure I understand the code since it refers to 32, 64, 80,
and 128bit floats, whereas Emacs normally only uses 64bit floats.


        Stefan


> Thanks,
> Jay


> /** Code:*/
> #include <config.h>
> #include "lisp.h"

> #if HAVE_X_WINDOWS
> #include "xterm.h"
> #endif

> #if STDC_HEADERS
> #include <float.h>
> #endif

> /* if HAVE_DOUBLE_SIZEOF is not defined, then 
>    try to define it from float.h, when possible, the 
>    values in table below are considered.

>    Otherwise you may define it from the make command line if you 
>    dare!

>   +--------------------+----------------------------------------+
>   |HAVE_DOUBLE_SIZEOF  |meaning                                 |
>   +--------------------+----------------------------------------+
>   |0                   |binary IEEE754 is not used by Lisp      |
>   |                    |floats, or Lisp floats are not supported|
>   +--------------------+----------------------------------------+
>   |4                   |32bit binary IEEE is used by Lisp floats|
>   |                    |                                        |
>   +--------------------+----------------------------------------+
>   |8                   |64bit binary IEEE is used by Lisp floats|
>   +--------------------+----------------------------------------+
>   |10                  |79bit binary IEEE is used by Lisp floats|
>   +--------------------+----------------------------------------+
>   |16                  |128bit binary IEEE is used by Lisp      |
>   |                    |floats                                  |
>   +--------------------+----------------------------------------+
>  */
   
> #ifndef HAVE_DOUBLE_SIZEOF 
> #  if !defined LISP_FLOAT_TYPE \
>   || !defined FLT_RADIX || !defined DBL_MANT_DIG \
>   || !defined DBL_MIN_EXP || !defined DBL_MAX_EXP \
>   || FLT_RADIX != 2
> #    define HAVE_DOUBLE_SIZEOF 0
> #  elif (DBL_MANT_DIG == 1+23 \
>                && DBL_MIN_EXP == 1+-126 \
>                && DBL_MAX_EXP == 1+127)
> #    define HAVE_DOUBLE_SIZEOF 4
> #  elif (DBL_MANT_DIG == 1+52 \
>                && DBL_MIN_EXP == 1+-1022 \
>                && DBL_MAX_EXP == 1+1023)
> #    define HAVE_DOUBLE_SIZEOF 8
> #  elif (DBL_MANT_DIG == 1+64 \
>                && DBL_MIN_EXP == 1+-16382 \
>                && DBL_MAX_EXP == 1+16383)
> #    define HAVE_DOUBLE_SIZEOF 10
> #  elif (DBL_MANT_DIG == 1+112 \
>                && DBL_MIN_EXP == 1+-16382 \
>                && DBL_MAX_EXP == 1+16383)
> #    define HAVE_DOUBLE_SIZEOF 16
> #  else
> #    define HAVE_DOUBLE_SIZEOF 0
> #  endif

> #endif

> /* end of to be moved to config.h */

> /* #########################################################################*/
> /* Essential macro definitions for definition of IEEE754 format */
> /* =========================================================================*/

> #define HAVE_BINARY_IEEE754 1
> #ifdef HAVE_DOUBLE_SIZEOF 
> #  if HAVE_DOUBLE_SIZEOF == 4
> #    define IEEE754_EXP_SIZE   8
> #    define IEEE754_MANT_SIZE  23/* no phantom*/
> #  elif HAVE_DOUBLE_SIZEOF == 8
> #    define IEEE754_EXP_SIZE   11
> #    define IEEE754_MANT_SIZE  52/* no phantom*/
> #  elif HAVE_DOUBLE_SIZEOF == 10
> #    define IEEE754_EXP_SIZE   15
> #    define IEEE754_MANT_SIZE  64/* no phantom*/
> #  elif HAVE_DOUBLE_SIZEOF == 16
> #    define IEEE754_EXP_SIZE   15
> #    define IEEE754_MANT_SIZE  112/* no phantom*/
> #  else
> #    undef HAVE_BINARY_IEEE754
> #  endif
> #else
> #  error incomplete config, missing HAVE_DOUBLE_SIZEOF 
> #endif

> /* #########################################################################*/
> /* Lisp variables support in C */
> /* =========================================================================*/
> int        math_lisp_float_binary_ieee754_conformance;
> #ifdef HAVE_BINARY_IEEE754
> EMACS_INT  math_binary_ieee754_float_max_exp,
>   math_binary_ieee754_float_min_exp,
>   math_binary_ieee754_mant_size;
> #endif

> /* #########################################################################*/
> /* Other macro definitions for handling of IEEE754 binary format */
> /* =========================================================================*/
> #ifdef HAVE_BINARY_IEEE754
> #  define MATH_DOUBLE_INT_COUNT ((sizeof(double)+sizeof(unsigned 
> int)-1)/sizeof(unsigned int))

> #  define MATH_SIGN_POS (IEEE754_EXP_SIZE+IEEE754_MANT_SIZE \
>                                          
> -(MATH_DOUBLE_INT_COUNT-1)*sizeof(unsigned int)*8)
> #  define MATH_EXP_MASK ((1<<IEEE754_EXP_SIZE)-1)
> #  define MATH_EXP_BIAS ((1<<(IEEE754_EXP_SIZE-1))-1)
> #  define MATH_MANT_LIST_SIZE ((IEEE754_MANT_SIZE+1 + VALBITS-2)/(VALBITS-1))

> #  ifndef WORDS_BIG_ENDIAN
> /* little endian case */
> #    define MATH_DOUBLE_INT_LSB 0
> #    define MATH_DOUBLE_INT_MSB (MATH_DOUBLE_INT_COUNT-1)
> #    define MATH_DOUBLE_INT_INC 1
> #  else
> /* big endian case */
> #    define MATH_DOUBLE_INT_LSB (MATH_DOUBLE_INT_COUNT-1)
> #    define MATH_DOUBLE_INT_MSB 0
> #    define MATH_DOUBLE_INT_INC -1
> #  endif
> typedef union union_double_map_ot
> {
>   unsigned int m_aui[MATH_DOUBLE_INT_COUNT];
>   double       m_f;
> }
> double_map_ot;
> #endif


> /* #########################################################################*/
> /* Private functions */
> /* =========================================================================*/
> #ifdef HAVE_BINARY_IEEE754

> /* Produce a Not-a-Number. `sign' is the sign bit and shall be either 0 or 1. 
>  */
> static Lisp_Object math_make_nan(int sign)
> {
>   double_map_ot double_map_o;
>   int           i;
>   for(i = 0; i < MATH_DOUBLE_INT_COUNT; ++i)
>       double_map_o.m_aui[i] = ~0;

>   if(!sign)
>       double_map_o.m_aui[MATH_DOUBLE_INT_MSB] &= ~(1<<MATH_SIGN_POS);

>   return make_float(double_map_o.m_f);
> }

> /* fill a float starting from MSB*/
> static double math_make_float(int bitfield_count,unsigned int const 
> bit_val_aui[],int const bit_count_ai[])
> {
>   double_map_ot double_map_o;
>   int j = 0;
>   int i = MATH_DOUBLE_INT_MSB;
>   unsigned int bit_val;
>   int bit_count;
>   int remain_count;
>   int to_fill;
>   unsigned int word = 0;
  
>   remain_count = 0;
>   for(j = 0; j < bitfield_count; ++j)
>       remain_count += bit_count_ai[j];
>   to_fill = MATH_SIGN_POS+1;
>   bit_val   = bit_val_aui[j = 0];
>   bit_count = bit_count_ai[j++];

          
>   for(;;)
>       {
>         if(to_fill >= bit_count)
>               {
>                 bit_val <<= (to_fill - bit_count);
>                 word |= bit_val;
>                 remain_count -= bit_count;
>                 to_fill -= bit_count;
>                 bit_count = 0;
>               }
>         else
>               {
>                 word |= bit_val>>(bit_count - to_fill);
>                 bit_count -= to_fill;
>                 remain_count -= to_fill;
>                 /* erase all bits more significant than bit_count in 
>                        bit_val */
>                 to_fill = sizeof(unsigned int)*8 - bit_count;
>                 bit_val <<= to_fill;
>                 bit_val >>= to_fill;
>                 to_fill = 0;
>               }
>         if(remain_count == 0)
>               goto fill_complete;
>         else if(remain_count < 0)
>               {
>                 word <<= -remain_count;
>                 goto fill_complete;
>               }
>         else
>               {
>                 if(to_fill == 0)
>                       {
>                         double_map_o.m_aui[i] = word;
>                         word = 0;
>                         i -= MATH_DOUBLE_INT_INC;
>                         to_fill = sizeof(unsigned int)*8;
>                       }
>                 if(bit_count == 0)
>                       {
>                         bit_val   = bit_val_aui[j];
>                         bit_count = bit_count_ai[j++];
>                       }
>               }
>       }
>  fill_complete:
>   do
>       {
>         double_map_o.m_aui[i] = word;
>         word = 0;
>       }
>   while((i -= MATH_DOUBLE_INT_INC) != 
>               (MATH_DOUBLE_INT_LSB-MATH_DOUBLE_INT_INC));

>   /* printf("math_make_float: %ud 
> %ud\n",double_map_o.m_aui[1],double_map_o.m_aui[0]);*/
>   return double_map_o.m_f;
> }

> /* sign bit shall be 0 or 1 */
> static Lisp_Object math_make_infty(int sign)
> {
>   unsigned int bit_val_aui[] = {sign, MATH_EXP_MASK};
>   static int const bit_count_ai[] = { 1, IEEE754_EXP_SIZE};

>   return make_float(math_make_float(2,bit_val_aui,bit_count_ai));
> }

> /* make infinitly small float number, with sign (this is same as signed 0.0 
> in 
>    IEEE754) */ 
> static Lisp_Object math_make_epsilon(int is_neg)
> {
>   double_map_ot double_map_o;
>   int i = MATH_DOUBLE_INT_LSB;
>   double_map_o.m_aui[i] = 0;
>   while(i != MATH_DOUBLE_INT_MSB)
>       {
>         i += MATH_DOUBLE_INT_INC;
>         double_map_o.m_aui[i] = 0;
>       }
>   if(is_neg)
>       {
>         double_map_o.m_aui[i] |= (1<<MATH_SIGN_POS);
>       }
>   return make_float(double_map_o.m_f);
> }

> /* return sign bit of x */
> static int math_get_sign(double x)
> {
>   double_map_ot double_map_o;
>   double_map_o.m_f = x;
>   return double_map_o.m_aui[MATH_DOUBLE_INT_MSB] >> MATH_SIGN_POS;
> }


> /* return biased exp value of x */
> static unsigned int math_get_exp(double_map_ot dm_o)
> {
>   int i = MATH_DOUBLE_INT_MSB;
>   unsigned int word;
>   int remain_count;
>   int to_fill = sizeof(unsigned int)*8 - MATH_SIGN_POS;
>   word = (dm_o.m_aui[i] << to_fill)>>to_fill;
>   remain_count = IEEE754_EXP_SIZE-MATH_SIGN_POS;
>   if(remain_count < 0)
>       word >>= -remain_count;
>   else while(remain_count > 0)
>       {
>         to_fill = remain_count;
>         if(to_fill > sizeof(unsigned int)*8)
>               to_fill = sizeof(unsigned int)*8;
>         remain_count -= to_fill;
>         if(remain_count < 0)
>               {
>                 to_fill += remain_count;
>                 remain_count = 0;
>               }
>         word <<= to_fill;
>         i += MATH_DOUBLE_INT_INC;
>         word |= dm_o.m_aui[i] >> to_fill;
>       }
>   return word;
> }

> static unsigned int math_rotate_left(unsigned int what,int how_much)
> {
>   what = (what<<how_much)|(what>>(sizeof(unsigned int)*8-how_much));
>   return what;
> }
> static unsigned int math_rotate_right(unsigned int what,int how_much)
> {
>   what = (what>>how_much)|(what<<(sizeof(unsigned int)*8-how_much));
>   return what;
> }
> /* could be optimized if there is some intrinsic function */
> #define MATH_ROTATE_LEFT(W,H) W = math_rotate_left(W,H)
> #define MATH_ROTATE_RIGHT(W,H) W = math_rotate_right(W,H)

> /* set biased exp value of x */
> static void math_set_exp(double_map_ot* dm_po,unsigned int word)
> {
>   int i = MATH_DOUBLE_INT_MSB;
>   int remain_count = IEEE754_EXP_SIZE;
>   int pos = MATH_SIGN_POS;
>   int to_fill;

>   for(;;)
>       {
>         unsigned int bit_val = dm_po->m_aui[i];
>         to_fill = pos;
>         remain_count -= to_fill;
>         if(remain_count < 0)
>               {
>                 to_fill += remain_count;
>                 remain_count = 0;
>               }
>         /* bring exponent bits to lsb position */
>         pos -= to_fill;
>         MATH_ROTATE_RIGHT(bit_val,pos);
>         /* erase exponent bits */
>         bit_val >>= to_fill;
>         bit_val <<= to_fill;
>         /* fill exponent bits */
>         bit_val |= word>> remain_count;
>         if(remain_count != 0)
>               {
>                 /* no need to rotate left, if we get here this means that  
>                        pos = 0*/
dm_po-> m_aui[i] = bit_val;
>                 i += MATH_DOUBLE_INT_INC;
>                 pos  = sizeof(unsigned int)*8;
>                 /* erase all bits in word that are more significant than 
>                  the remain_count least significant bits*/
>                 to_fill = pos-remain_count;
>                 word <<= to_fill;
>                 word >>= to_fill;

>               }
>         else
>               {
>                 /* bring exponent bits back to orginal position */
>                 MATH_ROTATE_LEFT(bit_val,pos);
dm_po-> m_aui[i] = bit_val;
>                 break;
>               }
>       }
  
> }
> #endif

> /* #########################################################################*/
> /* public functions (visible from Lisp) */
> /* =========================================================================*/
> #ifdef HAVE_BINARY_IEEE754
> DEFUN ("math-floatnum-binary-ieee754-formatter", 
>          Fmath_floatnum_binary_ieee754_formatter, 
>          Smath_floatnum_binary_ieee754_formatter, 
>          4, 4, 
>          0,
>        doc: /* Formats floating point number into the binary IEEE754 floating
> point format used by native Lisp floating point numbers.

> When this number is not special, its value is :

>       SIGN*MANTISSA * 2^EXPONENT 

> where SIGN, MANTISSA and EXPONENT are numbers derived from
> IS-NEG, MANT and EXP arguments.

> 1. IS-NEG is t or nil, t for negative, and nil otherwise. SIGN is
> -1 for negative, and +1 otherwise.

> 2. IS-SPECIAL is `inf', `nan', `zero', `epsilon' or `nil'
> repectively for infinitly big, Not-a-Number, zero, infinitely
> small, and a normal non-zerofloating point number.

> 3. When IS-SPECIAL is non nil, MANT shall be nil, otherwise MANT shall be a
> non empty list of cons cells representing mantissa bits from most
> significant downto least significant:

>  (N_k . B_k) (N_{k-1} . B_{k-1}) ...  (N_1 . B_1) (N_0 . B_0)

> where N_i and B_i are integers such that

> for all i from 0 to k,
>    0 <B_i, and
>    0 =< N_i < 2^B_i

> and 2^{(B_k)-1} =< N_k
    
> and such that the MANTISSA can be written as 

> \(N_0 + 2^B_0 *
>  \(N_1 + 2^B_1 *
>     \( ...
>        \(N_{k-1}+2^B_{k-1}*N_k\) ...\)\)\)
> * 2^-(B_0+B_1+...+B_{k-1}+B_k - 1)

> That is to say the mantissa can be obtained by concatenating k binary words
> N_0 though N_k with respective bit counts B_0 through B_k, where N_k is the
> most significant part first in list, and N_0 the least significant one, last
> in the list, and by placing point just after the most significant bit.

> 4. EXP is an integer giving the exponent. Bias is not applied. In
> other words EXPONENT is directly equal to EXP.*/)
>        (is_neg,is_special,mant,exp)
>        Lisp_Object is_neg;
>        Lisp_Object is_special;
>        Lisp_Object mant;
>        Lisp_Object exp;
> {
>   register Lisp_Object val = Qnil;
>   struct gcpro gcpro1, gcpro2, gcpro3, gcpro4;

>   GCPRO4 (is_neg,is_special,mant,exp);
>   CHECK_SYMBOL(is_neg);
>   CHECK_SYMBOL(is_special);
>   CHECK_NUMBER(exp);

>   /* voir (elisp) Writing Emacs Primitives*/
>   if(NILP(is_special))
>       {
>         double_map_ot double_map_o;
>         int exp_val = XINT(exp);
>         unsigned int word = 0;
>         int i = MATH_DOUBLE_INT_MSB;
>         /* add one, because of phantom bit */
>         int to_fill =  1+MATH_SIGN_POS;
>         int bit_count;
>         unsigned int bit_val;
>         Lisp_Object elt;


>         CHECK_LIST_CONS(mant,mant);
          
>         /* remove phantom bit from first element of mantissa */
>         elt = XCAR(mant);
>         mant = XCDR(mant);
>         CHECK_NATNUM_CAR(elt);
>         CHECK_NATNUM_CDR(elt);
>         bit_val = XINT((XCAR (elt)));
>         bit_count = XINT((XCDR (elt)));
>         if(bit_count > 1)
>           {
>             --bit_count;
>             {
>               int to_erase = sizeof(unsigned int)*8-bit_count;
>               bit_val <<= to_erase;
>               bit_val >>= to_erase;
>             }
>             mant = Fcons(
>                          Fcons(make_number(bit_val),make_number(bit_count)),
>                          mant);
>           }
>         else if(bit_count != 1)
>           error("invalid value for (cdar mant), shall be a positive integer");
>         /* bias the exponent, and push it to mant list*/
>         exp_val += MATH_EXP_BIAS;
>         if(exp_val & ~MATH_EXP_MASK)
>               error("Invalid exponent value");
>         mant = 
> Fcons(Fcons(make_number(exp_val),make_number(IEEE754_EXP_SIZE)),mant);

>         /* push sign to mant list */
>         bit_val = NILP(is_neg)?0:1;
>         mant = Fcons(Fcons(make_number(bit_val),make_number(1)),mant);
          
>         /* now read out mant list into double_map_o */
>         do
>               {
>                 elt = XCAR(mant);
>                 mant = XCDR(mant);
>                 CHECK_NATNUM_CAR(elt);
>                 CHECK_NATNUM_CDR(elt);
>                 bit_val = XINT((XCAR (elt)));
>                 bit_count = XINT((XCDR (elt)));
>                 do
>                       {
>                         if(to_fill >= bit_count)
>                               {
>                                 bit_val <<= (to_fill - bit_count);
>                                 word |= bit_val;
>                                 to_fill -= bit_count;
>                                 bit_count = 0;
>                               }
>                         else
>                               {
>                                 word |= bit_val>>(bit_count - to_fill);
>                                 bit_count -= to_fill;
>                                 /* erase all bits more significant than 
> bit_count in 
>                                    bit_val */
>                                 to_fill = sizeof(unsigned int)*8 - bit_count;
>                                 bit_val <<= to_fill;
>                                 bit_val >>= to_fill;
>                                 to_fill = 0;
>                               }
>                         if(to_fill == 0)
>                               {
>                                 if(i == MATH_DOUBLE_INT_LSB)
>                                   goto mant_complete;
>                                 else
>                                       to_fill = sizeof(unsigned int)*8;
>                                 double_map_o.m_aui[i] = word;
>                                 word = 0;
>                                 i -= MATH_DOUBLE_INT_INC;
>                               }
>                       }
>                 while(bit_count != 0);

>               }
>         while (CONSP (mant));
> mant_complete:
>         do
>           {
>             double_map_o.m_aui[i] = word;
>             word = 0;
>           }
>         while((i -= MATH_DOUBLE_INT_INC) != 
>               (MATH_DOUBLE_INT_LSB-MATH_DOUBLE_INT_INC));

          
>         val = make_float(double_map_o.m_f);

>       }
>   else
>       {
>         char const* is_special_symbol_name = 
>               (char const*)( SDATA (SYMBOL_NAME (is_special)));

>         if(!EQ(mant,Qnil))
>           {
>             error("mant shall be `nil' when is-special is non `nil'");
>           }
>         if(strncmp (is_special_symbol_name,
>                                 "zero", sizeof ("zero") - 1) == 0)
>               {
>                 val = make_float(0.0);
>               }
>         else if(strncmp (is_special_symbol_name,
>                                 "epsilon", sizeof ("epsilon") - 1) == 0)
>               {
>                 val = math_make_epsilon(!NILP(is_neg));
>               }
>         else if(strncmp (is_special_symbol_name,
>                                 "nan", sizeof ("nan") - 1) == 0)
>               {
>                 val = math_make_nan(!NILP(is_neg));
>               }
>         else if(strncmp (is_special_symbol_name,
>                                 "infty", sizeof ("infty") - 1) == 0)
>               {
>                 val = math_make_infty(!NILP(is_neg));
>               }
>         else
>               error("unexpected symbol for is-special: 
> %s",is_special_symbol_name);
>       }

>   UNGCPRO;
>   return val;
> }

> #define MATH_IS_NEG     ret[0]
> #define MATH_IS_SPECIAL ret[1]
> #define MATH_MANT       ret[2]
> #define MATH_EXP        ret[3]
> DEFUN ("math-floatnum-binary-ieee754-expander",
>          Fmath_floatnum_binary_ieee754_expander,
>          Smath_floatnum_binary_ieee754_expander,
>          1,1,
>          0,
>          doc: /*
> Analyses a floating point number X conforming to the binary
> IEEE754 format used by Lisp floatingin point numbers, and returns
> a list \(IS-NEG IS-SPECIAL MANT EXP) where

> IS-NEG IS-SPECIAL MANT and EXP have the same meaning as in
> function `math-floatnum-binary-ieee754-formatter'.*/)
>        (x)
>          Lisp_Object x;
> {
>   register Lisp_Object val = Qnil;
>   Lisp_Object ret[4];
>   double_map_ot double_map_o;
>   int exp_val,sign,i,j,bit_count,remain_count,filled,to_fill;
>   unsigned int word,bit_val,pos;
>   int          bit_count_ai[2+MATH_MANT_LIST_SIZE];
>   unsigned int bit_val_aui[2+MATH_MANT_LIST_SIZE];
>   Lisp_Object  mant_list[MATH_MANT_LIST_SIZE];
>   int          is_mant_zero;

>   CHECK_FLOAT(x);
>   double_map_o.m_f = XFLOAT_DATA(x);

>   bit_count_ai[0] = 1; /* sign bit size */
>   bit_count_ai[1] = IEEE754_EXP_SIZE; /* sign bit size */
>   remain_count = IEEE754_MANT_SIZE;

>   i = 2;
>   to_fill = (VALBITS-1)-1; /* first word shall have one spare bit for phantom 
> */
>   do
>       {
>         to_fill = remain_count > to_fill ? to_fill: remain_count;
>         bit_count_ai[i++] = to_fill;
>         remain_count-= to_fill;
>         to_fill = VALBITS-1;
>       }
>   while(remain_count != 0);

>   /* analyse*/
>   remain_count = IEEE754_EXP_SIZE+IEEE754_MANT_SIZE+1;
>   i = MATH_DOUBLE_INT_MSB;
>   word = 0;
>   filled = 0;
>   j = 0;
>   pos = 1+MATH_SIGN_POS;
>   do
>       {

>         bit_val = double_map_o.m_aui[i];
>         i -= MATH_DOUBLE_INT_INC;
>         if(remain_count > pos)
>               {
>                 bit_count = pos;
>                 remain_count -= pos;
>               }
>         else
>               {
>                 bit_val >>= pos - remain_count;
>                 bit_count = remain_count;
>                 remain_count = 0;
>               }

>         do
>               {
>                 /* compute how many bits can be filled into word */
>                 to_fill = bit_count;
>                 if(to_fill + filled > bit_count_ai[j])
>                       to_fill = bit_count_ai[j] - filled;
                          
>                 /* fill those bits into word, bit_count becomes the remaining
>                        number of bits in bit_val */
>                 word <<= to_fill;
>                 bit_count -= to_fill;
>                 word |= bit_val>>bit_count;
>                 filled += to_fill;

>                 /* erase all bits of bit_val that are more significant than
>                        the bit_count least significant ones */
>                 to_fill = sizeof(unsigned int)*8-bit_count;
>                 bit_val <<= to_fill;
>                 bit_val >>= to_fill;

>                 /* if there is no more room in word, then take the next item
>                        in mant_list */
>                 if(filled == bit_count_ai[j])
>                       {
>                         bit_val_aui[j++] = word;
>                         filled = 0;
>                         word   = 0;
>                       }

>               }
>         while(bit_count != 0);

>         pos = sizeof(unsigned int)*8;
                  
>       }
>   while(remain_count != 0);
>   /* end of analyse */

>   is_mant_zero = 1;
>   for(i = 2; i < 2+MATH_MANT_LIST_SIZE; ++i)
>       {
>         if(bit_val_aui[i] != 0)
>               {
>                 is_mant_zero = 0;
>                 break;
>               }
>       }

>   /* check special values */ 
>   if(bit_val_aui[1] == 0 && is_mant_zero)
>       {
>         MATH_IS_SPECIAL = intern("zero");
>         MATH_MANT   = Qnil;
>         MATH_EXP    = make_number(0);
>         MATH_IS_NEG = bit_val_aui[0]?Qt:Qnil;
>       }
>   else if(bit_val_aui[1] == MATH_EXP_MASK)
>       {
>         /* Infinity or NaN */
>         if(is_mant_zero)
>               {
>                 MATH_IS_SPECIAL = intern("infty");
>               }
>         else
>               {
>                 MATH_IS_SPECIAL = intern("nan");
>               }
>         MATH_MANT   = Qnil;
>         MATH_EXP    = make_number(0);
>         MATH_IS_NEG = bit_val_aui[0]?Qt:Qnil;
>       }
>   else
>       {
>         /* normal number */
>         MATH_IS_SPECIAL = Qnil;

>         /* restore phantom */
>         bit_val_aui[2] |= 1<<bit_count_ai[2];
>         ++bit_count_ai[2];

>         for(i = 0; i < MATH_MANT_LIST_SIZE; ++i)
>               mant_list[i] = Fcons(make_number(bit_val_aui[i+2]),
>                                                        
> make_number(bit_count_ai[i+2]));

>         MATH_MANT   = Flist(MATH_MANT_LIST_SIZE,mant_list);
>         MATH_EXP    = make_number((int)bit_val_aui[1] - MATH_EXP_BIAS);
>         MATH_IS_NEG = bit_val_aui[0]?Qt:Qnil;

>       }
  
>   val = Flist(4,ret);

>   return val;
> }
> #endif

> /* #########################################################################*/
> /* Syms of Calc */
> /* =========================================================================*/
> void syms_of_calc(void)
> {
> #ifdef HAVE_BINARY_IEEE754
>   defsubr (&Smath_floatnum_binary_ieee754_expander);
>   defsubr (&Smath_floatnum_binary_ieee754_formatter);
> #endif

>   DEFVAR_BOOL ("math-lisp-float-binary-ieee754-conformance", 
>                          &math_lisp_float_binary_ieee754_conformance,
>              doc: /* t if Lisp floats are conformant to one of the binary
>                                  ieee754 format, nil otherwise */);
>  math_lisp_float_binary_ieee754_conformance = 
> #ifdef HAVE_BINARY_IEEE754
>    1
> #else
>    0
> #endif
>    ;
>   XSYMBOL (intern ("math-lisp-float-binary-ieee754-conformance"))->constant = 
> 1;

> #ifdef HAVE_BINARY_IEEE754
>   DEFVAR_INT ("math-binary-ieee754-float-max-exp", 
>                          &math_binary_ieee754_float_max_exp,
>              doc: /* The largest exponent value for a lisp float, defined 
> only on machines where lisp float are using one of the binary IEEE754 
> formats.  */);
>   math_binary_ieee754_float_max_exp = MATH_EXP_MASK-1-MATH_EXP_BIAS;
>   XSYMBOL (intern ("math-binary-ieee754-float-max-exp"))->constant = 1;

>   DEFVAR_INT ("math-binary-ieee754-float-min-exp", 
>                          &math_binary_ieee754_float_min_exp,
>              doc: /* The smalled exponent value for a lisp float, defined 
> only on machines where lisp float are using one of the binary IEEE754 
> formats.  */);
>   math_binary_ieee754_float_min_exp = 1-MATH_EXP_BIAS;
>   XSYMBOL (intern ("math-binary-ieee754-float-min-exp"))->constant = 1;

>   DEFVAR_INT ("math-binary-ieee754-mant-size", 
>                          &math_binary_ieee754_mant_size,
>              doc: /* The number of bits in a binary IEEE754 mantissa, 
> including
> the phantom bit. That is to say, the number of significant bits. */);
>   math_binary_ieee754_mant_size = IEEE754_MANT_SIZE+1;/* add 1 for phantom*/
>   XSYMBOL (intern ("math-binary-ieee754-mant-size"))->constant = 1;

> #endif
  
> }


> /** mathfloat.c ends here*/





reply via email to

[Prev in Thread] Current Thread [Next in Thread]