[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / e_exp.c

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, see <http://www.gnu.org/licenses/>.
 */
/***************************************************************************/
/*  MODULE_NAME:uexp.c                                                     */
/*                                                                         */
/*  FUNCTION:uexp                                                          */
/*           exp1                                                          */
/*                                                                         */
/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
/*              mpa.c mpexp.x slowexp.c                                    */
/*                                                                         */
/* An ultimate exp routine. Given an IEEE double machine number x          */
/* it computes the correctly rounded (to nearest) value of e^x             */
/* Assumption: Machine arithmetic operations are performed in              */
/* round to nearest mode of IEEE 754 standard.                             */
/*                                                                         */
/***************************************************************************/

#include <math.h>
#include "endian.h"
#include "uexp.h"
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
#include <math_private.h>
#include <fenv.h>
#include <float.h>

#ifndef SECTION
# define SECTION
#endif

double __slowexp (double);

/* An ultimate exp routine. Given an IEEE double machine number x it computes
   the correctly rounded (to nearest) value of e^x.  */
double
SECTION
__ieee754_exp (double x)
{
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp = {{0, 0}};
  int4 i, j, m, n, ex;
  double retval;

  {
    SET_RESTORE_ROUND (FE_TONEAREST);

    junk1.x = x;
    m = junk1.i[HIGH_HALF];
    n = m & hugeint;

    if (n > smallint && n < bigint)
      {
	y = x * log2e.x + three51.x;
	bexp = y - three51.x;	/*  multiply the result by 2**bexp        */

	junk1.x = y;

	eps = bexp * ln_two2.x;	/* x = bexp*ln(2) + t - eps               */
	t = x - bexp * ln_two1.x;

	y = t + three33.x;
	base = y - three33.x;	/* t rounded to a multiple of 2**-18      */
	junk2.x = y;
	del = (t - base) - eps;	/*  x = bexp*ln(2) + base + del           */
	eps = del + del * del * (p3.x * del + p2.x);

	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;

	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
	j = (junk2.i[LOW_HALF] & 511) << 1;

	al = coar.x[i] * fine.x[j];
	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	       + coar.x[i + 1] * fine.x[j + 1]);

	rem = (bet + bet * eps) + al * eps;
	res = al + rem;
	cor = (al - res) + rem;
	if (res == (res + cor * err_0))
	  {
	    retval = res * binexp.x;
	    goto ret;
	  }
	else
	  {
	    retval = __slowexp (x);
	    goto ret;
	  }			/*if error is over bound */
      }

    if (n <= smallint)
      {
	retval = 1.0;
	goto ret;
      }

    if (n >= badint)
      {
	if (n > infint)
	  {
	    retval = x + x;
	    goto ret;
	  }			/* x is NaN */
	if (n < infint)
	  {
	    if (x > 0)
	      goto ret_huge;
	    else
	      goto ret_tiny;
	  }
	/* x is finite,  cause either overflow or underflow  */
	if (junk1.i[LOW_HALF] != 0)
	  {
	    retval = x + x;
	    goto ret;
	  }			/*  x is NaN  */
	retval = (x > 0) ? inf.x : zero;	/* |x| = inf;  return either inf or 0 */
	goto ret;
      }

    y = x * log2e.x + three51.x;
    bexp = y - three51.x;
    junk1.x = y;
    eps = bexp * ln_two2.x;
    t = x - bexp * ln_two1.x;
    y = t + three33.x;
    base = y - three33.x;
    junk2.x = y;
    del = (t - base) - eps;
    eps = del + del * del * (p3.x * del + p2.x);
    i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
    j = (junk2.i[LOW_HALF] & 511) << 1;
    al = coar.x[i] * fine.x[j];
    bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	   + coar.x[i + 1] * fine.x[j + 1]);
    rem = (bet + bet * eps) + al * eps;
    res = al + rem;
    cor = (al - res) + rem;
    if (m >> 31)
      {
	ex = junk1.i[LOW_HALF];
	if (res < 1.0)
	  {
	    res += res;
	    cor += cor;
	    ex -= 1;
	  }
	if (ex >= -1022)
	  {
	    binexp.i[HIGH_HALF] = (1023 + ex) << 20;
	    if (res == (res + cor * err_0))
	      {
		retval = res * binexp.x;
		goto ret;
	      }
	    else
	      {
		retval = __slowexp (x);
		goto check_uflow_ret;
	      }			/*if error is over bound */
	  }
	ex = -(1022 + ex);
	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
	res *= binexp.x;
	cor *= binexp.x;
	eps = 1.0000000001 + err_0 * binexp.x;
	t = 1.0 + res;
	y = ((1.0 - t) + res) + cor;
	res = t + y;
	cor = (t - res) + y;
	if (res == (res + eps * cor))
	  {
	    binexp.i[HIGH_HALF] = 0x00100000;
	    retval = (res - 1.0) * binexp.x;
	    goto check_uflow_ret;
	  }
	else
	  {
	    retval = __slowexp (x);
	    goto check_uflow_ret;
	  }			/*   if error is over bound    */
      check_uflow_ret:
	if (retval < DBL_MIN)
	  {
#if FLT_EVAL_METHOD != 0
	    volatile
#endif
	      double force_underflow = tiny * tiny;
	    math_force_eval (force_underflow);
	  }
	if (retval == 0)
	  goto ret_tiny;
	goto ret;
      }
    else
      {
	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
	if (res == (res + cor * err_0))
	  retval = res * binexp.x * t256.x;
	else
	  retval = __slowexp (x);
	if (__isinf (retval))
	  goto ret_huge;
	else
	  goto ret;
      }
  }
