Subversion Repositories tendra.SVN

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 7u83 1 
/*
2 
    		 Crown Copyright (c) 1997
3 
 
4 
    This TenDRA(r) Computer Program is subject to Copyright
5 
    owned by the United Kingdom Secretary of State for Defence
6 
    acting through the Defence Evaluation and Research Agency
7 
    (DERA).  It is made available to Recipients with a
8 
    royalty-free licence for its use, reproduction, transfer
9 
    to other parties and amendment for any purpose not excluding
10 
    product development provided that any such use et cetera
11 
    shall be deemed to be acceptance of the following conditions:-
12 
 
13 
        (1) Its Recipients shall ensure that this Notice is
14 
        reproduced upon any copies or amended versions of it;
15 
 
16 
        (2) Any amended version of it shall be clearly marked to
17 
        show both the nature of and the organisation responsible
18 
        for the relevant amendment or amendments;
19 
 
20 
        (3) Its onward transfer from a recipient to another
21 
        party shall be deemed to be that party's acceptance of
22 
        these conditions;
23 
 
24 
        (4) DERA gives no warranty or assurance as to its
25 
        quality or suitability for any purpose and DERA accepts
26 
        no liability whatsoever in relation to any use to which
27 
        it may be put.
28 
*/
29 
 
30 
 
31 
/**********************************************************************
32 
$Author: release $
33 
$Date: 1998/01/17 15:55:47 $
34 
$Revision: 1.1.1.1 $
35 
$Log: flpt.c,v $
36 
 * Revision 1.1.1.1  1998/01/17  15:55:47  release
37 
 * First version to be checked into rolling release.
38 
 *
39 
 * Revision 1.4  1996/01/10  14:58:48  currie
40 
 * BIGEND var params chars & shorts
41 
 *
42 
 * Revision 1.3  1995/09/11  13:58:37  currie
43 
 * gcc pedantry
44 
 *
45 
 * Revision 1.2  1995/08/09  08:59:57  currie
46 
 * round bug
47 
 *
48 
 * Revision 1.1  1995/04/06  10:44:05  currie
49 
 * Initial revision
50 
 *
51 
***********************************************************************/
52 
 
53 
 
54 
 
55 
/***********************************************************************
56 
                                flpt.c
57 
 
58 
  Routines for handling the internal floating point representation.
59 
 
60 
 ***********************************************************************/
61 
#include "config.h"
62 
#include <limits.h>
63 
#include "common_types.h"
64 
#include "flpttypes.h"
65 
#include "exp.h"
66 
#include "xalloc.h"
67 
#include "szs_als.h"
68 
#include "messages_c.h"
69 
#include "basicread.h"
70 
#include "installglob.h"
71 
 
72 
 
73 
#include "flpt.h"
74 
 
75 
/* MACROS */
76 
 
77 
#define MAX_USEFUL_DECEXP 10000
78 
	/* plenty big enough for IEEE and VAX */
79 
#define initial_flpts 50
80 
#define extra_flpts 50
81 
 
82 
 
83 
 
84 
/***********************************************************************
85 
          ---------- Floating point emulator code ----------
86 
 NB - lines which assume the use of an ASCII character set are accompanied
87 
      by the comment ASCII
88 
 ***********************************************************************/
89 
 
90 
typedef struct _dbl {
91 
				/* Double precision type - used for
92 
				   internal working */
93 
  Fdig  mant[2 * MANT_SIZE];
94 
  int   sign;
95 
  int  exp;
96 
}                   dbl;
97 
 
98 
/* The current rounding rule. This may be one of the following :
99 
     R2ZERO - round towards zero.
100 
     R2PINF - round towards positive infinity.
101 
     R2NINF - round towards negative infinity.
102 
     R2NEAR - round to nearest integer.
103 
 It should be set using the "set_round" function.
104 
*/
105 
 
106 
/* IDENTITIES */
107 
 
108 
static Fdig bitround [16] = {0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80,
109 
			     0x100, 0x200, 0x400, 0x800, 0x1000, 0x2000,
110 
			     0x4000, 0x8000};
111 
static Fdig bitmask [17] = {0xffff, 0xfffe, 0xfffc, 0xfff8,
112 
			    0xfff0, 0xffe0, 0xffc0, 0xff80,
113 
			    0xff00, 0xfe00, 0xfc00, 0xf800,
114 
			    0xf000, 0xe000, 0xc000, 0x8000, 0x0};
115 
 
116 
static flt powers[16] = {
117 
	{{10, 0, 0, 0, 0, 0, 0, 0},1, 0},
118 
	{{100, 0, 0, 0, 0, 0, 0, 0},1, 0},
119 
	{{1000, 0, 0, 0, 0, 0, 0, 0},1, 0},
120 
	{{10000, 0, 0, 0, 0, 0, 0, 0},1, 0},
121 
	{{1, 34464, 0, 0, 0, 0, 0, 0},1, 1},
122 
	{{15, 16960, 0, 0, 0, 0, 0, 0},1, 1},
123 
	{{152, 38528, 0, 0, 0, 0, 0, 0},1, 1},
124 
	{{1525, 57600, 0, 0, 0, 0, 0, 0},1, 1},
125 
	{{15258, 51712, 0, 0, 0, 0, 0, 0},1, 1},
126 
	{{2, 21515, 58368, 0, 0, 0, 0, 0},1, 2},
127 
	{{23, 18550, 59392, 0, 0, 0, 0, 0},1, 2},
128 
	{{232, 54437, 4096, 0, 0, 0, 0, 0},1, 2},
129 
	{{2328, 20082, 40960, 0, 0, 0, 0, 0},1, 2},
130 
	{{23283, 4218, 16384, 0, 0, 0, 0, 0},1, 2},
131 
	{{3, 36222, 42182, 32768, 0, 0, 0, 0},1, 3},
132 
	{{35, 34546, 28609, 0, 0, 0, 0, 0},1, 3}
133 
};
134 
 
135 
static int two_powers[] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024,
136 
			    2048, 4096, 8192, 16384, 32768};
137 
static int  round_type = R2NEAR;
138 
 
139 
/* VARIABLES */
140 
/* All variables initialised */
141 
 
142 
int  tot_flpts;		/* total number of floating point numbers
143 
				*/
144 
int  flptfree;			/* floating point free list */
145 
int  flpt_left;		/* number of floating point numbers left
146 
				*/
147 
 
148 
flt * flptnos;			/* the extendable array of floating point
149 
				   numbers */
150 
 
151 
int fzero_no;
152 
int fone_no;
153 
 
154 
 
155 
 
156 
/* PROCEDURES */
157 
 
158 
flpt new_flpt PROTO_S ((void));
159 
 
160 
/* ----- internal functions ----- */
161 
 
162 
 
163 
/* "f1" and "f2" are two legal single
164 
   precision numbers. The result is -1, 0
165 
   or +1 depending on whether the
166 
   magnitude of "f1" is less than, equal
167 
   to or greater than that of "f2" */
168 
static int  mag_cmp
169 
    PROTO_N ( (f1, f2) )
170 
    PROTO_T ( flt f1 X flt f2 )
171 
{
172 
  int   i;
173 
  if (f1.sign == 0) {
174 
    if (f2.sign == 0) {
175 
      return (0);
176 
    }
177 
    return (-1);
178 
  }
179 
  if (f2.sign == 0) {
180 
    return (1);
181 
  }
182 
  if (f1.exp != f2.exp) {
183 
    return ((f1.exp < f2.exp) ? -1 : 1);
184 
  }
185 
  for (i = 0; i < MANT_SIZE; i++) {
186 
    if ((f1.mant[i]) != (f2.mant[i])) {
187 
      return (((f1.mant[i]) < (f2.mant[i])) ? -1 : 1);
188 
    }
189 
  }
190 
  return (0);
191 
}
192 
 
193 
 
194 
/* "d1" and "d2" are two legal double
195 
   precision numbers. The result is -1, 0
196 
   or +1 depending on whether the
197 
   magnitude of "d1" is less than, equal
198 
   to or greater than that of "d2" */
199 
static int  dbl_mag_cmp
200 
    PROTO_N ( (d1, d2) )
201 
    PROTO_T ( dbl d1 X dbl d2 )
