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

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib/] [newlib/] [libc/] [stdlib/] [mprec.c] - Blame information for rev 39

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

Line No. Rev Author Line
1 39 lampret
/****************************************************************
2
 *
3
 * The author of this software is David M. Gay.
4
 *
5
 * Copyright (c) 1991 by AT&T.
6
 *
7
 * Permission to use, copy, modify, and distribute this software for any
8
 * purpose without fee is hereby granted, provided that this entire notice
9
 * is included in all copies of any software which is or includes a copy
10
 * or modification of this software and in all copies of the supporting
11
 * documentation for such software.
12
 *
13
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17
 *
18
 ***************************************************************/
19
 
20
/* Please send bug reports to
21
        David M. Gay
22
        AT&T Bell Laboratories, Room 2C-463
23
        600 Mountain Avenue
24
        Murray Hill, NJ 07974-2070
25
        U.S.A.
26
        dmg@research.att.com or research!dmg
27
 */
28
 
29
/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
30
 *
31
 * This strtod returns a nearest machine number to the input decimal
32
 * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
33
 * broken by the IEEE round-even rule.  Otherwise ties are broken by
34
 * biased rounding (add half and chop).
35
 *
36
 * Inspired loosely by William D. Clinger's paper "How to Read Floating
37
 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
38
 *
39
 * Modifications:
40
 *
41
 *      1. We only require IEEE, IBM, or VAX double-precision
42
 *              arithmetic (not IEEE double-extended).
43
 *      2. We get by with floating-point arithmetic in a case that
44
 *              Clinger missed -- when we're computing d * 10^n
45
 *              for a small integer d and the integer n is not too
46
 *              much larger than 22 (the maximum integer k for which
47
 *              we can represent 10^k exactly), we may be able to
48
 *              compute (d*10^k) * 10^(e-k) with just one roundoff.
49
 *      3. Rather than a bit-at-a-time adjustment of the binary
50
 *              result in the hard case, we use floating-point
51
 *              arithmetic to determine the adjustment to within
52
 *              one bit; only in really hard cases do we need to
53
 *              compute a second residual.
54
 *      4. Because of 3., we don't need a large table of powers of 10
55
 *              for ten-to-e (just some small tables, e.g. of 10^k
56
 *              for 0 <= k <= 22).
57
 */
58
 
59
/*
60
 * #define IEEE_8087 for IEEE-arithmetic machines where the least
61
 *      significant byte has the lowest address.
62
 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
63
 *      significant byte has the lowest address.
64
 * #define Sudden_Underflow for IEEE-format machines without gradual
65
 *      underflow (i.e., that flush to zero on underflow).
66
 * #define IBM for IBM mainframe-style floating-point arithmetic.
67
 * #define VAX for VAX-style floating-point arithmetic.
68
 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
69
 * #define No_leftright to omit left-right logic in fast floating-point
70
 *      computation of dtoa.
71
 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
72
 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
73
 *      that use extended-precision instructions to compute rounded
74
 *      products and quotients) with IBM.
75
 * #define ROUND_BIASED for IEEE-format with biased rounding.
76
 * #define Inaccurate_Divide for IEEE-format with correctly rounded
77
 *      products but inaccurate quotients, e.g., for Intel i860.
78
 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
79
 *      integer arithmetic.  Whether this speeds things up or slows things
80
 *      down depends on the machine and the number being converted.
81
 */
82
 
83
#include <_ansi.h>
84
#include <stdlib.h>
85
#include <string.h>
86
#include <reent.h>
87
#include "mprec.h"
88
 
89
/* reent.c knows this value */
90
#define _Kmax 15
91
 
