OpenCores
URL https://opencores.org/ocsvn/eco32/eco32/trunk

Subversion Repositories eco32

[/] [eco32/] [trunk/] [fp/] [implementation/] [ultimate2/] [fp.c] - Blame information for rev 153

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 15 hellwig 
/*
2 
 * fp.c -- floating point arithmetic
3 
 */
4 
 
5 
 
6 
#include "fp.h"
7 
 
8 
 
9 
#define sign_bit ((unsigned) 0x80000000)
10 
 
11 
#define ROUND_OFF 1
12 
#define ROUND_UP 2
13 
#define ROUND_DOWN 3
14 
#define ROUND_NEAR 4
15 
 
16 
#define X_BIT (1<<8)
17 
#define Z_BIT (1<<9)
18 
#define U_BIT (1<<10)
19 
#define O_BIT (1<<11)
20 
#define I_BIT (1<<12)
21 
#define W_BIT (1<<13)
22 
#define V_BIT (1<<14)
23 
#define D_BIT (1<<15)
24 
#define E_BIT (1<<18)
25 
 
26 
#define zero_exponent (-1000)
27 
 
28 
#define bignum_prec 157
29 
 
30 
#define magic_offset 2112
31 
#define origin 37
32 
 
33 
#define buf0 (buf+8)
34 
#define buf_max (buf+777)
35 
 
36 
#include <stdio.h>
37 
#include <string.h>
38 
#include <ctype.h>
39 
 
40 
typedef struct {
41 
  int a;
42 
  int b;
43 
  tetra dat[bignum_prec];
44 
} bignum;
45 
 
46 
octa zero_octa;
47 
octa neg_one = { -1, -1 };
48 
octa inf_octa = { 0x7ff00000, 0 };
49 
octa standard_NaN = { 0x7ff80000, 0 };
50 
octa aux;
51 
bool overflow;
52 
 
53 
int cur_round;
54 
 
55 
int exceptions;
56 
 
57 
octa val;
58 
char *next_char;
59 
 
60 
static char buf[785] = "00000000";
61 
 
62 
 
63 
static void bignum_times_ten (bignum *);
64 
static void bignum_dec (bignum *, bignum *, tetra);
65 
static int bignum_compare (bignum *, bignum *);
66 
static void bignum_double (bignum *);
67 
 
68 
 
69 
octa oplus(octa y, octa z) {
70 
  octa x;
71 
 
72 
  x.h = y.h + z.h;
73 
  x.l = y.l + z.l;
74 
  if (x.l < y.l)
75 
    x.h++;
76 
  return x;
77 
}
78 
 
79 
 
80 
octa ominus(octa y, octa z) {
81 
  octa x;
82 
 
83 
  x.h = y.h - z.h;
84 
  x.l = y.l - z.l;
85 
  if (x.l > y.l)
86 
    x.h--;
87 
  return x;
88 
}
89 
 
90 
 
91 
octa incr(octa y, int delta) {
92 
  octa x;
93 
 
94 
  x.h = y.h;
95 
  x.l = y.l + delta;
96 
  if (delta >= 0)
97 
    {
98 
      if (x.l < y.l)
99 
        x.h++;
100 
    }
101 
  else if (x.l > y.l)
102 
    x.h--;
103 
  return x;
104 
}
105 
 
106 
 
107 
octa shift_left(octa y, int s) {
108 
  while (s >= 32)
109 
    y.h = y.l, y.l = 0, s -= 32;
110 
  if (s)
111 
    {
112 
      tetra yhl = y.h << s, ylh = y.l >> (32 - s);
113 
      y.h = yhl + ylh;
114 
      y.l <<= s;
115 
    }
116 
  return y;
117 
}
118 
 
119 
 
120 
octa shift_right(octa y, int s, int u) {
121 
  while (s >= 32)
122 
    y.l = y.h, y.h = (u ? 0 : -(y.h >> 31)), s -= 32;
123 
  if (s)
124 
    {
125 
      tetra yhl = y.h << (32 - s), ylh = y.l >> s;
126 
      y.h = (u ? 0 : (-(y.h >> 31)) << (32 - s)) + (y.h >> s);
127 
      y.l = yhl + ylh;
128 
    }
129 
  return y;
130 
}
131 
 
132 
 
133 
octa omult(octa y, octa z) {
134 
  int i, j, k;
135 
  tetra u[4], v[4], w[8];
136 
  tetra t;
137 
  octa acc;
138 
 
139 
  u[3] = y.h >> 16, u[2] = y.h & 0xffff, u[1] = y.l >> 16, u[0] =
140 
    y.l & 0xffff;
141 
  v[3] = z.h >> 16, v[2] = z.h & 0xffff, v[1] = z.l >> 16, v[0] =
142 
    z.l & 0xffff;
143 
 
144 
  for (j = 0; j < 4; j++)
145 
    w[j] = 0;
146 
  for (j = 0; j < 4; j++)
147 
    if (!v[j])
148 
      w[j + 4] = 0;
149 
    else
150 
      {
151 
        for (i = k = 0; i < 4; i++)
152 
          {
153 
            t = u[i] * v[j] + w[i + j] + k;
154 
            w[i + j] = t & 0xffff, k = t >> 16;
155 
          }
156 
        w[j + 4] = k;
157 
      }
158 
 
159 
  aux.h = (w[7] << 16) + w[6], aux.l = (w[5] << 16) + w[4];
160 
  acc.h = (w[3] << 16) + w[2], acc.l = (w[1] << 16) + w[0];
161 
 
162 
  return acc;
163 
}
164 
 
165 
 
166 
octa odiv(octa x, octa y, octa z) {
167 
  int i, j, k, n, d;
168 
  tetra u[8], v[4], q[4], mask, qhat, rhat, vh, vmh;
169 
  tetra t;
170 
  octa acc;
171 
 
172 
  if (x.h > z.h || (x.h == z.h && x.l >= z.l))
173 
    {
174 
      aux = y;
175 
      return x;
176 
    }
177 
 
178 
  u[7] = x.h >> 16, u[6] = x.h & 0xffff, u[5] = x.l >> 16, u[4] =
179 
    x.l & 0xffff;
180 
  u[3] = y.h >> 16, u[2] = y.h & 0xffff, u[1] = y.l >> 16, u[0] =
181 
    y.l & 0xffff;
182 
  v[3] = z.h >> 16, v[2] = z.h & 0xffff, v[1] = z.l >> 16, v[0] =
183 
    z.l & 0xffff;
184 
 
185 
  for (n = 4; v[n - 1] == 0; n--);
186 
 
187 
  vh = v[n - 1];
188 
  for (d = 0; vh < 0x8000; d++, vh <<= 1);
189 
  for (j = k = 0; j < n + 4; j++)
190 
    {
191 
      t = (u[j] << d) + k;
192 
      u[j] = t & 0xffff, k = t >> 16;
193 
    }
194 
  for (j = k = 0; j < n; j++)
195 
    {
196 
      t = (v[j] << d) + k;
197 
      v[j] = t & 0xffff, k = t >> 16;
198 
    }
199 
  vh = v[n - 1];
200 
  vmh = (n > 1 ? v[n - 2] : 0);
201 
 
202 
  for (j = 3; j >= 0; j--)
203 
 
204 
    {
205 
 
206 
      t = (u[j + n] << 16) + u[j + n - 1];
207 
      qhat = t / vh, rhat = t - vh * qhat;
208 
      if (n > 1)
209 
        while (qhat == 0x10000 || qhat * vmh > (rhat << 16) + u[j + n - 2])
210 
          {
211 
            qhat--, rhat += vh;
212 
            if (rhat >= 0x10000)
213 
              break;
214 
          }
215 
 
216 
      for (i = k = 0; i < n; i++)
217 
        {
218 
          t = u[i + j] + 0xffff0000 - k - qhat * v[i];
219 
          u[i + j] = t & 0xffff, k = 0xffff - (t >> 16);
220 
        }
221 
 
222 
      if (u[j + n] != k)
223 
        {
224 
          qhat--;
225 
          for (i = k = 0; i < n; i++)
226 
            {
227 
              t = u[i + j] + v[i] + k;
228 
              u[i + j] = t & 0xffff, k = t >> 16;
229 
            }
230 
        }
231 
 
232 
      q[j] = qhat;
233 
    }
234 
 
235 
  mask = (1 << d) - 1;
236 
  for (j = 3; j >= n; j--)
237 
    u[j] = 0;
238 
  for (k = 0; j >= 0; j--)
239 
    {
240 
      t = (k << 16) + u[j];
241 
      u[j] = t >> d, k = t & mask;
242 
    }
243 
 
244 
  acc.h = (q[3] << 16) + q[2], acc.l = (q[1] << 16) + q[0];
245 
  aux.h = (u[3] << 16) + u[2], aux.l = (u[1] << 16) + u[0];
246 
 
247 
  return acc;
248 
}
249 
 