202 
{
203 
  int   i;
204 
  if (d1.sign == 0) {
205 
    if (d2.sign == 0) {
206 
      return (0);
207 
    }
208 
    return (-1);
209 
  }
210 
  if (d2.sign == 0) {
211 
    return (1);
212 
  }
213 
  if (d1.exp != d2.exp) {
214 
    return ((d1.exp < d2.exp) ? -1 : 1);
215 
  }
216 
  for (i = 0; i < (2 * MANT_SIZE); i++) {
217 
    if ((d1.mant[i]) != (d2.mant[i])) {
218 
      return (((d1.mant[i]) < (d2.mant[i])) ? -1 : 1);
219 
    }
220 
  }
221 
  return (0);
222 
}
223 
 
224 
 
225 
/* "f" is a pointer to a legal single
226 
   precision number. "c" is a character
227 
   value from 0 to FBASE-1. On return, the
228 
   single precision number pointed to by
229 
   "f" will have been modified as if "c"
230 
   were the digit following the least
231 
   significant digit of the number pointed
232 
   to by "f" and the number has been
233 
   rounded according to the current
234 
   rounding rule */
235 
static  void flt_int_round
236 
    PROTO_N ( (f, ic) )
237 
    PROTO_T ( flt * f X unsigned int ic )
238 
{
239 
  int   i;
240 
  switch (round_type) {
241 
    default:
242 
    case R2ZERO:
243 
      return;
244 
    case R2PINF:
245 
      if (((f -> sign) != 1) || (ic == 0))
246 
	return;
247 
      break;
248 
    case R2NINF:
249 
      if (((f -> sign) != -1) || (ic == 0))
250 
	return;
251 
      break;
252 
    case R2NEAR:
253 
      if (ic < (unsigned int)(FBASE/2))
254 
	return;
255 
      break;
256 
  }
257 
  ic = 1;
258 
  for (i = (MANT_SIZE - 1); i >= 0; i--) {
259 
    (f -> mant)[i] = (unsigned short)
260 
			((ic += (unsigned int)((f -> mant)[i])) % FBASE);
261 
	/* CAST:jmf: unsigned short since FBASE <= MAX_USHORT+1 */
262 
    if ((ic /= FBASE) == 0) {
263 
      return;
264 
    }
265 
  }
266 
  (f -> exp)++;
267 
  (f -> mant)[0] = 1;
268 
  for (i = 1; i < MANT_SIZE; i++) {
269 
    (f -> mant)[i] = 0;
270 
  }
271 
}
272 
 
273 
 
274 
/* "d" is a legal double precision number.
275 
   "f" is a pointer to a single precision
276 
   number. On return, the number pointed
277 
   to by "f" contains the value of "d" but
278 
   to less precision. The current rounding
279 
   rule is used to truncate the mantissa
280 
*/
281 
static  void dbl2float
282 
    PROTO_N ( (d, f) )
283 
    PROTO_T ( dbl d X flt * f )
284 
{
285 
  int   i;
286 
  (f -> sign) = d.sign;
287 
  (f -> exp) = d.exp;
288 
  for (i = 0; i < MANT_SIZE; i++) {
289 
    (f -> mant)[i] = d.mant[i];
290 
  }
291 
  flt_int_round (f, (unsigned int)d.mant[MANT_SIZE]);
292 
}
293 
 
294 
 
295 
/* "f" is a legal single precision number.
296 
   "d" is a pointer to a double precision
297 
   number. On return, the number pointed
298 
   to by "d" will contain the value of "f"
299 
   */
300 
static  void flt2double
301 
    PROTO_N ( (f, d) )
302 
    PROTO_T ( flt f X dbl * d )
303 
{
304 
  int   i;
305 
  (d -> sign) = f.sign;
306 
  (d -> exp) = f.exp;
307 
  for (i = 0; i < MANT_SIZE; i++) {
308 
    (d -> mant)[i] = f.mant[i];
309 
  }
310 
  for (i = MANT_SIZE; i < (2 * MANT_SIZE); i++) {
311 
    (d -> mant)[i] = 0;
312 
  }
313 
}
314 
 
315 
 
316 
/* "d" is a pointer to a legal double
317 
   precision number. "dist" is a int with
318 
   value greater than zero. On return, the
319 
   number pointed to by "d" keeps the same
320 
   value ( although some loss of precision
321 
   may occur ) but is no longer in
322 
   normalised form. The mantissa is
323 
   shifted right "dist" places, and the
324 
   exponent is adjusted to keep the value
325 
   the same */
326 
static  void shift_right
327 
    PROTO_N ( (d, dist) )
328 
    PROTO_T ( dbl * d X int dist )
329 
{
330 
  int   i;
331 
  int  j = (2 * MANT_SIZE) - 1 - dist;
332 
  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--, j--) {
333 
    (d -> mant)[i] = (unsigned short)((j >= 0) ? ((d -> mant)[j]) : 0);
334 
  }
335 
  (d -> exp) += dist;
336 
}
337 
 
338 
 
339 
/* "d1" and "d2" are double precision
340 
   numbers, with the same exponent and
341 
   opposite (non-zero) signs. "res" is a
342 
   pointer to a single precision number,
343 
   whose exponent has been set to the same
344 
   value as "d1" and "d2". On return, the
345 
   number pointed to by "res" is the
346 
   result of adding "d1" to "d2" and
347 
   rounding it by the current rounding
348 
   rule. The return value is either OKAY
349 
   or EXP2BIG ( in which case the value of
350 
   the number pointed to by "res" is
351 
   undefined ) */
352 
static int  sub_mantissas
353 
    PROTO_N ( (d1, d2, res) )
354 
    PROTO_T ( dbl d1 X dbl d2 X flt * res )
355 
{
356 
  dbl * dp1 = &d1, *dp2 = &d2, rdbl;
357 
  unsigned int  c = 0;
358 
  int   i,
359 
        cmp = dbl_mag_cmp (d1, d2);
360 
  if (cmp == 0) {
361 
    (res -> sign) = 0;
362 
    (res -> exp) = 0;
363 
    for (i = 0; i < MANT_SIZE; i++) {
364 
      (res -> mant)[i] = 0;
365 
    }
366 
    return (OKAY);
367 
  }
368 
  else
369 
    if (cmp == -1) {
370 
      dp1 = &d2;
371 
      dp2 = &d1;
372 
    }
373 
  rdbl.exp = (res -> exp);
374 
  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {
375 
    int   s = (int)(((dp1 -> mant)[i]) - ((dp2 -> mant)[i]) -
376 
            c);
377 
    c = 0;
378 
    if (s < 0) {
379 
      c = 1;
380 
      s += FBASE;
381 
    }
382 
    rdbl.mant[i] = (unsigned short)s;	/* CAST:jmf: */
383 
  }
384 
  rdbl.sign = (dp1 -> sign);
385 
  while ((rdbl.mant[0]) == 0) {
386 
    for (i = 1; i < (2 * MANT_SIZE); i++) {
387 
      rdbl.mant[i - 1] = rdbl.mant[i];
388 
    }
389 
    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;
390 
    if ((--rdbl.exp) < E_MIN) {
391 
      return (EXP2BIG);
392 
    }
393 
  }
394 
  dbl2float (rdbl, res);
395 
  return (OKAY);
396 
}
397 
 
398 
 
399 
 
400 
/* "d1" and "d2" are double precision
401 
   numbers, with the same exponent and
402 
   non-zero signs. "res" is a pointer to a
403 
   single precision number which has the
404 
   same exponent as "d1" and "d2". On
405 
   return, the number pointed to by "res"
406 
   will be the result of adding "d1" to
407 
   "d2" and rounding it according to the
408 
   current rounding rule. The return value
409 
   will be either OKAY or EXP2BIG ( in
410 
   which case the value of the number
411 
   pointed to by "res" is undefined ) */
412 
static int  add_mantissas
413 
    PROTO_N ( (d1, d2, res) )
414 
    PROTO_T ( dbl d1 X dbl d2 X flt * res )
415 
{
416 
  int   i;
417 
  unsigned int  c;
418 
  if (d1.sign == d2.sign) {
419 
    dbl rdbl;
420 
    (res -> sign) = d1.sign;
421 
    flt2double (*res, &rdbl);
422 
    c = 0;
423 
    for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {
424 
      c += (unsigned int)(d1.mant[i] + d2.mant[i]);
425 
      rdbl.mant[i] = (unsigned short)(c % FBASE);	/* CAST:jmf: */
426 
      c /= FBASE;
427 
    }
428 
    if (c) {
429 
      for (i = (2 * MANT_SIZE) - 1; i > 0; i--) {
430 
	rdbl.mant[i] = rdbl.mant[i - 1];
431 
      }
432 
      rdbl.mant[0] = (unsigned short)c;	/* CAST:jmf: */
433 
      if (((++rdbl.exp) >= E_MAX) ||
434 
	  (rdbl.exp <= E_MIN)) {
435 
	return (EXP2BIG);
436 
      }
437 
    }
438 
    dbl2float (rdbl, res);
439 
    return (OKAY);
440 
  }
441 
  return (sub_mantissas (d1, d2, res));
442 
}
443 
 