92
_Bigint *
93
_DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k)
94
{
95
  int x;
96
  _Bigint *rv ;
97
 
98
  if (ptr->_freelist == NULL)
99
    {
100
      /* Allocate a list of pointers to the mprec objects */
101
      ptr->_freelist = (struct _Bigint **) _calloc_r (ptr,
102
                                                      sizeof (struct _Bigint *),
103
                                                      _Kmax + 1);
104
      if (ptr->_freelist == NULL)
105
        {
106
          return NULL;
107
        }
108
    }
109
 
110
  if (rv = ptr->_freelist[k])
111
    {
112
      ptr->_freelist[k] = rv->_next;
113
    }
114
  else
115
    {
116
      x = 1 << k;
117
      /* Allocate an mprec Bigint and stick in in the freelist */
118
      rv = (_Bigint *) _calloc_r (ptr,
119
                                  1,
120
                                  sizeof (_Bigint) +
121
                                  (x-1) * sizeof(rv->_x));
122
      if (rv == NULL) return NULL;
123
      rv->_k = k;
124
      rv->_maxwds = x;
125
    }
126
  rv->_sign = rv->_wds = 0;
127
  return rv;
128
}
129
 
130
void
131
_DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v)
132
{
133
  if (v)
134
    {
135
      v->_next = ptr->_freelist[v->_k];
136
      ptr->_freelist[v->_k] = v;
137
    }
138
}
139
 
140
_Bigint *
141
_DEFUN (multadd, (ptr, b, m, a),
142
        struct _reent *ptr _AND
143
        _Bigint * b _AND
144
        int m _AND
145
        int a)
146
{
147
  int i, wds;
148
  unsigned long *x, y;
149
#ifdef Pack_32
150
  unsigned long xi, z;
151
#endif
152
  _Bigint *b1;
153
 
154
  wds = b->_wds;
155
  x = b->_x;
156
  i = 0;
157
  do
158
    {
159
#ifdef Pack_32
160
      xi = *x;
161
      y = (xi & 0xffff) * m + a;
162
      z = (xi >> 16) * m + (y >> 16);
163
      a = (int) (z >> 16);
164
      *x++ = (z << 16) + (y & 0xffff);
165
#else
166
      y = *x * m + a;
167
      a = (int) (y >> 16);
168
      *x++ = y & 0xffff;
169
#endif
170
    }
171
  while (++i < wds);
172
  if (a)
173
    {
174
      if (wds >= b->_maxwds)
175
        {
176
          b1 = Balloc (ptr, b->_k + 1);
177
          Bcopy (b1, b);
178
          Bfree (ptr, b);
179
          b = b1;
180
        }
181
      b->_x[wds++] = a;
182
      b->_wds = wds;
183
    }
184
  return b;
185
}
186
 
187
_Bigint *
188
_DEFUN (s2b, (ptr, s, nd0, nd, y9),
189
        struct _reent * ptr _AND
190
        _CONST char *s _AND
191
        int nd0 _AND
192
        int nd _AND
193
        unsigned long y9)
194
{
195
  _Bigint *b;
196
  int i, k;
197
  long x, y;
198
 
199
  x = (nd + 8) / 9;
200
  for (k = 0, y = 1; x > y; y <<= 1, k++);
201
#ifdef Pack_32
202
  b = Balloc (ptr, k);
203
  b->_x[0] = y9;
204
  b->_wds = 1;
205
#else
206
  b = Balloc (ptr, k + 1);
207
  b->_x[0] = y9 & 0xffff;
208
  b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
209
#endif
210
 
211
  i = 9;
212
  if (9 < nd0)
213
    {
214
      s += 9;
215
      do
216
        b = multadd (ptr, b, 10, *s++ - '0');
217
      while (++i < nd0);
218
      s++;
219
    }
220
  else
221
    s += 10;
222
  for (; i < nd; i++)
223
    b = multadd (ptr, b, 10, *s++ - '0');
224
  return b;
225
}
226
 
227
int
228
_DEFUN (hi0bits,
229
        (x), register unsigned long x)
230
{
231
  register int k = 0;
232
 
233
  if (!(x & 0xffff0000))
234
    {
235
      k = 16;
236
      x <<= 16;
237
    }
238
  if (!(x & 0xff000000))
239
    {
240
      k += 8;
241
      x <<= 8;
242
    }
243
  if (!(x & 0xf0000000))
244
    {
245
      k += 4;
246
      x <<= 4;
247
    }
248
  if (!(x & 0xc0000000))
249
    {
250
      k += 2;
251
      x <<= 2;
252
    }
253
  if (!(x & 0x80000000))
254
    {
255
      k++;
256
      if (!(x & 0x40000000))
257
        return 32;
258
    }
259
  return k;
260
}
261
 
