Subversion Repositories tendra.SVN

/*

 * Copyright (c) 2002-2006 The TenDRA Project <http://www.tendra.org/>.

 * All rights reserved.

 *

 * Redistribution and use in source and binary forms, with or without

 * modification, are permitted provided that the following conditions are met:

 *

 * 1. Redistributions of source code must retain the above copyright notice,

 *    this list of conditions and the following disclaimer.

 * 2. Redistributions in binary form must reproduce the above copyright notice,

 *    this list of conditions and the following disclaimer in the documentation

 *    and/or other materials provided with the distribution.

 * 3. Neither the name of The TenDRA Project nor the names of its contributors

 *    may be used to endorse or promote products derived from this software

 *    without specific, prior written permission.

 *

 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS

 * IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,

 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR

 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR

 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

 * EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;

 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,

 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR

 * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF

 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 *

 * $Id$

 */

/*

    		 Crown Copyright (c) 1997

    This TenDRA(r) Computer Program is subject to Copyright

    owned by the United Kingdom Secretary of State for Defence

    acting through the Defence Evaluation and Research Agency

    (DERA).  It is made available to Recipients with a

    royalty-free licence for its use, reproduction, transfer

    to other parties and amendment for any purpose not excluding

    product development provided that any such use et cetera

    shall be deemed to be acceptance of the following conditions:-

        (1) Its Recipients shall ensure that this Notice is

        reproduced upon any copies or amended versions of it;

        (2) Any amended version of it shall be clearly marked to

        show both the nature of and the organisation responsible

        for the relevant amendment or amendments;

        (3) Its onward transfer from a recipient to another

        party shall be deemed to be that party's acceptance of

        these conditions;

        (4) DERA gives no warranty or assurance as to its

        quality or suitability for any purpose and DERA accepts

        no liability whatsoever in relation to any use to which

        it may be put.

*/

/**********************************************************************

$Author: release $

$Date: 1998/01/17 15:55:47 $

$Revision: 1.1.1.1 $

$Log: flpt.c,v $

 * Revision 1.1.1.1  1998/01/17  15:55:47  release

 * First version to be checked into rolling release.

 *

 * Revision 1.4  1996/01/10  14:58:48  currie

 * BIGEND var params chars & shorts

 *

 * Revision 1.3  1995/09/11  13:58:37  currie

 * gcc pedantry

 *

 * Revision 1.2  1995/08/09  08:59:57  currie

 * round bug

 *

 * Revision 1.1  1995/04/06  10:44:05  currie

 * Initial revision

 *

***********************************************************************/

/***********************************************************************

                                flpt.c

  Routines for handling the internal floating point representation.

 ***********************************************************************/

#include "config.h"

#include <limits.h>

#include "common_types.h"

#include "flpttypes.h"

#include "exp.h"

#include "xalloc.h"

#include "szs_als.h"

#include "messages_c.h"

#include "basicread.h"

#include "installglob.h"

#include "flpt.h"

/* MACROS */

#define MAX_USEFUL_DECEXP 10000

/* plenty big enough for IEEE and VAX */

#define initial_flpts 50

#define extra_flpts 50

/***********************************************************************

          ---------- Floating point emulator code ----------

 NB - lines which assume the use of an ASCII character set are accompanied

      by the comment ASCII

 ***********************************************************************/

typedef struct _dbl {

  /* Double precision type - used for internal working */

  Fdig mant[2 * MANT_SIZE];

  int  sign;

  int  exp;

} dbl;

/* The current rounding rule. This may be one of the following :

     R2ZERO - round towards zero.

     R2PINF - round towards positive infinity.

     R2NINF - round towards negative infinity.

     R2NEAR - round to nearest integer.

 It should be set using the "set_round" function.

*/

/* IDENTITIES */

static Fdig bitround[16] = {0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80,

			     0x100, 0x200, 0x400, 0x800, 0x1000, 0x2000,

			     0x4000, 0x8000};

static Fdig bitmask[17] = {0xffff, 0xfffe, 0xfffc, 0xfff8,

			    0xfff0, 0xffe0, 0xffc0, 0xff80,

			    0xff00, 0xfe00, 0xfc00, 0xf800,

			    0xf000, 0xe000, 0xc000, 0x8000, 0x0};

