[thirdparty/gcc.git] / gcc / sreal.c

/* Simple data type for positive real numbers for the GNU compiler.
   Copyright (C) 2002, 2003 Free Software Foundation, Inc.

This file is part of GCC.

GCC is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free
Software Foundation; either version 2, or (at your option) any later
version.

GCC is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with GCC; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA.  */

/* This library supports positive real numbers and 0;
   inf and nan are NOT supported.
   It is written to be simple and fast.

   Value of sreal is
	x = sig * 2 ^ exp
   where
	sig = significant
	  (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS)
	exp = exponent

   One HOST_WIDE_INT is used for the significant on 64-bit (and more than
   64-bit) machines,
   otherwise two HOST_WIDE_INTs are used for the significant.
   Only a half of significant bits is used (in normalized sreals) so that we do
   not have problems with overflow, for example when c->sig = a->sig * b->sig.
   So the precision for 64-bit and 32-bit machines is 32-bit.

   Invariant: The numbers are normalized before and after each call of sreal_*.

   Normalized sreals:
   All numbers (except zero) meet following conditions:
	 SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG
	-SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP

   If the number would be too large, it is set to upper bounds of these
   conditions.

   If the number is zero or would be too small it meets following conditions:
	sig == 0 && exp == -SREAL_MAX_EXP
*/

#include "config.h"
#include "system.h"
#include "coretypes.h"
#include "tm.h"
#include "sreal.h"

static inline void copy (sreal *, sreal *);
static inline void shift_right (sreal *, int);
static void normalize (sreal *);

/* Print the content of struct sreal.  */

void
dump_sreal (FILE *file, sreal *x)
{
#if SREAL_PART_BITS < 32
  fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
	   HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
	   x->sig_hi, x->sig_lo, x->exp);
#else
  fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
#endif
}

/* Copy the sreal number.  */

static inline void
copy (sreal *r, sreal *a)
{
#if SREAL_PART_BITS < 32
  r->sig_lo = a->sig_lo;
  r->sig_hi = a->sig_hi;
#else
  r->sig = a->sig;
#endif
  r->exp = a->exp;
}

/* Shift X right by S bits.  Needed: 0 < S <= SREAL_BITS.
   When the most significant bit shifted out is 1, add 1 to X (rounding).  */

static inline void
shift_right (sreal *x, int s)
{
#ifdef ENABLE_CHECKING
  if (s <= 0 || s > SREAL_BITS)
    abort ();
  if (x->exp + s > SREAL_MAX_EXP)
    {
      /* Exponent should never be so large because shift_right is used only by
	 sreal_add and sreal_sub ant thus the number cannot be shifted out from
	 exponent range.  */
      abort ();
    }
#endif

  x->exp += s;

#if SREAL_PART_BITS < 32
  if (s > SREAL_PART_BITS)
    {
      s -= SREAL_PART_BITS;
      x->sig_hi += (uhwi) 1 << (s - 1);
      x->sig_lo = x->sig_hi >> s;
      x->sig_hi = 0;
    }
  else
    {
      x->sig_lo += (uhwi) 1 << (s - 1);
      if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
	{
	  x->sig_hi++;
	  x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
	}
      x->sig_lo >>= s;
      x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
      x->sig_hi >>= s;
    }
#else
  x->sig += (uhwi) 1 << (s - 1);
  x->sig >>= s;
#endif
}

/* Normalize *X.  */