262
int
263
_DEFUN (lo0bits, (y), unsigned long *y)
264
{
265
  register int k;
266
  register unsigned long x = *y;
267
 
268
  if (x & 7)
269
    {
270
      if (x & 1)
271
        return 0;
272
      if (x & 2)
273
        {
274
          *y = x >> 1;
275
          return 1;
276
        }
277
      *y = x >> 2;
278
      return 2;
279
    }
280
  k = 0;
281
  if (!(x & 0xffff))
282
    {
283
      k = 16;
284
      x >>= 16;
285
    }
286
  if (!(x & 0xff))
287
    {
288
      k += 8;
289
      x >>= 8;
290
    }
291
  if (!(x & 0xf))
292
    {
293
      k += 4;
294
      x >>= 4;
295
    }
296
  if (!(x & 0x3))
297
    {
298
      k += 2;
299
      x >>= 2;
300
    }
301
  if (!(x & 1))
302
    {
303
      k++;
304
      x >>= 1;
305
      if (!x & 1)
306
        return 32;
307
    }
308
  *y = x;
309
  return k;
310
}
311
 
312
_Bigint *
313
_DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i)
314
{
315
  _Bigint *b;
316
 
317
  b = Balloc (ptr, 1);
318
  b->_x[0] = i;
319
  b->_wds = 1;
320
  return b;
321
}
322
 
323
_Bigint *
324
_DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b)
325
{
326
  _Bigint *c;
327
  int k, wa, wb, wc;
328
  unsigned long carry, y, z;
329
  unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
330
#ifdef Pack_32
331
  unsigned long z2;
332
#endif
333
 
334
  if (a->_wds < b->_wds)
335
    {
336
      c = a;
337
      a = b;
338
      b = c;
339
    }
340
  k = a->_k;
341
  wa = a->_wds;
342
  wb = b->_wds;
343
  wc = wa + wb;
344
  if (wc > a->_maxwds)
345
    k++;
346
  c = Balloc (ptr, k);
347
  for (x = c->_x, xa = x + wc; x < xa; x++)
348
    *x = 0;
349
  xa = a->_x;
350
  xae = xa + wa;
351
  xb = b->_x;
352
  xbe = xb + wb;
353
  xc0 = c->_x;
354
#ifdef Pack_32
355
  for (; xb < xbe; xb++, xc0++)
356
    {
357
      if (y = *xb & 0xffff)
358
        {
359
          x = xa;
360
          xc = xc0;
361
          carry = 0;
362
          do
363
            {
364
              z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
365
              carry = z >> 16;
366
              z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
367
              carry = z2 >> 16;
368
              Storeinc (xc, z2, z);
369
            }
370
          while (x < xae);
371
          *xc = carry;
372
        }
373
      if (y = *xb >> 16)
374
        {
375
          x = xa;
376
          xc = xc0;
377
          carry = 0;
378
          z2 = *xc;
379
          do
380
            {
381
              z = (*x & 0xffff) * y + (*xc >> 16) + carry;
382
              carry = z >> 16;
383
              Storeinc (xc, z, z2);
384
              z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
385
              carry = z2 >> 16;
386
            }
387
          while (x < xae);
388
          *xc = z2;
389
        }
390
    }
391
#else
392
  for (; xb < xbe; xc0++)
393
    {
394
      if (y = *xb++)
395
        {
396
          x = xa;
397
          xc = xc0;
398
          carry = 0;
399
          do
400
            {
401
              z = *x++ * y + *xc + carry;
402
              carry = z >> 16;
403
              *xc++ = z & 0xffff;
404
            }
405
          while (x < xae);
406
          *xc = carry;
407
        }
408
    }
409
#endif
410
  for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
411
  c->_wds = wc;
412
  return c;
413
}
414
 