444 
 
445 
/* ----- "public" emulator functions ----- */
446 
 
447 
/* "f" is a pointer to a single precision
448 
   number. On return, the number pointed
449 
   to by "f" will have the value zero */
450 
void flt_zero
451 
    PROTO_N ( (f) )
452 
    PROTO_T ( flt * f )
453 
{
454 
  int   i;
455 
  (f -> sign) = 0;
456 
  (f -> exp) = 0;
457 
  for (i = 0; i < MANT_SIZE; i++) {
458 
    (f -> mant)[i] = 0;
459 
  }
460 
}
461 
 
462 
 
463 
/* "f" is a legal single precision number.
464 
   "res" is a pointer to a single
465 
   precision number. On return, the number
466 
   pointed to by "res" is the same as "f"
467 
*/
468 
void flt_copy
469 
    PROTO_N ( (f, res) )
470 
    PROTO_T ( flt f X flt *res )
471 
{
472 
  int   i;
473 
  (res -> sign) = f.sign;
474 
  (res -> exp) = f.exp;
475 
  for (i = 0; i < MANT_SIZE; i++) {
476 
    (res -> mant)[i] = f.mant[i];
477 
  }
478 
  return;
479 
}
480 
 
481 
 
482 
/* "f1" and "f2" are legal single
483 
   precision numbers. "res" is a pointer
484 
   to a single precision number. On
485 
   return, if the result is OKAY, then the
486 
   number pointed to by "res" will contain
487 
   the value of adding "f1" to "f2". If
488 
   the result is EXP2BIG, then the value
489 
   of the number pointed to by "res" is
490 
   undefined */
491 
int   flt_add
492 
    PROTO_N ( (f1, f2, res) )
493 
    PROTO_T ( flt f1 X flt f2 X flt *res )
494 
{
495 
  dbl d1, d2;
496 
  if (f1.sign == 0) {
497 
    flt_copy (f2, res);
498 
    return (OKAY);
499 
  }
500 
  else
501 
    if (f2.sign == 0) {
502 
      flt_copy (f1, res);
503 
      return (OKAY);
504 
    }
505 
  flt2double (f1, &d1);
506 
  flt2double (f2, &d2);
507 
  if (d1.exp < d2.exp) {
508 
    shift_right (&d1, d2.exp - d1.exp);
509 
  }
510 
  else
511 
    if (d1.exp > d2.exp) {
512 
      shift_right (&d2, d1.exp - d2.exp);
513 
    }
514 
  (res -> exp) = d1.exp;
515 
  return (add_mantissas (d1, d2, res));
516 
}
517 
 
518 
 
519 
/* "f1" and "f2" are legal single
520 
   precision numbers. "res" is a pointer
521 
   to a single precision number. On
522 
   return, if the result is OKAY, then the
523 
   number pointed to by "res" will contain
524 
   the value of subtracting "f2" from
525 
   "f1". If the result is EXP2BIG, then
526 
   the value of the number pointed to by
527 
   "res" is undefined */
528 
int   flt_sub
529 
    PROTO_N ( (f1, f2, res) )
530 
    PROTO_T ( flt f1 X flt f2 X flt *res )
531 
{
532 
  f2.sign = -f2.sign;
533 
  return (flt_add (f1, f2, res));
534 
}
535 
 
536 
#if FBASE == 10
537 
 
538 
/* "f1" and "f2" are legal single
539 
   precision numbers. "res" is a pointer
540 
   to a single precision number. On
541 
   return, if the result is OKAY, then the
542 
   number pointed to by "res" will contain
543 
   the value of multiplying "f1" by "f2".
544 
   If the result is EXP2BIG, then the
545 
   value of the number pointed to by "res"
546 
   is undefined */
547 
int   flt_mul
548 
    PROTO_N ( (f1, f2, res) )
549 
    PROTO_T ( flt f1 X flt f2 X flt *res )
550 
{
551 
  dbl rdbl;
552 
  int   i,
553 
        j;
554 
  unsigned int c = 0;
555 
  unsigned int acc[2 * MANT_SIZE];
556 
  if ((f1.sign == 0) || (f2.sign == 0)) {
557 
    flt_zero (res);
558 
    return (OKAY);
559 
  }
560 
  rdbl.sign = ((f1.sign == f2.sign) ? 1 : -1);
561 
  rdbl.exp = (f1.exp + f2.exp + 1);
562 
  if ((rdbl.exp >= E_MAX) ||
563 
      (rdbl.exp <= E_MIN)) {
564 
    return (EXP2BIG);
565 
  }
566 
  for (i = 0; i < (2 * MANT_SIZE); i++) {
567 
    acc[i] = 0;
568 
  }
569 
  for (i = MANT_SIZE - 1; i >= 0; i--) {
570 
    for (j = MANT_SIZE - 1; j >= 0; j--) {
571 
      acc[i + j + 1] += (f1.mant[i] * f2.mant[j]);
572 
    }
573 
  }
574 
  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {
575 
    c += acc[i];
576 
    rdbl.mant[i] = (c % FBASE);
577 
    c /= FBASE;
578 
  }
579 
  while (rdbl.mant[0] == 0) {
580 
    for (i = 1; i < (2 * MANT_SIZE); i++) {
581 
      rdbl.mant[i - 1] = rdbl.mant[i];
582 
    }
583 
    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;
584 
    if ((--rdbl.exp) <= E_MIN) {
585 
      return (EXP2BIG);
586 
    }
587 
  }
588 
  dbl2float (rdbl, res);
589 
  return (OKAY);
590 
}
591 
#else
592 
/* "f1" and "f2" are legal single
593 
   precision numbers. "res" is a pointer
594 
   to a single precision number. On
595 
   return, if the result is OKAY, then the
596 
   number pointed to by "res" will contain
597 
   the value of multiplying "f1" by "f2".
598 
   If the result is EXP2BIG, then the
599 
   value of the number pointed to by "res"
600 
   is undefined */
601 
int   flt_mul
602 
    PROTO_N ( (f1, f2, res) )
603 
    PROTO_T ( flt f1 X flt f2 X flt *res )
604 
{
605 
  dbl rdbl;
606 
  int   i,
607 
        j,
608 
	k;
609 
  unsigned int acc[2 * MANT_SIZE];
610 
  if ((f1.sign == 0) || (f2.sign == 0)) {
611 
    flt_zero (res);
612 
    return (OKAY);
613 
  }
614 
  rdbl.sign = ((f1.sign == f2.sign) ? 1 : -1);
615 
  rdbl.exp = (f1.exp + f2.exp + 1);
616 
  if ((rdbl.exp >= E_MAX) ||
617 
      (rdbl.exp <= E_MIN)) {
618 
    return (EXP2BIG);
619 
  }
620 
  for (i = 0; i < (2 * MANT_SIZE); i++) {
621 
    acc[i] = 0;
622 
  }
623 
  SET(acc);
624 
 
625 
  for (i = MANT_SIZE - 1; i >= 0; i--) {
626 
    for (j = MANT_SIZE - 1; j >= 0; j--) {
627 
      unsigned int temp = (unsigned int)(f1.mant[i] * f2.mant[j]);
628 
      unsigned int atl;
629 
      for (k = i + j + 1; temp != 0; --k) {
630 
	atl = acc[k] + (temp & 0xffff);
631 
	acc[k] = atl & 0xffff;
632 
	temp = (atl >> 16) + (temp >> 16);
633 
      };
634 
    };
635 
  };
636 
 
637 
  for (i = (2 * MANT_SIZE) - 1; i >= 0; i--) {
638 
    rdbl.mant[i] = (unsigned short)(acc[i]);	/* CAST:jmf: */
639 
  };
640 
 
641 
  while (rdbl.mant[0] == 0) {
642 
    for (i = 1; i < (2 * MANT_SIZE); i++) {
643 
      rdbl.mant[i - 1] = rdbl.mant[i];
644 
    }
645 
    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;
646 
    if ((--rdbl.exp) <= E_MIN) {
647 
      return (EXP2BIG);
648 
    }
649 
  };
650 
 
651 
  dbl2float (rdbl, res);
652 
  return (OKAY);
653 
}
654 
#endif
655 
 
656 
#if FBASE == 10
657 
 