static void
normalize (sreal *x)
{
#if SREAL_PART_BITS < 32
  int shift;
  HOST_WIDE_INT mask;

  if (x->sig_lo == 0 && x->sig_hi == 0)
    {
      x->exp = -SREAL_MAX_EXP;
    }
  else if (x->sig_hi < SREAL_MIN_SIG)
    {
      if (x->sig_hi == 0)
	{
	  /* Move lower part of significant to higher part.  */
	  x->sig_hi = x->sig_lo;
	  x->sig_lo = 0;
	  x->exp -= SREAL_PART_BITS;
	}
      shift = 0;
      while (x->sig_hi < SREAL_MIN_SIG)
	{
	  x->sig_hi <<= 1;
	  x->exp--;
	  shift++;
	}
      /* Check underflow.  */
      if (x->exp < -SREAL_MAX_EXP)
	{
	  x->exp = -SREAL_MAX_EXP;
	  x->sig_hi = 0;
	  x->sig_lo = 0;
	}
      else if (shift)
	{
	  mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
	  x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
	  x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
	}
    }
  else if (x->sig_hi > SREAL_MAX_SIG)
    {
      unsigned HOST_WIDE_INT tmp = x->sig_hi;

      /* Find out how many bits will be shifted.  */
      shift = 0;
      do
	{
	  tmp >>= 1;
	  shift++;
	}
      while (tmp > SREAL_MAX_SIG);

      /* Round the number.  */
      x->sig_lo += (uhwi) 1 << (shift - 1);

      x->sig_lo >>= shift;
      x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
		    << (SREAL_PART_BITS - shift));
      x->sig_hi >>= shift;
      x->exp += shift;
      if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
	{
	  x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
	  x->sig_hi++;
	  if (x->sig_hi > SREAL_MAX_SIG)
	    {
	      /* x->sig_hi was SREAL_MAX_SIG before increment
		 so now last bit is zero.  */
	      x->sig_hi >>= 1;
	      x->sig_lo >>= 1;
	      x->exp++;
	    }
	}

      /* Check overflow.  */
      if (x->exp > SREAL_MAX_EXP)
	{
	  x->exp = SREAL_MAX_EXP;
	  x->sig_hi = SREAL_MAX_SIG;
	  x->sig_lo = SREAL_MAX_SIG;
	}
    }
#else
  if (x->sig == 0)
    {
      x->exp = -SREAL_MAX_EXP;
    }
  else if (x->sig < SREAL_MIN_SIG)
    {
      do
	{
	  x->sig <<= 1;
	  x->exp--;
	}
      while (x->sig < SREAL_MIN_SIG);

      /* Check underflow.  */
      if (x->exp < -SREAL_MAX_EXP)
	{
	  x->exp = -SREAL_MAX_EXP;
	  x->sig = 0;
	}
    }
  else if (x->sig > SREAL_MAX_SIG)
    {
      int last_bit;
      do
	{
	  last_bit = x->sig & 1;
	  x->sig >>= 1;
	  x->exp++;
	}
      while (x->sig > SREAL_MAX_SIG);

      /* Round the number.  */
      x->sig += last_bit;
      if (x->sig > SREAL_MAX_SIG)
	{
	  x->sig >>= 1;
	  x->exp++;
	}

      /* Check overflow.  */
      if (x->exp > SREAL_MAX_EXP)
	{
	  x->exp = SREAL_MAX_EXP;
	  x->sig = SREAL_MAX_SIG;
	}
    }
#endif
}

/* Set *R to SIG * 2 ^ EXP.  Return R.  */

sreal *
sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
{
#if SREAL_PART_BITS < 32
  r->sig_lo = 0;
  r->sig_hi = sig;
  r->exp = exp - 16;
#else
  r->sig = sig;
  r->exp = exp;
#endif
  normalize (r);
  return r;
}

/* Return integer value of *R.  */

HOST_WIDE_INT
sreal_to_int (sreal *r)
{
#if SREAL_PART_BITS < 32
  if (r->exp <= -SREAL_BITS)
    return 0;
  if (r->exp >= 0)
    return MAX_HOST_WIDE_INT;
  return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
#else
  if (r->exp <= -SREAL_BITS)
    return 0;
  if (r->exp >= SREAL_PART_BITS)
    return MAX_HOST_WIDE_INT;
  if (r->exp > 0)
    return r->sig << r->exp;
  if (r->exp < 0)
    return r->sig >> -r->exp;
  return r->sig;
#endif
}

/* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B.  */

int
sreal_compare (sreal *a, sreal *b)
{
  if (a->exp > b->exp)
    return 1;
  if (a->exp < b->exp)
    return -1;
#if SREAL_PART_BITS < 32
  if (a->sig_hi > b->sig_hi)
    return 1;
  if (a->sig_hi < b->sig_hi)
    return -1;
  if (a->sig_lo > b->sig_lo)
    return 1;
  if (a->sig_lo < b->sig_lo)
    return -1;
#else
  if (a->sig > b->sig)
    return 1;
  if (a->sig < b->sig)
    return -1;
#endif
  return 0;
}

/* *R = *A + *B.  Return R.  */

sreal *
sreal_add (sreal *r, sreal *a, sreal *b)
{
  int dexp;
  sreal tmp;
  sreal *bb;

  if (sreal_compare (a, b) < 0)
    {
      sreal *swap;
      swap = a;
      a = b;
      b = swap;
    }

  dexp = a->exp - b->exp;
  r->exp = a->exp;
  if (dexp > SREAL_BITS)
    {
#if SREAL_PART_BITS < 32
      r->sig_hi = a->sig_hi;
      r->sig_lo = a->sig_lo;
#else
      r->sig = a->sig;
#endif
      return r;
    }

  if (dexp == 0)
    bb = b;
  else
    {
      copy (&tmp, b);
      shift_right (&tmp, dexp);
      bb = &tmp;
    }

#if SREAL_PART_BITS < 32
  r->sig_hi = a->sig_hi + bb->sig_hi;
  r->sig_lo = a->sig_lo + bb->sig_lo;
  if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
    {
      r->sig_hi++;
      r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
    }
#else
  r->sig = a->sig + bb->sig;
#endif
  normalize (r);
  return r;
}

/* *R = *A - *B.  Return R.  */

sreal *
sreal_sub (sreal *r, sreal *a, sreal *b)
{
  int dexp;
  sreal tmp;
  sreal *bb;

  if (sreal_compare (a, b) < 0)
    {
      abort ();
    }

  dexp = a->exp - b->exp;
  r->exp = a->exp;
  if (dexp > SREAL_BITS)
    {
#if SREAL_PART_BITS < 32
      r->sig_hi = a->sig_hi;
      r->sig_lo = a->sig_lo;
#else
      r->sig = a->sig;
#endif
      return r;
    }
  if (dexp == 0)
    bb = b;
  else
    {
      copy (&tmp, b);
      shift_right (&tmp, dexp);
      bb = &tmp;
    }

#if SREAL_PART_BITS < 32
  if (a->sig_lo < bb->sig_lo)
    {
      r->sig_hi = a->sig_hi - bb->sig_hi - 1;
      r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
    }
  else
    {
      r->sig_hi = a->sig_hi - bb->sig_hi;
      r->sig_lo = a->sig_lo - bb->sig_lo;
    }
#else
  r->sig = a->sig - bb->sig;
#endif
  normalize (r);
  return r;
}

/* *R = *A * *B.  Return R.  */

sreal *
sreal_mul (sreal *r, sreal *a, sreal *b)
{
#if SREAL_PART_BITS < 32
  if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
    {
      r->sig_lo = 0;
      r->sig_hi = 0;
      r->exp = -SREAL_MAX_EXP;
    }
  else
    {
      unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
      if (sreal_compare (a, b) < 0)
	{
	  sreal *swap;
	  swap = a;
	  a = b;
	  b = swap;
	}

      r->exp = a->exp + b->exp + SREAL_PART_BITS;

      tmp1 = a->sig_lo * b->sig_lo;
      tmp2 = a->sig_lo * b->sig_hi;
      tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);

      r->sig_hi = a->sig_hi * b->sig_hi;
      r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
      tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
      tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
      tmp1 = tmp2 + tmp3;

      r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
      r->sig_hi += tmp1 >> SREAL_PART_BITS;

      normalize (r);
    }