250 
 
251 
octa
252 
fpack (f, e, s, r)
253 
     octa f;
254 
     int e;
255 
     char s;
256 
     int r;
257 
{
258 
  octa o;
259 
  if (e > 0x7fd)
260 
    e = 0x7ff, o = zero_octa;
261 
  else
262 
    {
263 
      if (e < 0)
264 
        {
265 
          if (e < -54)
266 
            o.h = 0, o.l = 1;
267 
          else
268 
            {
269 
              octa oo;
270 
              o = shift_right (f, -e, 1);
271 
              oo = shift_left (o, -e);
272 
              if (oo.l != f.l || oo.h != f.h)
273 
                o.l |= 1;
274 
 
275 
            }
276 
          e = 0;
277 
        }
278 
      else
279 
        o = f;
280 
    }
281 
 
282 
  if (o.l & 3)
283 
    exceptions |= X_BIT;
284 
  switch (r)
285 
    {
286 
    case ROUND_DOWN:
287 
      if (s == '-')
288 
        o = incr (o, 3);
289 
      break;
290 
    case ROUND_UP:
291 
      if (s != '-')
292 
        o = incr (o, 3);
293 
    case ROUND_OFF:
294 
      break;
295 
    case ROUND_NEAR:
296 
      o = incr (o, o.l & 4 ? 2 : 1);
297 
      break;
298 
    }
299 
  o = shift_right (o, 2, 1);
300 
  o.h += e << 20;
301 
  if (o.h >= 0x7ff00000)
302 
    exceptions |= O_BIT + X_BIT;
303 
  else if (o.h < 0x100000)
304 
    exceptions |= U_BIT;
305 
  if (s == '-')
306 
    o.h |= sign_bit;
307 
  return o;
308 
 
309 
  ;
310 
}
311 
 
312 
tetra
313 
sfpack (f, e, s, r)
314 
     octa f;
315 
     int e;
316 
     char s;
317 
     int r;
318 
{
319 
  tetra o;
320 
  if (e > 0x47d)
321 
    e = 0x47f, o = 0;
322 
  else
323 
    {
324 
      o = shift_left (f, 3).h;
325 
      if (f.l & 0x1fffffff)
326 
        o |= 1;
327 
      if (e < 0x380)
328 
        {
329 
          if (e < 0x380 - 25)
330 
            o = 1;
331 
          else
332 
            {
333 
              tetra o0, oo;
334 
              o0 = o;
335 
              o = o >> (0x380 - e);
336 
              oo = o << (0x380 - e);
337 
              if (oo != o0)
338 
                o |= 1;
339 
 
340 
            }
341 
          e = 0x380;
342 
        }
343 
    }
344 
 
345 
  if (o & 3)
346 
    exceptions |= X_BIT;
347 
  switch (r)
348 
    {
349 
    case ROUND_DOWN:
350 
      if (s == '-')
351 
        o += 3;
352 
      break;
353 
    case ROUND_UP:
354 
      if (s != '-')
355 
        o += 3;
356 
    case ROUND_OFF:
357 
      break;
358 
    case ROUND_NEAR:
359 
      o += (o & 4 ? 2 : 1);
360 
      break;
361 
    }
362 
  o = o >> 2;
363 
  o += (e - 0x380) << 23;
364 
  if (o >= 0x7f800000)
365 
    exceptions |= O_BIT + X_BIT;
366 
  else if (o < 0x100000)
367 
    exceptions |= U_BIT;
368 
  if (s == '-')
369 
    o |= sign_bit;
370 
  return o;
371 
 
372 
  ;
373 
}
374 
 
375 
ftype
376 
funpack (x, f, e, s)
377 
     octa x;
378 
     octa *f;
379 
     int *e;
380 
     char *s;
381 
{
382 
  int ee;
383 
  exceptions = 0;
384 
  *s = (x.h & sign_bit ? '-' : '+');
385 
  *f = shift_left (x, 2);
386 
  f->h &= 0x3fffff;
387 
  ee = (x.h >> 20) & 0x7ff;
388 
  if (ee)
389 
    {
390 
      *e = ee - 1;
391 
      f->h |= 0x400000;
392 
      return (ee < 0x7ff ? num : f->h == 0x400000 && !f->l ? inf : nan);
393 
    }
394 
  if (!x.l && !f->h)
395 
    {
396 
      *e = zero_exponent;
397 
      return zro;
398 
    }
399 
  do
400 
    {
401 
      ee--;
402 
      *f = shift_left (*f, 1);
403 
    }
404 
  while (!(f->h & 0x400000));
405 
  *e = ee;
406 
  return num;
407 
}
408 
 
409 
ftype
410 
sfunpack (x, f, e, s)
411 
     tetra x;
412 
     octa *f;
413 
     int *e;
414 
     char *s;
415 
{
416 
  int ee;
417 
  exceptions = 0;
418 
  *s = (x & sign_bit ? '-' : '+');
419 
  f->h = (x >> 1) & 0x3fffff, f->l = x << 31;
420 
  ee = (x >> 23) & 0xff;
421 
  if (ee)
422 
    {
423 
      *e = ee + 0x380 - 1;
424 
      f->h |= 0x400000;
425 
      return (ee < 0xff ? num : (x & 0x7fffffff) == 0x7f800000 ? inf : nan);
426 
    }
427 
  if (!(x & 0x7fffffff))
428 
    {
429 
      *e = zero_exponent;
430 
      return zro;
431 
    }
432 
  do
433 
    {
434 
      ee--;
435 
      *f = shift_left (*f, 1);
436 
    }
437 
  while (!(f->h & 0x400000));
438 
  *e = ee + 0x380;
439 
  return num;
440 
}
441 
 