658 
/* "f1" and "f2" are legal single
659 
   precision numbers, and "f2" is
660 
   non-zero. "res" is a pointer to a
661 
   single precision number. On return, if
662 
   the result is OKAY, then the number
663 
   pointed to by "res" contains the result
664 
   of dividing "f1" by "f2". If the result
665 
   is EXP2BIG or DIVBY0 then the contents
666 
   of the number pointed to by "res" are
667 
   undefined */
668 
int   flt_div
669 
    PROTO_N ( (f1, f2, res) )
670 
    PROTO_T ( flt f1 X flt f2 X flt *res )
671 
{
672 
  dbl rdbl;
673 
  int   i = 0;
674 
  flt f3;
675 
  if (f2.sign == 0) {
676 
    return (DIVBY0);
677 
  }
678 
  if (f1.sign == 0) {
679 
    flt_zero (res);
680 
    return (OKAY);
681 
  }
682 
  rdbl.sign = ((f1.sign == f2.sign) ? 1 : -1);
683 
  rdbl.exp = (f1.exp - f2.exp);
684 
  if ((rdbl.exp >= E_MAX) ||
685 
      (rdbl.exp <= E_MIN)) {
686 
    return (EXP2BIG);
687 
  }
688 
  f1.sign = 1;
689 
  f2.sign = 1;
690 
  f1.exp = 0;
691 
  f2.exp = 0;
692 
  flt_copy (f1, &f3);
693 
  while (i < (2 * MANT_SIZE)) {
694 
    int   count = -1;
695 
    while (f3.sign != -1) {
696 
      if (flt_sub (f3, f2, &f3) != OKAY) {
697 
	return (EXP2BIG);
698 
      }
699 
      count++;
700 
    }
701 
    if (flt_add (f3, f2, &f3) != OKAY) {
702 
      return (EXP2BIG);
703 
    }
704 
    rdbl.mant[i++] = count;
705 
    if (f3.sign == 0) {
706 
      break;
707 
    }
708 
    if ((--f2.exp) <= E_MIN) {
709 
      return (EXP2BIG);
710 
    }
711 
  }
712 
  while (i < (2 * MANT_SIZE)) {
713 
    rdbl.mant[i++] = 0;
714 
  }
715 
  while (rdbl.mant[0] == 0) {
716 
    for (i = 1; i < (2 * MANT_SIZE); i++) {
717 
      rdbl.mant[i - 1] = rdbl.mant[i];
718 
    }
719 
    rdbl.mant[(2 * MANT_SIZE) - 1] = 0;
720 
    if ((--rdbl.exp) <= E_MIN) {
721 
      return (EXP2BIG);
722 
    }
723 
  }
724 
  dbl2float (rdbl, res);
725 
  return (OKAY);
726 
}
727 
 
728 
#else
729 
 
730 
int   flt_div
731 
    PROTO_N ( (f1, f2, res) )
732 
    PROTO_T ( flt f1 X flt f2 X flt *res )
733 
{
734 
  Fdig a1[MANT_SIZE+1];
735 
  Fdig a2[MANT_SIZE+1];
736 
  Fdig r[MANT_SIZE+1];
737 
  int bit_diff = 0;
738 
  int i;
739 
  int t, s;
740 
  unsigned int c = 0;
741 
  int bit, bitpos;
742 
  int sg;
743 
  int keep_on;
744 
  int final_expt = f1.exp - f2.exp;
745 
  int final_shift = 0;
746 
  unsigned int k;
747 
 
748 
  if (f2.sign == 0) {
749 
    return (DIVBY0);
750 
  };
751 
 
752 
  if (f1.sign == 0) {
753 
    flt_zero (res);
754 
    return (OKAY);
755 
  };
756 
 
757 
  for (i = 0; i < MANT_SIZE; ++i) {
758 
    a1[i] = f1.mant[i];
759 
    a2[i] = f2.mant[i];
760 
    r[i] = 0;
761 
  };
762 
  a1[MANT_SIZE] = 0;
763 
  a2[MANT_SIZE] = 0;
764 
  r[MANT_SIZE] = 0;
765 
 
766 
	/* Shift first digit of a1 or a2 right until
767 
	   a1[0] < a2[0] and a1[0] >= 2*a2[0]
768 
	   Count the places in bit_diff.
769 
	*/
770 
  t = (int)(a1[0]);
771 
  s = (int)(a2[0]);
772 
  if (t >= s) {
773 
    for (; t >= s; ++bit_diff)
774 
	t >>= 1;
775 
  }
776 
  else {
777 
    for (; t < s; --bit_diff)
778 
	s >>= 1;
779 
    ++bit_diff;
780 
  };
781 
 
782 
 
783 
	/* Shift a1 bit_diff places right if bit_diff positive.
784 
	   Shift a2 -bit_diff places right if bit_diff negative.
785 
	*/
786 
  if (bit_diff > 0) {
787 
    for (i = 0; i < (MANT_SIZE+1); ++i) {
788 
      c = (unsigned int)(a1[i] + (c << 16));
789 
      a1[i] = (unsigned short)(c >> bit_diff);	/* CAST:jmf: */
790 
      c &= (unsigned int)((1 << (bit_diff+1)) -1);
791 
    };
792 
  }
793 
  else
794 
  if (bit_diff < 0) {
795 
    for (i = 0; i < (MANT_SIZE+1); ++i) {
796 
      c = (unsigned int)(a2[i] + (c << 16));
797 
      a2[i] = (unsigned short)(c >> -bit_diff);	/* CAST:jmf: */
798 
      c &= (unsigned int)((1 << (-bit_diff+1)) -1);
799 
    };
800 
  };
801 
 
802 
	/* do the division */
803 
  bit = 0;	/* current bit of result */
804 
  bitpos = 0;	/* current digit of result */
805 
  sg = 1;	/* 1 means subtract, 0 means add */
806 
  keep_on = 1;
807 
  while (keep_on) {
808 
    c = 0;
809 
	/* subtract or add */
810 
    if (sg) {
811 
      for (i = MANT_SIZE; i >= 0; --i) {
812 
	c = (unsigned int)(a1[i] - a2[i] + c);
813 
	a1[i] = (unsigned short)(c & 0xffff);	/* CAST:jmf: */
814 
	if (c & 0x10000)
815 
	  c = (unsigned int)(-1);
816 
	else
817 
	  c = 0;
818 
      };
819 
      sg = (c) ? 0 : 1;
820 
    }
821 
    else {
822 
      for (i = MANT_SIZE; i >= 0; --i) {
823 
	c = (unsigned int)(a1[i] + a2[i] + c);
824 
	a1[i] = (unsigned short)(c & 0xffff);	/* CAST:jmf: */
825 
	if (c & 0x10000)
826 
	  c = 1;
827 
	else
828 
	  c = 0;
829 
      };
830 
      sg = (c) ? 1 : 0;
831 
    };
832 
    if (sg) {
833 
      r[bitpos] = (unsigned short)((int)r[bitpos] | (1 << bit));
834 
		/* CAST:jmf: */
835 
    };
836 
    if (bit == 0) {
837 
      bit = 15;
838 
      ++bitpos;
839 
    }
840 
    else
841 
      --bit;
842 
    keep_on = 0;
843 
	/* shift f2 right one place, if zero keep_on = 0 */
844 
    c = 0;
845 
    for (i = 0; i < (MANT_SIZE+1); ++i) {
846 
      c = (unsigned int)(a2[i] + (c << 16));
847 
      if (c)
848 
	keep_on = 1;
849 
      a2[i] = (unsigned short)(c >> 1);	/* CAST:jmf: */
850 
      c &= 1;
851 
    };
852 
  };
853 
 
854 
	/* correct line-up of r */
855 
  if (bit_diff > 0) {
856 
    final_shift = bit_diff;
857 
  }
858 
  else {
859 
    --final_expt;
860 
    if (bit_diff > -16)
861 
      final_shift = 16 + bit_diff;
862 
  };
863 
  k = (unsigned int)((((int)r[MANT_SIZE] << final_shift) & 0x8000) >> 15);
864 
	/* (int) coercion OK because r is shorter */
865 
  k = (unsigned int)(k + (r[MANT_SIZE] >> (16 - final_shift)));
866 
  for (i = MANT_SIZE-1; i >= 0; --i) {
867 
    k = (unsigned int)((r[i] << final_shift) + k);
868 
    res ->mant[i] = (unsigned short)(k & 0xffff);;	/* CAST:jmf: */
869 
    k >>= 16;
870 
  };
871 
 
872 
  if (res->mant[0] == 0) {
873 
    for (i = 0; i < MANT_SIZE-1; ++i)
874 
      res->mant[i] = res->mant[i+1];
875 
      res->mant[MANT_SIZE-1] = 0;
876 
    --final_expt;
877 
  };
878 
 
879 
  res -> exp = final_expt;
880 
  res -> sign = (f1.sign == f2.sign) ? 1 : -1;
881 
  return OKAY;
882 
}
883 
 