ret:
  return retval;

 ret_huge:
  return hhuge * hhuge;

 ret_tiny:
  return tiny * tiny;
}
#ifndef __ieee754_exp
strong_alias (__ieee754_exp, __exp_finite)
#endif

/* Compute e^(x+xx).  The routine also receives bound of error of previous
   calculation.  If after computing exp the error exceeds the allowed bounds,
   the routine returns a non-positive number.  Otherwise it returns the
   computed result, which is always positive.  */
double
SECTION
__exp1 (double x, double xx, double error)
{
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp = {{0, 0}};
  int4 i, j, m, n, ex;

  junk1.x = x;
  m = junk1.i[HIGH_HALF];
  n = m & hugeint;		/* no sign */

  if (n > smallint && n < bigint)
    {
      y = x * log2e.x + three51.x;
      bexp = y - three51.x;	/*  multiply the result by 2**bexp        */

      junk1.x = y;

      eps = bexp * ln_two2.x;	/* x = bexp*ln(2) + t - eps               */
      t = x - bexp * ln_two1.x;

      y = t + three33.x;
      base = y - three33.x;	/* t rounded to a multiple of 2**-18      */
      junk2.x = y;
      del = (t - base) + (xx - eps);	/*  x = bexp*ln(2) + base + del      */
      eps = del + del * del * (p3.x * del + p2.x);

      binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;

      i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
      j = (junk2.i[LOW_HALF] & 511) << 1;

      al = coar.x[i] * fine.x[j];
      bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	     + coar.x[i + 1] * fine.x[j + 1]);

      rem = (bet + bet * eps) + al * eps;
      res = al + rem;
      cor = (al - res) + rem;
      if (res == (res + cor * (1.0 + error + err_1)))
	return res * binexp.x;
      else
	return -10.0;
    }

  if (n <= smallint)
    return 1.0;			/*  if x->0 e^x=1 */

  if (n >= badint)
    {
      if (n > infint)
	return (zero / zero);	/* x is NaN,  return invalid */
      if (n < infint)
	return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
      /* x is finite,  cause either overflow or underflow  */
      if (junk1.i[LOW_HALF] != 0)
	return (zero / zero);	/*  x is NaN  */
      return ((x > 0) ? inf.x : zero);	/* |x| = inf;  return either inf or 0 */
    }

  y = x * log2e.x + three51.x;
  bexp = y - three51.x;
  junk1.x = y;
  eps = bexp * ln_two2.x;
  t = x - bexp * ln_two1.x;
  y = t + three33.x;
  base = y - three33.x;
  junk2.x = y;
  del = (t - base) + (xx - eps);
  eps = del + del * del * (p3.x * del + p2.x);
  i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
  j = (junk2.i[LOW_HALF] & 511) << 1;
  al = coar.x[i] * fine.x[j];
  bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	 + coar.x[i + 1] * fine.x[j + 1]);
  rem = (bet + bet * eps) + al * eps;
  res = al + rem;
  cor = (al - res) + rem;
  if (m >> 31)
    {
      ex = junk1.i[LOW_HALF];
      if (res < 1.0)
	{
	  res += res;
	  cor += cor;
	  ex -= 1;
	}
      if (ex >= -1022)
	{
	  binexp.i[HIGH_HALF] = (1023 + ex) << 20;
	  if (res == (res + cor * (1.0 + error + err_1)))
	    return res * binexp.x;
	  else
	    return -10.0;
	}
      ex = -(1022 + ex);
      binexp.i[HIGH_HALF] = (1023 - ex) << 20;
      res *= binexp.x;
      cor *= binexp.x;
      eps = 1.00000000001 + (error + err_1) * binexp.x;
      t = 1.0 + res;
      y = ((1.0 - t) + res) + cor;
      res = t + y;
      cor = (t - res) + y;
      if (res == (res + eps * cor))
	{
	  binexp.i[HIGH_HALF] = 0x00100000;
	  return (res - 1.0) * binexp.x;
	}
      else
	return -10.0;
    }
  else
    {
      binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
      if (res == (res + cor * (1.0 + error + err_1)))
	return res * binexp.x * t256.x;
      else
	return -10.0;
    }
}
Commit	Line	Data
e4d82761 UD	1	/*
e4d82761 UD	2	* IBM Accurate Mathematical Library
aeb25823	3	* written by International Business Machines Corp.
b168057a	4	* Copyright (C) 2001-2015 Free Software Foundation, Inc.
e4d82761 UD	5	*
	6	* This program is free software; you can redistribute it and/or modify
	7	* it under the terms of the GNU Lesser General Public License as published by
cc7375ce	8	* the Free Software Foundation; either version 2.1 of the License, or
e4d82761 UD	9	* (at your option) any later version.
	10	*
	11	* This program is distributed in the hope that it will be useful,
	12	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c6c6dd48	14	* GNU Lesser General Public License for more details.
e4d82761 UD	15	*
e4d82761 UD	16	* You should have received a copy of the GNU Lesser General Public License
59ba27a6	17	* along with this program; if not, see <http://www.gnu.org/licenses/>.
e4d82761 UD	18	*/
	19	/***************************************************************************/
	20	/* MODULE_NAME:uexp.c */
	21	/* */
	22	/* FUNCTION:uexp */
	23	/* exp1 */
	24	/* */
	25	/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */
	26	/* mpa.c mpexp.x slowexp.c */
	27	/* */
	28	/* An ultimate exp routine. Given an IEEE double machine number x */
	29	/* it computes the correctly rounded (to nearest) value of e^x */
	30	/* Assumption: Machine arithmetic operations are performed in */
	31	/* round to nearest mode of IEEE 754 standard. */
	32	/* */
	33	/***************************************************************************/
	34