#else
  if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
    {
      r->sig = 0;
      r->exp = -SREAL_MAX_EXP;
    }
  else
    {
      r->sig = a->sig * b->sig;
      r->exp = a->exp + b->exp;
      normalize (r);
    }
#endif
  return r;
}

/* *R = *A / *B.  Return R.  */

sreal *
sreal_div (sreal *r, sreal *a, sreal *b)
{
#if SREAL_PART_BITS < 32
  unsigned HOST_WIDE_INT tmp, tmp1, tmp2;

  if (b->sig_hi < SREAL_MIN_SIG)
    {
      abort ();
    }
  else if (a->sig_hi < SREAL_MIN_SIG)
    {
      r->sig_hi = 0;
      r->sig_lo = 0;
      r->exp = -SREAL_MAX_EXP;
    }
  else
    {
      /* Since division by the whole number is pretty ugly to write
	 we are dividing by first 3/4 of bits of number.  */

      tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
      tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
	      + (b->sig_lo >> (SREAL_PART_BITS / 2)));
      if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
	tmp2++;

      r->sig_lo = 0;
      tmp = tmp1 / tmp2;
      tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
      r->sig_hi = tmp << SREAL_PART_BITS;

      tmp = tmp1 / tmp2;
      tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
      r->sig_hi += tmp << (SREAL_PART_BITS / 2);

      tmp = tmp1 / tmp2;
      r->sig_hi += tmp;

      r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
      normalize (r);
    }
#else
  if (b->sig == 0)
    {
      abort ();
    }
  else
    {
      r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
      r->exp = a->exp - b->exp - SREAL_PART_BITS;
      normalize (r);
    }
#endif
  return r;
}
Commit	Line	Data
ac5e69da JZ	1	/* Simple data type for positive real numbers for the GNU compiler.
	2	Copyright (C) 2002, 2003 Free Software Foundation, Inc.
	3
	4	This file is part of GCC.
	5
	6	GCC is free software; you can redistribute it and/or modify it under
	7	the terms of the GNU General Public License as published by the Free
	8	Software Foundation; either version 2, or (at your option) any later
	9	version.
	10
	11	GCC is distributed in the hope that it will be useful, but WITHOUT ANY
	12	WARRANTY; without even the implied warranty of MERCHANTABILITY or
	13	FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
	14	for more details.
	15
	16	You should have received a copy of the GNU General Public License
	17	along with GCC; see the file COPYING. If not, write to the Free
	18	Software Foundation, 59 Temple Place - Suite 330, Boston, MA
	19	02111-1307, USA. */
	20
	21	/* This library supports positive real numbers and 0;
	22	inf and nan are NOT supported.
	23	It is written to be simple and fast.
	24
	25	Value of sreal is
	26	x = sig * 2 ^ exp
46c5ad27	27	where
ac5e69da JZ	28	sig = significant
	29	(for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS)
	30	exp = exponent
	31
	32	One HOST_WIDE_INT is used for the significant on 64-bit (and more than
	33	64-bit) machines,
	34	otherwise two HOST_WIDE_INTs are used for the significant.
	35	Only a half of significant bits is used (in normalized sreals) so that we do
	36	not have problems with overflow, for example when c->sig = a->sig * b->sig.
e0bb17a8	37	So the precision for 64-bit and 32-bit machines is 32-bit.
46c5ad27	38
ac5e69da JZ	39	Invariant: The numbers are normalized before and after each call of sreal_*.
	40
	41	Normalized sreals:
	42	All numbers (except zero) meet following conditions:
	43	SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG
46c5ad27	44	-SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP
ac5e69da JZ	45
	46	If the number would be too large, it is set to upper bounds of these
	47	conditions.
	48
	49	If the number is zero or would be too small it meets following conditions:
	50	sig == 0 && exp == -SREAL_MAX_EXP
	51	*/
	52
	53	#include "config.h"
	54	#include "system.h"
	55	#include "coretypes.h"
	56	#include "tm.h"
	57	#include "sreal.h"
	58