884 
#endif
885 
 
886 
/* "f1" and "f2" are legal single
887 
   precision numbers. The return value
888 
   will be -1, 0 or 1 depending whether
889 
   "f1" is less than, equal to or greater
890 
   than "f2" */
891 
int   flt_cmp
892 
    PROTO_N ( (f1, f2) )
893 
    PROTO_T ( flt f1 X flt f2 )
894 
{
895 
  int   ret;
896 
  if (f1.sign < f2.sign) {
897 
    return (-1);
898 
  }
899 
  if (f1.sign > f2.sign) {
900 
    return (1);
901 
  }
902 
  ret = mag_cmp (f1, f2);
903 
  return ((f1.sign == -1) ? -ret : ret);
904 
}
905 
 
906 
 
907 
/* "f" is a legal single precision number.
908 
   "res" is a pointer to a single
909 
   precision number. On return, the number
910 
   pointed to by "res" will be the integer
911 
   value of "f" rounded using the current
912 
   rounding rule */
913 
void flt_round
914 
    PROTO_N ( (f, res) )
915 
    PROTO_T ( flt f X flt *res )
916 
{
917 
  int   i;
918 
  unsigned int  ex;
919 
  if (f.exp < -1) {
920 
    flt_zero (res);
921 
    return;
922 
  }
923 
  flt_copy (f, res);
924 
  if ((f.exp + 1) >= MANT_SIZE) {
925 
    return;
926 
  }
927 
  ex = ((res -> mant)[f.exp + 1]);
928 
  for (i = f.exp + 1; i < MANT_SIZE; i++) {
929 
    (res -> mant)[i] = 0;
930 
  }
931 
  if (f.exp == -1) {
932 
    (res -> sign) = 0;
933 
    (res -> exp) = 0;
934 
  }
935 
  switch (round_type) {
936 
    default:
937 
    case R2ZERO:
938 
      return;
939 
    case R2PINF:
940 
      if ((f.sign != 1) || (ex == 0))
941 
	return;
942 
      break;
943 
    case R2NINF:
944 
      if ((f.sign != -1) || (ex == 0))
945 
	return;
946 
      break;
947 
    case R2NEAR:
948 
      if (ex < (unsigned int)(FBASE/2))
949 
	return;
950 
      break;
951 
  }
952 
  ex = 1;
953 
  for (i = f.exp; i >= 0; i--) {
954 
    (res -> mant)[i] =
955 
	(unsigned short)((ex = ex +
956 
			(unsigned int)((res -> mant)[i])) % FBASE);
957 
		/* CAST:jmf */
958 
    if ((ex /= FBASE) == 0) {
959 
      return;
960 
    }
961 
  }
962 
  if (f.exp != -1) {
963 
    (res -> exp)++;
964 
  }
965 
  else {
966 
    (res -> sign) = f.sign;
967 
  }
968 
  (res -> mant)[0] = 1;
969 
  for (i = 1; i < MANT_SIZE; i++) {
970 
    (res -> mant)[i] = 0;
971 
  }
972 
}
973 
 
974 
 
975 
/* "f" is a legal single precision number.
976 
   "res" is a pointer to a single
977 
   precision number. On return, the number
978 
   pointed to by "res" will be the integer
979 
   value of "f" rounded towards zero */
980 
void flt_trunc
981 
    PROTO_N ( (f, res) )
982 
    PROTO_T ( flt f X flt *res )
983 
{
984 
  int   i;
985 
  if (f.exp < 0) {
986 
    flt_zero (res);
987 
    return;
988 
  }
989 
  flt_copy (f, res);
990 
  for (i = f.exp + 1; i < MANT_SIZE; i++) {
991 
    (res -> mant)[i] = 0;
992 
  }
993 
}
994 
 
995 
 
996 
/* "s" is a pointer to an array of characters.
997 
   "f" is a pointer to a single precision number.
998 
   "ret" is NULL, or a pointer to a pointer to a character.
999 
   On return, if the return value is OKAY, then the number pointed to
1000 
   by "f"
1001 
   is the floating point number represented in the string "s". If "ret" was
1002 
   not NULL, then it will point to a pointer to the next character
1003 
   in "s" not
1004 
   used in the number. If the return value is SYNTAX or EXP2BIG, then the
1005 
   value of the number pointed to by "f" is undefined. In this case,
1006 
   if "ret"
1007 
   was not NULL, then it will point to a pointer to the start of
1008 
   the string.
1009 
   Leading white space in the string will be ignored. Numbers have the
1010 
   following form :
1011 
       [+-]?(([0-9]+(\.[0-9]*)?)|([0-9]*\.[0-9]+))([Ee][+-]?[0-9]+)?
1012 
*/
1013 
 
1014 
 
1015 
#if FBASE == 10
1016 
 
1017 
int   str2flt
1018 
    PROTO_N ( (s, f, r) )
1019 
    PROTO_T ( char *s X flt * f X char **r )
1020 
{
1021 
  int   i = 0,
1022 
        ids = -1,
1023 
        rounded = 0,
1024 
        mant_empty = 1,
1025 
        zero = 1;
1026 
  (f -> sign) = 1;
1027 
  (f -> exp) = 0;
1028 
  if (r) {
1029 
    *r = s;
1030 
  }
1031 
  while ((*s) && (' ' == (*s))) {
1032 
    s++;
1033 
  }
1034 
  if (*s == '-') {
1035 
    (f -> sign) = -1;
1036 
    s++;
1037 
  }
1038 
  else
1039 
    if (*s == '+') {
1040 
      s++;
1041 
    }
1042 
  while (*s == '0') {
1043 
    mant_empty = 0;
1044 
    s++;
1045 
  }
1046 
  while (*s >= '0' && *s <= '9') {
1047 
    mant_empty = 0;
1048 
    zero = 0;
1049 
    if (i < MANT_SIZE) {
1050 
      (f -> mant)[i++] = ((*s) - '0');/* ASCII */
1051 
    }
1052 
    else
1053 
      if (!rounded) {
1054 
	flt_int_round (f, (*s) - '0');
1055 
	rounded = 1;		/* ASCII */
1056 
      }
1057 
    if (ids >= E_MAX) {
1058 
      return (EXP2BIG);
1059 
    }
1060 
    ids++;
1061 
    s++;
1062 
  }
1063 
  if (*s == '.') {
1064 
    s++;
1065 
    if (zero) {
1066 
      while (*s == '0') {
1067 
	if (ids <= E_MIN) {
1068 
	  return (EXP2BIG);
1069 
	}
1070 
	ids--;
1071 
	s++;
1072 
	mant_empty = 0;
1073 
      }
1074 
    }
1075 
    while (*s >= '0' && *s <= '9') {
1076 
      if (i < MANT_SIZE) {
1077 
	(f -> mant)[i++] = ((*s) - '0');/* ASCII */
1078 
      }
1079 
      else
1080 
	if (!rounded) {
1081 
	  flt_int_round (f, (*s) - '0');
1082 
	  rounded = 1;		/* ASCII */
1083 
	}
1084 
      s++;
1085 
      zero = 0;
1086 
      mant_empty = 0;
1087 
    }
1088 
  }
1089 
  while (i < MANT_SIZE) {
1090 
    (f -> mant)[i++] = 0;
1091 
  }
1092 
  if ((*s == 'E') || (*s == 'e')) {
1093 
    int   e_sign = 1,
1094 
          exp_empty = 1;
1095 
    int  e = 0;
1096 
    if (*++s == '-') {
1097 
      e_sign = -1;
1098 
      s++;
1099 
    }
1100 
    else
1101 
      if (*s == '+') {
1102 
	s++;
1103 
      }
1104 
    while ((*s) && (*s >= '0' && *s <= '9')) {
1105 
      if (e >= (E_MAX / 10)) {
1106 
	return (EXP2BIG);
1107 
      }
1108 
      e = e * 10 + ((*s++) - '0');
1109 
      exp_empty = 0;		/* ASCII */
1110 
    }
1111 
    if (exp_empty) {
1112 
      return (SYNTAX);
1113 
    }
1114 
    e *= e_sign;
1115 
    (f -> exp) += (e + ids);
1116 
  }
1117 
  else {
1118 
    (f -> exp) += ids;
1119 
  }
1120 
  if (((f -> exp) >= E_MAX) ||
1121 
      ((f -> exp) <= E_MIN)) {
1122 
    return (EXP2BIG);
1123 
  }
1124 
  if (zero) {
1125 
    (f -> sign) = 0;
1126 
    (f -> exp) = 0;
1127 
  }
1128 
  if (mant_empty) {
1129 
    return (SYNTAX);
1130 
  }
1131 
  if (r) {
1132 
    *r = s;
1133 
  }
1134 
  return (OKAY);
1135 
}
1136 
 