415
_Bigint *
416
_DEFUN (pow5mult,
417
        (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
418
{
419
  _Bigint *b1, *p5, *p51;
420
  int i;
421
  static _CONST int p05[3] = {5, 25, 125};
422
 
423
  if (i = k & 3)
424
    b = multadd (ptr, b, p05[i - 1], 0);
425
 
426
  if (!(k >>= 2))
427
    return b;
428
  if (!(p5 = ptr->_p5s))
429
    {
430
      /* first time */
431
      p5 = ptr->_p5s = i2b (ptr, 625);
432
      p5->_next = 0;
433
    }
434
  for (;;)
435
    {
436
      if (k & 1)
437
        {
438
          b1 = mult (ptr, b, p5);
439
          Bfree (ptr, b);
440
          b = b1;
441
        }
442
      if (!(k >>= 1))
443
        break;
444
      if (!(p51 = p5->_next))
445
        {
446
          p51 = p5->_next = mult (ptr, p5, p5);
447
          p51->_next = 0;
448
        }
449
      p5 = p51;
450
    }
451
  return b;
452
}
453
 
454
_Bigint *
455
_DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
456
{
457
  int i, k1, n, n1;
458
  _Bigint *b1;
459
  unsigned long *x, *x1, *xe, z;
460
 
461
#ifdef Pack_32
462
  n = k >> 5;
463
#else
464
  n = k >> 4;
465
#endif
466
  k1 = b->_k;
467
  n1 = n + b->_wds + 1;
468
  for (i = b->_maxwds; n1 > i; i <<= 1)
469
    k1++;
470
  b1 = Balloc (ptr, k1);
471
  x1 = b1->_x;
472
  for (i = 0; i < n; i++)
473
    *x1++ = 0;
474
  x = b->_x;
475
  xe = x + b->_wds;
476
#ifdef Pack_32
477
  if (k &= 0x1f)
478
    {
479
      k1 = 32 - k;
480
      z = 0;
481
      do
482
        {
483
          *x1++ = *x << k | z;
484
          z = *x++ >> k1;
485
        }
486
      while (x < xe);
487
      if (*x1 = z)
488
        ++n1;
489
    }
490
#else
491
  if (k &= 0xf)
492
    {
493
      k1 = 16 - k;
494
      z = 0;
495
      do
496
        {
497
          *x1++ = *x << k & 0xffff | z;
498
          z = *x++ >> k1;
499
        }
500
      while (x < xe);
501
      if (*x1 = z)
502
        ++n1;
503
    }
504
#endif
505
  else
506
    do
507
      *x1++ = *x++;
508
    while (x < xe);
509
  b1->_wds = n1 - 1;
510
  Bfree (ptr, b);
511
  return b1;
512
}
513
 
514
int
515
_DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b)
516
{
517
  unsigned long *xa, *xa0, *xb, *xb0;
518
  int i, j;
519
 
520
  i = a->_wds;
521
  j = b->_wds;
522
#ifdef DEBUG
523
  if (i > 1 && !a->_x[i - 1])
524
    Bug ("cmp called with a->_x[a->_wds-1] == 0");
525
  if (j > 1 && !b->_x[j - 1])
526
    Bug ("cmp called with b->_x[b->_wds-1] == 0");
527
#endif
528
  if (i -= j)
529
    return i;
530
  xa0 = a->_x;
531
  xa = xa0 + j;
532
  xb0 = b->_x;
533
  xb = xb0 + j;
534
  for (;;)
535
    {
536
      if (*--xa != *--xb)
537
        return *xa < *xb ? -1 : 1;
538
      if (xa <= xa0)
539
        break;
540
    }
541
  return 0;
542
}
543
 
544
_Bigint *
545
_DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND
546
        _Bigint * a _AND _Bigint * b)