442 
octa
443 
load_sf (z)
444 
     tetra z;
445 
{
446 
  octa f, x;
447 
  int e;
448 
  char s;
449 
  ftype t;
450 
  t = sfunpack (z, &f, &e, &s);
451 
  switch (t)
452 
    {
453 
    case zro:
454 
      x = zero_octa;
455 
      break;
456 
    case num:
457 
      return fpack (f, e, s, ROUND_OFF);
458 
    case inf:
459 
      x = inf_octa;
460 
      break;
461 
    case nan:
462 
      x = shift_right (f, 2, 1);
463 
      x.h |= 0x7ff00000;
464 
      break;
465 
    }
466 
  if (s == '-')
467 
    x.h |= sign_bit;
468 
  return x;
469 
}
470 
 
471 
tetra
472 
store_sf (x)
473 
     octa x;
474 
{
475 
  octa f;
476 
  tetra z;
477 
  int e;
478 
  char s;
479 
  ftype t;
480 
  t = funpack (x, &f, &e, &s);
481 
  switch (t)
482 
    {
483 
    case zro:
484 
      z = 0;
485 
      break;
486 
    case num:
487 
      return sfpack (f, e, s, cur_round);
488 
    case inf:
489 
      z = 0x7f800000;
490 
      break;
491 
    case nan:
492 
      if (!(f.h & 0x200000))
493 
        {
494 
          f.h |= 0x200000;
495 
          exceptions |= I_BIT;
496 
        }
497 
      z = 0x7f800000 | (f.h << 1) | (f.l >> 31);
498 
      break;
499 
    }
500 
  if (s == '-')
501 
    z |= sign_bit;
502 
  return z;
503 
}
504 
 
505 
octa
506 
fmult (y, z)
507 
     octa y, z;
508 
{
509 
  ftype yt, zt;
510 
  int ye, ze;
511 
  char ys, zs;
512 
  octa x, xf, yf, zf;
513 
  int xe;
514 
  char xs;
515 
  yt = funpack (y, &yf, &ye, &ys);
516 
  zt = funpack (z, &zf, &ze, &zs);
517 
  xs = ys + zs - '+';
518 
  switch (4 * yt + zt)
519 
    {
520 
 
521 
    case 4 * nan + nan:
522 
      if (!(y.h & 0x80000))
523 
        exceptions |= I_BIT;
524 
    case 4 * zro + nan:
525 
    case 4 * num + nan:
526 
    case 4 * inf + nan:
527 
      if (!(z.h & 0x80000))
528 
        exceptions |= I_BIT, z.h |= 0x80000;
529 
      return z;
530 
    case 4 * nan + zro:
531 
    case 4 * nan + num:
532 
    case 4 * nan + inf:
533 
      if (!(y.h & 0x80000))
534 
        exceptions |= I_BIT, y.h |= 0x80000;
535 
      return y;
536 
 
537 
      ;
538 
    case 4 * zro + zro:
539 
    case 4 * zro + num:
540 
    case 4 * num + zro:
541 
      x = zero_octa;
542 
      break;
543 
    case 4 * num + inf:
544 
    case 4 * inf + num:
545 
    case 4 * inf + inf:
546 
      x = inf_octa;
547 
      break;
548 
    case 4 * zro + inf:
549 
    case 4 * inf + zro:
550 
      x = standard_NaN;
551 
      exceptions |= I_BIT;
552 
      break;
553 
    case 4 * num + num:
554 
 
555 
      xe = ye + ze - 0x3fd;
556 
      x = omult (yf, shift_left (zf, 9));
557 
      if (aux.h >= 0x400000)
558 
        xf = aux;
559 
      else
560 
        xf = shift_left (aux, 1), xe--;
561 
      if (x.h || x.l)
562 
        xf.l |= 1;
563 
      return fpack (xf, xe, xs, cur_round);
564 
 
565 
      ;
566 
    }
567 
  if (xs == '-')
568 
    x.h |= sign_bit;
569 
  return x;
570 
}
571 
 
572 
octa
573 
fdivide (y, z)
574 
     octa y, z;
575 
{
576 
  ftype yt, zt;
577 
  int ye, ze;
578 
  char ys, zs;
579 
  octa x, xf, yf, zf;
580 
  int xe;
581 
  char xs;
582 
  yt = funpack (y, &yf, &ye, &ys);
583 
  zt = funpack (z, &zf, &ze, &zs);
584 
  xs = ys + zs - '+';
585 
  switch (4 * yt + zt)
586 
    {
587 
 
588 
    case 4 * nan + nan:
589 
      if (!(y.h & 0x80000))
590 
        exceptions |= I_BIT;
591 
    case 4 * zro + nan:
592 
    case 4 * num + nan:
593 
    case 4 * inf + nan:
594 
      if (!(z.h & 0x80000))
595 
        exceptions |= I_BIT, z.h |= 0x80000;
596 
      return z;
597 
    case 4 * nan + zro:
598 
    case 4 * nan + num:
599 
    case 4 * nan + inf:
600 
      if (!(y.h & 0x80000))
601 
        exceptions |= I_BIT, y.h |= 0x80000;
602 
      return y;
603 
 
604 
      ;
605 
    case 4 * zro + inf:
606 
    case 4 * zro + num:
607 
    case 4 * num + inf:
608 
      x = zero_octa;
609 
      break;
610 
    case 4 * num + zro:
611 
      exceptions |= Z_BIT;
612 
    case 4 * inf + num:
613 
    case 4 * inf + zro:
614 
      x = inf_octa;
615 
      break;
616 
    case 4 * zro + zro:
617 
    case 4 * inf + inf:
618 
      x = standard_NaN;
619 
      exceptions |= I_BIT;
620 
      break;
621 
    case 4 * num + num:
622 
 
623 
      xe = ye - ze + 0x3fd;
624 
      xf = odiv (yf, zero_octa, shift_left (zf, 9));
625 
      if (xf.h >= 0x800000)
626 
        {
627 
          aux.l |= xf.l & 1;
628 
          xf = shift_right (xf, 1, 1);
629 
          xe++;
630 
        }
631 
      if (aux.h || aux.l)
632 
        xf.l |= 1;
633 
      return fpack (xf, xe, xs, cur_round);
634 
 
635 
      ;
636 
    }
637 
  if (xs == '-')
638 
    x.h |= sign_bit;
639 
  return x;
640 
}
641 
 
642 
octa
643 
fplus (y, z)
644 
     octa y, z;