1137 
#endif
1138 
 
1139 
 
1140 
/***********************************************************************
1141 
          ========== interface to  TDF translator ==========
1142 
 ***********************************************************************/
1143 
 
1144 
/* memory allocation */
1145 
 
1146 
 
1147 
void init_flpt
1148 
    PROTO_Z ()
1149 
{
1150 
  /* initialise */
1151 
  int  i;
1152 
  flt * fzr;
1153 
  flt * forf;
1154 
  flptnos = (flt *) xcalloc (initial_flpts, sizeof (flt));
1155 
  tot_flpts = initial_flpts;
1156 
  for (i = 1; i < tot_flpts; ++i)
1157 
    flptnos[i].exp = i - 1;
1158 
  flptfree = tot_flpts - 1;
1159 
  flpt_left = tot_flpts;
1160 
 
1161 
  fzero_no = new_flpt();
1162 
  fone_no = new_flpt();
1163 
  fzr = &flptnos[fzero_no];
1164 
  fzr -> sign = 0;
1165 
  fzr -> exp = 0;
1166 
  forf = &flptnos[fone_no];
1167 
  forf -> sign = 1;
1168 
  forf -> exp = 0;
1169 
  for (i = 0; i < MANT_SIZE; i++) {
1170 
    (fzr -> mant)[i] = 0;
1171 
    (forf -> mant)[i] = 0;
1172 
  };
1173 
  (forf -> mant)[0] = 1;
1174 
 
1175 
  return;
1176 
}
1177 
 
1178 
void more_flpts
1179 
    PROTO_Z ()
1180 
{
1181 
  /* extend the floating point array */
1182 
  int  i;
1183 
  tot_flpts += extra_flpts;
1184 
  flptnos = (flt *) xrealloc ((void*)(CH flptnos),
1185 
  				TOSZ (tot_flpts * (int)sizeof (flt)));
1186 
  for (i = tot_flpts - extra_flpts + 1; i < tot_flpts; ++i)
1187 
    flptnos[i].exp = i - 1;
1188 
  flptfree = tot_flpts - 1;
1189 
  flpt_left = extra_flpts;
1190 
  return;
1191 
}
1192 
 
1193 
flpt new_flpt
1194 
    PROTO_Z ()
1195 
{
1196 
  /* allocate space for a new floating point number */
1197 
  flpt r = flptfree;
1198 
  flptfree = flptnos[r].exp;
1199 
  if (--flpt_left == 0)
1200 
    more_flpts ();
1201 
  return (r);
1202 
}
1203 
 
1204 
void flpt_ret
1205 
    PROTO_N ( (f) )
1206 
    PROTO_T ( flpt f )
1207 
{
1208 
  /* return a floating point number to free */
1209 
  flptnos[f].exp = flptfree;
1210 
  flptfree = f;
1211 
  ++flpt_left;
1212 
  return;
1213 
}
1214 
 
1215 
 
1216 
/***********************************************************************/
1217 
 
1218 
 
1219 
/* do the test testno on and and b deliver
1220 
   1 if true, 0 otherwise */
1221 
int cmpflpt
1222 
    PROTO_N ( (a, b, testno) )
1223 
    PROTO_T ( flpt a X flpt b X int testno )
1224 
{
1225 
  int   res = flt_cmp (flptnos[a], flptnos[b]);
1226 
 
1227 
  switch (testno) {
1228 
    case 1:
1229 
      return (res == 1);
1230 
    case 2:
1231 
      return (res != -1);
1232 
    case 3:
1233 
      return (res == -1);
1234 
    case 4:
1235 
      return (res != 1);
1236 
    case 5:
1237 
      return (res == 0);
1238 
    default:
1239 
      return (res != 0);
1240 
  };
1241 
}
1242 
 
1243 
 
1244 
#if FBASE == 10
1245 
 
1246 
flpt floatrep
1247 
    PROTO_N ( (n) )
1248 
    PROTO_T ( int n )
1249 
{
1250 
  flpt res = new_flpt ();
1251 
  char *pr = intchars (n);
1252 
	/* this is wrong for unsigned integers */
1253 
  flt fr;
1254 
  IGNORE str2flt (pr, &fr, (char **) 0);
1255 
  flptnos[res] = fr;
1256 
  return (res);
1257 
}
1258 
 
1259 
flpt floatrep_unsigned
1260 
    PROTO_N ( (n) )
1261 
    PROTO_T ( unsigned int n )
1262 
{
1263 
  failer("floatrep_unsigned not used");
1264 
  return 0;
1265 
}
1266 
 
1267 
int flpt_bits
1268 
    PROTO_N ( (fv) )
1269 
    PROTO_T ( floating_variety fv )
1270 
{
1271 
  return 0;
1272 
}
1273 
 
1274 
void flpt_round
1275 
    PROTO_N ( (round_t, posn, res) )
1276 
    PROTO_T ( int round_t X int posn X flt * res )
1277 
{
1278 
  return;
1279 
}
1280 
 
1281 
#else
1282 
 
1283 
static flpt floatrep_aux
1284 
    PROTO_N ( (n, sign) )
1285 
    PROTO_T ( int n X int sign )
1286 
{
1287 
  flpt res = new_flpt ();
1288 
  flt fr;
1289 
  unsigned int mask = 0xffff;
1290 
  int i;
1291 
  int supp = 1;
1292 
  int index = 0;
1293 
 
1294 
  flt_zero(&fr);
1295 
  if (n == 0)  {
1296 
    flptnos[res] = fr;
1297 
    return (res);
1298 
  };
1299 
 
1300 
  fr.sign = sign;
1301 
 
1302 
  for (i = 1; i >= 0; --i)  {
1303 
    int t = (int)((n >> (i*16)) & (int)mask);
1304 
    if (supp && t)  {
1305 
      supp = 0;
1306 
      fr.exp = i;
1307 
    };
1308 
    if (!supp) {
1309 
      fr.mant[index++] = (unsigned short)t;	/* CAST:jmf: */
1310 
    };
1311 
  };
1312 
 
1313 
  flptnos[res] = fr;
1314 
  return (res);
1315 
}
1316 
 
1317 
flpt floatrep
1318 
    PROTO_N ( (n) )
1319 
    PROTO_T ( int n )
1320 
{
1321 
  if (n < 0)
1322 
    return floatrep_aux(-n, -1);
1323 
  return floatrep_aux(n, 1);
1324 
}
1325 
 
1326 
flpt floatrep_unsigned
1327 
    PROTO_N ( (n) )
1328 
    PROTO_T ( unsigned int n )
1329 
{
1330 
  return floatrep_aux((int)n, 1);
1331 
}
1332 
 
1333 
void flpt_newdig
1334 
    PROTO_N ( (dig, res, base) )
1335 
    PROTO_T ( unsigned int dig X flt * res X int base )
1336 
{
1337 
  int i;
1338 
  unsigned int c = 0;
1339 
 
1340 
  if (dig != 0)
1341 
    res -> sign = 1;
1342 
 
1343 
  i = MANT_SIZE - 1;
1344 
 
1345 
  for ( ; i >= 0; --i)  {
1346 
    c = (unsigned int)(((int)res -> mant[i] * base) + (int)c);
1347 
    res -> mant[i] = (unsigned short)(c % FBASE);	/* CAST:jmf: */
1348 
    c = c / FBASE;
1349 
  };
1350 
  if (c) {
1351 
    ++res -> exp;
1352 
    i = res -> exp;
1353 
    if (i >= MANT_SIZE)
1354 
      i = MANT_SIZE - 1;
1355 
    for ( ; i >= 0; --i)
1356 
      res -> mant[i] = res -> mant[i-1];
1357 
    res -> mant[0] = (unsigned short)c;	/* CAST:jmf: */
1358 
  };
1359 
 
1360 
  if (res -> exp >= MANT_SIZE)
1361 
    return;
1362 
 
1363 
  c = dig;
1364 
  for (i = res -> exp; c && i >= 0; --i)  {
1365 
    c = (unsigned int)(res -> mant[i] + c);
1366 
    res -> mant[i] = (unsigned short)(c % FBASE);	/* CAST:jmf: */
1367 
    c = c / FBASE;
1368 
  };
1369 
  if (c) {
1370 
    ++res -> exp;
1371 
    i = res -> exp;
1372 
    if (i >= MANT_SIZE)
1373 
      i = MANT_SIZE - 1;
1374 
    for ( ; i >= 0; --i)
1375 
      res -> mant[i] = res -> mant[i-1];
1376 
    res -> mant[0] = (unsigned short)c;	/* CAST:jmf: */
1377 
  };
1378 
  return;
1379 
}
1380 
 
