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

Subversion Repositories eco32

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

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

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

powered by: WebSVN 2.1.0

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