f3426898	35	#include <math.h>
e4d82761	36	#include "endian.h"
e859d1d9	37	#include "uexp.h"
e4d82761 UD	38	#include "mydefs.h"
	39	#include "MathLib.h"
	40	#include "uexp.tbl"
1ed0291c	41	#include <math_private.h>
28afd92d	42	#include <fenv.h>
749008ff	43	#include <float.h>
e859d1d9	44
31d3cc00 UD	45	#ifndef SECTION
	46	# define SECTION
	47	#endif
	48
e7b2d1dd	49	double __slowexp (double);
e4d82761	50
e7b2d1dd SP	51	/* An ultimate exp routine. Given an IEEE double machine number x it computes
e7b2d1dd SP	52	the correctly rounded (to nearest) value of e^x. */
31d3cc00 UD	53	double
31d3cc00 UD	54	SECTION
e7b2d1dd SP	55	__ieee754_exp (double x)
e7b2d1dd SP	56	{
e4d82761	57	double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
e7b2d1dd SP	58	mynumber junk1, junk2, binexp = {{0, 0}};
e7b2d1dd SP	59	int4 i, j, m, n, ex;
28afd92d JM	60	double retval;
28afd92d JM	61
b376a11a JM	62	{
	63	SET_RESTORE_ROUND (FE_TONEAREST);
	64
	65	junk1.x = x;
	66	m = junk1.i[HIGH_HALF];
	67	n = m & hugeint;
	68
	69	if (n > smallint && n < bigint)
	70	{
	71	y = x * log2e.x + three51.x;
	72	bexp = y - three51.x; /* multiply the result by 2*bexp /
	73
	74	junk1.x = y;
	75
	76	eps = bexp * ln_two2.x; /* x = bexpln(2) + t - eps /
	77	t = x - bexp * ln_two1.x;
	78
	79	y = t + three33.x;
	80	base = y - three33.x; /* t rounded to a multiple of 2*-18 /
	81	junk2.x = y;
	82	del = (t - base) - eps; /* x = bexpln(2) + base + del /
	83	eps = del + del * del * (p3.x * del + p2.x);
	84
	85	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
	86
	87	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
	88	j = (junk2.i[LOW_HALF] & 511) << 1;
	89
	90	al = coar.x[i] * fine.x[j];
	91	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	92	+ coar.x[i + 1] * fine.x[j + 1]);
	93
	94	rem = (bet + bet * eps) + al * eps;
	95	res = al + rem;
	96	cor = (al - res) + rem;
	97	if (res == (res + cor * err_0))
	98	{
	99	retval = res * binexp.x;
	100	goto ret;
	101	}
	102	else
	103	{
	104	retval = __slowexp (x);
	105	goto ret;
	106	} /if error is over bound /
	107	}
	108
	109	if (n <= smallint)
	110	{
	111	retval = 1.0;
	112	goto ret;
	113	}
	114
	115	if (n >= badint)
	116	{
	117	if (n > infint)
	118	{
	119	retval = x + x;
	120	goto ret;
	121	} /* x is NaN */
	122	if (n < infint)
	123	{
	124	if (x > 0)
	125	goto ret_huge;
126	else
127	goto ret_tiny;
128	}
129	/* x is finite, cause either overflow or underflow */
130	if (junk1.i[LOW_HALF] != 0)
131	{
132	retval = x + x;
133	goto ret;
134	} /* x is NaN */
135	retval = (x > 0) ? inf.x : zero; /* \|x\| = inf; return either inf or 0 */
136	goto ret;
137	}
138
139	y = x * log2e.x + three51.x;
140	bexp = y - three51.x;
141	junk1.x = y;
142	eps = bexp * ln_two2.x;
143	t = x - bexp * ln_two1.x;
144	y = t + three33.x;
145	base = y - three33.x;
146	junk2.x = y;
147	del = (t - base) - eps;
148	eps = del + del * del * (p3.x * del + p2.x);
149	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
150	j = (junk2.i[LOW_HALF] & 511) << 1;
151	al = coar.x[i] * fine.x[j];
152	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
153	+ coar.x[i + 1] * fine.x[j + 1]);
154	rem = (bet + bet * eps) + al * eps;
155	res = al + rem;
156	cor = (al - res) + rem;
157	if (m >> 31)
158	{
159	ex = junk1.i[LOW_HALF];
160	if (res < 1.0)
161	{
162	res += res;
163	cor += cor;
164	ex -= 1;
165	}
166	if (ex >= -1022)
167	{
168	binexp.i[HIGH_HALF] = (1023 + ex) << 20;
169	if (res == (res + cor * err_0))
170	{
171	retval = res * binexp.x;
172	goto ret;
173	}
174	else
175	{
176	retval = __slowexp (x);
177	goto check_uflow_ret;
178	} /if error is over bound /
179	}
180	ex = -(1022 + ex);
181	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
182	res *= binexp.x;
183	cor *= binexp.x;
184	eps = 1.0000000001 + err_0 * binexp.x;
185	t = 1.0 + res;
186	y = ((1.0 - t) + res) + cor;
187	res = t + y;
188	cor = (t - res) + y;
189	if (res == (res + eps * cor))
190	{
191	binexp.i[HIGH_HALF] = 0x00100000;
192	retval = (res - 1.0) * binexp.x;
193	goto check_uflow_ret;
194	}
195	else
196	{
197	retval = __slowexp (x);
198	goto check_uflow_ret;
199	} /* if error is over bound */
200	check_uflow_ret:
201	if (retval < DBL_MIN)
202	{
749008ff	203	#if FLT_EVAL_METHOD != 0
b376a11a	204	volatile
749008ff	205	#endif
b376a11a JM	206	double force_underflow = tiny * tiny;
	207	math_force_eval (force_underflow);
	208	}
	209	if (retval == 0)
	210	goto ret_tiny;
	211	goto ret;
	212	}
	213	else
	214	{
	215	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
	216	if (res == (res + cor * err_0))
e7b2d1dd	217	retval = res * binexp.x * t256.x;
b376a11a	218	else
e7b2d1dd	219	retval = __slowexp (x);
b376a11a JM	220	if (__isinf (retval))
	221	goto ret_huge;
	222	else
e7b2d1dd	223	goto ret;
b376a11a JM	224	}
b376a11a JM	225	}
e7b2d1dd	226	ret:
28afd92d	227	return retval;
b376a11a JM	228
	229	ret_huge:
	230	return hhuge * hhuge;
	231
	232	ret_tiny:
	233	return tiny * tiny;