645 
{
646 
  ftype yt, zt;
647 
  int ye, ze;
648 
  char ys, zs;
649 
  octa x, xf, yf, zf;
650 
  int xe, d;
651 
  char xs;
652 
  yt = funpack (y, &yf, &ye, &ys);
653 
  zt = funpack (z, &zf, &ze, &zs);
654 
  switch (4 * yt + zt)
655 
    {
656 
 
657 
    case 4 * nan + nan:
658 
      if (!(y.h & 0x80000))
659 
        exceptions |= I_BIT;
660 
    case 4 * zro + nan:
661 
    case 4 * num + nan:
662 
    case 4 * inf + nan:
663 
      if (!(z.h & 0x80000))
664 
        exceptions |= I_BIT, z.h |= 0x80000;
665 
      return z;
666 
    case 4 * nan + zro:
667 
    case 4 * nan + num:
668 
    case 4 * nan + inf:
669 
      if (!(y.h & 0x80000))
670 
        exceptions |= I_BIT, y.h |= 0x80000;
671 
      return y;
672 
 
673 
      ;
674 
    case 4 * zro + num:
675 
      return fpack (zf, ze, zs, ROUND_OFF);
676 
      break;
677 
    case 4 * num + zro:
678 
      return fpack (yf, ye, ys, ROUND_OFF);
679 
      break;
680 
    case 4 * inf + inf:
681 
      if (ys != zs)
682 
        {
683 
          exceptions |= I_BIT;
684 
          x = standard_NaN;
685 
          xs = zs;
686 
          break;
687 
        }
688 
    case 4 * num + inf:
689 
    case 4 * zro + inf:
690 
      x = inf_octa;
691 
      xs = zs;
692 
      break;
693 
    case 4 * inf + num:
694 
    case 4 * inf + zro:
695 
      x = inf_octa;
696 
      xs = ys;
697 
      break;
698 
    case 4 * num + num:
699 
      if (y.h != (z.h ^ 0x80000000) || y.l != z.l)
700 
 
701 
        {
702 
          octa o, oo;
703 
          if (ye < ze
704 
              || (ye == ze && (yf.h < zf.h || (yf.h == zf.h && yf.l < zf.l))))
705 
 
706 
            {
707 
              o = yf, yf = zf, zf = o;
708 
              d = ye, ye = ze, ze = d;
709 
              d = ys, ys = zs, zs = d;
710 
            }
711 
 
712 
          ;
713 
          d = ye - ze;
714 
          xs = ys, xe = ye;
715 
          if (d)
716 
 
717 
            {
718 
              if (d <= 2)
719 
                zf = shift_right (zf, d, 1);
720 
              else if (d > 53)
721 
                zf.h = 0, zf.l = 1;
722 
              else
723 
                {
724 
                  if (ys != zs)
725 
                    d--, xe--, yf = shift_left (yf, 1);
726 
                  o = zf;
727 
                  zf = shift_right (o, d, 1);
728 
                  oo = shift_left (zf, d);
729 
                  if (oo.l != o.l || oo.h != o.h)
730 
                    zf.l |= 1;
731 
                }
732 
            }
733 
 
734 
          ;
735 
          if (ys == zs)
736 
            {
737 
              xf = oplus (yf, zf);
738 
              if (xf.h >= 0x800000)
739 
                xe++, d = xf.l & 1, xf = shift_right (xf, 1, 1), xf.l |= d;
740 
            }
741 
          else
742 
            {
743 
              xf = ominus (yf, zf);
744 
              if (xf.h >= 0x800000)
745 
                xe++, d = xf.l & 1, xf = shift_right (xf, 1, 1), xf.l |= d;
746 
              else
747 
                while (xf.h < 0x400000)
748 
                  xe--, xf = shift_left (xf, 1);
749 
            }
750 
          return fpack (xf, xe, xs, cur_round);
751 
        }
752 
 
753 
      ;
754 
    case 4 * zro + zro:
755 
      x = zero_octa;
756 
      xs = (ys == zs ? ys : cur_round == ROUND_DOWN ? '-' : '+');
757 
      break;
758 
    }
759 
  if (xs == '-')
760 
    x.h |= sign_bit;
761 
  return x;
762 
}
763 
 
764 
int
765 
fepscomp (y, z, e, s)
766 
     octa y, z, e;
767 
     int s;
768 
{
769 
  octa yf, zf, ef, o, oo;
770 
  int ye, ze, ee;
771 
  char ys, zs, es;
772 
  int yt, zt, et, d;
773 
  et = funpack (e, &ef, &ee, &es);
774 
  if (es == '-')
775 
    return 2;
776 
  switch (et)
777 
    {
778 
    case nan:
779 
      return 2;
780 
    case inf:
781 
      ee = 10000;
782 
    case num:
783 
    case zro:
784 
      break;
785 
    }
786 
  yt = funpack (y, &yf, &ye, &ys);
787 
  zt = funpack (z, &zf, &ze, &zs);
788 
  switch (4 * yt + zt)
789 
    {
790 
    case 4 * nan + nan:
791 
    case 4 * nan + inf:
792 
    case 4 * nan + num:
793 
    case 4 * nan + zro:
794 
    case 4 * inf + nan:
795 
    case 4 * num + nan:
796 
    case 4 * zro + nan:
797 
      return 2;
798 
    case 4 * inf + inf:
799 
      return (ys == zs || ee >= 1023);
800 
    case 4 * inf + num:
801 
    case 4 * inf + zro:
802 
    case 4 * num + inf:
803 
    case 4 * zro + inf:
804 
      return (s && ee >= 1022);
805 
    case 4 * zro + zro:
806 
      return 1;
807 
    case 4 * zro + num:
808 
    case 4 * num + zro:
809 
      if (!s)
810 
        return 0;
811 
    case 4 * num + num:
812 
      break;
813 
    }
814 
 
815 
 
816 
  if (ye < 0 && yt != zro)
817 
    yf = shift_left (y, 2), ye = 0;
818 
  if (ze < 0 && zt != zro)
819 
    zf = shift_left (z, 2), ze = 0;
820 
 
821 
  ;
822 
  if (ye < ze || (ye == ze && (yf.h < zf.h || (yf.h == zf.h && yf.l < zf.l))))
823 
 
824 
    {
825 
      o = yf, yf = zf, zf = o;
826 
      d = ye, ye = ze, ze = d;
827 
      d = ys, ys = zs, zs = d;
828 
    }
829 
 
830 
  ;
831 
  if (ze == zero_exponent)
832 
    ze = ye;
833 
  d = ye - ze;
834 
  if (!s)
835 
    ee -= d;
836 
  if (ee >= 1023)
837 
    return 1;
838 
 
839 
  if (d > 54)
840 
    o = zero_octa, oo = zf;
841 
  else
842 
    o = shift_right (zf, d, 1), oo = shift_left (o, d);
843 
  if (oo.h != zf.h || oo.l != zf.l)
844 
    {
845 
      if (ee < 1020)
846 
        return 0;
847 
      o = incr (o, ys == zs ? 0 : 1);
848 
    }
849 
  o = (ys == zs ? ominus (yf, o) : oplus (yf, o));
850 
 
851 
  ;
852 
  if (!o.h && !o.l)
853 
    return 1;
854 
  if (ee < 968)
855 
    return 0;
856 
  if (ee >= 1021)
857 
    ef = shift_left (ef, ee - 1021);
858 
  else
859 
    ef = shift_right (ef, 1021 - ee, 1);
860 
  return o.h < ef.h || (o.h == ef.h && o.l <= ef.l);
861 
 
862 
  ;
863 
}
864 
 
865 
void
866 
print_float (x)
867 
     octa x;