static flt powers[16] = {

	{{10, 0, 0, 0, 0, 0, 0, 0},1, 0},

	{{100, 0, 0, 0, 0, 0, 0, 0},1, 0},

	{{1000, 0, 0, 0, 0, 0, 0, 0},1, 0},

	{{10000, 0, 0, 0, 0, 0, 0, 0},1, 0},

	{{1, 34464, 0, 0, 0, 0, 0, 0},1, 1},

	{{15, 16960, 0, 0, 0, 0, 0, 0},1, 1},

	{{152, 38528, 0, 0, 0, 0, 0, 0},1, 1},

	{{1525, 57600, 0, 0, 0, 0, 0, 0},1, 1},

	{{15258, 51712, 0, 0, 0, 0, 0, 0},1, 1},

	{{2, 21515, 58368, 0, 0, 0, 0, 0},1, 2},

	{{23, 18550, 59392, 0, 0, 0, 0, 0},1, 2},

	{{232, 54437, 4096, 0, 0, 0, 0, 0},1, 2},

	{{2328, 20082, 40960, 0, 0, 0, 0, 0},1, 2},

	{{23283, 4218, 16384, 0, 0, 0, 0, 0},1, 2},

	{{3, 36222, 42182, 32768, 0, 0, 0, 0},1, 3},

	{{35, 34546, 28609, 0, 0, 0, 0, 0},1, 3}

};

static int two_powers[] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024,

			    2048, 4096, 8192, 16384, 32768};

static int round_type = R2NEAR;

/* VARIABLES */

/* All variables initialised */

int tot_flpts;		/* total number of floating point numbers */

int flptfree;		/* floating point free list */

int flpt_left;		/* number of floating point numbers left */

flt *flptnos;		/* the extendable array of floating point numbers */

int fzero_no;

int fone_no;

/* PROCEDURES */

flpt new_flpt(void);

/* ----- internal functions ----- */

/* "f1" and "f2" are two legal single precision numbers. The result is -1, 0 or

 * +1 depending on whether the magnitude of "f1" is less than, equal to or

 * greater than that of "f2" */

static int

mag_cmp(flt f1, flt f2)

{

  int i;

  if (f1.sign == 0) {

    if (f2.sign == 0) {

      return(0);

    }

    return(-1);

  }

  if (f2.sign == 0) {

    return(1);

  }

  if (f1.exp != f2.exp) {

    return((f1.exp < f2.exp) ? -1 : 1);

  }

  for (i = 0; i < MANT_SIZE; i++) {

    if ((f1.mant[i]) != (f2.mant[i])) {

      return(((f1.mant[i]) < (f2.mant[i])) ? -1 : 1);

    }

  }

  return(0);

}

/* "d1" and "d2" are two legal double precision numbers. The result is -1, 0 or

 * +1 depending on whether the magnitude of "d1" is less than, equal to or

 * greater than that of "d2" */

static int

dbl_mag_cmp (dbl d1, dbl d2)

{

  int i;

  if (d1.sign == 0) {

    if (d2.sign == 0) {

      return(0);

    }

    return(-1);

  }

  if (d2.sign == 0) {

    return(1);

  }

  if (d1.exp != d2.exp) {

    return((d1.exp < d2.exp) ? -1 : 1);

  }

  for (i = 0; i < (2 * MANT_SIZE); i++) {

    if ((d1.mant[i]) != (d2.mant[i])) {

      return(((d1.mant[i]) < (d2.mant[i]))? -1 : 1);

    }

  }

  return(0);

}

/* "f" is a pointer to a legal single precision number. "c" is a character

 * value from 0 to FBASE-1. On return, the single precision number pointed to

 * by "f" will have been modified as if "c" were the digit following the least

 * significant digit of the number pointed to by "f" and the number has been

 * rounded according to the current rounding rule */

static void

flt_int_round(flt *f, unsigned int ic)

{

  int i;

  switch (round_type) {

  default:

  case R2ZERO:

    return;

  case R2PINF:

    if (((f->sign) != 1) || (ic == 0))

      return;

    break;

  case R2NINF:

    if (((f->sign) != -1) || (ic == 0))

      return;

    break;

  case R2NEAR:

    if (ic < (unsigned int)(FBASE / 2)) {

      return;

    }

    break;

  }

  ic = 1;

  for (i = (MANT_SIZE - 1); i >= 0; i--) {

    (f->mant)[i] = (unsigned short)((ic += (unsigned int)((f->mant)[i])) % FBASE);

    /* CAST:jmf: unsigned short since FBASE <= MAX_USHORT+1 */

    if ((ic /= FBASE) == 0) {

      return;

    }

  }

  (f->exp)++;

  (f->mant)[0] = 1;

  for (i = 1; i < MANT_SIZE; i++) {

    (f->mant)[i] = 0;

  }

}

