URL https://opencores.org/ocsvn/scarts/scarts/trunk

Subversion Repositories scarts

[/] [scarts/] [trunk/] [toolchain/] [scarts-gcc/] [gcc-4.1.1/] [gcc/] [sreal.c] - Blame information for rev 12

Details | Compare with Previous | View Log


/* Simple data type for positive real numbers for the GNU compiler.
   Copyright (C) 2002, 2003, 2004 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, 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, 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)
{
  gcc_assert (s > 0);
  gcc_assert (s <= SREAL_BITS);
  /* 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.  */
  gcc_assert (x->exp + s <= SREAL_MAX_EXP);
 
  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;
 
  gcc_assert (sreal_compare (a, b) >= 0);
 
  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;
 
  gcc_assert (b->sig_hi >= SREAL_MIN_SIG);
  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
  gcc_assert (b->sig != 0);
  r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
  r->exp = a->exp - b->exp - SREAL_PART_BITS;
  normalize (r);
#endif
  return r;
}

Browse

Tools

Subversion Repositories scarts

[/] [scarts/] [trunk/] [toolchain/] [scarts-gcc/] [gcc-4.1.1/] [gcc/] [sreal.c] - Blame information for rev 12

Line No.	Rev	Author	Line
1	12	jlechner	`/* Simple data type for positive real numbers for the GNU compiler.`
2			`Copyright (C) 2002, 2003, 2004 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, 51 Franklin Street, Fifth Floor, Boston, MA`
19			`02110-1301, 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`
27			`where`
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.`
37			`So the precision for 64-bit and 32-bit machines is 32-bit.`
38
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`
44			`-SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP`
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
59			`static inline void copy (sreal , sreal );`
60			`static inline void shift_right (sreal *, int);`
61			`static void normalize (sreal *);`
62
63			`/* Print the content of struct sreal. */`
64
65			`void`
66			`dump_sreal (FILE file, sreal x)`
67			`{`
68			`#if SREAL_PART_BITS < 32`
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);`
72			`#else`
73			`fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);`
74			`#endif`
75			`}`
76
77			`/* Copy the sreal number. */`
78
79			`static inline void`
80			`copy (sreal r, sreal a)`
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`
95			`shift_right (sreal *x, int s)`
96			`{`
97			`gcc_assert (s > 0);`
98			`gcc_assert (s <= SREAL_BITS);`
99			`/* Exponent should never be so large because shift_right is used only by`
100			`sreal_add and sreal_sub ant thus the number cannot be shifted out from`
101			`exponent range. */`
102			`gcc_assert (x->exp + s <= SREAL_MAX_EXP);`
103
104			`x->exp += s;`
105
106			`#if SREAL_PART_BITS < 32`
107			`if (s > SREAL_PART_BITS)`
108			`{`
109			`s -= SREAL_PART_BITS;`
110			`x->sig_hi += (uhwi) 1 << (s - 1);`
111			`x->sig_lo = x->sig_hi >> s;`
112			`x->sig_hi = 0;`
113			`}`
114			`else`
115			`{`
116			`x->sig_lo += (uhwi) 1 << (s - 1);`
117			`if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))`
118			`{`
119			`x->sig_hi++;`
120			`x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;`
121			`}`
122			`x->sig_lo >>= s;`
123			`x->sig_lo \|= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);`
124			`x->sig_hi >>= s;`
125			`}`
126			`#else`
127			`x->sig += (uhwi) 1 << (s - 1);`
128			`x->sig >>= s;`
129			`#endif`
130			`}`
131
132			`/* Normalize X. /`
133
134			`static void`
135			`normalize (sreal *x)`
136			`{`
137			`#if SREAL_PART_BITS < 32`
138			`int shift;`
139			`HOST_WIDE_INT mask;`
140
141			`if (x->sig_lo == 0 && x->sig_hi == 0)`
142			`{`
143			`x->exp = -SREAL_MAX_EXP;`
144			`}`
145			`else if (x->sig_hi < SREAL_MIN_SIG)`
146			`{`
147			`if (x->sig_hi == 0)`
148			`{`
149			`/* Move lower part of significant to higher part. */`
150			`x->sig_hi = x->sig_lo;`
151			`x->sig_lo = 0;`
152			`x->exp -= SREAL_PART_BITS;`
153			`}`
154			`shift = 0;`
155			`while (x->sig_hi < SREAL_MIN_SIG)`
156			`{`
157			`x->sig_hi <<= 1;`
158			`x->exp--;`
159			`shift++;`
160			`}`
161			`/* Check underflow. */`
162			`if (x->exp < -SREAL_MAX_EXP)`
163			`{`
164			`x->exp = -SREAL_MAX_EXP;`
165			`x->sig_hi = 0;`
166			`x->sig_lo = 0;`
167			`}`
168			`else if (shift)`
169			`{`
170			`mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));`
171			`x->sig_hi \|= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);`
172			`x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);`
173			`}`
174			`}`
175			`else if (x->sig_hi > SREAL_MAX_SIG)`
176			`{`
177			`unsigned HOST_WIDE_INT tmp = x->sig_hi;`
178
179			`/* Find out how many bits will be shifted. */`
180			`shift = 0;`
181			`do`
182			`{`
183			`tmp >>= 1;`
184			`shift++;`
185			`}`
186			`while (tmp > SREAL_MAX_SIG);`
187
188			`/* Round the number. */`
189			`x->sig_lo += (uhwi) 1 << (shift - 1);`
190
191			`x->sig_lo >>= shift;`
192			`x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))`
193			`<< (SREAL_PART_BITS - shift));`
194			`x->sig_hi >>= shift;`
195			`x->exp += shift;`
196			`if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))`
197			`{`
198			`x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;`
199			`x->sig_hi++;`
200			`if (x->sig_hi > SREAL_MAX_SIG)`
201			`{`
202			`/* x->sig_hi was SREAL_MAX_SIG before increment`
203			`so now last bit is zero. */`
204			`x->sig_hi >>= 1;`
205			`x->sig_lo >>= 1;`
206			`x->exp++;`
207			`}`
208			`}`
209
210			`/* Check overflow. */`
211			`if (x->exp > SREAL_MAX_EXP)`
212			`{`
213			`x->exp = SREAL_MAX_EXP;`
214			`x->sig_hi = SREAL_MAX_SIG;`
215			`x->sig_lo = SREAL_MAX_SIG;`
216			`}`
217			`}`
218			`#else`
219			`if (x->sig == 0)`
220			`{`
221			`x->exp = -SREAL_MAX_EXP;`
222			`}`
223			`else if (x->sig < SREAL_MIN_SIG)`
224			`{`
225			`do`
226			`{`
227			`x->sig <<= 1;`
228			`x->exp--;`
229			`}`
230			`while (x->sig < SREAL_MIN_SIG);`
231
232			`/* Check underflow. */`
233			`if (x->exp < -SREAL_MAX_EXP)`
234			`{`
235			`x->exp = -SREAL_MAX_EXP;`
236			`x->sig = 0;`
237			`}`
238			`}`
239			`else if (x->sig > SREAL_MAX_SIG)`
240			`{`
241			`int last_bit;`
242			`do`
243			`{`
244			`last_bit = x->sig & 1;`
245			`x->sig >>= 1;`
246			`x->exp++;`
247			`}`
248			`while (x->sig > SREAL_MAX_SIG);`
249
250			`/* Round the number. */`
251			`x->sig += last_bit;`
252			`if (x->sig > SREAL_MAX_SIG)`
253			`{`
254			`x->sig >>= 1;`
255			`x->exp++;`
256			`}`
257
258			`/* Check overflow. */`
259			`if (x->exp > SREAL_MAX_EXP)`
260			`{`
261			`x->exp = SREAL_MAX_EXP;`
262			`x->sig = SREAL_MAX_SIG;`
263			`}`
264			`}`
265			`#endif`
266			`}`
267
268			`/* Set R to SIG 2 ^ EXP. Return R. */`
269
270			`sreal *`
271			`sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)`
272			`{`
273			`#if SREAL_PART_BITS < 32`
274			`r->sig_lo = 0;`
275			`r->sig_hi = sig;`
276			`r->exp = exp - 16;`
277			`#else`
278			`r->sig = sig;`
279			`r->exp = exp;`
280			`#endif`
281			`normalize (r);`
282			`return r;`
283			`}`
284
285			`/* Return integer value of R. /`
286
287			`HOST_WIDE_INT`
288			`sreal_to_int (sreal *r)`
289			`{`
290			`#if SREAL_PART_BITS < 32`
291			`if (r->exp <= -SREAL_BITS)`
292			`return 0;`
293			`if (r->exp >= 0)`
294			`return MAX_HOST_WIDE_INT;`
295			`return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;`
296			`#else`
297			`if (r->exp <= -SREAL_BITS)`
298			`return 0;`
299			`if (r->exp >= SREAL_PART_BITS)`
300			`return MAX_HOST_WIDE_INT;`
301			`if (r->exp > 0)`
302			`return r->sig << r->exp;`
303			`if (r->exp < 0)`
304			`return r->sig >> -r->exp;`
305			`return r->sig;`
306			`#endif`
307			`}`
308
309			`/* Compare A and B. Return -1 if A < B, 1 if A > B and 0 if A == B. */`
310
311			`int`
312			`sreal_compare (sreal a, sreal b)`
313			`{`
314			`if (a->exp > b->exp)`
315			`return 1;`
316			`if (a->exp < b->exp)`
317			`return -1;`
318			`#if SREAL_PART_BITS < 32`
319			`if (a->sig_hi > b->sig_hi)`
320			`return 1;`
321			`if (a->sig_hi < b->sig_hi)`
322			`return -1;`
323			`if (a->sig_lo > b->sig_lo)`
324			`return 1;`
325			`if (a->sig_lo < b->sig_lo)`
326			`return -1;`
327			`#else`
328			`if (a->sig > b->sig)`
329			`return 1;`
330			`if (a->sig < b->sig)`
331			`return -1;`
332			`#endif`
333			`return 0;`
334			`}`
335
336			`/* R = A + B. Return R. /`
337
338			`sreal *`
339			`sreal_add (sreal r, sreal a, sreal *b)`
340			`{`
341			`int dexp;`
342			`sreal tmp;`
343			`sreal *bb;`
344
345			`if (sreal_compare (a, b) < 0)`
346			`{`
347			`sreal *swap;`
348			`swap = a;`
349			`a = b;`
350			`b = swap;`
351			`}`
352
353			`dexp = a->exp - b->exp;`
354			`r->exp = a->exp;`
355			`if (dexp > SREAL_BITS)`
356			`{`
357			`#if SREAL_PART_BITS < 32`
358			`r->sig_hi = a->sig_hi;`
359			`r->sig_lo = a->sig_lo;`
360			`#else`
361			`r->sig = a->sig;`
362			`#endif`
363			`return r;`
364			`}`
365
366			`if (dexp == 0)`
367			`bb = b;`
368			`else`
369			`{`
370			`copy (&tmp, b);`
371			`shift_right (&tmp, dexp);`
372			`bb = &tmp;`
373			`}`
374
375			`#if SREAL_PART_BITS < 32`
376			`r->sig_hi = a->sig_hi + bb->sig_hi;`
377			`r->sig_lo = a->sig_lo + bb->sig_lo;`
378			`if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))`
379			`{`
380			`r->sig_hi++;`
381			`r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;`
382			`}`
383			`#else`
384			`r->sig = a->sig + bb->sig;`
385			`#endif`
386			`normalize (r);`
387			`return r;`
388			`}`
389
390			`/* R = A - B. Return R. /`
391
392			`sreal *`
393			`sreal_sub (sreal r, sreal a, sreal *b)`
394			`{`
395			`int dexp;`
396			`sreal tmp;`
397			`sreal *bb;`
398
399			`gcc_assert (sreal_compare (a, b) >= 0);`
400
401			`dexp = a->exp - b->exp;`
402			`r->exp = a->exp;`
403			`if (dexp > SREAL_BITS)`
404			`{`
405			`#if SREAL_PART_BITS < 32`
406			`r->sig_hi = a->sig_hi;`
407			`r->sig_lo = a->sig_lo;`
408			`#else`
409			`r->sig = a->sig;`
410			`#endif`
411			`return r;`
412			`}`
413			`if (dexp == 0)`
414			`bb = b;`
415			`else`
416			`{`
417			`copy (&tmp, b);`
418			`shift_right (&tmp, dexp);`
419			`bb = &tmp;`
420			`}`
421
422			`#if SREAL_PART_BITS < 32`
423			`if (a->sig_lo < bb->sig_lo)`
424			`{`
425			`r->sig_hi = a->sig_hi - bb->sig_hi - 1;`
426			`r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;`
427			`}`
428			`else`
429			`{`
430			`r->sig_hi = a->sig_hi - bb->sig_hi;`
431			`r->sig_lo = a->sig_lo - bb->sig_lo;`
432			`}`
433			`#else`
434			`r->sig = a->sig - bb->sig;`
435			`#endif`
436			`normalize (r);`
437			`return r;`
438			`}`
439
440			`/* R = A * B. Return R. /`
441
442			`sreal *`
443			`sreal_mul (sreal r, sreal a, sreal *b)`
444			`{`
445			`#if SREAL_PART_BITS < 32`
446			`if (a->sig_hi < SREAL_MIN_SIG \|\| b->sig_hi < SREAL_MIN_SIG)`
447			`{`
448			`r->sig_lo = 0;`
449			`r->sig_hi = 0;`
450			`r->exp = -SREAL_MAX_EXP;`
451			`}`
452			`else`
453			`{`
454			`unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;`
455			`if (sreal_compare (a, b) < 0)`
456			`{`
457			`sreal *swap;`
458			`swap = a;`
459			`a = b;`
460			`b = swap;`
461			`}`
462
463			`r->exp = a->exp + b->exp + SREAL_PART_BITS;`
464
465			`tmp1 = a->sig_lo * b->sig_lo;`
466			`tmp2 = a->sig_lo * b->sig_hi;`
467			`tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);`
468
469			`r->sig_hi = a->sig_hi * b->sig_hi;`
470			`r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);`
471			`tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;`
472			`tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;`
473			`tmp1 = tmp2 + tmp3;`
474
475			`r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);`
476			`r->sig_hi += tmp1 >> SREAL_PART_BITS;`
477
478			`normalize (r);`
479			`}`
480			`#else`
481			`if (a->sig < SREAL_MIN_SIG \|\| b->sig < SREAL_MIN_SIG)`
482			`{`
483			`r->sig = 0;`
484			`r->exp = -SREAL_MAX_EXP;`
485			`}`
486			`else`
487			`{`
488			`r->sig = a->sig * b->sig;`
489			`r->exp = a->exp + b->exp;`
490			`normalize (r);`
491			`}`
492			`#endif`
493			`return r;`
494			`}`
495
496			`/* R = A / B. Return R. /`
497
498			`sreal *`
499			`sreal_div (sreal r, sreal a, sreal *b)`
500			`{`
501			`#if SREAL_PART_BITS < 32`
502			`unsigned HOST_WIDE_INT tmp, tmp1, tmp2;`
503
504			`gcc_assert (b->sig_hi >= SREAL_MIN_SIG);`
505			`if (a->sig_hi < SREAL_MIN_SIG)`
506			`{`
507			`r->sig_hi = 0;`
508			`r->sig_lo = 0;`
509			`r->exp = -SREAL_MAX_EXP;`
510			`}`
511			`else`
512			`{`
513			`/* Since division by the whole number is pretty ugly to write`
514			`we are dividing by first 3/4 of bits of number. */`
515
516			`tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;`
517			`tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))`
518			`+ (b->sig_lo >> (SREAL_PART_BITS / 2)));`
519			`if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))`
520			`tmp2++;`
521
522			`r->sig_lo = 0;`
523			`tmp = tmp1 / tmp2;`
524			`tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);`
525			`r->sig_hi = tmp << SREAL_PART_BITS;`
526
527			`tmp = tmp1 / tmp2;`
528			`tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);`
529			`r->sig_hi += tmp << (SREAL_PART_BITS / 2);`
530
531			`tmp = tmp1 / tmp2;`
532			`r->sig_hi += tmp;`
533
534			`r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;`
535			`normalize (r);`
536			`}`
537			`#else`
538			`gcc_assert (b->sig != 0);`
539			`r->sig = (a->sig << SREAL_PART_BITS) / b->sig;`
540			`r->exp = a->exp - b->exp - SREAL_PART_BITS;`
541			`normalize (r);`
542			`#endif`
543			`return r;`
544			`}`