e4d82761	234	}
af968f62	235	#ifndef __ieee754_exp
bcf01e6d	236	strong_alias (__ieee754_exp, __exp_finite)
af968f62	237	#endif
e4d82761	238
e7b2d1dd SP	239	/* Compute e^(x+xx). The routine also receives bound of error of previous
	240	calculation. If after computing exp the error exceeds the allowed bounds,
	241	the routine returns a non-positive number. Otherwise it returns the
	242	computed result, which is always positive. */
31d3cc00 UD	243	double
31d3cc00 UD	244	SECTION
e7b2d1dd SP	245	__exp1 (double x, double xx, double error)
e7b2d1dd SP	246	{
e4d82761	247	double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
e7b2d1dd SP	248	mynumber junk1, junk2, binexp = {{0, 0}};
e7b2d1dd SP	249	int4 i, j, m, n, ex;
e4d82761 UD	250
	251	junk1.x = x;
	252	m = junk1.i[HIGH_HALF];
e7b2d1dd	253	n = m & hugeint; /* no sign */
e4d82761	254
e7b2d1dd SP	255	if (n > smallint && n < bigint)
	256	{
	257	y = x * log2e.x + three51.x;
	258	bexp = y - three51.x; /* multiply the result by 2*bexp /
e4d82761	259
e7b2d1dd	260	junk1.x = y;
e4d82761	261
e7b2d1dd SP	262	eps = bexp * ln_two2.x; /* x = bexpln(2) + t - eps /
e7b2d1dd SP	263	t = x - bexp * ln_two1.x;
e4d82761	264
e7b2d1dd SP	265	y = t + three33.x;
	266	base = y - three33.x; /* t rounded to a multiple of 2*-18 /
	267	junk2.x = y;
	268	del = (t - base) + (xx - eps); /* x = bexpln(2) + base + del /
	269	eps = del + del * del * (p3.x * del + p2.x);
e4d82761	270
e7b2d1dd	271	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
e4d82761	272
e7b2d1dd SP	273	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
e7b2d1dd SP	274	j = (junk2.i[LOW_HALF] & 511) << 1;
e4d82761	275
e7b2d1dd SP	276	al = coar.x[i] * fine.x[j];
	277	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	278	+ coar.x[i + 1] * fine.x[j + 1]);
e4d82761	279
e7b2d1dd SP	280	rem = (bet + bet * eps) + al * eps;
	281	res = al + rem;
	282	cor = (al - res) + rem;
	283	if (res == (res + cor * (1.0 + error + err_1)))
	284	return res * binexp.x;
	285	else
	286	return -10.0;
	287	}
e4d82761	288
e7b2d1dd SP	289	if (n <= smallint)
	290	return 1.0; /* if x->0 e^x=1 */
	291
	292	if (n >= badint)
	293	{
	294	if (n > infint)
	295	return (zero / zero); /* x is NaN, return invalid */
	296	if (n < infint)
	297	return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
	298	/* x is finite, cause either overflow or underflow */
	299	if (junk1.i[LOW_HALF] != 0)
	300	return (zero / zero); /* x is NaN */
	301	return ((x > 0) ? inf.x : zero); /* \|x\| = inf; return either inf or 0 */
	302	}
e4d82761	303
e7b2d1dd	304	y = x * log2e.x + three51.x;
e4d82761 UD	305	bexp = y - three51.x;
e4d82761 UD	306	junk1.x = y;
e7b2d1dd SP	307	eps = bexp * ln_two2.x;
e7b2d1dd SP	308	t = x - bexp * ln_two1.x;
e4d82761 UD	309	y = t + three33.x;
	310	base = y - three33.x;
	311	junk2.x = y;
e7b2d1dd SP	312	del = (t - base) + (xx - eps);
	313	eps = del + del * del * (p3.x * del + p2.x);
	314	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
	315	j = (junk2.i[LOW_HALF] & 511) << 1;
	316	al = coar.x[i] * fine.x[j];
	317	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
	318	+ coar.x[i + 1] * fine.x[j + 1]);
	319	rem = (bet + bet * eps) + al * eps;
e4d82761 UD	320	res = al + rem;
e4d82761 UD	321	cor = (al - res) + rem;
e7b2d1dd SP	322	if (m >> 31)
	323	{
	324	ex = junk1.i[LOW_HALF];
	325	if (res < 1.0)
	326	{
	327	res += res;
	328	cor += cor;
	329	ex -= 1;
	330	}
	331	if (ex >= -1022)
	332	{
	333	binexp.i[HIGH_HALF] = (1023 + ex) << 20;
	334	if (res == (res + cor * (1.0 + error + err_1)))
	335	return res * binexp.x;
	336	else
	337	return -10.0;
	338	}
	339	ex = -(1022 + ex);
	340	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
	341	res *= binexp.x;
	342	cor *= binexp.x;
	343	eps = 1.00000000001 + (error + err_1) * binexp.x;
	344	t = 1.0 + res;
	345	y = ((1.0 - t) + res) + cor;
	346	res = t + y;
	347	cor = (t - res) + y;
	348	if (res == (res + eps * cor))
	349	{
	350	binexp.i[HIGH_HALF] = 0x00100000;
	351	return (res - 1.0) * binexp.x;
	352	}
	353	else
	354	return -10.0;
	355	}
	356	else
	357	{
	358	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
	359	if (res == (res + cor * (1.0 + error + err_1)))
	360	return res * binexp.x * t256.x;
	361	else
	362	return -10.0;
e4d82761	363	}
f7eac6eb	364	}