868 
{
869 
 
870 
  octa f, g;
871 
  int e;
872 
  int j, k;
873 
 
874 
  bignum ff, gg;
875 
  bignum tt;
876 
  char s[18];
877 
  char *p;
878 
 
879 
  ;
880 
  if (x.h & sign_bit)
881 
    printf ("-");
882 
 
883 
  f = shift_left (x, 1);
884 
  e = f.h >> 21;
885 
  f.h &= 0x1fffff;
886 
  if (!f.h && !f.l)
887 
 
888 
    {
889 
      if (!e)
890 
        {
891 
          printf ("0.");
892 
          return;
893 
        }
894 
      if (e == 0x7ff)
895 
        {
896 
          printf ("Inf");
897 
          return;
898 
        }
899 
      e--;
900 
      f.h = 0x3fffff, f.l = 0xffffffff;
901 
      g.h = 0x400000, g.l = 2;
902 
    }
903 
 
904 
 
905 
  else
906 
    {
907 
      g = incr (f, 1);
908 
      f = incr (f, -1);
909 
      if (!e)
910 
        e = 1;
911 
      else if (e == 0x7ff)
912 
        {
913 
          printf ("NaN");
914 
          if (g.h == 0x100000 && g.l == 1)
915 
            return;
916 
          e = 0x3ff;
917 
        }
918 
      else
919 
        f.h |= 0x200000, g.h |= 0x200000;
920 
    }
921 
 
922 
  ;
923 
 
924 
  k = (magic_offset - e) / 28;
925 
  ff.dat[k - 1] =
926 
    shift_right (f, magic_offset + 28 - e - 28 * k, 1).l & 0xfffffff;
927 
  gg.dat[k - 1] =
928 
    shift_right (g, magic_offset + 28 - e - 28 * k, 1).l & 0xfffffff;
929 
  ff.dat[k] = shift_right (f, magic_offset - e - 28 * k, 1).l & 0xfffffff;
930 
  gg.dat[k] = shift_right (g, magic_offset - e - 28 * k, 1).l & 0xfffffff;
931 
  ff.dat[k + 1] =
932 
    shift_left (f, e + 28 * k - (magic_offset - 28)).l & 0xfffffff;
933 
  gg.dat[k + 1] =
934 
    shift_left (g, e + 28 * k - (magic_offset - 28)).l & 0xfffffff;
935 
  ff.a = (ff.dat[k - 1] ? k - 1 : k);
936 
  ff.b = (ff.dat[k + 1] ? k + 1 : k);
937 
  gg.a = (gg.dat[k - 1] ? k - 1 : k);
938 
  gg.b = (gg.dat[k + 1] ? k + 1 : k);
939 
 
940 
  ;
941 
 
942 
  if (e > 0x401)
943 
 
944 
    {
945 
      int open = x.l & 1;
946 
      tt.dat[origin] = 10;
947 
      tt.a = tt.b = origin;
948 
      for (e = 1; bignum_compare (&gg, &tt) >= open; e++)
949 
        bignum_times_ten (&tt);
950 
      p = s;
951 
      while (1)
952 
        {
953 
          bignum_times_ten (&ff);
954 
          bignum_times_ten (&gg);
955 
          for (j = '0'; bignum_compare (&ff, &tt) >= 0; j++)
956 
            bignum_dec (&ff, &tt, 0x10000000), bignum_dec (&gg, &tt,
957 
                                                           0x10000000);
958 
          if (bignum_compare (&gg, &tt) >= open)
959 
            break;
960 
          *p++ = j;
961 
          if (ff.a == bignum_prec - 1 && !open)
962 
            goto done;
963 
        }
964 
      for (k = j; bignum_compare (&gg, &tt) >= open; k++)
965 
        bignum_dec (&gg, &tt, 0x10000000);
966 
      *p++ = (j + 1 + k) >> 1;
967 
    done:;
968 
    }
969 
 
970 
 
971 
  else
972 
    {
973 
      if (ff.a > origin)
974 
        ff.dat[origin] = 0;
975 
      for (e = 1, p = s; gg.a > origin || ff.dat[origin] == gg.dat[origin];)
976 
        {
977 
          if (gg.a > origin)
978 
            e--;
979 
          else
980 
            *p++ = ff.dat[origin] + '0', ff.dat[origin] = 0, gg.dat[origin] =
981 
              0;
982 
          bignum_times_ten (&ff);
983 
          bignum_times_ten (&gg);
984 
        }
985 
      *p++ = ((ff.dat[origin] + 1 + gg.dat[origin]) >> 1) + '0';
986 
    }
987 
  *p = '\0';
988 
 
989 
  ;
990 
 
991 
  if (e > 17 || e < (int) strlen (s) - 17)
992 
    printf ("%c%s%se%d", s[0], (s[1] ? "." : ""), s + 1, e - 1);
993 
  else if (e < 0)
994 
    printf (".%0*d%s", -e, 0, s);
995 
  else if (strlen (s) >= e)
996 
    printf ("%.*s.%s", e, s, s + e);
997 
  else
998 
    printf ("%s%0*d.", s, e - (int) strlen (s), 0);
999 
 
1000 
  ;
1001 
}
1002 
 
1003 
static void
1004 
bignum_times_ten (f)
1005 
     bignum *f;
1006 
{
1007 
  tetra *p, *q;
1008 
  tetra x, carry;
1009 
  for (p = &f->dat[f->b], q = &f->dat[f->a], carry = 0; p >= q; p--)
1010 
    {
1011 
      x = *p * 10 + carry;
1012 
      *p = x & 0xfffffff;
1013 
      carry = x >> 28;
1014 
    }
1015 
  *p = carry;
1016 
  if (carry)
1017 
    f->a--;
1018 
  if (f->dat[f->b] == 0 && f->b > f->a)
1019 
    f->b--;
1020 
}
1021 
 
1022 
static int
1023 
bignum_compare (f, g)
1024 
     bignum *f, *g;
1025 
{
1026 
  tetra *p, *pp, *q, *qq;
1027 
  if (f->a != g->a)
1028 
    return f->a > g->a ? -1 : 1;
1029 
  pp = &f->dat[f->b], qq = &g->dat[g->b];
1030 
  for (p = &f->dat[f->a], q = &g->dat[g->a]; p <= pp; p++, q++)
1031 
    {
1032 
      if (*p != *q)
1033 
        return *p < *q ? -1 : 1;
1034 
      if (q == qq)
1035 
        return p < pp;
1036 
    }
1037 
  return -1;
1038 
}
1039 
 
1040 
static void
1041 
bignum_dec (f, g, r)
1042 
     bignum *f, *g;
1043 
     tetra r;
1044 
{
1045 
  tetra *p, *q, *qq;
1046 
  int x, borrow;
1047 
  while (g->b > f->b)
1048 
    f->dat[++f->b] = 0;
1049 
  qq = &g->dat[g->a];
1050 
  for (p = &f->dat[g->b], q = &g->dat[g->b], borrow = 0; q >= qq; p--, q--)
1051 
    {
1052 
      x = *p - *q - borrow;
1053 
      if (x >= 0)
1054 
        borrow = 0, *p = x;
1055 
      else
1056 
        borrow = 1, *p = x + r;
1057 
    }
1058 
  for (; borrow; p--)
1059 
    if (*p)
1060 
      borrow = 0, *p = *p - 1;
1061 
    else
1062 
      *p = r - 1;
1063 
  while (f->dat[f->a] == 0)
1064 
    {
1065 
      if (f->a == f->b)
1066 
        {
1067 
          f->a = f->b = bignum_prec - 1, f->dat[bignum_prec - 1] = 0;
1068 
          return;
1069 
        }
1070 
      f->a++;
1071 
    }
1072 
  while (f->dat[f->b] == 0)
1073 
    f->b--;
1074 
}
1075 
 
