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

Subversion Repositories eco32

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

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
/*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.