547
{
548
  _Bigint *c;
549
  int i, wa, wb;
550
  long borrow, y;               /* We need signed shifts here. */
551
  unsigned long *xa, *xae, *xb, *xbe, *xc;
552
#ifdef Pack_32
553
  long z;
554
#endif
555
 
556
  i = cmp (a, b);
557
  if (!i)
558
    {
559
      c = Balloc (ptr, 0);
560
      c->_wds = 1;
561
      c->_x[0] = 0;
562
      return c;
563
    }
564
  if (i < 0)
565
    {
566
      c = a;
567
      a = b;
568
      b = c;
569
      i = 1;
570
    }
571
  else
572
    i = 0;
573
  c = Balloc (ptr, a->_k);
574
  c->_sign = i;
575
  wa = a->_wds;
576
  xa = a->_x;
577
  xae = xa + wa;
578
  wb = b->_wds;
579
  xb = b->_x;
580
  xbe = xb + wb;
581
  xc = c->_x;
582
  borrow = 0;
583
#ifdef Pack_32
584
  do
585
    {
586
      y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
587
      borrow = y >> 16;
588
      Sign_Extend (borrow, y);
589
      z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
590
      borrow = z >> 16;
591
      Sign_Extend (borrow, z);
592
      Storeinc (xc, z, y);
593
    }
594
  while (xb < xbe);
595
  while (xa < xae)
596
    {
597
      y = (*xa & 0xffff) + borrow;
598
      borrow = y >> 16;
599
      Sign_Extend (borrow, y);
600
      z = (*xa++ >> 16) + borrow;
601
      borrow = z >> 16;
602
      Sign_Extend (borrow, z);
603
      Storeinc (xc, z, y);
604
    }
605
#else
606
  do
607
    {
608
      y = *xa++ - *xb++ + borrow;
609
      borrow = y >> 16;
610
      Sign_Extend (borrow, y);
611
      *xc++ = y & 0xffff;
612
    }
613
  while (xb < xbe);
614
  while (xa < xae)
615
    {
616
      y = *xa++ + borrow;
617
      borrow = y >> 16;
618
      Sign_Extend (borrow, y);
619
      *xc++ = y & 0xffff;
620
    }
621
#endif
622
  while (!*--xc)
623
    wa--;
624
  c->_wds = wa;
625
  return c;
626
}
627
 
628
double
629
_DEFUN (ulp, (_x), double _x)
630
{
631
  union double_union x, a;
632
  register long L;
633
 
634
  x.d = _x;
635
 
636
  L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
637
#ifndef Sudden_Underflow
638
  if (L > 0)
639
    {
640
#endif
641
#ifdef IBM
642
      L |= Exp_msk1 >> 4;
643
#endif
644
      word0 (a) = L;
645
#ifndef _DOUBLE_IS_32BITS
646
      word1 (a) = 0;
647
#endif
648
 
649
#ifndef Sudden_Underflow
650
    }
651
  else
652
    {
653
      L = -L >> Exp_shift;
654
      if (L < Exp_shift)
655
        {
656
          word0 (a) = 0x80000 >> L;
657
#ifndef _DOUBLE_IS_32BITS
658
          word1 (a) = 0;
659
#endif
660
        }
661
      else
662
        {
663
          word0 (a) = 0;
664
          L -= Exp_shift;
665
#ifndef _DOUBLE_IS_32BITS
666
          word1 (a) = L >= 31 ? 1 : 1 << 31 - L;
667
#endif
668
        }
669
    }
670
#endif
671
  return a.d;
672
}
673
 
674
double
675
_DEFUN (b2d, (a, e),
676
        _Bigint * a _AND int *e)