1076 
int
1077 
scan_const (s)
1078 
     char *s;
1079 
{
1080 
 
1081 
  char *p, *q;
1082 
  bool NaN;
1083 
  int sign;
1084 
 
1085 
 
1086 
  char *dec_pt;
1087 
  int exp;
1088 
  int zeros;
1089 
 
1090 
 
1091 
  int k, x;
1092 
  char *pp;
1093 
  bignum ff, tt;
1094 
 
1095 
  ;
1096 
  val.h = val.l = 0;
1097 
  p = s;
1098 
  if (*p == '+' || *p == '-')
1099 
    sign = *p++;
1100 
  else
1101 
    sign = '+';
1102 
  if (strncmp (p, "NaN", 3) == 0)
1103 
    NaN = true, p += 3;
1104 
  else
1105 
    NaN = false;
1106 
  if ((isdigit (*p) && !NaN) || (*p == '.' && isdigit (*(p + 1))))
1107 
 
1108 
    {
1109 
      for (q = buf0, dec_pt = (char *) 0; isdigit (*p); p++)
1110 
        {
1111 
          val = oplus (val, shift_left (val, 2));
1112 
          val = incr (shift_left (val, 1), *p - '0');
1113 
          if (q > buf0 || *p != '0') {
1114 
            if (q < buf_max)
1115 
              *q++ = *p;
1116 
            else if (*(q - 1) == '0')
1117 
              *(q - 1) = *p;
1118 
          }
1119 
        }
1120 
      if (NaN)
1121 
        *q++ = '1';
1122 
      if (*p == '.')
1123 
 
1124 
        {
1125 
          dec_pt = q;
1126 
          p++;
1127 
          for (zeros = 0; isdigit (*p); p++)
1128 
            if (*p == '0' && q == buf0)
1129 
              zeros++;
1130 
            else if (q < buf_max)
1131 
              *q++ = *p;
1132 
            else if (*(q - 1) == '0')
1133 
              *(q - 1) = *p;
1134 
        }
1135 
 
1136 
      ;
1137 
      next_char = p;
1138 
      if (*p == 'e' && !NaN)
1139 
 
1140 
        {
1141 
          char exp_sign;
1142 
          p++;
1143 
          if (*p == '+' || *p == '-')
1144 
            exp_sign = *p++;
1145 
          else
1146 
            exp_sign = '+';
1147 
          if (isdigit (*p))
1148 
            {
1149 
              for (exp = *p++ - '0'; isdigit (*p); p++)
1150 
                if (exp < 1000)
1151 
                  exp = 10 * exp + *p - '0';
1152 
              if (!dec_pt)
1153 
                dec_pt = q, zeros = 0;
1154 
              if (exp_sign == '-')
1155 
                exp = -exp;
1156 
              next_char = p;
1157 
            }
1158 
        }
1159 
 
1160 
 
1161 
      else
1162 
        exp = 0;
1163 
      if (dec_pt)
1164 
 
1165 
        {
1166 
 
1167 
          x = buf + 341 + zeros - dec_pt - exp;
1168 
          if (q == buf0 || x >= 1413)
1169 
            {
1170 
              exp = -99999;
1171 
              goto packit;
1172 
            }
1173 
          if (x < 0)
1174 
            {
1175 
            make_it_infinite:exp = 99999;
1176 
              goto packit;
1177 
            }
1178 
          ff.a = x / 9;
1179 
          for (p = q; p < q + 8; p++)
1180 
            *p = '0';
1181 
          q = q - 1 - (q + 341 + zeros - dec_pt - exp) % 9;
1182 
          for (p = buf0 - x % 9, k = ff.a; p <= q && k <= 156; p += 9, k++)
1183 
 
1184 
            {
1185 
              for (x = *p - '0', pp = p + 1; pp < p + 9; pp++)
1186 
                x = 10 * x + *pp - '0';
1187 
              ff.dat[k] = x;
1188 
            }
1189 
 
1190 
          ;
1191 
          ff.b = k - 1;
1192 
          for (x = 0; p <= q; p += 9)
1193 
            if (strncmp (p, "000000000", 9) != 0)
1194 
              x = 1;
1195 
          ff.dat[156] += x;
1196 
 
1197 
          while (ff.dat[ff.b] == 0)
1198 
            ff.b--;
1199 
 
1200 
          ;
1201 
 
1202 
          val = zero_octa;
1203 
          if (ff.a > 36)
1204 
            {
1205 
              for (exp = 0x3fe; ff.a > 36; exp--)
1206 
                bignum_double (&ff);
1207 
              for (k = 54; k; k--)
1208 
                {
1209 
                  if (ff.dat[36])
1210 
                    {
1211 
                      if (k >= 32)
1212 
                        val.h |= 1 << (k - 32);
1213 
                      else
1214 
                        val.l |= 1 << k;
1215 
                      ff.dat[36] = 0;
1216 
                      if (ff.b == 36)
1217 
                        break;
1218 
                    }
1219 
                  bignum_double (&ff);
1220 
                }
1221 
            }
1222 
          else
1223 
            {
1224 
              tt.a = tt.b = 36, tt.dat[36] = 2;
1225 
              for (exp = 0x3fe; bignum_compare (&ff, &tt) >= 0; exp++)
1226 
                bignum_double (&tt);
1227 
              for (k = 54; k; k--)
1228 
                {
1229 
                  bignum_double (&ff);
1230 
                  if (bignum_compare (&ff, &tt) >= 0)
1231 
                    {
1232 
                      if (k >= 32)
1233 
                        val.h |= 1 << (k - 32);
1234 
                      else
1235 
                        val.l |= 1 << k;
1236 
                      bignum_dec (&ff, &tt, 1000000000);
1237 
                      if (ff.a == bignum_prec - 1)
1238 
                        break;
1239 
                    }
1240 
                }
1241 
            }
1242 
          if (k == 0)
1243 
            val.l |= 1;
1244 
 
1245 
          ;
1246 
        packit:
1247 
 
1248 
          val = fpack (val, exp, sign, ROUND_NEAR);
1249 
          if (NaN)
1250 
            {
1251 
              if ((val.h & 0x7fffffff) == 0x40000000)
1252 
                val.h |= 0x7fffffff, val.l = 0xffffffff;
1253 
              else if ((val.h & 0x7fffffff) == 0x3ff00000 && !val.l)
1254 
                val.h |= 0x40000000, val.l = 1;
1255 
              else
1256 
                val.h |= 0x40000000;
1257 
            }
1258 
 
1259 
          ;
1260 
          return 1;
1261 
        }
1262 
 
1263 
      ;
1264 
      if (sign == '-')
1265 
        val = ominus (zero_octa, val);
1266 
      return 0;
1267 
    }
1268 
 
1269 
  ;
1270 
  if (NaN)
1271 
 
1272 
    {
1273 
      next_char = p;
1274 
      val.h = 0x600000, exp = 0x3fe;
1275 
      goto packit;
1276 
    }
1277 
 
1278 
  ;
1279 
  if (strncmp (p, "Inf", 3) == 0)
1280 
 
1281 
    {
1282 
      next_char = p + 3;
1283 
      goto make_it_infinite;
1284 
    }
1285 
 
1286 
  ;
1287 
  next_char = s;
1288 
  return -1;
1289 
}
1290 
 
