Code for converting between Elisp and Calc floats

Code for converting between Elisp and Calc floats

[Jay Belanger]

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.

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*/

Re: Code for converting between Elisp and Calc floats

[James Cloos]

>>>>> "Jay" == Jay Belanger <jay.p.belanger@gmail.com> writes:

Jay>   +--------------------+----------------------------------------+
Jay>   |10                  |79bit binary IEEE is used by Lisp floats|
Jay>   +--------------------+----------------------------------------+

Jay> #  elif (DBL_MANT_DIG == 1+64 \
Jay> 		 && DBL_MIN_EXP == 1+-16382 \
Jay> 		 && DBL_MAX_EXP == 1+16383)
Jay> #    define HAVE_DOUBLE_SIZEOF 10
Jay> #  elif (DBL_MANT_DIG == 1+112 \
Jay> 		 && DBL_MIN_EXP == 1+-16382 \
Jay> 		 && DBL_MAX_EXP == 1+16383)
Jay> #    define HAVE_DOUBLE_SIZEOF 16
Jay> #  else
Jay> #    define HAVE_DOUBLE_SIZEOF 0
Jay> #  endif

As there is no implied bit in the i387's 80-bit float; I suspect
that the DBL_MANT_DIG == 1+64 will never match.

Also -- and I've no idea whether this is relevant -- the 128 bit long
doubles used on IBM's mainframes (and perhaps elsewhere) are the sum
of two 64 bit doubles; they have more significant bits, but their
exponents are no larger than a double.

-JimC
-- 
James Cloos <cloos@jhcloos.com>         OpenPGP: 1024D/ED7DAEA6

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

> 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*/

Re: Code for converting between Elisp and Calc floats

[Richard Stallman]