677
{
678
  unsigned long *xa, *xa0, w, y, z;
679
  int k;
680
  union double_union d;
681
#ifdef VAX
682
  unsigned long d0, d1;
683
#else
684
#define d0 word0(d)
685
#define d1 word1(d)
686
#endif
687
 
688
  xa0 = a->_x;
689
  xa = xa0 + a->_wds;
690
  y = *--xa;
691
#ifdef DEBUG
692
  if (!y)
693
    Bug ("zero y in b2d");
694
#endif
695
  k = hi0bits (y);
696
  *e = 32 - k;
697
#ifdef Pack_32
698
  if (k < Ebits)
699
    {
700
      d0 = Exp_1 | y >> Ebits - k;
701
      w = xa > xa0 ? *--xa : 0;
702
#ifndef _DOUBLE_IS_32BITS
703
      d1 = y << (32 - Ebits) + k | w >> Ebits - k;
704
#endif
705
      goto ret_d;
706
    }
707
  z = xa > xa0 ? *--xa : 0;
708
  if (k -= Ebits)
709
    {
710
      d0 = Exp_1 | y << k | z >> 32 - k;
711
      y = xa > xa0 ? *--xa : 0;
712
#ifndef _DOUBLE_IS_32BITS
713
      d1 = z << k | y >> 32 - k;
714
#endif
715
    }
716
  else
717
    {
718
      d0 = Exp_1 | y;
719
#ifndef _DOUBLE_IS_32BITS
720
      d1 = z;
721
#endif
722
    }
723
#else
724
  if (k < Ebits + 16)
725
    {
726
      z = xa > xa0 ? *--xa : 0;
727
      d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
728
      w = xa > xa0 ? *--xa : 0;
729
      y = xa > xa0 ? *--xa : 0;
730
      d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
731
      goto ret_d;
732
    }
733
  z = xa > xa0 ? *--xa : 0;
734
  w = xa > xa0 ? *--xa : 0;
735
  k -= Ebits + 16;
736
  d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
737
  y = xa > xa0 ? *--xa : 0;
738
  d1 = w << k + 16 | y << k;
739
#endif
740
ret_d:
741
#ifdef VAX
742
  word0 (d) = d0 >> 16 | d0 << 16;
743
  word1 (d) = d1 >> 16 | d1 << 16;
744
#else
745
#undef d0
746
#undef d1
747
#endif
748
  return d.d;
749
}
750
 
751
_Bigint *
752
_DEFUN (d2b,
753
        (ptr, _d, e, bits),
754
        struct _reent * ptr _AND
755
        double _d _AND
756
        int *e _AND
757
        int *bits)
758
 
759
{
760
  union double_union d;
761
  _Bigint *b;
762
  int de, i, k;
763
  unsigned long *x, y, z;
764
#ifdef VAX
765
  unsigned long d0, d1;
766
  d.d = _d;
767
  d0 = word0 (d) >> 16 | word0 (d) << 16;
768
  d1 = word1 (d) >> 16 | word1 (d) << 16;
769
#else
770
#define d0 word0(d)
771
#define d1 word1(d)
772
#endif
773
 
774
#ifdef Pack_32
775
  b = Balloc (ptr, 1);
776
#else
777
  b = Balloc (ptr, 2);
778
#endif
779
  x = b->_x;
780
 
781
  z = d0 & Frac_mask;
782
  d0 &= 0x7fffffff;             /* clear sign bit, which we ignore */
783
#ifdef Sudden_Underflow
784
  de = (int) (d0 >> Exp_shift);
785
#ifndef IBM
786
  z |= Exp_msk11;
787
#endif
788
#else
789
  if (de = (int) (d0 >> Exp_shift))
790
    z |= Exp_msk1;
791
#endif
792
#ifdef Pack_32
793
#ifndef _DOUBLE_IS_32BITS
794
  if (y = d1)
795
    {
796
      if (k = lo0bits (&y))
797
        {
798
          x[0] = y | z << 32 - k;
799
          z >>= k;
800
        }
801
      else
802
        x[0] = y;
803
      i = b->_wds = (x[1] = z) ? 2 : 1;
804
    }
805
  else
806
#endif
807
    {
808
#ifdef DEBUG
809
      if (!z)
810
        Bug ("Zero passed to d2b");
811
#endif
812
      k = lo0bits (&z);
813
      x[0] = z;
814
      i = b->_wds = 1;
815
#ifndef _DOUBLE_IS_32BITS
816
      k += 32;
817
#endif
818
    }
819
#else
820
  if (y = d1)
821
    {
822
      if (k = lo0bits (&y))
823
        if (k >= 16)
824
          {
825
            x[0] = y | z << 32 - k & 0xffff;
826
            x[1] = z >> k - 16 & 0xffff;
827
            x[2] = z >> k;
828
            i = 2;
829
          }
830
        else
831
          {
832
            x[0] = y & 0xffff;
833
            x[1] = y >> 16 | z << 16 - k & 0xffff;
834
            x[2] = z >> k & 0xffff;
835
            x[3] = z >> k + 16;
836
            i = 3;
837
          }
838
      else
839
        {
840
          x[0] = y & 0xffff;
841
          x[1] = y >> 16;
842
          x[2] = z & 0xffff;
843
          x[3] = z >> 16;
844
          i = 3;
845
        }
846
    }
847
  else
848
    {
849
#ifdef DEBUG
850
      if (!z)
851
        Bug ("Zero passed to d2b");
852
#endif
853
      k = lo0bits (&z);
854
      if (k >= 16)
855
        {
856
          x[0] = z;
857
          i = 0;
858
        }
859
      else
860
        {
861
          x[0] = z & 0xffff;
862
          x[1] = z >> 16;
863
          i = 1;
864
        }
865
      k += 32;
866
    }
867
  while (!x[i])
868
    --i;
869
  b->_wds = i + 1;
870
#endif
871
#ifndef Sudden_Underflow
872
  if (de)
873
    {
874
#endif
875
#ifdef IBM
876
      *e = (de - Bias - (P - 1) << 2) + k;
877
      *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
878
#else
879
      *e = de - Bias - (P - 1) + k;
880
      *bits = P - k;
881
#endif
882
#ifndef Sudden_Underflow
883
    }
884
  else
885
    {
886
      *e = de - Bias - (P - 1) + 1 + k;
887
#ifdef Pack_32
888
      *bits = 32 * i - hi0bits (x[i - 1]);
889
#else
890
      *bits = (i + 2) * 16 - hi0bits (x[i]);
891
#endif
892
    }
893
#endif
894
  return b;
895
}
896
#undef d0
897
#undef d1
898
 