1291 
static void
1292 
bignum_double (f)
1293 
     bignum *f;
1294 
{
1295 
  tetra *p, *q;
1296 
  int x, carry;
1297 
  for (p = &f->dat[f->b], q = &f->dat[f->a], carry = 0; p >= q; p--)
1298 
    {
1299 
      x = *p + *p + carry;
1300 
      if (x >= 1000000000)
1301 
        carry = 1, *p = x - 1000000000;
1302 
      else
1303 
        carry = 0, *p = x;
1304 
    }
1305 
  *p = carry;
1306 
  if (carry)
1307 
    f->a--;
1308 
  if (f->dat[f->b] == 0 && f->b > f->a)
1309 
    f->b--;
1310 
}
1311 
 
1312 
int
1313 
fcomp (y, z)
1314 
     octa y, z;
1315 
{
1316 
  ftype yt, zt;
1317 
  int ye, ze;
1318 
  char ys, zs;
1319 
  octa yf, zf;
1320 
  int x;
1321 
  yt = funpack (y, &yf, &ye, &ys);
1322 
  zt = funpack (z, &zf, &ze, &zs);
1323 
  switch (4 * yt + zt)
1324 
    {
1325 
    case 4 * nan + nan:
1326 
    case 4 * zro + nan:
1327 
    case 4 * num + nan:
1328 
    case 4 * inf + nan:
1329 
    case 4 * nan + zro:
1330 
    case 4 * nan + num:
1331 
    case 4 * nan + inf:
1332 
      return 2;
1333 
    case 4 * zro + zro:
1334 
      return 0;
1335 
    case 4 * zro + num:
1336 
    case 4 * num + zro:
1337 
    case 4 * zro + inf:
1338 
    case 4 * inf + zro:
1339 
    case 4 * num + num:
1340 
    case 4 * num + inf:
1341 
    case 4 * inf + num:
1342 
    case 4 * inf + inf:
1343 
      if (ys != zs)
1344 
        x = 1;
1345 
      else if (y.h > z.h)
1346 
        x = 1;
1347 
      else if (y.h < z.h)
1348 
        x = -1;
1349 
      else if (y.l > z.l)
1350 
        x = 1;
1351 
      else if (y.l < z.l)
1352 
        x = -1;
1353 
      else
1354 
        return 0;
1355 
      break;
1356 
    }
1357 
  return (ys == '-' ? -x : x);
1358 
}
1359 
 
1360 
octa
1361 
fintegerize (z, r)
1362 
     octa z;
1363 
     int r;
1364 
{
1365 
  ftype zt;
1366 
  int ze;
1367 
  char zs;
1368 
  octa xf, zf;
1369 
  zt = funpack (z, &zf, &ze, &zs);
1370 
  if (!r)
1371 
    r = cur_round;
1372 
  switch (zt)
1373 
    {
1374 
    case nan:
1375 
      if (!(z.h & 0x80000))
1376 
        {
1377 
          exceptions |= I_BIT;
1378 
          z.h |= 0x80000;
1379 
        }
1380 
    case inf:
1381 
    case zro:
1382 
      return z;
1383 
    case num:
1384 
 
1385 
      if (ze >= 1074)
1386 
        return fpack (zf, ze, zs, ROUND_OFF);
1387 
      if (ze <= 1020)
1388 
        xf.h = 0, xf.l = 1;
1389 
      else
1390 
        {
1391 
          octa oo;
1392 
          xf = shift_right (zf, 1074 - ze, 1);
1393 
          oo = shift_left (xf, 1074 - ze);
1394 
          if (oo.l != zf.l || oo.h != zf.h)
1395 
            xf.l |= 1;
1396 
 
1397 
        }
1398 
      switch (r)
1399 
        {
1400 
        case ROUND_DOWN:
1401 
          if (zs == '-')
1402 
            xf = incr (xf, 3);
1403 
          break;
1404 
        case ROUND_UP:
1405 
          if (zs != '-')
1406 
            xf = incr (xf, 3);
1407 
        case ROUND_OFF:
1408 
          break;
1409 
        case ROUND_NEAR:
1410 
          xf = incr (xf, xf.l & 4 ? 2 : 1);
1411 
          break;
1412 
        }
1413 
      xf.l &= 0xfffffffc;
1414 
      if (ze >= 1022)
1415 
        return fpack (shift_left (xf, 1074 - ze), ze, zs, ROUND_OFF);
1416 
      if (xf.l)
1417 
        xf.h = 0x3ff00000, xf.l = 0;
1418 
      if (zs == '-')
1419 
        xf.h |= sign_bit;
1420 
      return xf;
1421 
 
1422 
      ;
1423 
    }
1424 
  /* never reached */
1425 
  return zero_octa;
1426 
}
1427 
 
1428 
octa
1429 
fixit (z, r)
1430 
     octa z;
1431 
     int r;
1432 
{
1433 
  ftype zt;
1434 
  int ze;
1435 
  char zs;
1436 
  octa zf, o;
1437 
  zt = funpack (z, &zf, &ze, &zs);
1438 
  if (!r)
1439 
    r = cur_round;
1440 
  switch (zt)
1441 
    {
1442 
    case nan:
1443 
    case inf:
1444 
      exceptions |= I_BIT;
1445 
      return z;
1446 
    case zro:
1447 
      return zero_octa;
1448 
    case num:
1449 
      if (funpack (fintegerize (z, r), &zf, &ze, &zs) == zro)
1450 
        return zero_octa;
1451 
      if (ze <= 1076)
1452 
        o = shift_right (zf, 1076 - ze, 1);
1453 
      else
1454 
        {
1455 
          if (ze > 1085 || (ze == 1085 && (zf.h > 0x400000 ||
1456 
                                           (zf.h == 0x400000
1457 
                                            && (zf.l || zs != '-')))))
1458 
            exceptions |= W_BIT;
1459 
          if (ze >= 1140)
1460 
            return zero_octa;
1461 
          o = shift_left (zf, ze - 1076);
1462 
        }
1463 
      return (zs == '-' ? ominus (zero_octa, o) : o);
1464 
    }
1465 
  /* never reached */
1466 
  return zero_octa;
1467 
}
1468 
 
1469 
octa
1470 
floatit (z, r, u, p)
1471 
     octa z;
1472 
     int r;
1473 
     int u;
1474 
     int p;
1475 
{
1476 
  int e;
1477 
  char s;
1478 
  int t;
1479 
  exceptions = 0;
1480 
  if (!z.h && !z.l)
1481 
    return zero_octa;
1482 
  if (!r)
1483 
    r = cur_round;
1484 
  if (!u && (z.h & sign_bit))
1485 
    s = '-', z = ominus (zero_octa, z);
1486 
  else
1487 
    s = '+';
1488 
  e = 1076;
1489 
  while (z.h < 0x400000)
1490 
    e--, z = shift_left (z, 1);
1491 
  while (z.h >= 0x800000)
1492 
    {
1493 
      e++;
1494 
      t = z.l & 1;
1495 
      z = shift_right (z, 1, 1);
1496 
      z.l |= t;
1497 
    }
1498 
  if (p)
1499 
 
1500 
    {
1501 
      int ex;
1502 
      tetra t;
1503 
      t = sfpack (z, e, s, r);
1504 
      ex = exceptions;
1505 
      sfunpack (t, &z, &e, &s);
1506 
      exceptions = ex;
1507 
    }
1508 
 
1509 
  ;
1510 
  return fpack (z, e, s, r);
1511 
}
1512 
 