1381 
 
1382 
void flpt_scale
1383 
    PROTO_N ( (expt, res, base) )
1384 
    PROTO_T ( int expt X flt * res X int base )
1385 
{
1386 
  flt ft;
1387 
 
1388 
  if (base == 10) {
1389 
    if (expt > 0) {
1390 
      if (expt > MAX_USEFUL_DECEXP) {
1391 
  	failer(BIG_FLPT);
1392 
	exit(EXIT_FAILURE);
1393 
  	/* UNREACHED */
1394 
      };
1395 
      while (expt > 16) {
1396 
        IGNORE flt_mul(*res, powers[15], &ft); /* cannot fail */
1397 
        flt_copy(ft, res);
1398 
        expt -= 16;
1399 
      };
1400 
      IGNORE flt_mul(*res, powers[expt-1], &ft); /* cannot fail */
1401 
      flt_copy(ft, res);
1402 
    }
1403 
    else
1404 
    if (expt < 0) {
1405 
      if (expt < - MAX_USEFUL_DECEXP) {
1406 
	  flt_zero(res);
1407 
	  return;
1408 
      };
1409 
      while (expt < -16) {
1410 
        IGNORE flt_div(*res, powers[15], &ft); /* cannot fail */
1411 
        flt_copy(ft, res);
1412 
        expt += 16;
1413 
      };
1414 
      IGNORE flt_div(*res, powers[-1 - expt], &ft); /* cannot fail */
1415 
      flt_copy(ft, res);
1416 
    };
1417 
  }
1418 
  else  {
1419 
    if (base == 4)
1420 
      expt += expt;
1421 
    if (base == 8)
1422 
      expt *= 3;
1423 
    if (base == 16)
1424 
      expt *= 4;
1425 
 
1426 
    if (expt > 0) {
1427 
      res->exp += (expt / 16);
1428 
      expt = expt % 16;
1429 
      if (expt != 0)
1430 
        flpt_newdig((unsigned int)0, res, two_powers[expt]);
1431 
      return;
1432 
    }
1433 
    else
1434 
    if (expt < 0) {
1435 
      int temp = ((-expt) / 16);
1436 
      expt = (-expt) % 16;
1437 
      if (expt == 0)
1438 
	res->exp -= temp;
1439 
      else  {
1440 
        flpt_newdig((unsigned int)0, res, two_powers[16 - expt]);
1441 
	res->exp -= (temp+1);
1442 
      };
1443 
    };
1444 
  };
1445 
  return;
1446 
}
1447 
 
1448 
/* posn is the number of bits to be left correct after the
1449 
   rounding operation */
1450 
void flpt_round
1451 
    PROTO_N ( (round_t, posn, res) )
1452 
    PROTO_T ( int round_t X int posn X flt * res )
1453 
{
1454 
  unsigned int c = res -> mant[0];
1455 
  unsigned int mask;
1456 
  int ndig0 = 0;
1457 
  int i, bitpos;
1458 
  int digpos = 0;
1459 
  int bits_discarded = 0;
1460 
 
1461 
  if (res -> sign == 0)
1462 
    return;
1463 
 
1464 
  for (mask = FBASE - 1; mask & c; mask <<= 1)
1465 
   ++ndig0;
1466 
 
1467 
  bitpos = ndig0 - posn;
1468 
 
1469 
  while (bitpos < 1)  {
1470 
    bitpos += FBITS;
1471 
    ++digpos;
1472 
  };
1473 
 
1474 
  --bitpos;
1475 
 
1476 
  /* digpos holds the number of the FBASE digit to be rounded,
1477 
     bitpos holds the number of the bit within that digit to have
1478 
     one added
1479 
  */
1480 
 
1481 
  for (i = digpos + 1; i < MANT_SIZE; ++i)  {
1482 
    if (res -> mant[i]) {
1483 
      bits_discarded = 1;
1484 
      res -> mant[i] = 0;
1485 
    };
1486 
  };
1487 
 
1488 
  switch (round_t) {
1489 
    default:
1490 
    case R2ZERO:
1491 
      res -> mant[digpos] =
1492 
	(unsigned short)(res -> mant[digpos] & bitmask[bitpos+1]);
1493 
      return;
1494 
    case R2PINF:
1495 
      if (res -> sign != 1 ||
1496 
	    (!bits_discarded && ((res -> mant[digpos] & bitmask[bitpos+1])
1497 
		 == (int)res -> mant[digpos]))) {
1498 
        res -> mant[digpos] =
1499 
		(unsigned short)(res -> mant[digpos] & bitmask[bitpos+1]);
1500 
 
1501 
        return;
1502 
      };
1503 
      res -> mant[digpos] =
1504 
	(unsigned short)(res -> mant[digpos] | bitround[bitpos]);
1505 
      break;
1506 
    case R2NINF:
1507 
      if (res -> sign == 1 ||
1508 
	    (!bits_discarded && ((res -> mant[digpos] & bitmask[bitpos+1])
1509 
		 == (int)res -> mant[digpos]))) {
1510 
        res -> mant[digpos] =
1511 
	(unsigned short)((int)res -> mant[digpos] & (int)bitmask[bitpos+1]);
1512 
 
1513 
        return;
1514 
      };
1515 
      res -> mant[digpos] =
1516 
	(unsigned short)((int)res -> mant[digpos] | (int)bitround[bitpos]);
1517 
      break;
1518 
    case R2NEAR:
1519 
      break;
1520 
   };
1521 
 
1522 
   res -> mant[digpos] =
1523 
	(unsigned short)((int)res -> mant[digpos] & (int)bitmask[bitpos]);
1524 
 
1525 
   c = bitround[bitpos];
1526 
    for (i = digpos; c && i >= 0; --i)  {
1527 
        c = (unsigned int)((int)res -> mant[i] + (int)c);
1528 
        res -> mant[i] = (unsigned short)(c % FBASE);	/* CAST:jmf: */
1529 
        c = c / FBASE;
1530 
    };
1531 
   if (c) {
1532 
        ++res -> exp;
1533 
        i = res -> exp;
1534 
        if (i >= MANT_SIZE)
1535 
          i = MANT_SIZE - 1;
1536 
        for ( ; i >= 0; --i)
1537 
          res -> mant[i] = res -> mant[i-1];
1538 
        res -> mant[0] = (unsigned short)c;
1539 
   };
1540 
  res -> mant[digpos+(int)c] =
1541 
	(unsigned short)((int)res -> mant[digpos+(int)c] &
1542 
				(int)bitmask[bitpos+1]);
1543 
  return;
1544 
}
1545 
 
1546 
int flpt_bits
1547 
    PROTO_N ( (fv) )
1548 
    PROTO_T ( floating_variety fv )
1549 
{
1550 
  switch (fv)
1551 
   {
1552 
     case 0: return FLOAT_BITS;
1553 
				/* FLOAT_BITS is defined in config.h
1554 
				   24 for IEEE */
1555 
     case 1: return DOUBLE_BITS;
1556 
				/* FLOAT_BITS is defined in config.h
1557 
				   53 for IEEE */
1558 
     case 2: return LDOUBLE_BITS;
1559 
				/* FLOAT_BITS is defined in config.h
1560 
				   64 for IEEE */
1561 
   };
1562 
  return 0;
1563 
}
1564 
 
1565 
int flpt_round_to_integer
1566 
    PROTO_N ( (rndmd, f) )
1567 
    PROTO_T ( int rndmd X flt * f )