46c5ad27 AJ	59	static inline void copy (sreal , sreal );
	60	static inline void shift_right (sreal *, int);
	61	static void normalize (sreal *);
ac5e69da JZ	62
	63	/* Print the content of struct sreal. */
	64
	65	void
46c5ad27	66	dump_sreal (FILE file, sreal x)
ac5e69da JZ	67	{
ac5e69da JZ	68	#if SREAL_PART_BITS < 32
90ff44cf KG	69	fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
	70	HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
	71	x->sig_hi, x->sig_lo, x->exp);
ac5e69da	72	#else
90ff44cf	73	fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
ac5e69da JZ	74	#endif
	75	}
	76
	77	/* Copy the sreal number. */
	78
	79	static inline void
46c5ad27	80	copy (sreal r, sreal a)
ac5e69da JZ	81	{
	82	#if SREAL_PART_BITS < 32
	83	r->sig_lo = a->sig_lo;
	84	r->sig_hi = a->sig_hi;
	85	#else
	86	r->sig = a->sig;
	87	#endif
	88	r->exp = a->exp;
	89	}
	90
	91	/* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
	92	When the most significant bit shifted out is 1, add 1 to X (rounding). */
	93
	94	static inline void
46c5ad27	95	shift_right (sreal *x, int s)
ac5e69da JZ	96	{
	97	#ifdef ENABLE_CHECKING
	98	if (s <= 0 \|\| s > SREAL_BITS)
	99	abort ();
	100	if (x->exp + s > SREAL_MAX_EXP)
	101	{
	102	/* Exponent should never be so large because shift_right is used only by
	103	sreal_add and sreal_sub ant thus the number cannot be shifted out from
	104	exponent range. */
	105	abort ();
	106	}
	107	#endif
	108
	109	x->exp += s;
	110
	111	#if SREAL_PART_BITS < 32
	112	if (s > SREAL_PART_BITS)
	113	{
	114	s -= SREAL_PART_BITS;
	115	x->sig_hi += (uhwi) 1 << (s - 1);
	116	x->sig_lo = x->sig_hi >> s;
	117	x->sig_hi = 0;
	118	}
	119	else
	120	{
	121	x->sig_lo += (uhwi) 1 << (s - 1);
	122	if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
	123	{
	124	x->sig_hi++;
	125	x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
	126	}
	127	x->sig_lo >>= s;
	128	x->sig_lo \|= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
	129	x->sig_hi >>= s;
	130	}
	131	#else
	132	x->sig += (uhwi) 1 << (s - 1);
	133	x->sig >>= s;
	134	#endif
	135	}
	136
	137	/* Normalize X. /
	138
	139	static void
46c5ad27	140	normalize (sreal *x)
ac5e69da JZ	141	{
	142	#if SREAL_PART_BITS < 32
	143	int shift;
	144	HOST_WIDE_INT mask;
46c5ad27	145
ac5e69da JZ	146	if (x->sig_lo == 0 && x->sig_hi == 0)
	147	{
	148	x->exp = -SREAL_MAX_EXP;
	149	}
	150	else if (x->sig_hi < SREAL_MIN_SIG)
	151	{
	152	if (x->sig_hi == 0)
	153	{
	154	/* Move lower part of significant to higher part. */
	155	x->sig_hi = x->sig_lo;
	156	x->sig_lo = 0;
	157	x->exp -= SREAL_PART_BITS;
	158	}
	159	shift = 0;
	160	while (x->sig_hi < SREAL_MIN_SIG)
	161	{
	162	x->sig_hi <<= 1;
	163	x->exp--;
	164	shift++;
	165	}
	166	/* Check underflow. */
	167	if (x->exp < -SREAL_MAX_EXP)
	168	{
	169	x->exp = -SREAL_MAX_EXP;
	170	x->sig_hi = 0;
	171	x->sig_lo = 0;
	172	}
	173	else if (shift)
	174	{
	175	mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
	176	x->sig_hi \|= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
	177	x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
	178	}
	179	}
	180	else if (x->sig_hi > SREAL_MAX_SIG)
	181	{
	182	unsigned HOST_WIDE_INT tmp = x->sig_hi;
	183
	184	/* Find out how many bits will be shifted. */
	185	shift = 0;
	186	do
	187	{
	188	tmp >>= 1;
	189	shift++;
	190	}
	191	while (tmp > SREAL_MAX_SIG);
	192
	193	/* Round the number. */
	194	x->sig_lo += (uhwi) 1 << (shift - 1);
	195
	196	x->sig_lo >>= shift;
	197	x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
	198	<< (SREAL_PART_BITS - shift));
	199	x->sig_hi >>= shift;
	200	x->exp += shift;
	201	if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
	202	{
	203	x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
	204	x->sig_hi++;
	205	if (x->sig_hi > SREAL_MAX_SIG)
	206	{
	207	/* x->sig_hi was SREAL_MAX_SIG before increment
	208	so now last bit is zero. */
	209	x->sig_hi >>= 1;
210	x->sig_lo >>= 1;
211	x->exp++;
212	}
213	}
214
215	/* Check overflow. */
216	if (x->exp > SREAL_MAX_EXP)
217	{
218	x->exp = SREAL_MAX_EXP;
219	x->sig_hi = SREAL_MAX_SIG;
220	x->sig_lo = SREAL_MAX_SIG;
221	}
222	}
223	#else
224	if (x->sig == 0)
225	{
226	x->exp = -SREAL_MAX_EXP;
227	}
228	else if (x->sig < SREAL_MIN_SIG)
229	{
230	do
231	{
232	x->sig <<= 1;
233	x->exp--;
234	}
235	while (x->sig < SREAL_MIN_SIG);
236
237	/* Check underflow. */
238	if (x->exp < -SREAL_MAX_EXP)
239	{
240	x->exp = -SREAL_MAX_EXP;
241	x->sig = 0;
242	}
243	}
244	else if (x->sig > SREAL_MAX_SIG)
245	{
246	int last_bit;
247	do
248	{
249	last_bit = x->sig & 1;
250	x->sig >>= 1;
251	x->exp++;
252	}
253	while (x->sig > SREAL_MAX_SIG);
254
255	/* Round the number. */
256	x->sig += last_bit;
257	if (x->sig > SREAL_MAX_SIG)
258	{
259	x->sig >>= 1;
260	x->exp++;
261	}
262
263	/* Check overflow. */
264	if (x->exp > SREAL_MAX_EXP)
265	{
266	x->exp = SREAL_MAX_EXP;
267	x->sig = SREAL_MAX_SIG;
268	}
269	}
270	#endif
271	}
272
273	/* Set R to SIG 2 ^ EXP. Return R. */
274
275	sreal *
46c5ad27	276	sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
ac5e69da JZ	277	{
	278	#if SREAL_PART_BITS < 32
	279	r->sig_lo = 0;
	280	r->sig_hi = sig;
	281	r->exp = exp - 16;
	282	#else
	283	r->sig = sig;
	284	r->exp = exp;
	285	#endif
	286	normalize (r);
	287	return r;
	288	}
	289
	290	/* Return integer value of R. /
	291
	292	HOST_WIDE_INT
46c5ad27	293	sreal_to_int (sreal *r)
ac5e69da JZ	294	{
	295	#if SREAL_PART_BITS < 32
	296	if (r->exp <= -SREAL_BITS)
	297	return 0;
	298	if (r->exp >= 0)
	299	return MAX_HOST_WIDE_INT;
	300	return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
	301	#else
	302	if (r->exp <= -SREAL_BITS)
	303	return 0;
	304	if (r->exp >= SREAL_PART_BITS)
	305	return MAX_HOST_WIDE_INT;
	306	if (r->exp > 0)
	307	return r->sig << r->exp;
	308	if (r->exp < 0)
	309	return r->sig >> -r->exp;
	310	return r->sig;
	311	#endif
	312	}
	313
	314	/* Compare A and B. Return -1 if A < B, 1 if A > B and 0 if A == B. */
	315
	316	int