/* "d" is a legal double precision number.  "f" is a pointer to a single

 * precision number. On return, the number pointed to by "f" contains the value

 * of "d" but to less precision. The current rounding rule is used to truncate

 * the mantissa */

static void

dbl2float(dbl d, flt *f)

{

  int i;

  (f->sign) = d.sign;

  (f->exp) = d.exp;

  for (i = 0; i < MANT_SIZE; i++) {

    (f->mant)[i] = d.mant[i];

  }

  flt_int_round(f, (unsigned int)d.mant[MANT_SIZE]);

}

/* "f" is a legal single precision number.  "d" is a pointer to a double

 * precision number. On return, the number pointed to by "d" will contain the

 * value of "f" */

static void

flt2double(flt f, dbl *d)

{

  int i;

  (d->sign) = f.sign;

  (d->exp) = f.exp;

  for (i = 0; i < MANT_SIZE; i++) {

    (d->mant)[i] = f.mant[i];

  }

  for (i = MANT_SIZE; i < (2 * MANT_SIZE); i++) {

    (d->mant)[i] = 0;

  }

}

/* "d" is a pointer to a legal double precision number. "dist" is a int with

 * value greater than zero. On return, the number pointed to by "d" keeps the

 * same value ( although some loss of precision may occur ) but is no longer in

 * normalised form. The mantissa is shifted right "dist" places, and the

 * exponent is adjusted to keep the value the same */

static void

shift_right(dbl *d, int dist)

{

  int i;

  int j = (2 * MANT_SIZE) - 1 - dist;

  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--, j--) {

    (d->mant)[i] = (unsigned short)((j >= 0) ? ((d->mant)[j]) : 0);

  }

  (d->exp) += dist;

}

/* "d1" and "d2" are double precision numbers, with the same exponent and

 * opposite (non-zero) signs. "res" is a pointer to a single precision number,

 * whose exponent has been set to the same value as "d1" and "d2". On return,

 * the number pointed to by "res" is the result of adding "d1" to "d2" and

 * rounding it by the current rounding rule. The return value is either OKAY or

 * EXP2BIG ( in which case the value of the number pointed to by "res" is

 * undefined ) */

static int

sub_mantissas(dbl d1, dbl d2, flt *res)

{

  dbl *dp1 = &d1, *dp2 = &d2, rdbl;

  unsigned int  c = 0;

  int i, cmp = dbl_mag_cmp(d1, d2);

  if (cmp == 0) {

    (res->sign) = 0;

    (res->exp) = 0;

    for (i = 0; i < MANT_SIZE; i++) {

      (res->mant)[i] = 0;

    }

    return(OKAY);

  } else if (cmp == -1) {

    dp1 = &d2;

    dp2 = &d1;

  }

  rdbl.exp = (res->exp);

  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {

    int s = (int)(((dp1->mant)[i]) - ((dp2->mant)[i]) - c);

    c = 0;

    if (s < 0) {

      c = 1;

      s += FBASE;

    }

    rdbl.mant[i] = (unsigned short)s;	/* CAST:jmf: */

  }

  rdbl.sign = (dp1->sign);

  while ((rdbl.mant[0]) == 0) {

    for (i = 1; i < (2 * MANT_SIZE); i++) {

      rdbl.mant[i - 1] = rdbl.mant[i];

    }

    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;

    if ((--rdbl.exp) < E_MIN) {

      return(EXP2BIG);

    }

  }

  dbl2float(rdbl, res);

  return(OKAY);

}

/* "d1" and "d2" are double precision numbers, with the same exponent and

 * non-zero signs. "res" is a pointer to a single precision number which has

 * the same exponent as "d1" and "d2". On return, the number pointed to by

 * "res" will be the result of adding "d1" to "d2" and rounding it according to

 * the current rounding rule. The return value will be either OKAY or EXP2BIG (

 * in which case the value of the number pointed to by "res" is undefined ) */

static int

add_mantissas(dbl d1, dbl d2, flt * res)