1568 
{
1569 
  unsigned int c = f -> mant[0];
1570 
  unsigned int mask;
1571 
  int ndig0 = 0;
1572 
  int intbits;
1573 
  int res;
1574 
 
1575 
  if (f -> sign == 0)
1576 
    return 0;
1577 
 
1578 
  for (mask = FBASE - 1; mask & c; mask <<= 1)
1579 
   ++ndig0;
1580 
  intbits = ndig0 + (f -> exp * FBITS);
1581 
 
1582 
  if (intbits >= 0)
1583 
    flpt_round(rndmd, intbits, f);
1584 
  else {
1585 
    switch(rndmd) EXHAUSTIVE {
1586 
      case R2ZERO:
1587 
        break;
1588 
      case R2PINF:
1589 
	if (f->sign == 1)
1590 
	  flt_copy(flptnos[fone_no], f);
1591 
	else
1592 
	  flt_zero(f);
1593 
        break;
1594 
      case R2NINF:
1595 
	if (f->sign == -1) {
1596 
	  flt_copy(flptnos[fone_no], f);
1597 
	  f->sign = -1;
1598 
	}
1599 
	else
1600 
	  flt_zero(f);
1601 
        break;
1602 
      case R2NEAR:
1603 
        break;
1604 
    };
1605 
  };
1606 
 
1607 
 
1608 
/*
1609 
#if has64bits
1610 
  if (f -> exp > 3) {
1611 
    int ij;
1612 
    f -> mant[0] = f -> mant[f -> exp - 3];
1613 
    f -> mant[1] = f -> mant[f -> exp - 2];
1614 
    f -> mant[2] = f -> mant[f -> exp - 1];
1615 
    f -> mant[3] = f -> mant[f -> exp];
1616 
    for (ij = 4; ij <= f -> exp; ++ij) f -> mant[ij] = 0;
1617 
    f -> exp = 3;
1618 
  };
1619 
#else
1620 
  if (f -> exp > 1) {
1621 
    f -> mant[0] = f -> mant[f -> exp - 1];
1622 
    f -> mant[1] = f -> mant[f -> exp];
1623 
    f -> exp = 1;
1624 
  };
1625 
#endif
1626 
*/
1627 
  if (f -> exp == 1)
1628 
    res = ((int)f -> mant[1] + ((int)f -> mant[0] << 16));
1629 
  else
1630 
    res = (int)(f -> mant[0]);
1631 
  return res * f -> sign;  /* this is only valid for 32-bit results */
1632 
 
1633 
}
1634 
 
1635 
/* convert flt* to IEEE representation in ints.
1636 
   sw is 0 for single, 1 for double,
1637 
   2 for extended */
1638 
r2l real2longs_IEEE
1639 
    PROTO_N ( (fp, sw) )
1640 
    PROTO_T ( flt * fp X int sw )
1641 
{
1642 
  r2l res;
1643 
  flt f;
1644 
  unsigned int c = fp -> mant[0];
1645 
  unsigned int mask;
1646 
  int ndig0 = 0;
1647 
  int expt;
1648 
  int bias, expt_size, precision;
1649 
  unsigned int sig1, sig2 = 0, sig3 = 0, sig4 = 0;
1650 
  int i, index;
1651 
  res.i1 = 0;
1652 
  res.i2 = 0;
1653 
  res.i3 = 0;
1654 
  res.i4 = 0;
1655 
 
1656 
  if (fp -> sign == 0) {
1657 
    return res;
1658 
  };
1659 
 
1660 
  f = *fp;
1661 
 
1662 
  switch (sw) EXHAUSTIVE {
1663 
    case 0:
1664 
      bias = 127;
1665 
      expt_size = 8;
1666 
      precision = 24;
1667 
      break;
1668 
    case 1:
1669 
      bias = 1023;
1670 
      expt_size = 11;
1671 
      precision = 53;
1672 
      break;
1673 
    case 2:
1674 
      bias = 16383;
1675 
      expt_size = 15;
1676 
      precision = LDOUBLE_BITS;
1677 
      break;
1678 
  };
1679 
 
1680 
    for (mask = FBASE - 1; mask & c; mask <<= 1)
1681 
     ++ndig0;
1682 
    expt = ndig0 + (f.exp * FBITS) - 1;
1683 
 
1684 
    if (expt < (1-bias-precision)) return res; /* 0 - underflowed */
1685 
 
1686 
    if (expt > bias) {
1687 
      if (flpt_const_overflow_fail) {
1688 
        failer(BIG_FLPT);
1689 
	exit(EXIT_FAILURE);
1690 
      };
1691 
      switch (sw) {
1692 
	case 0:
1693 
	  if (f.sign == -1)
1694 
	    res.i1 = 0x80000000;
1695 
	  res.i1 += 0x7f800000;
1696 
	  return res;
1697 
	case 1:
1698 
	  if (f.sign == -1)
1699 
	    res.i2 = 0x80000000;
1700 
	  res.i2 += 0x7ff00000;
1701 
	  return res;
1702 
	case 2:
1703 
#if is80x86
1704 
	if (f.sign == -1)
1705 
	  res.i3 = 0x8000;
1706 
	res.i3 += 0x7fff;
1707 
	res.i2 = 0x80000000;
1708 
	return res;
1709 
#else
1710 
#if issparc || ishppa
1711 
	if (f.sign == -1)
1712 
	  res.i4 = 0x8000;
1713 
	res.i4 += 0x7fff;
1714 
	return res;
1715 
#else
1716 
      failer("long double not implemented");
1717 
      return res;
1718 
#endif
1719 
#endif
1720 
      };
1721 
    };
1722 
 
1723 
    expt += bias;
1724 
    i = precision - ndig0;;
1725 
    sig1 = f.mant[0];
1726 
    index = 1;
1727 
    while (i >= 16)  {
1728 
      sig4 <<= 16;
1729 
      sig4 += (sig3 >> 16);
1730 
      sig3 <<= 16;
1731 
      sig3 += (sig2 >> 16);
1732 
      sig2 <<= 16;
1733 
      sig2 += (sig1 >> 16);
1734 
      sig1 <<= 16;
1735 
      sig1 = (unsigned int)(sig1 + f.mant[index++]);
1736 
      i -= 16;
1737 
    };
1738 
    if (i != 0) {
1739 
      sig4 <<= i;
1740 
      sig4 += (sig3 >> (32-i));
1741 
      sig3 <<= i;
1742 
      sig3 += (sig2 >> (32-i));
1743 
      sig2 <<= i;
1744 
      sig2 += (sig1 >> (32-i));
1745 
      sig1 <<= i;
1746 
      sig1 = (unsigned int)(sig1 + (f.mant[index] >> (16-i)));
1747 
    };
1748 
 
1749 
    if (expt < 1) {
1750 
      int places = 1 - expt;
1751 
      while (places >= 32) {
1752 
	sig1 = sig2;
1753 
        sig2 = sig3;
1754 
	sig3 = sig4;
1755 
	sig4 = 0;
1756 
	places -= 32;
1757 
      };
1758 
      if (places > 0) {
1759 
        sig1 = (sig1 >> places) + (sig2 << (32 - places));
1760 
        sig2 = (sig2 >> places) + (sig3 << (32 - places));
1761 
        sig3 = (sig3 >> places) + (sig4 << (32 - places));
1762 
        sig4 >>= places;
1763 
      };
1764 
      expt = 0;
1765 
    }
1766 
    else
1767 
      expt &= ((1 << (expt_size+1)) -1);
1768 
 
1769 
    switch (sw) {
1770 
      case 0:
1771 
	if (f.sign == -1)
1772 
	  res.i1 = 0x80000000;
1773 
	res.i1 += (expt << 23);
1774 
	res.i1 += (int)(sig1 & 0x007fffff);
1775 
	break;
1776 
      case 1:
1777 
	if (f.sign == -1)
1778 
	  res.i2 = 0x80000000;
1779 
	res.i2 += (expt << 20);
1780 
	res.i2 += (int)(sig2 & 0xfffff);
1781 
	res.i1 = (int)sig1;
1782 
	break;
1783 
      case 2:
1784 
#if is80x86
1785 
	res.i1 = (int)sig1;
1786 
	res.i2 = (int)sig2;
1787 
	if (f.sign == -1)
1788 
	  res.i3 = 0x8000;
1789 
	res.i3 += expt;
1790 
	UNUSED(sig3);
1791 
	UNUSED(sig4);
1792 
	break;
1793 
#else
1794 
#if issparc || ishppa
1795 
	if (f.sign == -1)
1796 
	  res.i4 = 0x80000000;
1797 
	res.i4 += (expt << 16);
1798 
	res.i4 += (int)(sig4 & 0xffff);
1799 
	res.i3 = (int)sig3;
1800 
	res.i2 = (int)sig2;
1801 
	res.i1 = (int)sig1;
1802 
	break;
1803 
#else
1804 
      failer("long double not implemented");
1805 
      return res;
1806 
#endif
1807 
#endif
1808 
    };
1809 
 
1810 
  return res;
1811 
}
1812 
 
1813 
 
1814 
#endif
1815