46c5ad27	317	sreal_compare (sreal a, sreal b)
ac5e69da JZ	318	{
	319	if (a->exp > b->exp)
	320	return 1;
	321	if (a->exp < b->exp)
	322	return -1;
	323	#if SREAL_PART_BITS < 32
	324	if (a->sig_hi > b->sig_hi)
	325	return 1;
	326	if (a->sig_hi < b->sig_hi)
	327	return -1;
	328	if (a->sig_lo > b->sig_lo)
	329	return 1;
	330	if (a->sig_lo < b->sig_lo)
	331	return -1;
	332	#else
	333	if (a->sig > b->sig)
	334	return 1;
	335	if (a->sig < b->sig)
	336	return -1;
	337	#endif
	338	return 0;
	339	}
	340
	341	/* R = A + B. Return R. /
	342
	343	sreal *
46c5ad27	344	sreal_add (sreal r, sreal a, sreal *b)
ac5e69da JZ	345	{
	346	int dexp;
	347	sreal tmp;
	348	sreal *bb;
	349
	350	if (sreal_compare (a, b) < 0)
	351	{
	352	sreal *swap;
	353	swap = a;
	354	a = b;
	355	b = swap;
	356	}
	357
	358	dexp = a->exp - b->exp;
	359	r->exp = a->exp;
	360	if (dexp > SREAL_BITS)
	361	{
	362	#if SREAL_PART_BITS < 32
	363	r->sig_hi = a->sig_hi;
	364	r->sig_lo = a->sig_lo;
	365	#else
	366	r->sig = a->sig;
	367	#endif
	368	return r;
	369	}
	370
	371	if (dexp == 0)
	372	bb = b;
	373	else
	374	{
	375	copy (&tmp, b);
	376	shift_right (&tmp, dexp);
	377	bb = &tmp;
	378	}
	379
	380	#if SREAL_PART_BITS < 32
	381	r->sig_hi = a->sig_hi + bb->sig_hi;
	382	r->sig_lo = a->sig_lo + bb->sig_lo;
	383	if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
	384	{
	385	r->sig_hi++;
	386	r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
	387	}
	388	#else
	389	r->sig = a->sig + bb->sig;
	390	#endif
	391	normalize (r);
	392	return r;
	393	}
	394
	395	/* R = A - B. Return R. /
	396
	397	sreal *
46c5ad27	398	sreal_sub (sreal r, sreal a, sreal *b)
ac5e69da JZ	399	{
	400	int dexp;
	401	sreal tmp;
	402	sreal *bb;
	403
	404	if (sreal_compare (a, b) < 0)
	405	{
	406	abort ();
	407	}
	408
	409	dexp = a->exp - b->exp;
	410	r->exp = a->exp;
	411	if (dexp > SREAL_BITS)
	412	{
	413	#if SREAL_PART_BITS < 32
	414	r->sig_hi = a->sig_hi;
	415	r->sig_lo = a->sig_lo;
	416	#else
	417	r->sig = a->sig;
	418	#endif
	419	return r;
	420	}
	421	if (dexp == 0)
	422	bb = b;
	423	else
	424	{
	425	copy (&tmp, b);
	426	shift_right (&tmp, dexp);
	427	bb = &tmp;
	428	}
	429
	430	#if SREAL_PART_BITS < 32
	431	if (a->sig_lo < bb->sig_lo)
	432	{
	433	r->sig_hi = a->sig_hi - bb->sig_hi - 1;
	434	r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
	435	}
	436	else
	437	{
	438	r->sig_hi = a->sig_hi - bb->sig_hi;
	439	r->sig_lo = a->sig_lo - bb->sig_lo;
	440	}
	441	#else
	442	r->sig = a->sig - bb->sig;
	443	#endif
	444	normalize (r);
	445	return r;
	446	}
	447
	448	/* R = A * B. Return R. /
	449
	450	sreal *
46c5ad27	451	sreal_mul (sreal r, sreal a, sreal *b)
ac5e69da JZ	452	{
	453	#if SREAL_PART_BITS < 32
	454	if (a->sig_hi < SREAL_MIN_SIG \|\| b->sig_hi < SREAL_MIN_SIG)
	455	{
	456	r->sig_lo = 0;
	457	r->sig_hi = 0;
	458	r->exp = -SREAL_MAX_EXP;
	459	}
	460	else
	461	{
	462	unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
	463	if (sreal_compare (a, b) < 0)
	464	{
	465	sreal *swap;
	466	swap = a;
	467	a = b;
	468	b = swap;
	469	}
	470
	471	r->exp = a->exp + b->exp + SREAL_PART_BITS;
	472
	473	tmp1 = a->sig_lo * b->sig_lo;
	474	tmp2 = a->sig_lo * b->sig_hi;
	475	tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
	476
	477	r->sig_hi = a->sig_hi * b->sig_hi;
	478	r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
	479	tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
	480	tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
	481	tmp1 = tmp2 + tmp3;
	482
	483	r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
	484	r->sig_hi += tmp1 >> SREAL_PART_BITS;
	485
	486	normalize (r);
	487	}
	488	#else
	489	if (a->sig < SREAL_MIN_SIG \|\| b->sig < SREAL_MIN_SIG)
	490	{
	491	r->sig = 0;
	492	r->exp = -SREAL_MAX_EXP;
	493	}
	494	else
	495	{
	496	r->sig = a->sig * b->sig;
	497	r->exp = a->exp + b->exp;
	498	normalize (r);
	499	}
	500	#endif
	501	return r;
	502	}
	503
	504	/* R = A / B. Return R. /
	505
	506	sreal *
46c5ad27	507	sreal_div (sreal r, sreal a, sreal *b)
ac5e69da JZ	508	{
	509	#if SREAL_PART_BITS < 32
	510	unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
	511
	512	if (b->sig_hi < SREAL_MIN_SIG)
	513	{
	514	abort ();
	515	}
	516	else if (a->sig_hi < SREAL_MIN_SIG)
	517	{
	518	r->sig_hi = 0;
	519	r->sig_lo = 0;
	520	r->exp = -SREAL_MAX_EXP;
	521	}
	522	else
	523	{
	524	/* Since division by the whole number is pretty ugly to write
	525	we are dividing by first 3/4 of bits of number. */
	526
	527	tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
	528	tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
	529	+ (b->sig_lo >> (SREAL_PART_BITS / 2)));
	530	if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
	531	tmp2++;
	532
	533	r->sig_lo = 0;
	534	tmp = tmp1 / tmp2;
	535	tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
	536	r->sig_hi = tmp << SREAL_PART_BITS;
	537
	538	tmp = tmp1 / tmp2;
	539	tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
	540	r->sig_hi += tmp << (SREAL_PART_BITS / 2);
	541
	542	tmp = tmp1 / tmp2;
	543	r->sig_hi += tmp;
	544
	545	r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
	546	normalize (r);
	547	}
	548	#else
	549	if (b->sig == 0)
	550	{
	551	abort ();
	552	}
	553	else
	554	{
	555	r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
	556	r->exp = a->exp - b->exp - SREAL_PART_BITS;
	557	normalize (r);
	558	}
	559	#endif
	560	return r;
	561	}