{

  int i;

  unsigned int c;

  if (d1.sign == d2.sign) {

    dbl rdbl;

    (res->sign) = d1.sign;

    flt2double(*res, &rdbl);

    c = 0;

    for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {

      c += (unsigned int)(d1.mant[i] + d2.mant[i]);

      rdbl.mant[i] = (unsigned short)(c % FBASE);	/* CAST:jmf: */

      c /= FBASE;

    }

    if (c) {

      for (i = (2 * MANT_SIZE) - 1; i > 0; i--) {

	rdbl.mant[i] = rdbl.mant[i - 1];

      }

      rdbl.mant[0] = (unsigned short)c;	/* CAST:jmf: */

      if (((++rdbl.exp) >= E_MAX) || (rdbl.exp <= E_MIN)) {

	return(EXP2BIG);

      }

    }

    dbl2float(rdbl, res);

    return(OKAY);

  }

  return(sub_mantissas(d1, d2, res));

}

/* ----- "public" emulator functions ----- */

/* "f" is a pointer to a single precision number. On return, the number pointed

 * to by "f" will have the value zero */

void

flt_zero(flt *f)

{

  int i;

  (f->sign) = 0;

  (f->exp) = 0;

  for (i = 0; i < MANT_SIZE; i++) {

    (f->mant)[i] = 0;

  }

}

/* "f" is a legal single precision number.  "res" is a pointer to a single

 * precision number. On return, the number pointed to by "res" is the same as

 * "f" */

void

flt_copy(flt f, flt *res)

{

  int i;

  (res->sign) = f.sign;

  (res->exp) = f.exp;

  for (i = 0; i < MANT_SIZE; i++) {

    (res->mant)[i] = f.mant[i];

  }

  return;

}

/* "f1" and "f2" are legal single precision numbers. "res" is a pointer to a

 * single precision number. On return, if the result is OKAY, then the number

 * pointed to by "res" will contain the value of adding "f1" to "f2". If the

 * result is EXP2BIG, then the value of the number pointed to by "res" is

 * undefined */

int

flt_add(flt f1, flt f2, flt *res)

{

  dbl d1, d2;

  if (f1.sign == 0) {

    flt_copy(f2, res);

    return(OKAY);

  } else if (f2.sign == 0) {

    flt_copy(f1, res);

    return(OKAY);

  }

  flt2double(f1, &d1);

  flt2double(f2, &d2);

  if (d1.exp < d2.exp) {

    shift_right(&d1, d2.exp - d1.exp);

  } else if (d1.exp > d2.exp) {

    shift_right(&d2, d1.exp - d2.exp);

  }

  (res->exp) = d1.exp;

  return(add_mantissas(d1, d2, res));

}

/* "f1" and "f2" are legal single precision numbers. "res" is a pointer to a

 * single precision number. On return, if the result is OKAY, then the number

 * pointed to by "res" will contain the value of subtracting "f2" from "f1". If

 * the result is EXP2BIG, then the value of the number pointed to by "res" is

 * undefined */

int

flt_sub(flt f1, flt f2, flt *res)

{

  f2.sign = -f2.sign;

  return(flt_add(f1, f2, res));

}

#if FBASE == 10

/* "f1" and "f2" are legal single precision numbers. "res" is a pointer to a

 * single precision number. On return, if the result is OKAY, then the number

 * pointed to by "res" will contain the value of multiplying "f1" by "f2".  If

 * the result is EXP2BIG, then the value of the number pointed to by "res" is

 * undefined */

int

flt_mul(flt f1, flt f2, flt *res)

{

  dbl rdbl;

  int i, j;

  unsigned int c = 0;

  unsigned int acc[2 * MANT_SIZE];

  if ((f1.sign == 0) || (f2.sign == 0)) {

    flt_zero(res);

    return(OKAY);

  }

  rdbl.sign = ((f1.sign == f2.sign)? 1 : -1);

  rdbl.exp = (f1.exp + f2.exp + 1);

  if ((rdbl.exp >= E_MAX) || (rdbl.exp <= E_MIN)) {

    return(EXP2BIG);

  }

  for (i = 0; i < (2 * MANT_SIZE); i++) {

    acc[i] = 0;

  }

  for (i = MANT_SIZE - 1; i >= 0; i--) {

    for (j = MANT_SIZE - 1; j >= 0; j--) {

      acc[i + j + 1] += (f1.mant[i] * f2.mant[j]);

    }

  }

  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {

    c += acc[i];

    rdbl.mant[i] = (c % FBASE);

    c /= FBASE;

  }

  while (rdbl.mant[0] == 0) {

    for (i = 1; i < (2 * MANT_SIZE); i++) {

      rdbl.mant[i - 1] = rdbl.mant[i];

    }

    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;

    if ((--rdbl.exp) <= E_MIN) {

      return(EXP2BIG);

    }

  }

  dbl2float(rdbl, res);

  return(OKAY);

}

#else