899
double
900
_DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b)
901
 
902
{
903
  union double_union da, db;
904
  int k, ka, kb;
905
 
906
  da.d = b2d (a, &ka);
907
  db.d = b2d (b, &kb);
908
#ifdef Pack_32
909
  k = ka - kb + 32 * (a->_wds - b->_wds);
910
#else
911
  k = ka - kb + 16 * (a->_wds - b->_wds);
912
#endif
913
#ifdef IBM
914
  if (k > 0)
915
    {
916
      word0 (da) += (k >> 2) * Exp_msk1;
917
      if (k &= 3)
918
        da.d *= 1 << k;
919
    }
920
  else
921
    {
922
      k = -k;
923
      word0 (db) += (k >> 2) * Exp_msk1;
924
      if (k &= 3)
925
        db.d *= 1 << k;
926
    }
927
#else
928
  if (k > 0)
929
    word0 (da) += k * Exp_msk1;
930
  else
931
    {
932
      k = -k;
933
      word0 (db) += k * Exp_msk1;
934
    }
935
#endif
936
  return da.d / db.d;
937
}
938
 
939
 
940
_CONST double
941
  tens[] =
942
{
943
  1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
944
  1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
945
  1e20, 1e21, 1e22, 1e23, 1e24
946
 
947
};
948
 
949
#if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
950
_CONST double bigtens[] =
951
{1e16, 1e32, 1e64, 1e128, 1e256};
952
 
953
_CONST double tinytens[] =
954
{1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
955
#else
956
_CONST double bigtens[] =
957
{1e16, 1e32};
958
 
959
_CONST double tinytens[] =
960
{1e-16, 1e-32};
961
#endif
962
 
963
 
964
double
965
_DEFUN (_mprec_log10, (dig),
966
        int dig)
967
{
968
  double v = 1.0;
969
  if (dig < 24)
970
    return tens[dig];
971
  while (dig > 0)
972
    {
973
      v *= 10;
974
      dig--;
975
    }
976
  return v;
977
}

powered by: WebSVN 2.1.0

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