1513 
octa
1514 
froot (z, r)
1515 
     octa z;
1516 
     int r;
1517 
{
1518 
  ftype zt;
1519 
  int ze;
1520 
  char zs;
1521 
  octa x, xf, rf, zf;
1522 
  int xe, k;
1523 
  if (!r)
1524 
    r = cur_round;
1525 
  zt = funpack (z, &zf, &ze, &zs);
1526 
  if (zs == '-' && zt != zro)
1527 
    exceptions |= I_BIT, x = standard_NaN;
1528 
  else
1529 
    switch (zt)
1530 
      {
1531 
      case nan:
1532 
        if (!(z.h & 0x80000))
1533 
          exceptions |= I_BIT, z.h |= 0x80000;
1534 
        return z;
1535 
      case inf:
1536 
      case zro:
1537 
        x = z;
1538 
        break;
1539 
      case num:
1540 
 
1541 
        xf.h = 0, xf.l = 2;
1542 
        xe = (ze + 0x3fe) >> 1;
1543 
        if (ze & 1)
1544 
          zf = shift_left (zf, 1);
1545 
        rf.h = 0, rf.l = (zf.h >> 22) - 1;
1546 
        for (k = 53; k; k--)
1547 
          {
1548 
            rf = shift_left (rf, 2);
1549 
            xf = shift_left (xf, 1);
1550 
            if (k >= 43)
1551 
              rf = incr (rf, (zf.h >> (2 * (k - 43))) & 3);
1552 
            else if (k >= 27)
1553 
              rf = incr (rf, (zf.l >> (2 * (k - 27))) & 3);
1554 
            if ((rf.l > xf.l && rf.h >= xf.h) || rf.h > xf.h)
1555 
              {
1556 
                xf.l++;
1557 
                rf = ominus (rf, xf);
1558 
                xf.l++;
1559 
              }
1560 
          }
1561 
        if (rf.h || rf.l)
1562 
          xf.l++;
1563 
        return fpack (xf, xe, '+', r);
1564 
 
1565 
        ;
1566 
      }
1567 
  if (zs == '-')
1568 
    x.h |= sign_bit;
1569 
  return x;
1570 
}
1571 
 
1572 
 
1573 
octa fremstep(octa y, octa z, int delta) {
1574 
  ftype yt, zt;
1575 
  int ye, ze;
1576 
  char xs, ys, zs;
1577 
  octa x, xf, yf, zf;
1578 
  int xe, thresh, odd;
1579 
 
1580 
  yt = funpack (y, &yf, &ye, &ys);
1581 
  zt = funpack (z, &zf, &ze, &zs);
1582 
  switch (4 * yt + zt)
1583 
    {
1584 
 
1585 
    case 4 * nan + nan:
1586 
      if (!(y.h & 0x80000))
1587 
        exceptions |= I_BIT;
1588 
    case 4 * zro + nan:
1589 
    case 4 * num + nan:
1590 
    case 4 * inf + nan:
1591 
      if (!(z.h & 0x80000))
1592 
        exceptions |= I_BIT, z.h |= 0x80000;
1593 
      return z;
1594 
    case 4 * nan + zro:
1595 
    case 4 * nan + num:
1596 
    case 4 * nan + inf:
1597 
      if (!(y.h & 0x80000))
1598 
        exceptions |= I_BIT, y.h |= 0x80000;
1599 
      return y;
1600 
 
1601 
    case 4 * zro + zro:
1602 
    case 4 * num + zro:
1603 
    case 4 * inf + zro:
1604 
    case 4 * inf + num:
1605 
    case 4 * inf + inf:
1606 
      x = standard_NaN;
1607 
      exceptions |= I_BIT;
1608 
      break;
1609 
    case 4 * zro + num:
1610 
    case 4 * zro + inf:
1611 
    case 4 * num + inf:
1612 
      return y;
1613 
    case 4 * num + num:
1614 
 
1615 
      odd = 0;
1616 
      thresh = ye - delta;
1617 
      if (thresh < ze)
1618 
        thresh = ze;
1619 
      while (ye >= thresh)
1620 
 
1621 
        {
1622 
          if (yf.h == zf.h && yf.l == zf.l)
1623 
            goto zero_out;
1624 
          if (yf.h < zf.h || (yf.h == zf.h && yf.l < zf.l))
1625 
            {
1626 
              if (ye == ze)
1627 
                goto try_complement;
1628 
              ye--, yf = shift_left (yf, 1);
1629 
            }
1630 
          yf = ominus (yf, zf);
1631 
          if (ye == ze)
1632 
            odd = 1;
1633 
          while (yf.h < 0x400000)
1634 
            ye--, yf = shift_left (yf, 1);
1635 
        }
1636 
 
1637 
      if (ye >= ze)
1638 
        {
1639 
          exceptions |= E_BIT;
1640 
          return fpack (yf, ye, ys, ROUND_OFF);
1641 
        }
1642 
      if (ye < ze - 1)
1643 
        return fpack (yf, ye, ys, ROUND_OFF);
1644 
      yf = shift_right (yf, 1, 1);
1645 
    try_complement:xf = ominus (zf, yf), xe = ze, xs = '+' + '-' - ys;
1646 
      if (xf.h > yf.h
1647 
          || (xf.h == yf.h && (xf.l > yf.l || (xf.l == yf.l && !odd))))
1648 
        xf = yf, xs = ys;
1649 
      while (xf.h < 0x400000)
1650 
        xe--, xf = shift_left (xf, 1);
1651 
      return fpack (xf, xe, xs, ROUND_OFF);
1652 
 
1653 
    zero_out:x = zero_octa;
1654 
    }
1655 
  if (ys == '-')
1656 
    x.h |= sign_bit;
1657 
  return x;
1658 
}
1659 
 
1660 
 
1661 
/**************************************************************/
1662 
 
1663 
 
1664 
static octa freg[32];
1665 
 
1666 
 
1667 
void addFloat(int rd, int rs1, int rs2) {
1668 
  freg[rd] = fplus(freg[rs1], freg[rs2]);
1669 
}
1670 
 
1671 
 
1672 
void subFloat(int rd, int rs1, int rs2) {
1673 
  octa aux;
1674 
 
1675 
  aux = freg[rs2];
1676 
  if (fcomp(aux, zero_octa) != 2) {
1677 
    aux.h ^= sign_bit;
1678 
  }
1679 
  freg[rd] = fplus(freg[rs1], aux);
1680 
}
1681 
 
1682 
 
1683 
void mulFloat(int rd, int rs1, int rs2) {
1684 
  freg[rd] = fmult(freg[rs1], freg[rs2]);
1685 
}
1686 
 
1687 
 
1688 
void divFloat(int rd, int rs1, int rs2) {
1689 
  freg[rd] = fdivide(freg[rs1], freg[rs2]);
1690 
}
1691 
 
1692 
 
1693 
/**************************************************************/
1694 
 
1695 
 
1696 
int main(void) {
1697 
  return 0;
1698 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.