/* "f1" and "f2" are legal single precision numbers. "res" is a pointer to a

 * single precision number. On return, if the result is OKAY, then the number

 * pointed to by "res" will contain the value of multiplying "f1" by "f2".  If

 * the result is EXP2BIG, then the value of the number pointed to by "res" is

 * undefined */

int

flt_mul(flt f1, flt f2, flt *res)

{

  dbl rdbl;

  int i, j, k;

  unsigned int acc[2 * MANT_SIZE];

  if ((f1.sign == 0) || (f2.sign == 0)) {

    flt_zero(res);

    return(OKAY);

  }

  rdbl.sign = ((f1.sign == f2.sign)? 1 : -1);

  rdbl.exp = (f1.exp + f2.exp + 1);

  if ((rdbl.exp >= E_MAX) || (rdbl.exp <= E_MIN)) {

    return(EXP2BIG);

  }

  for (i = 0; i < (2 * MANT_SIZE); i++) {

    acc[i] = 0;

  }

  SET(acc);

  for (i = MANT_SIZE - 1; i >= 0; i--) {

    for (j = MANT_SIZE - 1; j >= 0; j--) {

      unsigned int temp = (unsigned int)(f1.mant[i] * f2.mant[j]);

      unsigned int atl;

      for (k = i + j + 1; temp != 0; --k) {

	atl = acc[k] + (temp & 0xffff);

	acc[k] = atl & 0xffff;

	temp = (atl >> 16) + (temp >> 16);

      }

    }

  }

  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {

    rdbl.mant[i] = (unsigned short)(acc[i]);	/* CAST:jmf: */

  }

  while (rdbl.mant[0] == 0) {

    for (i = 1; i < (2 * MANT_SIZE); i++) {

      rdbl.mant[i - 1] = rdbl.mant[i];

    }

    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;

    if ((--rdbl.exp) <= E_MIN) {

      return(EXP2BIG);

    }

  }

  dbl2float(rdbl, res);

  return(OKAY);

}

#endif

#if FBASE == 10

/* "f1" and "f2" are legal single precision numbers, and "f2" is non-zero.

 * "res" is a pointer to a single precision number. On return, if the result is

 * OKAY, then the number pointed to by "res" contains the result of dividing

 * "f1" by "f2". If the result is EXP2BIG or DIVBY0 then the contents of the

 * number pointed to by "res" are undefined */

int

flt_div(flt f1, flt f2, flt *res)

{

  dbl rdbl;

  int i = 0;

  flt f3;

  if (f2.sign == 0) {

    return(DIVBY0);

  }

  if (f1.sign == 0) {

    flt_zero(res);

    return(OKAY);

  }

  rdbl.sign = ((f1.sign == f2.sign) ? 1 : -1);

  rdbl.exp = (f1.exp - f2.exp);

  if ((rdbl.exp >= E_MAX) || (rdbl.exp <= E_MIN)) {

    return(EXP2BIG);

  }

  f1.sign = 1;

  f2.sign = 1;

  f1.exp = 0;

  f2.exp = 0;

  flt_copy(f1, &f3);

  while (i < (2 * MANT_SIZE)) {

    int count = -1;

    while (f3.sign != -1) {

      if (flt_sub(f3, f2, &f3) != OKAY) {

	return(EXP2BIG);

      }

      count++;

    }

    if (flt_add(f3, f2, &f3) != OKAY) {

      return(EXP2BIG);

    }

    rdbl.mant[i++] = count;

    if (f3.sign == 0) {

      break;

    }

    if ((--f2.exp) <= E_MIN) {

      return(EXP2BIG);

    }

  }

  while (i < (2 * MANT_SIZE)) {

    rdbl.mant[i++] = 0;

  }

  while (rdbl.mant[0] == 0) {

    for (i = 1; i < (2 * MANT_SIZE); i++) {

      rdbl.mant[i - 1] = rdbl.mant[i];

    }

    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;

    if ((--rdbl.exp) <= E_MIN) {

      return(EXP2BIG);

    }

  }

  dbl2float(rdbl, res);

  return(OKAY);

}

#else

int

flt_div(flt f1, flt f2, flt *res)

