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

Subversion Repositories eco32

[/] [eco32/] [trunk/] [fp/] [implementation/] [ultimate/] [fp.c] - Blame information for rev 15

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

powered by: WebSVN 2.1.0

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