That code needs a lot more comments to explain what's happening.
Also, would you please format the code in accord with the Emacs coding
style?

  /* voir (elisp) Writing Emacs Primitives*/
  if(NILP(is_special))
	{
	  double_map_ot double_map_o;
	  int exp_val = XINT(exp);

should be

  /* voir (elisp) Writing Emacs Primitives*/
  if (NILP (is_special))
    {
      double_map_ot double_map_o;
      int exp_val = XINT(exp);


And why does it say "voir"?  Is that French?  Is it a typo for `void'? 
If it is a typo for `void', I don't understand the meaning.

Re: Code for converting between Elisp and Calc floats

[David Kastrup]

Richard Stallman <rms@gnu.org> writes:

> That code needs a lot more comments to explain what's happening.
> Also, would you please format the code in accord with the Emacs coding
> style?
>
>   /* voir (elisp) Writing Emacs Primitives*/
>   if(NILP(is_special))
> 	{
> 	  double_map_ot double_map_o;
> 	  int exp_val = XINT(exp);
>
> should be
>
>   /* voir (elisp) Writing Emacs Primitives*/
>   if (NILP (is_special))
>     {
>       double_map_ot double_map_o;
>       int exp_val = XINT(exp);
>
>
> And why does it say "voir"?  Is that French?  Is it a typo for `void'? 
> If it is a typo for `void', I don't understand the meaning.

It's basically "See (elisp) Writing Emacs Primitives"

-- 
David Kastrup

Re: Code for converting between Elisp and Calc floats

[Richard Stallman]

    > And why does it say "voir"?  Is that French?  Is it a typo for `void'? 
    > If it is a typo for `void', I don't understand the meaning.

    It's basically "See (elisp) Writing Emacs Primitives"

Comments in Emacs should be written in English.
I can guess what "voir" means here, since I speak French,
but most Emacs hackers don't speak French.

Re: Code for converting between Elisp and Calc floats

[Jay Belanger]

Stefan Monnier <monnier@iro.umontreal.ca> writes:

> 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...?

Vincent explained this to me, but I expect he'll come along and do a
better job explaining than I can.

> E.g. an "obvious" alternative would be to add builtin "big floats" via
> libgmp or some such.

Using libgmp in Emacs has been discussed in the past, but as I recall
hasn't had much support.

Jay

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

>> 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...?
> Vincent explained this to me, but I expect he'll come along and do a
> better job explaining than I can.

Is he reading this thread?

>> E.g. an "obvious" alternative would be to add builtin "big floats" via
>> libgmp or some such.
> Using libgmp in Emacs has been discussed in the past, but as I recall
> hasn't had much support.

It hasn't had much support indeed, mostly for lack of interest
in bignums.  But if we decide bignums deserve some C code, then maybe
libgmp would be a good idea.  Then again, without knowing the whys and
hows, it's hard to tell.


        Stefan

Re: Code for converting between Elisp and Calc floats

[David Kastrup]

Stefan Monnier <monnier@iro.umontreal.ca> writes:

>>> 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...?
>> Vincent explained this to me, but I expect he'll come along and do a
>> better job explaining than I can.
>
> Is he reading this thread?
>
>>> E.g. an "obvious" alternative would be to add builtin "big floats" via
>>> libgmp or some such.
>> Using libgmp in Emacs has been discussed in the past, but as I recall
>> hasn't had much support.
>
> It hasn't had much support indeed, mostly for lack of interest
> in bignums.  But if we decide bignums deserve some C code, then maybe
> libgmp would be a good idea.  Then again, without knowing the whys and
> hows, it's hard to tell.

There are few areas of importance.  One would be buffer_size.  Putting
them transparently in there would imply having them all over the place
in overlay positions, point and so on.  I am also not sure that bignums
of equal value compare as eq (floats don't) and there is likely code
around relying on it.  The buffer size problem would likely be less
radically addressed by adding another bit of tagging, and then use half
of the possible values for integers, giving a range of -2^30 .. 2^30-1.
That's what XEmacs does, and it goes a far way of extending the point of
nuisance.

Of course calc could profit if one ripped out all of its
multiple-precision list-based code (does not even use floats).  But
that would not require a transparent bignum type.  An on-demand type
would suffice.

-- 
David Kastrup

RE: Code for converting between Elisp and Calc floats

[Vincent Belaïche]

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

 	<87aazj2r7d.fsf_-_@vh213601.truman.edu>

 	<jwv7hunkkf4.fsf-monnier+emacs@gnu.org> <87ws2kmuez.fsf@gmail.com>

 	<jwvhbtoi9gj.fsf-monnier+emacs@gnu.org>
 
 <87hbto11d8.fsf@lola.goethe.zz>
MIME-Version: 1.0

[-- Attachment #2: Type: text/plain, Size: 2958 bytes --]


Hello,


Sorry for coming so late into this thread. Thanks Jay for introducing my contribution, and thanks to you all for the feedback. 

I fully agree with most of the comments: 

* code presentation (indentation, commenting) probably needs some improvement (actually it took some time before I got something working well, and now I realize that I've neglected this aspect)
* For your information "voir" just mean "see" in French, sorry for I forgot this in the code. Richard: if ever you get in touch with the French bureaucracy and have to fill in some form, then please remember that to fill an item by "void", just write "néant",
* I developped and tested the code for the 64bit little endian float implementation of emacs (on a Mingw Windows2000 platform). I have put some code so that porting to other type of platform is easier, however I cannot guarantee that it would work at all for those ports. I did not know that emacs had 64bits floats only in all platform (and it seems to me quite unwise to base some code on the assumption that it will always & everywhere be that way).

Now coming back to the contribution itself: of course we could discuss at length whether or not I did it the right way. However, we should first agree that there some need for having conversion function between Lisp floats and calc floats. Lisp float are the usual IEEE753 64bits radix 2 floats, while Calc floats are list based radix-10, configurable precision.

The need which I see is that I often use calc as a math library. Native lisp has a few math primitives, however, when it comes to matrices, polynomials, and suchlikes, only Calc can make it. Most calc users use calc only as a Desk calculator, not as a math library. Most often I use it as a library for some preprocessor code that I include inline within some LaTeX document to generate some calculated texted values (I intend to make this preprocessor code public some day when it is mature enough, but I can provide draft code to anyone interested in it). 

When you use Calc as a library, then you would like it to interact more easily with some other piece of code, or simply to write the code more simply when the precision you need is only that of a 64bit float: nobody will deny that it is simpler to write: 

0.5 

instead of 

'(float 5 -1)


However, currently calc do not accept numbers like 0.5 (those would be processed as if they were integers, and this  causes at best an error throw, and at worst a buggy code).

So if we come to the conclusion that indeed there is a need for Calc-to-Lisp + Lisp-to-Calc float conversion, then we can go into the next step of the discussion that is to discuss what is the best way to do it.

Regards,

    Vincent.
 		 	   		  
_________________________________________________________________
Nouveau ! Tout Windows débarque dans votre téléphone. Voir les Windows phone
http://clk.atdmt.com/FRM/go/175819071/direct/01/

[-- Attachment #3: Type: text/html, Size: 3179 bytes --]

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

> There are few areas of importance.  One would be buffer_size.  Putting

bignums for buffer size is a bad idea.  It would require an enormous
amount of code changes and a very significant slow down (buffer
positions are used extensively everywhere in the C code, including
performance critical code such as the rendering engine, and changing
them from scalar values to something else would slow things down
tremendously).

> The buffer size problem would likely be less radically addressed by
> adding another bit of tagging, and then use half of the possible
> values for integers, giving a range of -2^30 .. 2^30-1.

Currently we have -2^28..2^28-1 and we can easily extend it to
2^29..2^29-1 without any noticeable downside (I have been using such
a change for a while now).

> That's what XEmacs does, and it goes a far way of extending the point
> of nuisance.

Not far enough, AFAICT: if you regularl edit files somewhere between
512M and 1G, then it's quite likely that you'll also want to edit files
larger than 1G.  Currently in 64bit mode, you supposedly can edit files
larger but in practice when I tried it, I bumped into other bugs when
reaching the 2G limit (or maybe even earlier) IIRC.


        Stefan

Re: Code for converting between Elisp and Calc floats

[Chong Yidong]

Vincent Belaïche <vincent.b.1@hotmail.fr> writes:

> When you use Calc as a library, then you would like it to interact more easily
> with some other piece of code, or simply to write the code more simply when the
> precision you need is only that of a 64bit float: nobody will deny that it is
> simpler to write:
>
> 0.5
>
> instead of
>
> '(float 5 -1)

Could you discuss why this needs built-in support?  It should be
possible to write a Lisp function or macro that automatically converts
Lisp expressions based on Lisp floats into Lisp expressions based on
calc floats.

Re: Code for converting between Elisp and Calc floats

[Eli Zaretskii]

> From: Stefan Monnier <monnier@iro.umontreal.ca>
> Date: Sun, 25 Oct 2009 09:57:49 -0400
> Cc: emacs-devel@gnu.org
> 
> Currently in 64bit mode, you supposedly can edit files larger but in
> practice when I tried it, I bumped into other bugs when reaching the
> 2G limit (or maybe even earlier) IIRC.

There are lots of places that use int instead of EMACS_INT.  If we
want to get full 64-bit buffer size, Someone(TM) should hunt them down
and fix them.

RE: Code for converting between Elisp and Calc floats

[Vincent Belaïche]

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

 	<BLU104-W139663DC7482450EC1140284BB0@phx.gbl>
 

 <87d44bef7f.fsf@stupidchicken.com>
MIME-Version: 1.0

[-- Attachment #2: Type: text/plain, Size: 3638 bytes --]


Hello and thanks for the feedback,

First, I would like to clarify that I was not aware of the issue of buffer size and pointer representation and that this was not my intent. Up to now, I have not felt concerned by this issue and have little idea on the best way to solve it, but my first feeling from the guts is that some kind of bigger-integer with fixed dynamic, e.g. 64bit integer would be a better option than arbitrary dynamic true bignums. Bignums should be really reserved for some math or financial things were you need no limitation.

Then concerning the need for built-in support (point raised by cyd(string 64)stupidchicken.com), I fully agree that it is possible to do without it but as far as I made it, it was at the cost of some accuracy. 

First notes that the code that I submitted use built-in support only for a tiny part of the job that is kind-of bitwise access to  IEEE754 floats, the main part of the job (the computation from radix 10 to radix 2) is entirely done in Lisp. So the built-in function are quite generic ones and could be used by other libraries needing to generate specific IEEE754 floats.

Then also note that the Lisp code which I submitted actually contains both types of implementation: that using builtin support, and that not using it. The former is selected when the presence of builtin implementation is detected, while the other one should be selected otherwise (well I have not fully checked whether that switch works well).

The reason why builtin support is needed is that without it there are losts of cases when some number that should convert without any loss of information, will suffer from some rounding error. This is error is in the magnitude of the quantization accuracy of the number, so it should not be too harmful: however I felt that it is a pity when you can get something better not to have it. 

The reason for this issue may be simply that the Lisp-only code which I wrote is not good enough, I have no formal proof that to that extent using builtin is a must. So if you can take this code and achieve to make it all Lisp, please feel free to do that. For instance if you could write Lisp only function doing exactly the same job as the builtin functions, then this would be quite an acceptable way to me. I did not try that, because I felt that it was so considerably more simple to do this in C than in Lisp, that I did not think that it would be an issue if this part is done in C, once people agree on the need for the contribution, and on the algorithm used to make the conversion, along with the split between the Lisp only part and the C code.

Note that it took to me some time to make this work, the first attempt was a very simple Calc format to a string, followed by string to number, or vice versa. Then there was this lisp only implementation (that is still in the submitted code), and then only I made a novel implementation with the use of two builtin primitives.

   Vincent.



> From: cyd@stupidchicken.com
> To: vincent.b.1@hotmail.fr
> Date: Sun, 25 Oct 2009 10:11:16 -0400
> CC: emacs-devel@gnu.org
> Subject: Re: Code for converting between Elisp and Calc floats
>
[SNIP]

> Could you discuss why this needs built-in support?  It should be
> possible to write a Lisp function or macro that automatically converts
> Lisp expressions based on Lisp floats into Lisp expressions based on
> calc floats.
> 
> 
 		 	   		  
_________________________________________________________________
Nouveau Windows 7 : Simplifiez votre PC ! Trouvez le PC qui vous convient.
http://clk.atdmt.com/FRM/go/181574580/direct/01/

[-- Attachment #3: Type: text/html, Size: 3910 bytes --]

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

> First, I would like to clarify that I was not aware of the issue of
> buffer size and pointer representation and that this was not my
> intent.

Indeed it has nothing to do with the discussion at hand.

> Then concerning the need for built-in support (point raised by
> cyd(string 64)stupidchicken.com), I fully agree that it is possible to
> do without it but as far as I made it, it was at the cost of
> some accuracy.

That's interesting, so the main purpose is accuracy.

> First notes that the code that I submitted use built-in support only
> for a tiny part of the job that is kind-of bitwise access to  IEEE754
> floats, the main part of the job (the computation from radix 10 to
> radix 2) is entirely done in Lisp. So the built-in function are quite
> generic ones and could be used by other libraries needing to generate
> specific IEEE754 floats.

I didn't see this code, then (the code I saw only included C code).
Could someone describe the general skeleton of the whole code (plus
which part is in C, and what are the entry points)?

> The reason why builtin support is needed is that without it there are
> losts of cases when some number that should convert without any loss
> of information, will suffer from some rounding error. This is error is
> in the magnitude of the quantization accuracy of the number, so it
> should not be too harmful: however I felt that it is a pity when you
> can get something better not to have it. 

If you take (calc-in "3.5") as input rather than 3.5, then there is no
such issue, right?

> The reason for this issue may be simply that the Lisp-only code which
> I wrote is not good enough, I have no formal proof that to that extent
> using builtin is a must. So if you can take this code and achieve to
> make it all Lisp, please feel free to do that. For instance if you
> could write Lisp only function doing exactly the same job as the
> builtin functions, then this would be quite an acceptable way to
> me. I did not try that, because I felt that it was so considerably
> more simple to do this in C than in Lisp, that I did not think that it
> would be an issue if this part is done in C, once people agree on the
> need for the contribution, and on the algorithm used to make the
> conversion, along with the split between the Lisp only part and the
> C code.

I think I see what your two functions are doing, but I think their
naming and interface is a bit more complex than needed.  I'd suggest to
drop the reference to ieee754 from the name and doc (they should just
work for whatever floats are used internally, which in practice is just
what `double' uses in the C compiler used to compile Emacs).
Also, the mantissa should probably be represented as either a single
integer or a list of integers where each integer provides 16bits of data
(that's a format already used elsewhere in Emacs to represent "large
integers" such as for time).

You could rename them to something like construct-float and deconstruct-float.

I'd also rename the mantissa and exponent size info, or maybe even
consider providing it differently (in case its use is always linked to
calls to one of the previous 2 functions, maybe those functions could
return the relevant info.  E.g. make-float could return the part of the
mantissa it ignored).  Since you're using those constants, you clearly
know better how they're used, so I'll let you figure out whether there's
a better way to ptovide the same info.

Also, how is math-lisp-float-binary-ieee754-conformance used?

> Note that it took to me some time to make this work, the first attempt
> was a very simple Calc format to a string, followed by string to
> number, or vice versa.

What was the problem with this simple solution?


        Stefan

RE: Code for converting between Elisp and Calc floats

[Vincent Belaïche]

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

 	<BLU104-W164CDF6A8CEB6630D267D484BB0@phx.gbl>
 

 <jwvljiyykp3.fsf-monnier+emacs@gnu.org>
MIME-Version: 1.0

[-- Attachment #2.1: Type: text/plain, Size: 4403 bytes --]


Hello Stefan,

Thanks for your valuable feedback.

I attached the Lisp part of the contribution. The entry point is just the two functions math-numfloat and math-floatnum. The former converts from Lisp to Calc, and the latter the other way round.


> 

> If you take (calc-in "3.5") as input rather than 3.5, then there is no
> such issue, right?
>

I just tried it (with setting math-lisp-float-binary-ieee754-conformance to nil so that no builtin support is used)


(math-numfloat  3.5)

evals correctly to
(float 35 -1)

but 

(math-floatnum '(float 35 -1))

evals to:

3.4999999999968168

Note that by default Calc has a 12 digit accuracy (that is less than a double that has almost 16 digit accuracy). So, to that extent, this latter conversion does not betray the accuracy of Calc (the 1st false digit is after 12 good digits).

If I had used the version with builtin support, there would not be this issue  (to be 100% frank I did not checked it just now, but I am quite confident that (math-floatnum '(float 35 -1)) would have evaluated to 3.5).


[SNIP]

> Also, the mantissa should probably be represented as either a single

> integer or a list of integers where each integer provides 16bits of data
> (that's a format already used elsewhere in Emacs to represent "large
> integers" such as for time).
>

I did it this way because it allowed a simpler algorithm: in the conversion algorithm the integer size step is based on the number of decimal digits that Calc uses in representing big integers. The Calc big integers are not composed of 16bits integers (and actually there is not an integral number of bits, because these components are in radix 10). From that perspective it was easier if the number of bits of each component integer is passed along with the component integer in the interface.



> You could rename them to something like construct-float and deconstruct-float.
>

No problem to rename them, I fully agree that what name you propose is better.



> I'd also rename the mantissa and exponent size info, or maybe even
> consider providing it differently (in case its use is always linked to
> calls to one of the previous 2 functions, maybe those functions could
> return the relevant info.  E.g. make-float could return the part of the
> mantissa it ignored).  Since you're using those constants, you clearly
> know better how they're used, so I'll let you figure out whether there's
> a better way to ptovide the same info.
>


OK, this can be improved. Actually, my intention was not that these builtin functions were kind-of part of Calc

library API, these are really internal things that the user developing above Calc should not need.


> Also, how is math-lisp-float-binary-ieee754-conformance used?
>

This is the switch between the two types of implementations: the one using builtin support, and the one not needing builtin support.


Note that I can imagine that it is possible to implement the construct/deconstruct builtin functions just using lisp, this is after all just a matter of generating clean powers of 2. But, since this code would be so machine/system dependant, and since Lisp would be so more uselessly complex than C code,I thought that it was a better choice to make it in C.


> > Note that it took to me some time to make this work, the first attempt
> > was a very simple Calc format to a string, followed by string to
> > number, or vice versa.
> 
> What was the problem with this simple solution?
>

String formatting/reading may depend on many parameters, including locale settings. It seemed to me unwise to go this way. You could of course use some default locale, and make all the settings in accordance with the objective of making this conversion, but conceptually I did not like the idea that you relie on functions whose objective is human/machine *interfacing* in order to make *computation*. Interfacing and Inner guts should be cleanly separate in any good software. These kind of function may depend on 3rd parties libraries from the compiling chain, there would certainly be porting issues, ...


    Vincent.


 
>         Stefan
 		 	   		  
_________________________________________________________________
Nouveau Windows 7 : Trouvez le PC qui vous convient. En savoir plus.
http://clk.atdmt.com/FRM/go/181574580/direct/01/

[-- Attachment #2.2: Type: text/html, Size: 4850 bytes --]

[-- Attachment #3: mathfloat.el.gz --]
[-- Type: application/x-gzip, Size: 4372 bytes --]

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

>> Also, the mantissa should probably be represented as either a single

>> integer or a list of integers where each integer provides 16bits of data
>> (that's a format already used elsewhere in Emacs to represent "large
>> integers" such as for time).

> I did it this way because it allowed a simpler algorithm: in the
> conversion algorithm the integer size step is based on the number of
> decimal digits that Calc uses in representing big integers. The Calc
> big integers are not composed of 16bits integers (and actually there
> is not an integral number of bits, because these components are in
> radix 10). From that perspective it was easier if the number of bits
> of each component integer is passed along with the component integer
> in the interface.

I understand the code is simpler on the Lisp side, but I'm more
concerned about the C side.  Of course, if the Lisp side becomes
unmanageable it's also relevant.

>> I'd also rename the mantissa and exponent size info, or maybe even
>> consider providing it differently (in case its use is always linked to
>> calls to one of the previous 2 functions, maybe those functions could
>> return the relevant info.  E.g. make-float could return the part of the
>> mantissa it ignored).  Since you're using those constants, you clearly
>> know better how they're used, so I'll let you figure out whether there's
>> a better way to ptovide the same info.

> OK, this can be improved. Actually, my intention was not that these
> builtin functions were kind-of part of Calc library API, these are
> really internal things that the user developing above Calc should
> not need.

I understand they're not Calc-specific althought they'd currently only
be used by Calc.  That's good.  Still, Calc is the best sample-use we
have, so it's what we have to guess at what a good API should look like.

>> Also, how is math-lisp-float-binary-ieee754-conformance used?

> This is the switch between the two types of implementations: the one
> using builtin support, and the one not needing builtin support.

So (fboundp 'construct-float) would work as well?

> Note that I can imagine that it is possible to implement the
> construct/deconstruct builtin functions just using Lisp, this is after
> all just a matter of generating clean powers of 2.  But, since this
> code would be so machine/system dependant, and since Lisp would be so
> more uselessly complex than C code,I thought that it was a better
> choice to make it in C.

I think it fits better in C.  Basically, it would be the equivalent to
C's frexp.


        Stefan

RE: Code for converting between Elisp and Calc floats

[Vincent Belaïche]

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














Hello Stefan,


> [SNIP]

> 
> I understand the code is simpler on the Lisp side, but I'm more
> concerned about the C side.  Of course, if the Lisp side becomes
> unmanageable it's also relevant.
>

For the C code anyhow, even if the integer component word size was not passed along in the interface, you would have anyhow to manage strange word sizes. This is because the mantissa is 53 bits, and 53 is not a multiple of 16. Also the mantissa is not aligned on the double word boundary (the sign + exponent not being a muliple of 16 bits) and
you also have this phantom bit case to handle.

I agree however that the C could have been slightly simpler if the integer component word size was forced to be 16bits.

Maybe we could add one more layer in Lisp that would encapsulate the C function and add the word lengths to be 16 for those who would like to use the builtin function for other purpose, but with a simpler interface. This would basically achieve the interface that you envisionned as far as I can understand. Calc would however not use this additional layer of code.


> [SNIP] 
> > This is the switch between the two types of implementations: the one
> > using builtin support, and the one not needing builtin support.
> 
> So (fboundp 'construct-float) would work as well?
>

Yes, definitely. That could be done this way, the  reason why I used a specific variable is more historical than rational.


> [SNIP] 
> I think it fits better in C. Basically, it would be the equivalent to
> C's frexp.

C++' frexp and ldexp are not a direct alternative, because the significand is still a non integer number (between 0.5 and 1 for frexp). So anyhow you would have to build clean powers of two to convert from the frexp returned significand to Calc internal mantissa that is an integer. Also the dynamic of the mantissa is higher than that of a Lisp integer, so any how you have so stick/split bits together/apart. However by combining calls to frexp, ldexp and floor one could get around.


Frankly speaking, I must admit that I ignored the existance of frexp and ldexp standard functions. If I had known about them, maybe I would have tried to write the builtin functions with using only frexp ldexp and stuff like _isnan and suhclikes to handle special cases. That would have come at the cost of a slight additional processing cost, but this is completely neglectable in the case of Calc, due to the overhead of Lisp.

Another alternative would be that frexp and ldexp and suhlikes are made builtin functions, and that the construct-float, deconstruct-float functions are in Lisp using those builtins: that would indeed be my preference if everyone thinks that the submitted C code is too complex. 


I am not sure however that re-coding those functions with frexp, ldexp and floor, would be that significantly less complex to justify the effort, now that I have already done the job in a more primitive way.

     Vincent.



 		 	   		  
_________________________________________________________________
Nouveau ! Tout Windows débarque dans votre téléphone. Voir les Windows phone
http://clk.atdmt.com/FRM/go/175819071/direct/01/

[-- Attachment #2: Type: text/html, Size: 3742 bytes --]

Re: Code for converting between Elisp and Calc floats

[Stefan Monnier]

> For the C code anyhow, even if the integer component word size was not
> passed along in the interface, you would have anyhow to manage strange
> word sizes. This is because the mantissa is 53 bits, and 53 is not
> a multiple of 16. Also the mantissa is not aligned on the double word
> boundary (the sign + exponent not being a muliple of 16 bits) and you
> also have this phantom bit case to handle.

I think if you use frexp+floor, most of those issues won't matter much.

>> I think it fits better in C. Basically, it would be the equivalent to
>> C's frexp.
> C++' frexp and ldexp are not a direct alternative,

I understand.  The interface you proposed is closer to what we'd like
than what frexp has to offer.

> Another alternative would be that frexp and ldexp and suhlikes are
> made builtin functions, and that the construct-float,
> deconstruct-float functions are in Lisp using those builtins: that
> would indeed be my preference if everyone thinks that the submitted
> C code is too complex. 

Actually my preference would be to keep the functionality you proposed
(i.e. two functions (de)construct-float written in C), except force the
base to be 65536 (i.e. 16bits per integer), and rewrite the C code to
use standard functions like frexp, floor, isnan, and friends.


        Stefan

RE: Code for converting between Elisp and Calc floats

[Vincent Belaïche]

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

 <jwvocnte5aa.fsf-monnier+emacs@gnu.org>
MIME-Version: 1.0

[-- Attachment #2: Type: text/plain, Size: 789 bytes --]



> Actually my preference would be to keep the functionality you proposed
> (i.e. two functions (de)construct-float written in C), except force the
> base to be 65536 (i.e. 16bits per integer), and rewrite the C code to
> use standard functions like frexp, floor, isnan, and friends.
>

Then I need some fonction to map the list of words with non-always-16-bit width to a list of words with 16-bit width. I would do this in Lisp with some use of `lsh' and `+', so that the Lisp part of my contribution could still use the same interface.


Regards,
    Vincent.







 		 	   		  
_________________________________________________________________
Nouveau Windows 7 : Trouvez le PC qui vous convient. En savoir plus.
http://clk.atdmt.com/FRM/go/181574580/direct/01/

[-- Attachment #3: Type: text/html, Size: 1010 bytes --]