{

  Fdig a1[MANT_SIZE+1];

  Fdig a2[MANT_SIZE+1];

  Fdig r[MANT_SIZE+1];

  int bit_diff = 0;

  int i;

  int t, s;

  unsigned int c = 0;

  int bit, bitpos;

  int sg;

  int keep_on;

  int final_expt = f1.exp - f2.exp;

  int final_shift = 0;

  unsigned int k;

  if (f2.sign == 0) {

    return(DIVBY0);

  }

  if (f1.sign == 0) {

    flt_zero(res);

    return(OKAY);

  }

  for (i = 0; i < MANT_SIZE; ++i) {

    a1[i] = f1.mant[i];

    a2[i] = f2.mant[i];

    r[i] = 0;

  }

  a1[MANT_SIZE] = 0;

  a2[MANT_SIZE] = 0;

  r[MANT_SIZE] = 0;

  /* Shift first digit of a1 or a2 right until a1[0] < a2[0] and a1[0] >=

   * 2*a2[0] Count the places in bit_diff.  */

  t = (int)(a1[0]);

  s = (int)(a2[0]);

  if (t >= s) {

    for (; t >= s; ++bit_diff) {

      t >>= 1;

    }

  } else {

    for (; t < s; --bit_diff) {

      s >>= 1;

    }

    ++bit_diff;

  }

  /* Shift a1 bit_diff places right if bit_diff positive.

     Shift a2 -bit_diff places right if bit_diff negative.

   */

  if (bit_diff > 0) {

    for (i = 0; i < (MANT_SIZE + 1); ++i) {

      c = (unsigned int)(a1[i] + (c << 16));

      a1[i] = (unsigned short)(c >> bit_diff);	/* CAST:jmf: */

      c &= (unsigned int)((1 << (bit_diff + 1)) -1);

    }

  } else if (bit_diff < 0) {

    for (i = 0; i < (MANT_SIZE + 1); ++i) {

      c = (unsigned int)(a2[i] + (c << 16));

      a2[i] = (unsigned short)(c >> -bit_diff);	/* CAST:jmf: */

      c &= (unsigned int)((1 << (-bit_diff + 1)) -1);

    }

  }

  /* do the division */

  bit = 0;	/* current bit of result */

  bitpos = 0;	/* current digit of result */

  sg = 1;	/* 1 means subtract, 0 means add */

  keep_on = 1;

  while (keep_on) {

    c = 0;

    /* subtract or add */

    if (sg) {

      for (i = MANT_SIZE; i >= 0; --i) {

	c = (unsigned int)(a1[i] - a2[i] + c);

	a1[i] = (unsigned short)(c & 0xffff);	/* CAST:jmf: */

	if (c & 0x10000) {

	  c = (unsigned int)(-1);

	} else {

	  c = 0;

	}

      }

      sg = (c)? 0 : 1;

    } else {

      for (i = MANT_SIZE; i >= 0; --i) {

	c = (unsigned int)(a1[i] + a2[i] + c);

	a1[i] = (unsigned short)(c & 0xffff);	/* CAST:jmf: */

	if (c & 0x10000) {

	  c = 1;

	} else {

	  c = 0;

	}

      }

      sg = (c)? 1 : 0;

    }

    if (sg) {

      r[bitpos] = (unsigned short)((int)r[bitpos] | (1 << bit));

      /* CAST:jmf: */

    }

    if (bit == 0) {

      bit = 15;

      ++bitpos;

    } else {

      --bit;

    }

    keep_on = 0;

    /* shift f2 right one place, if zero keep_on = 0 */

    c = 0;

    for (i = 0; i < (MANT_SIZE + 1); ++i) {

      c = (unsigned int)(a2[i] + (c << 16));

      if (c) {

	keep_on = 1;

      }

      a2[i] = (unsigned short)(c >> 1);	/* CAST:jmf: */

      c &= 1;

    }

  }

  /* correct line-up of r */

  if (bit_diff > 0) {

    final_shift = bit_diff;

  } else {

    --final_expt;

    if (bit_diff > -16) {

      final_shift = 16 + bit_diff;

    }

  }

  k = (unsigned int)((((int)r[MANT_SIZE] << final_shift) & 0x8000) >> 15);

  /* (int) coercion OK because r is shorter */

  k = (unsigned int)(k + (r[MANT_SIZE] >> (16 - final_shift)));

  for (i = MANT_SIZE-1; i >= 0; --i) {

    k = (unsigned int)((r[i] << final_shift) + k);

    res ->mant[i] = (unsigned short)(k & 0xffff);;	/* CAST:jmf: */

    k >>= 16;

  }

  if (res->mant[0] == 0) {

    for (i = 0; i < MANT_SIZE - 1; ++i) {

      res->mant[i] = res->mant[i + 1];

    }

    res->mant[MANT_SIZE - 1] = 0;

    --final_expt;

  }