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

Subversion Repositories eco32

[/] [eco32/] [trunk/] [fp/] [implementation/] [arith/] [mmix-arith.w] - Blame information for rev 301

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

Line No. Rev Author Line
1 15 hellwig
% This file is part of the MMIXware package (c) Donald E Knuth 1999
2
@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
3
 
4
\def\title{MMIX-ARITH}
5
 
6
\def\MMIX{\.{MMIX}}
7
\def\MMIXAL{\.{MMIXAL}}
8
\def\Hex#1{\hbox{$^{\scriptscriptstyle\#}$\tt#1}} % experimental hex constant
9
\def\dts{\mathinner{\ldotp\ldotp}}
10
\def\<#1>{\hbox{$\langle\,$#1$\,\rangle$}}\let\is=\longrightarrow
11
\def\ff{\\{ff\kern-.05em}}
12
@s ff TeX
13
@s bool normal @q unreserve a C++ keyword @>
14
@s xor normal @q unreserve a C++ keyword @>
15
 
16
@* Introduction. The subroutines below are used to simulate 64-bit \MMIX\
17
arithmetic on an old-fashioned 32-bit computer---like the one the author
18
had when he wrote \MMIXAL\ and the first \MMIX\ simulators in 1998 and 1999.
19
All operations are fabricated from 32-bit arithmetic, including
20
a full implementation of the IEEE floating point standard,
21
assuming only that the \CEE/ compiler has a 32-bit unsigned integer type.
22
 
23
Some day 64-bit machines will be commonplace and the awkward manipulations of
24
the present program will look quite archaic. Interested readers who have such
25
computers will be able to convert the code to a pure 64-bit form without
26
difficulty, thereby obtaining much faster and simpler routines. Meanwhile,
27
however, we can simulate the future and hope for continued progress.
28
 
29
This program module has a simple structure, intended to make it
30
suitable for loading with \MMIX\ simulators and assemblers.
31
 
32
@c
33
#include 
34
#include 
35
#include 
36
@@;
37
typedef enum{@+false,true@+} bool;
38
@@;
39
@@;
40
@@;
41
@
42
 
43
@ Subroutines of this program are declared first with a prototype,
44
as in {\mc ANSI C}, then with an old-style \CEE/ function definition.
45
Here are some preprocessor commands that make this work correctly with both
46
new-style and old-style compilers.
47
@^prototypes for functions@>
48
 
49
@=
50
#ifdef __STDC__
51
#define ARGS(list) list
52
#else
53
#define ARGS(list) ()
54
#endif
55
 
56
@ The definition of type \&{tetra} should be changed, if necessary, so that
57
it represents an unsigned 32-bit integer.
58
@^system dependencies@>
59
 
60
@=
61
typedef unsigned int tetra;
62
 /* for systems conforming to the LP-64 data model */
63
typedef struct { tetra h,l;} octa; /* two tetrabytes make one octabyte */
64
 
65
@ @d sign_bit ((unsigned)0x80000000)
66
 
67
@=
68
octa zero_octa; /* |zero_octa.h=zero_octa.l=0| */
69
octa neg_one={-1,-1}; /* |neg_one.h=neg_one.l=-1| */
70
octa inf_octa={0x7ff00000,0}; /* floating point $+\infty$ */
71
octa standard_NaN={0x7ff80000,0}; /* floating point NaN(.5) */
72
octa aux; /* auxiliary output of a subroutine */
73
bool overflow; /* set by certain subroutines for signed arithmetic */
74
 
75
@ It's easy to add and subtract octabytes, if we aren't terribly
76
worried about speed.
77
 
78
@=
79
octa oplus @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
80
octa oplus(y,z) /* compute $y+z$ */
81
  octa y,z;
82
{@+ octa x;
83
  x.h=y.h+z.h;@+
84
  x.l=y.l+z.l;
85
  if (x.l
86
  return x;
87
}
88
@#
89
octa ominus @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
90
octa ominus(y,z) /* compute $y-z$ */
91
  octa y,z;
92
{@+ octa x;
93
  x.h=y.h-z.h;@+
94
  x.l=y.l-z.l;
95
  if (x.l>y.l) x.h--;
96
  return x;
97
}
98
 
99
@ In the following subroutine, |delta| is a signed quantity that is
100
assumed to fit in a signed tetrabyte.
101
 
102
@=
103
octa incr @,@,@[ARGS((octa,int))@];@+@t}\6{@>
104
octa incr(y,delta) /* compute $y+\delta$ */
105
  octa y;
106
  int delta;
107
{@+ octa x;
108
  x.h=y.h;@+ x.l=y.l+delta;
109
  if (delta>=0) {
110
    if (x.l
111
  }@+else if (x.l>y.l) x.h--;
112
  return x;
113
}
114
 
115
@ Left and right shifts are only a bit more difficult.
116
 
117
@=
118
octa shift_left @,@,@[ARGS((octa,int))@];@+@t}\6{@>
119
octa shift_left(y,s) /* shift left by $s$ bits, where $0\le s\le64$ */
120
  octa y;
121
  int s;
122
{
123
  while (s>=32) y.h=y.l,y.l=0,s-=32;
124
  if (s) {@+register tetra yhl=y.h<>(32-s);
125
    y.h=yhl+ylh;@+ y.l<<=s;
126
  }
127
  return y;
128
}
129
@#
130
octa shift_right @,@,@[ARGS((octa,int,int))@];@+@t}\6{@>
131
octa shift_right(y,s,u) /* shift right, arithmetically if $u=0$ */
132
  octa y;
133
  int s,u;
134
{
135
  while (s>=32) y.l=y.h, y.h=(u?0: -(y.h>>31)), s-=32;
136
  if (s) {@+register tetra yhl=y.h<<(32-s),ylh=y.l>>s;
137
    y.h=(u? 0:(-(y.h>>31))<<(32-s))+(y.h>>s);@+ y.l=yhl+ylh;
138
  }
139
  return y;
140
}
141
 
142
@* Multiplication. We need to multiply two unsigned 64-bit integers, obtaining
143
an unsigned 128-bit product. It is easy to do this on a 32-bit machine
144
by using Algorithm 4.3.1M of {\sl Seminumerical Algorithms}, with $b=2^{16}$.
145
@^multiprecision multiplication@>
146
 
147
The following subroutine returns the lower half of the product, and
148
puts the upper half into a global octabyte called |aux|.
149
 
150
@=
151
octa omult @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
152
octa omult(y,z)
153
  octa y,z;
154
{
155
  register int i,j,k;
156
  tetra u[4],v[4],w[8];
157
  register tetra t;
158
  octa acc;
159
  @;
160
  for (j=0;j<4;j++) w[j]=0;
161
  for (j=0;j<4;j++)
162
    if (!v[j]) w[j+4]=0;
163
    else {
164
      for (i=k=0;i<4;i++) {
165
        t=u[i]*v[j]+w[i+j]+k;
166
        w[i+j]=t&0xffff, k=t>>16;
167
      }
168
      w[j+4]=k;
169
    }
170
  @;
171
  return acc;
172
}
173
 
174
@ @=
175
extern octa aux; /* secondary output of subroutines with multiple outputs */
176
extern bool overflow;
177
 
178
@ @=
179
u[3]=y.h>>16, u[2]=y.h&0xffff, u[1]= y.l>>16, u[0]=y.l&0xffff;
180
v[3]=z.h>>16, v[2]=z.h&0xffff, v[1]= z.l>>16, v[0]=z.l&0xffff;
181
 
182
@ @=
183
aux.h=(w[7]<<16)+w[6], aux.l=(w[5]<<16)+w[4];
184
acc.h=(w[3]<<16)+w[2], acc.l=(w[1]<<16)+w[0];
185
 
186
@ Signed multiplication has the same lower half product as unsigned
187
multiplication. The signed upper half product is obtained with at most two
188
further subtractions, after which the result has overflowed if and only if
189
the upper half is unequal to 64 copies of the sign bit in the lower half.
190
 
191
@=
192
octa signed_omult @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
193
octa signed_omult(y,z)
194
  octa y,z;
195
{
196
  octa acc;
197
  acc=omult(y,z);
198
  if (y.h&sign_bit) aux=ominus(aux,z);
199
  if (z.h&sign_bit) aux=ominus(aux,y);
200
  overflow=(aux.h!=aux.l || (aux.h^(aux.h>>1)^(acc.h&sign_bit)));
201
  return acc;
202
}
203
 
204
@* Division. Long division of an unsigned 128-bit integer by an unsigned
205
64-bit integer is, of course, one of the most challenging routines
206
needed for \MMIX\ arithmetic. The following program, based on
207
Algorithm 4.3.1D of {\sl Seminumerical Algorithms}, computes
208
octabytes $q$ and $r$ such that $(2^{64}x+y)=qz+r$ and $0\le r
209
given octabytes $x$, $y$, and~$z$, assuming that $x
210
(If $x\ge z$, it simply sets $q=x$ and $r=y$.)
211
The quotient~$q$ is returned by the subroutine;
212
the remainder~$r$ is stored in |aux|.
213
@^multiprecision division@>
214
 
215
@=
216
octa odiv @,@,@[ARGS((octa,octa,octa))@];@+@t}\6{@>
217
octa odiv(x,y,z)
218
  octa x,y,z;
219
{
220
  register int i,j,k,n,d;
221
  tetra u[8],v[4],q[4],mask,qhat,rhat,vh,vmh;
222
  register tetra t;
223
  octa acc;
224
  @;
225
  @;
226
  @;
227
  @;
228
  for (j=3;j>=0;j--) @;
229
  @;
230
  @;
231
  return acc;
232
}
233
 
234
@ @=
235
if (x.h>z.h || (x.h==z.h && x.l>=z.l)) {
236
  aux=y;@+ return x;
237
}
238
 
239
@ @=
240
u[7]=x.h>>16, u[6]=x.h&0xffff, u[5]=x.l>>16, u[4]=x.l&0xffff;
241
u[3]=y.h>>16, u[2]=y.h&0xffff, u[1]=y.l>>16, u[0]=y.l&0xffff;
242
v[3]=z.h>>16, v[2]=z.h&0xffff, v[1]=z.l>>16, v[0]=z.l&0xffff;
243
 
244
@ @=
245
for (n=4;v[n-1]==0;n--);
246
 
247
@ We shift |u| and |v| left by |d| places, where |d| is chosen to
248
make $2^{15}\le v_{n-1}<2^{16}$.
249
 
250
@=
251
vh=v[n-1];
252
for (d=0;vh<0x8000;d++,vh<<=1);
253
for (j=k=0; j
254
  t=(u[j]<
255
  u[j]=t&0xffff, k=t>>16;
256
}
257
for (j=k=0; j
258
  t=(v[j]<
259
  v[j]=t&0xffff, k=t>>16;
260
}
261
vh=v[n-1];
262
vmh=(n>1? v[n-2]: 0);
263
 
264
@ @=
265
mask=(1<
266
for (j=3; j>=n; j--) u[j]=0;
267
for (k=0;j>=0;j--) {
268
  t=(k<<16)+u[j];
269
  u[j]=t>>d, k=t&mask;
270
}
271
 
272
@ @=
273
acc.h=(q[3]<<16)+q[2], acc.l=(q[1]<<16)+q[0];
274
aux.h=(u[3]<<16)+u[2], aux.l=(u[1]<<16)+u[0];
275
 
276
@ @=
277
{
278
  @;
279
  @;
280
  @;
281
  q[j]=qhat;
282
}
283
 
284
@ @=
285
t=(u[j+n]<<16)+u[j+n-1];
286
qhat=t/vh, rhat=t-vh*qhat;
287
if (n>1) while (qhat==0x10000 || qhat*vmh>(rhat<<16)+u[j+n-2]) {
288
  qhat--, rhat+=vh;
289
  if (rhat>=0x10000) break;
290
}
291
 
292
@ After this step, |u[j+n]| will either equal |k| or |k-1|. The
293
true value of~|u| would be obtained by subtracting~|k| from |u[j+n]|;
294
but we don't have to fuss over |u[j+n]|, because it won't be examined later.
295
 
296
@=
297
for (i=k=0; i
298
  t=u[i+j]+0xffff0000-k-qhat*v[i];
299
  u[i+j]=t&0xffff, k=0xffff-(t>>16);
300
}
301
 
302
@ The correction here occurs only rarely, but it can be necessary---for
303
example, when dividing the number \Hex{7fff800100000000} by \Hex{800080020005}.
304
 
305
@=
306
if (u[j+n]!=k) {
307
  qhat--;
308
  for (i=k=0; i
309
    t=u[i+j]+v[i]+k;
310
    u[i+j]=t&0xffff, k=t>>16;
311
  }
312
}
313
 
314
@ Signed division can be reduced to unsigned division in a tedious
315
but straightforward manner. We assume that the divisor isn't zero.
316
 
317
@=
318
octa signed_odiv @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
319
octa signed_odiv(y,z)
320
  octa y,z;
321
{
322
  octa yy,zz,q;
323
  register int sy,sz;
324
  if (y.h&sign_bit) sy=2, yy=ominus(zero_octa,y);
325
  else sy=0, yy=y;
326
  if (z.h&sign_bit) sz=1, zz=ominus(zero_octa,z);
327
  else sz=0, zz=z;
328
  q=odiv(zero_octa,yy,zz);
329
  overflow=false;
330
  switch (sy+sz) {
331
 case 2+1: aux=ominus(zero_octa,aux);
332
  if (q.h==sign_bit) overflow=true;
333
 case 0+0: return q;
334
 case 2+0:@+ if (aux.h || aux.l) aux=ominus(zz,aux);
335
  goto negate_q;
336
 case 0+1:@+ if (aux.h || aux.l) aux=ominus(aux,zz);
337
 negate_q:@+ if (aux.h || aux.l) return ominus(neg_one,q);
338
  else return ominus(zero_octa,q);
339
  }
340
}
341
 
342
@* Bit fiddling. The bitwise operators of \MMIX\ are fairly easy to
343
implement directly, but three of them occur often enough to deserve
344
packaging as subroutines.
345
 
346
@=
347
octa oand @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
348
octa oand(y,z) /* compute $y\land z$ */
349
  octa y,z;
350
{@+ octa x;
351
  x.h=y.h&z.h;@+ x.l=y.l&z.l;
352
  return x;
353
}
354
@#
355
octa oandn @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
356
octa oandn(y,z) /* compute $y\land\bar z$ */
357
  octa y,z;
358
{@+ octa x;
359
  x.h=y.h&~z.h;@+ x.l=y.l&~z.l;
360
  return x;
361
}
362
@#
363
octa oxor @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
364
octa oxor(y,z) /* compute $y\oplus z$ */
365
  octa y,z;
366
{@+ octa x;
367
  x.h=y.h^z.h;@+ x.l=y.l^z.l;
368
  return x;
369
}
370
 
371
@ Here's a fun way to count the number of bits in a tetrabyte.
372
[This classical trick is called the ``Gillies--Miller method
373
for sideways addition'' in {\sl The Preparation of Programs
374
for an Electronic Digital Computer\/} by Wilkes, Wheeler, and
375
Gill, second edition (Reading, Mass.:\ Addison--Wesley, 1957),
376
191--193. Some of the tricks used here were suggested by
377
Balbir Singh, Peter Rossmanith, and Stefan Schwoon.]
378
@^Gillies, Donald Bruce@>
379
@^Miller, Jeffrey Charles Percy@>
380
@^Wilkes, Maurice Vincent@>
381
@^Wheeler, David John@>
382
@^Gill, Stanley@>
383
@^Singh, Balbir@>
384
@^Rossmanith, Peter@>
385
@^Schwoon, Stefan@>
386
 
387
@=
388
int count_bits @,@,@[ARGS((tetra))@];@+@t}\6{@>
389
int count_bits(x)
390
  tetra x;
391
{
392
  register int xx=x;
393
  xx=xx-((xx>>1)&0x55555555);
394
  xx=(xx&0x33333333)+((xx>>2)&0x33333333);
395
  xx=(xx+(xx>>4))&0x0f0f0f0f;
396
  xx=xx+(xx>>8);
397
  return (xx+(xx>>16)) & 0xff;
398
}
399
 
400
@ To compute the nonnegative byte differences of two given tetrabytes,
401
we can carry out the following 20-step branchless computation:
402
 
403
@=
404
tetra byte_diff @,@,@[ARGS((tetra,tetra))@];@+@t}\6{@>
405
tetra byte_diff(y,z)
406
  tetra y,z;
407
{
408
  register tetra d=(y&0x00ff00ff)+0x01000100-(z&0x00ff00ff);
409
  register tetra m=d&0x01000100;
410
  register tetra x=d&(m-(m>>8));
411
  d=((y>>8)&0x00ff00ff)+0x01000100-((z>>8)&0x00ff00ff);
412
  m=d&0x01000100;
413
  return x+((d&(m-(m>>8)))<<8);
414
}
415
 
416
@ To compute the nonnegative wyde differences of two tetrabytes,
417
another trick leads to a 15-step branchless computation.
418
(Research problem: Can |count_bits|, |byte_diff|, or |wyde_diff| be done
419
with fewer operations?)
420
 
421
@=
422
tetra wyde_diff @,@,@[ARGS((tetra,tetra))@];@+@t}\6{@>
423
tetra wyde_diff(y,z)
424
  tetra y,z;
425
{
426
  register tetra a=((y>>16)-(z>>16))&0x10000;
427
  register tetra b=((y&0xffff)-(z&0xffff))&0x10000;
428
  return y-(z^((y^z)&(b-a-(b>>16))));
429
}
430
 
431
@ The last bitwise subroutine we need is the most interesting:
432
It implements \MMIX's \.{MOR} and \.{MXOR} operations.
433
 
434
@=
435
octa bool_mult @,@,@[ARGS((octa,octa,bool))@];@+@t}\6{@>
436
octa bool_mult(y,z,xor)
437
  octa y,z; /* the operands */
438
  bool xor; /* do we do xor instead of or? */
439
{
440
  octa o,x;
441
  register tetra a,b,c;
442
  register int k;
443
  for (k=0,o=y,x=zero_octa;o.h||o.l;k++,o=shift_right(o,8,1))
444
    if (o.l&0xff) {
445
      a=((z.h>>k)&0x01010101)*0xff;
446
      b=((z.l>>k)&0x01010101)*0xff;
447
      c=(o.l&0xff)*0x01010101;
448
      if (xor) x.h^=a&c, x.l^=b&c;
449
      else x.h|=a&c, x.l|=b&c;
450
    }
451
  return x;
452
}
453
 
454
@* Floating point packing and unpacking. Standard IEEE floating binary
455
numbers pack a sign, exponent, and fraction into a tetrabyte
456
or octabyte. In this section we consider basic subroutines that
457
convert between IEEE format and the separate unpacked components.
458
 
459
@d ROUND_OFF 1
460
@d ROUND_UP 2
461
@d ROUND_DOWN 3
462
@d ROUND_NEAR 4
463
 
464
@=
465
int cur_round; /* the current rounding mode */
466
 
467
@ The |fpack| routine takes an octabyte $f$, a raw exponent~$e$,
468
and a sign~|s|, and packs them
469
into the floating binary number that corresponds to
470
$\pm2^{e-1076}f$, using a given rounding mode.
471
The value of $f$ should satisfy $2^{54}\le f\le 2^{55}$.
472
 
473
Thus, for example, the floating binary number $+1.0=\Hex{3ff0000000000000}$
474
is obtained when $f=2^{54}$, $e=\Hex{3fe}$, and |s='+'|.
475
The raw exponent~$e$ is usually one less than
476
the final exponent value; the leading bit of~$f$ is essentially added
477
to the exponent. (This trick works nicely for subnormal numbers, when
478
$e<0$, or in cases where the value of $f$ is rounded upwards to $2^{55}$.)
479
 
480
Exceptional events are noted by oring appropriate bits into
481
the global variable |exceptions|. Special considerations apply to
482
underflow, which is not fully specified by Section 7.4 of the IEEE standard:
483
Implementations of the standard are free to choose between two definitions
484
of ``tininess'' and two definitions of ``accuracy loss.''
485
\MMIX\ determines tininess {\it after\/} rounding, hence a result with
486
$e<0$ is not necessarily tiny; \MMIX\ treats accuracy loss as equivalent
487
to inexactness. Thus, a result underflows if and only if
488
it is tiny and either (i)~it is inexact or (ii)~the underflow trap is enabled.
489
The |fpack| routine sets |U_BIT| in |exceptions| if and only if the result is
490
tiny, |X_BIT| if and only if the result is inexact.
491
@^underflow@>
492
 
493
@d X_BIT (1<<8) /* floating inexact */
494
@d Z_BIT (1<<9) /* floating division by zero */
495
@d U_BIT (1<<10) /* floating underflow */
496
@d O_BIT (1<<11) /* floating overflow */
497
@d I_BIT (1<<12) /* floating invalid operation */
498
@d W_BIT (1<<13) /* float-to-fix overflow */
499
@d V_BIT (1<<14) /* integer overflow */
500
@d D_BIT (1<<15) /* integer divide check */
501
@d E_BIT (1<<18) /* external (dynamic) trap bit */
502
 
503
@=
504
octa fpack @,@,@[ARGS((octa,int,char,int))@];@+@t}\6{@>
505
octa fpack(f,e,s,r)
506
  octa f; /* the normalized fraction part */
507
  int e; /* the raw exponent */
508
  char s; /* the sign */
509
  int r; /* the rounding mode */
510
{
511
  octa o;
512
  if (e>0x7fd) e=0x7ff, o=zero_octa;
513
  else {
514
    if (e<0) {
515
      if (e<-54) o.h=0, o.l=1;
516
      else {@+octa oo;
517
        o=shift_right(f,-e,1);
518
        oo=shift_left(o,-e);
519
        if (oo.l!=f.l || oo.h!=f.h) o.l |= 1; /* sticky bit */
520
@^sticky bit@>
521
      }
522
      e=0;
523
    }@+else o=f;
524
  }
525
  @;
526
}
527
 
528
@ @=
529
int exceptions; /* bits possibly destined for rA */
530
 
531
@ Everything falls together so nicely here, it's almost too good to be true!
532
 
533
@=
534
if (o.l&3) exceptions |= X_BIT;
535
switch (r) {
536
 case ROUND_DOWN:@+ if (s=='-') o=incr(o,3);@+break;
537
 case ROUND_UP:@+ if (s!='-') o=incr(o,3);
538
 case ROUND_OFF: break;
539
 case ROUND_NEAR: o=incr(o, o.l&4? 2: 1);@+break;
540
}
541
o = shift_right(o,2,1);
542
o.h += e<<20;
543
if (o.h>=0x7ff00000) exceptions |= O_BIT+X_BIT; /* overflow */
544
else if (o.h<0x100000) exceptions |= U_BIT; /* tininess */
545
if (s=='-') o.h |= sign_bit;
546
return o;
547
 
548
@ Similarly, |sfpack| packs a short float, from inputs
549
having the same conventions as |fpack|.
550
 
551
@=
552
tetra sfpack @,@,@[ARGS((octa,int,char,int))@];@+@t}\6{@>
553
tetra sfpack(f,e,s,r)
554
  octa f; /* the fraction part */
555
  int e; /* the raw exponent */
556
  char s; /* the sign */
557
  int r; /* the rounding mode */
558
{
559
  register tetra o;
560
  if (e>0x47d) e=0x47f, o=0;
561
  else {
562
    o=shift_left(f,3).h;
563
    if (f.l&0x1fffffff) o|=1;
564
    if (e<0x380) {
565
      if (e<0x380-25) o=1;
566
      else {@+register tetra o0,oo;
567
        o0 = o;
568
        o = o>>(0x380-e);
569
        oo = o<<(0x380-e);
570
        if (oo!=o0) o |= 1; /* sticky bit */
571
@^sticky bit@>
572
      }
573
      e=0x380;
574
    }
575
  }
576
  @;
577
}
578
 
579
@ @=
580
if (o&3) exceptions |= X_BIT;
581
switch (r) {
582
 case ROUND_DOWN:@+ if (s=='-') o+=3;@+break;
583
 case ROUND_UP:@+ if (s!='-') o+=3;
584
 case ROUND_OFF: break;
585
 case ROUND_NEAR: o+=(o&4? 2: 1);@+break;
586
}
587
o = o>>2;
588
o += (e-0x380)<<23;
589
if (o>=0x7f800000) exceptions |= O_BIT+X_BIT; /* overflow */
590
else if (o<0x100000) exceptions |= U_BIT; /* tininess */
591
if (s=='-') o |= sign_bit;
592
return o;
593
 
594
@ The |funpack| routine is, roughly speaking, the opposite of |fpack|.
595
It takes a given floating point number~$x$ and separates out its
596
fraction part~$f$, exponent~$e$, and sign~$s$. It clears |exceptions|
597
to zero. It returns the type of value found: |zro|, |num|, |inf|,
598
or |nan|. When it returns |num|,
599
it will have set $f$, $e$, and~$s$
600
to the values from which |fpack| would produce the original number~$x$
601
without exceptions.
602
 
603
@d zero_exponent (-1000) /* zero is assumed to have this exponent */
604
 
605
@=
606
typedef enum {@!zro,@!num,@!inf,@!nan}@+ftype;
607
 
608
@ @=
609
ftype funpack @,@,@[ARGS((octa,octa*,int*,char*))@];@+@t}\6{@>
610
ftype funpack(x,f,e,s)
611
  octa x; /* the given floating point value */
612
  octa *f; /* address where the fraction part should be stored */
613
  int *e; /* address where the exponent part should be stored */
614
  char *s; /* address where the sign should be stored */
615
{
616
  register int ee;
617
  exceptions=0;
618
  *s=(x.h&sign_bit? '-': '+');
619
  *f=shift_left(x,2);
620
  f->h &= 0x3fffff;
621
  ee=(x.h>>20)&0x7ff;
622
  if (ee) {
623
    *e=ee-1;
624
    f->h |= 0x400000;
625
    return (ee<0x7ff? num: f->h==0x400000 && !f->l? inf: nan);
626
  }
627
  if (!x.l && !f->h) {
628
    *e=zero_exponent;@+ return zro;
629
  }
630
  do {@+ ee--;@+ *f=shift_left(*f,1);@+} while (!(f->h&0x400000));
631
  *e=ee;@+ return num;
632
}
633
 
634
@ @=
635
ftype sfunpack @,@,@[ARGS((tetra,octa*,int*,char*))@];@+@t}\6{@>
636
ftype sfunpack(x,f,e,s)
637
  tetra x; /* the given floating point value */
638
  octa *f; /* address where the fraction part should be stored */
639
  int *e; /* address where the exponent part should be stored */
640
  char *s; /* address where the sign should be stored */
641
{
642
  register int ee;
643
  exceptions=0;
644
  *s=(x&sign_bit? '-': '+');
645
  f->h=(x>>1)&0x3fffff, f->l=x<<31;
646
  ee=(x>>23)&0xff;
647
  if (ee) {
648
    *e=ee+0x380-1;
649
    f->h |= 0x400000;
650
    return (ee<0xff? num: (x&0x7fffffff)==0x7f800000? inf: nan);
651
  }
652
  if (!(x&0x7fffffff)) {
653
    *e=zero_exponent;@+return zro;
654
  }
655
  do {@+ ee--;@+ *f=shift_left(*f,1);@+} while (!(f->h&0x400000));
656
  *e=ee+0x380;@+ return num;
657
}
658
 
659
@ Since \MMIX\ downplays 32-bit operations, it uses |sfpack| and |sfunpack|
660
only when loading and storing short floats, or when converting
661
from fixed point to floating point.
662
 
663
@=
664
octa load_sf @,@,@[ARGS((tetra))@];@+@t}\6{@>
665
octa load_sf(z)
666
  tetra z; /* 32 bits to be loaded into a 64-bit register */
667
{
668
  octa f,x;@+int e;@+char s;@+ftype t;
669
  t=sfunpack(z,&f,&e,&s);
670
  switch (t) {
671
 case zro: x=zero_octa;@+break;
672
 case num: return fpack(f,e,s,ROUND_OFF);
673
 case inf: x=inf_octa;@+break;
674
 case nan: x=shift_right(f,2,1);@+x.h|=0x7ff00000;@+break;
675
  }
676
  if (s=='-') x.h|=sign_bit;
677
  return x;
678
}
679
 
680
@ @=
681
tetra store_sf @,@,@[ARGS((octa))@];@+@t}\6{@>
682
tetra store_sf(x)
683
  octa x; /* 64 bits to be loaded into a 32-bit word */
684
{
685
  octa f;@+tetra z;@+int e;@+char s;@+ftype t;
686
  t=funpack(x,&f,&e,&s);
687
  switch (t) {
688
 case zro: z=0;@+break;
689
 case num: return sfpack(f,e,s,cur_round);
690
 case inf: z=0x7f800000;@+break;
691
 case nan:@+ if (!(f.h&0x200000)) {
692
      f.h|=0x200000;@+exceptions|=I_BIT; /* NaN was signaling */
693
    }
694
    z=0x7f800000|(f.h<<1)|(f.l>>31);@+break;
695
  }
696
  if (s=='-') z|=sign_bit;
697
  return z;
698
}
699
 
700
@* Floating multiplication and division.
701
The hardest fixed point operations were multiplication and division;
702
but these two operations are the {\it easiest\/} to implement in floating point
703
arithmetic, once their fixed point counterparts are available.
704
 
705
@=
706
octa fmult @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
707
octa fmult(y,z)
708
  octa y,z;
709
{
710
  ftype yt,zt;
711
  int ye,ze;
712
  char ys,zs;
713
  octa x,xf,yf,zf;
714
  register int xe;
715
  register char xs;
716
  yt=funpack(y,&yf,&ye,&ys);
717
  zt=funpack(z,&zf,&ze,&zs);
718
  xs=ys+zs-'+'; /* will be |'-'| when the result is negative */
719
  switch (4*yt+zt) {
720
 @t\4@>@;
721
 case 4*zro+zro: case 4*zro+num: case 4*num+zro: x=zero_octa;@+break;
722
 case 4*num+inf: case 4*inf+num: case 4*inf+inf: x=inf_octa;@+break;
723
 case 4*zro+inf: case 4*inf+zro: x=standard_NaN;
724
  exceptions|=I_BIT;@+break;
725
 case 4*num+num: @;
726
  }
727
  if (xs=='-') x.h|=sign_bit;
728
  return x;
729
}
730
 
731
@ @=
732
case 4*nan+nan:@+if (!(y.h&0x80000)) exceptions|=I_BIT; /* |y| is signaling */
733
case 4*zro+nan: case 4*num+nan: case 4*inf+nan:
734
  if (!(z.h&0x80000)) exceptions|=I_BIT, z.h|=0x80000;
735
  return z;
736
case 4*nan+zro: case 4*nan+num: case 4*nan+inf:
737
  if (!(y.h&0x80000)) exceptions|=I_BIT, y.h|=0x80000;
738
  return y;
739
 
740
@ @=
741
xe=ye+ze-0x3fd; /* the raw exponent */
742
x=omult(yf,shift_left(zf,9));
743
if (aux.h>=0x400000) xf=aux;
744
else xf=shift_left(aux,1), xe--;
745
if (x.h||x.l) xf.l|=1; /* adjust the sticky bit */
746
return fpack(xf,xe,xs,cur_round);
747
 
748
@ @=
749
octa fdivide @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
750
octa fdivide(y,z)
751
  octa y,z;
752
{
753
  ftype yt,zt;
754
  int ye,ze;
755
  char ys,zs;
756
  octa x,xf,yf,zf;
757
  register int xe;
758
  register char xs;
759
  yt=funpack(y,&yf,&ye,&ys);
760
  zt=funpack(z,&zf,&ze,&zs);
761
  xs=ys+zs-'+'; /* will be |'-'| when the result is negative */
762
  switch (4*yt+zt) {
763
 @t\4@>@;
764
 case 4*zro+inf: case 4*zro+num: case 4*num+inf: x=zero_octa;@+break;
765
 case 4*num+zro: exceptions|=Z_BIT;
766
 case 4*inf+num: case 4*inf+zro: x=inf_octa;@+break;
767
 case 4*zro+zro: case 4*inf+inf: x=standard_NaN;
768
  exceptions|=I_BIT;@+break;
769
 case 4*num+num: @;
770
  }
771
  if (xs=='-') x.h|=sign_bit;
772
  return x;
773
}
774
 
775
@ @=
776
xe=ye-ze+0x3fd; /* the raw exponent */
777
xf=odiv(yf,zero_octa,shift_left(zf,9));
778
if (xf.h>=0x800000) {
779
  aux.l|=xf.l&1;
780
  xf=shift_right(xf,1,1);
781
  xe++;
782
}
783
if (aux.h||aux.l) xf.l|=1; /* adjust the sticky bit */
784
return fpack(xf,xe,xs,cur_round);
785
 
786
@*Floating addition and subtraction. Now for the bread-and-butter
787
operation, the sum of two floating point numbers.
788
It is not terribly difficult, but many cases need to be handled carefully.
789
 
790
@=
791
octa fplus @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
792
octa fplus(y,z)
793
  octa y,z;
794
{
795
  ftype yt,zt;
796
  int ye,ze;
797
  char ys,zs;
798
  octa x,xf,yf,zf;
799
  register int xe,d;
800
  register char xs;
801
  yt=funpack(y,&yf,&ye,&ys);
802
  zt=funpack(z,&zf,&ze,&zs);
803
  switch (4*yt+zt) {
804
 @t\4@>@;
805
 case 4*zro+num: return fpack(zf,ze,zs,ROUND_OFF);@+break; /* may underflow */
806
 case 4*num+zro: return fpack(yf,ye,ys,ROUND_OFF);@+break; /* may underflow */
807
 case 4*inf+inf:@+if (ys!=zs) {
808
    exceptions|=I_BIT;@+x=standard_NaN;@+xs=zs;@+break;
809
  }
810
 case 4*num+inf: case 4*zro+inf: x=inf_octa;@+xs=zs;@+break;
811
 case 4*inf+num: case 4*inf+zro: x=inf_octa;@+xs=ys;@+break;
812
 case 4*num+num:@+ if (y.h!=(z.h^0x80000000) || y.l!=z.l)
813
   @;
814
 case 4*zro+zro: x=zero_octa;
815
  xs=(ys==zs? ys: cur_round==ROUND_DOWN? '-': '+');@+break;
816
  }
817
  if (xs=='-') x.h|=sign_bit;
818
  return x;
819
}
820
 
821
@ @=
822
{@+octa o,oo;
823
  if (ye
824
    @;
825
  d=ye-ze;
826
  xs=ys, xe=ye;
827
  if (d) @;
828
  if (ys==zs) {
829
    xf=oplus(yf,zf);
830
    if (xf.h>=0x800000) xe++, d=xf.l&1, xf=shift_right(xf,1,1), xf.l|=d;
831
  }@+else {
832
    xf=ominus(yf,zf);
833
    if (xf.h>=0x800000) xe++, d=xf.l&1, xf=shift_right(xf,1,1), xf.l|=d;
834
    else@+ while (xf.h<0x400000) xe--, xf=shift_left(xf,1);
835
  }
836
  return fpack(xf,xe,xs,cur_round);
837
}
838
 
839
@ @=
840
{
841
  o=yf, yf=zf, zf=o;
842
  d=ye, ye=ze, ze=d;
843
  d=ys, ys=zs, zs=d;
844
}
845
 
846
@ Proper rounding requires two bits to the right of the fraction delivered
847
to~|fpack|. The first is the true next bit of the result;
848
the other is a ``sticky'' bit, which is nonzero if any further bits of the
849
true result are nonzero. Sticky rounding to an integer takes
850
$x$ into the number $\lfloor x/2\rfloor+\lceil x/2\rceil$.
851
@^sticky bit@>
852
 
853
Some subtleties need to be observed here, in order to
854
prevent the sticky bit from being shifted left. If we did not
855
shift |yf| left~1 before shifting |zf| to the right, an incorrect
856
answer would be obtained in certain cases---for example, if
857
$|yf|=2^{54}$, $|zf|=2^{54}+2^{53}-1$, $d=52$.
858
 
859
@=
860
{
861
  if (d<=2) zf=shift_right(zf,d,1); /* exact result */
862
  else if (d>53) zf.h=0, zf.l=1; /* tricky but OK */
863
  else {
864
    if (ys!=zs) d--,xe--,yf=shift_left(yf,1);
865
    o=zf;
866
    zf=shift_right(o,d,1);
867
    oo=shift_left(zf,d);
868
    if (oo.l!=o.l || oo.h!=o.h) zf.l|=1;
869
  }
870
}
871
 
872
@ The comparison of floating point numbers with respect to $\epsilon$
873
shares some of the characteristics of floating point addition/subtraction.
874
In some ways it is simpler, and in other ways it is more difficult;
875
we might as well deal with it now. % anyways
876
 
877
Subroutine |fepscomp(y,z,e,s)| returns 2 if |y|, |z|, or |e| is a NaN
878
or |e| is negative. It returns 1 if |s=0| and $y\approx z\ (e)$ or if
879
|s!=0| and $y\sim z\ (e)$,
880
as defined in Section~4.2.2 of {\sl Seminumerical Algorithms\/};
881
otherwise it returns~0.
882
 
883
@=
884
int fepscomp @,@,@[ARGS((octa,octa,octa,int))@];@+@t}\6{@>
885
int fepscomp(y,z,e,s)
886
  octa y,z,e; /* the operands */
887
  int s; /* test similarity? */
888
{
889
  octa yf,zf,ef,o,oo;
890
  int ye,ze,ee;
891
  char ys,zs,es;
892
  register int yt,zt,et,d;
893
  et=funpack(e,&ef,&ee,&es);
894
  if (es=='-') return 2;
895
  switch (et) {
896
 case nan: return 2;
897
 case inf: ee=10000;
898
 case num: case zro: break;
899
  }
900
  yt=funpack(y,&yf,&ye,&ys);
901
  zt=funpack(z,&zf,&ze,&zs);
902
  switch (4*yt+zt) {
903
 case 4*nan+nan: case 4*nan+inf: case 4*nan+num: case 4*nan+zro:
904
 case 4*inf+nan: case 4*num+nan: case 4*zro+nan: return 2;
905
 case 4*inf+inf: return (ys==zs || ee>=1023);
906
 case 4*inf+num: case 4*inf+zro: case 4*num+inf: case 4*zro+inf:
907
   return (s && ee>=1022);
908
 case 4*zro+zro: return 1;
909
 case 4*zro+num: case 4*num+zro:@+ if (!s) return 0;
910
 case 4*num+num: break;
911
  }
912
  @;
913
}
914
 
915
@ The relation $y\approx z\ (\epsilon)$ reduces to
916
$y\sim z\ (\epsilon/2^d)$, if $d$~is the difference between the
917
larger and smaller exponents of $y$ and~$z$.
918
 
919
@=
920
@;
921
if (ye
922
  @;
923
if (ze==zero_exponent) ze=ye;
924
d=ye-ze;
925
if (!s) ee-=d;
926
if (ee>=1023) return 1; /* if $\epsilon\ge2$, $z\in N_\epsilon(y)$ */
927
@;
928
if (!o.h && !o.l) return 1;
929
if (ee<968) return 0; /* if $y\ne z$ and $\epsilon<2^{-54}$, $y\not\sim z$ */
930
if (ee>=1021) ef=shift_left(ef,ee-1021);
931
else ef=shift_right(ef,1021-ee,1);
932
return o.h
933
 
934
@ @=
935
if (ye<0 && yt!=zro) yf=shift_left(y,2), ye=0;
936
if (ze<0 && zt!=zro) zf=shift_left(z,2), ze=0;
937
 
938
@ At this point $y\sim z$ if and only if
939
$$|yf|+(-1)^{[ys=zs]}|zf|/2^d\le 2^{ee-1021}|ef|=2^{55}\epsilon.$$
940
We need to evaluate this relation without overstepping the bounds of
941
our simulated 64-bit registers.
942
 
943
When $d>2$, the difference of fraction parts might not fit exactly
944
in an octabyte;
945
in that case the numbers are not similar unless $\epsilon>3/8$,
946
and we replace the difference by the ceiling of the
947
true result. When $\epsilon<1/8$, our program essentially replaces
948
$2^{55}\epsilon$ by $\lfloor2^{55}\epsilon\rfloor$. These
949
truncations are not needed simultaneously. Therefore the logic
950
is justified by the facts that, if $n$ is an integer, we have
951
$x\le n$ if and only if $\lceil x\rceil\le n$;
952
$n\le x$ if and only if $n\le\lfloor x\rfloor$. (Notice that the
953
concept of ``sticky bit'' is {\it not\/} appropriate here.)
954
@^sticky bit@>
955
 
956
@=
957
if (d>54) o=zero_octa,oo=zf;
958
else o=shift_right(zf,d,1),oo=shift_left(o,d);
959
if (oo.h!=zf.h || oo.l!=zf.l) { /* truncated result, hence $d>2$ */
960
  if (ee<1020) return 0; /* difference is too large for similarity */
961
  o=incr(o,ys==zs? 0: 1); /* adjust for ceiling */
962
}
963
o=(ys==zs? ominus(yf,o): oplus(yf,o));
964
 
965
@*Floating point output conversion.
966
The |print_float| routine converts an octabyte to a floating decimal
967
representation that will be input as precisely the same value.
968
@^binary-to-decimal conversion@>
969
@^radix conversion@>
970
@^multiprecision conversion@>
971
 
972
@=
973
static void bignum_times_ten @,@,@[ARGS((bignum*))@];
974
static void bignum_dec @,@,@[ARGS((bignum*,bignum*,tetra))@];
975
static int bignum_compare @,@,@[ARGS((bignum*,bignum*))@];
976
void print_float @,@,@[ARGS((octa))@];@+@t}\6{@>
977
void print_float(x)
978
  octa x;
979
{
980
  @;
981
  if (x.h&sign_bit) printf("-");
982
  @
983
       fraction interval $[f\dts g]$ or $(f\dts g)$@>;
984
  @;
985
  @;
986
  @;
987
}
988
 
989
@ One way to visualize the problem being solved here is to consider
990
the vastly simpler case in which there are only 2-bit exponents
991
and 2-bit fractions. Then the sixteen possible 4-bit combinations
992
have the following interpretations:
993
$$\def\\{\;\dts\;}
994
\vbox{\halign{#\qquad&$#$\hfil\cr
995
0000&[0\\0.125]\cr
996
0001&(0.125\\0.375)\cr
997
0010&[0.375\\0.625]\cr
998
0011&(0.625\\0.875)\cr
999
0100&[0.875\\1.125]\cr
1000
0101&(1.125\\1.375)\cr
1001
0110&[1.375\\1.625]\cr
1002
0111&(1.625\\1.875)\cr
1003
1000&[1.875\\2.25]\cr
1004
1001&(2.25\\2.75)\cr
1005
1010&[2.75\\3.25]\cr
1006
1011&(3.25\\3.75)\cr
1007
1100&[3.75\\\infty]\cr
1008
1101&\rm NaN(0\\0.375)\cr
1009
1110&\rm NaN[0.375\\0.625]\cr
1010
1111&\rm NaN(0.625\\1)\cr}}$$
1011
Notice that the interval is closed, $[f\dts g]$, when the fraction part
1012
is even; it is open, $(f\dts g)$, when the fraction part is odd.
1013
The printed outputs for these sixteen values, if we actually were
1014
dealing with such short exponents and fractions, would be
1015
\.{0.}, \.{.2}, \.{.5}, \.{.7}, \.{1.}, \.{1.2}, \.{1.5}, \.{1.7},
1016
\.{2.}, \.{2.5}, \.{3.}, \.{3.5}, \.{Inf}, \.{NaN.2}, \.{NaN}, \.{NaN.8},
1017
respectively.
1018
 
1019
@=
1020
f=shift_left(x,1);
1021
e=f.h>>21;
1022
f.h&=0x1fffff;
1023
if (!f.h && !f.l) @@;
1024
else {
1025
  g=incr(f,1);
1026
  f=incr(f,-1);
1027
  if (!e) e=1; /* subnormal */
1028
  else if (e==0x7ff) {
1029
    printf("NaN");
1030
    if (g.h==0x100000 && g.l==1) return; /* the ``standard'' NaN */
1031
    e=0x3ff; /* extreme NaNs come out OK even without adjusting |f| or |g| */
1032
  }@+else f.h|=0x200000, g.h|=0x200000;
1033
}
1034
 
1035
@ @=
1036
octa f,g; /* lower and upper bounds on the fraction part */
1037
register int e; /* exponent part */
1038
register int j,k; /* all purpose indices */
1039
 
1040
@ The transition points between exponents correspond to powers of~2. At
1041
such points the interval extends only half as far to the left of that
1042
power of~2 as it does to the right. For example, in the 4-bit minifloat numbers
1043
considered above, case 1000 corresponds to the interval $[1.875\;\dts\;2.25]$.
1044
 
1045
@=
1046
{
1047
  if (!e) {
1048
    printf("0.");@+return;
1049
  }
1050
  if (e==0x7ff) {
1051
    printf("Inf");@+return;
1052
  }
1053
  e--;
1054
  f.h=0x3fffff, f.l=0xffffffff;
1055
  g.h=0x400000, g.l=2;
1056
}
1057
 
1058
@ We want to find the ``simplest'' value in the interval corresponding
1059
to the given number, in the sense that it has fewest significant
1060
digits when expressed in decimal notation. Thus, for example,
1061
if the floating point number can be described by a relatively
1062
short string such as `\.{.1}' or `\.{37e100}', we want to discover that
1063
representation.
1064
 
1065
The basic idea is to generate the decimal representations of the
1066
two endpoints of the interval, outputting the leading digits where
1067
both endpoints agree, then making a final decision at the first place where
1068
they disagree.
1069
 
1070
The ``simplest'' value is not always unique. For example, in the
1071
case of 4-bit minifloat numbers we could represent the bit pattern 0001 as
1072
either \.{.2} or \.{.3}, and we could represent 1001 in five equally short
1073
ways: \.{2.3} or \.{2.4} or \.{2.5} or \.{2.6} or \.{2.7}. The
1074
algorithm below tries to choose the middle possibility in such cases.
1075
 
1076
[A solution to the analogous problem for fixed-point representations,
1077
without the additional complication of round-to-even, was used by
1078
the author in the program for \TeX; see {\sl Beauty is Our Business\/}
1079
(Springer, 1990), 233--242.]
1080
@^Knuth, Donald Ervin@>
1081
 
1082
Suppose we are given two fractions $f$ and $g$, where $0\le f
1083
we want to compute the shortest decimal in the closed interval $[f\dts g]$.
1084
If $f=0$, we are done. Otherwise let $10f=d+f'$ and $10g=e+g'$, where
1085
$0\le f'<1$ and $0\le g'<1$. If $d
1086
any of the digits $d+1$, \dots,~$e$; otherwise we output the
1087
common digit $d=e$, and repeat the process on the fractions $0\le f'
1088
A similar procedure works with respect to the open interval $(f\dts g)$.
1089
 
1090
@ The program below carries out the stated algorithm by using multiprecision
1091
arithmetic on 77-place integers with 28 bits each. This choice
1092
facilitates multiplication by~10, and allows us to deal with the whole range of
1093
floating binary numbers using fixed point arithmetic. We keep track of
1094
the leading and trailing digit positions so that trivial operations on
1095
zeros are avoided.
1096
 
1097
If |f| points to a \&{bignum}, its radix-$2^{28}$ digits are
1098
|f->dat[0]| through |f->dat[76]|, from most significant to least significant.
1099
We assume that all digit positions are zero unless they lie in the
1100
subarray between indices |f->a| and |f->b|, inclusive.
1101
Furthermore, both |f->dat[f->a]| and |f->dat[f->b]| are nonzero,
1102
unless |f->a=f->b=bignum_prec-1|.
1103
 
1104
The \&{bignum} data type can be used with any radix less than
1105
$2^{32}$; we will use it later with radix~$10^9$. The |dat| array
1106
is made large enough to accommodate both applications.
1107
 
1108
@d bignum_prec 157 /* would be 77 if we cared only about |print_float| */
1109
 
1110
@=
1111
typedef struct {
1112
  int a; /* index of the most significant digit */
1113
  int b; /* index of the least significant digit; must be $\ge a$ */
1114
  tetra dat[bignum_prec]; /* the digits; undefined except between |a| and |b| */
1115
} bignum;
1116
 
1117
@ Here, for example, is how we go from $f$ to $10f$, assuming that
1118
overflow will not occur and that the radix is $2^{28}$:
1119
 
1120
@=
1121
static void bignum_times_ten(f)
1122
  bignum *f;
1123
{
1124
  register tetra *p,*q; register tetra x,carry;
1125
  for (p=&f->dat[f->b],q=&f->dat[f->a],carry=0; p>=q; p--) {
1126
    x=*p*10+carry;
1127
    *p=x&0xfffffff;
1128
    carry=x>>28;
1129
  }
1130
  *p=carry;
1131
  if (carry) f->a--;
1132
  if (f->dat[f->b]==0 && f->b>f->a) f->b--;
1133
}
1134
 
1135
@ And here is how we test whether $fg$, using any
1136
radix whatever:
1137
 
1138
@=
1139
static int bignum_compare(f,g)
1140
  bignum *f,*g;
1141
{
1142
  register tetra *p,*pp,*q,*qq;
1143
  if (f->a!=g->a) return f->a > g->a? -1: 1;
1144
  pp=&f->dat[f->b], qq=&g->dat[g->b];
1145
  for (p=&f->dat[f->a],q=&g->dat[g->a]; p<=pp; p++,q++) {
1146
    if (*p!=*q) return *p<*q? -1: 1;
1147
    if (q==qq) return p
1148
  }
1149
  return -1;
1150
}
1151
 
1152
@ The following subroutine subtracts $g$ from~$f$, assuming that
1153
$f\ge g>0$ and using a given radix.
1154
 
1155
@=
1156
static void bignum_dec(f,g,r)
1157
  bignum *f,*g;
1158
  tetra r; /* the radix */
1159
{
1160
  register tetra *p,*q,*qq;
1161
  register int x,borrow;
1162
  while (g->b>f->b) f->dat[++f->b]=0;
1163
  qq=&g->dat[g->a];
1164
  for (p=&f->dat[g->b],q=&g->dat[g->b],borrow=0;q>=qq;p--,q--) {
1165
    x=*p - *q - borrow;
1166
    if (x>=0) borrow=0, *p=x;
1167
    else borrow=1, *p=x+r;
1168
  }
1169
  for (;borrow;p--)
1170
    if (*p) borrow=0, *p=*p-1;
1171
    else *p=r-1;
1172
  while (f->dat[f->a]==0) {
1173
    if (f->a==f->b) { /* the result is zero */
1174
      f->a=f->b=bignum_prec-1, f->dat[bignum_prec-1]=0;
1175
      return;
1176
    }
1177
    f->a++;
1178
  }
1179
  while (f->dat[f->b]==0) f->b--;
1180
}
1181
 
1182
@ Armed with these subroutines, we are ready to solve the problem.
1183
The first task is to put the numbers into \&{bignum} form.
1184
If the exponent is |e|, the number destined for digit |dat[k]| will
1185
consist of the rightmost 28 bits of the given fraction after it has
1186
been shifted right $c-e-28k$ bits, for some constant~$c$.
1187
We choose $c$ so that,
1188
when $e$ has its maximum value \Hex{7ff}, the leading digit will
1189
go into position |dat[1]|, and so that when the number to be printed
1190
is exactly~1 the integer part of~$g$ will also be exactly~1.
1191
 
1192
@d magic_offset 2112 /* the constant $c$ that makes it work */
1193
@d origin 37 /* the radix point follows |dat[37]| */
1194
 
1195
@=
1196
k=(magic_offset-e)/28;
1197
ff.dat[k-1]=shift_right(f,magic_offset+28-e-28*k,1).l&0xfffffff;
1198
gg.dat[k-1]=shift_right(g,magic_offset+28-e-28*k,1).l&0xfffffff;
1199
ff.dat[k]=shift_right(f,magic_offset-e-28*k,1).l&0xfffffff;
1200
gg.dat[k]=shift_right(g,magic_offset-e-28*k,1).l&0xfffffff;
1201
ff.dat[k+1]=shift_left(f,e+28*k-(magic_offset-28)).l&0xfffffff;
1202
gg.dat[k+1]=shift_left(g,e+28*k-(magic_offset-28)).l&0xfffffff;
1203
ff.a=(ff.dat[k-1]? k-1: k);
1204
ff.b=(ff.dat[k+1]? k+1: k);
1205
gg.a=(gg.dat[k-1]? k-1: k);
1206
gg.b=(gg.dat[k+1]? k+1: k);
1207
 
1208
@ If $e$ is sufficiently small, the fractions $f$ and $g$ will be less than~1,
1209
and we can use the stated algorithm directly. Of course, if $e$ is
1210
extremely small, a lot of leading zeros need to be lopped off; in the
1211
worst case, we may have to multiply $f$ and~$g$ by~10 more than 300 times.
1212
But hey, we don't need to do that extremely often, and computers are
1213
pretty fast nowadays.
1214
 
1215
In the small-exponent case, the computation always terminates before
1216
$f$ becomes zero, because the interval endpoints are fractions with
1217
denominator $2^t$ for some $t>50$.
1218
 
1219
The invariant relations |ff.dat[ff.a]!=0| and |gg.dat[gg.a]!=0| are
1220
not maintained by the computation here, when |ff.a=origin| or |gg.a=origin|.
1221
But no harm is done, because |bignum_compare| is not used.
1222
 
1223
@=
1224
if (e>0x401) @@;
1225
else@+{ /* if |e<=0x401| we have |gg.a>=origin| and |gg.dat[origin]<=8| */
1226
  if (ff.a>origin) ff.dat[origin]=0;
1227
  for (e=1, p=s; gg.a>origin || ff.dat[origin]==gg.dat[origin]; ) {
1228
    if (gg.a>origin) e--;
1229
    else *p++=ff.dat[origin]+'0', ff.dat[origin]=0, gg.dat[origin]=0;
1230
    bignum_times_ten(&ff);
1231
    bignum_times_ten(&gg);
1232
  }
1233
  *p++=((ff.dat[origin]+1+gg.dat[origin])>>1)+'0'; /* the middle digit */
1234
}
1235
*p='\0'; /* terminate the string |s| */
1236
 
1237
@ When |e| is large, we use the stated algorithm by considering $f$ and
1238
$g$ to be fractions whose denominator is a power of~10.
1239
 
1240
An interesting case arises when the number to be converted is
1241
\Hex{44ada56a4b0835bf}, since the interval turns out to be
1242
$$ (69999999999999991611392\ \ \dts\ \ 70000000000000000000000).$$
1243
If this were a closed interval, we could simply give the answer
1244
\.{7e22}; but the number \.{7e22} actually corresponds to
1245
\Hex{44ada56a4b0835c0}
1246
because of the round-to-even rule. Therefore the correct answer is, say,
1247
\.{6.9999999999999995e22}. This example shows that we need a slightly
1248
different strategy in the case of open intervals; we cannot simply
1249
look at the first position in which the endpoints have different
1250
decimal digits. Therefore we change the invariant relation to $0\le f
1251
when open intervals are involved,
1252
and we do not terminate the process when $f=0$ or $g=1$.
1253
 
1254
@=
1255
{@+register int open=x.l&1;
1256
  tt.dat[origin]=10;
1257
  tt.a=tt.b=origin;
1258
  for (e=1;bignum_compare(&gg,&tt)>=open;e++)
1259
    bignum_times_ten(&tt);
1260
  p=s;
1261
  while (1) {
1262
    bignum_times_ten(&ff);
1263
    bignum_times_ten(&gg);
1264
    for (j='0';bignum_compare(&ff,&tt)>=0;j++)
1265
      bignum_dec(&ff,&tt,0x10000000),bignum_dec(&gg,&tt,0x10000000);
1266
    if (bignum_compare(&gg,&tt)>=open) break;
1267
    *p++=j;
1268
    if (ff.a==bignum_prec-1 && !open)
1269
      goto done; /* $f=0$ in a closed interval */
1270
  }
1271
  for (k=j;bignum_compare(&gg,&tt)>=open;k++) bignum_dec(&gg,&tt,0x10000000);
1272
  *p++=(j+1+k)>>1; /* the middle digit */
1273
 done:;
1274
}
1275
 
1276
@ The length of string~|s| will be at most 17. For if $f$ and $g$
1277
agree to 17 places, we have $g/f<1+10^{-16}$; but the
1278
ratio $g/f$ is always $\ge(1+2^{-52}+2^{-53})/(1+2^{-52}-2^{-53})
1279
>1+2\times10^{-16}$.
1280
 
1281
@=
1282
bignum ff,gg; /* fractions or numerators of fractions */
1283
bignum tt; /* power of ten (used as the denominator) */
1284
char s[18];
1285
register char *p;
1286
 
1287
@ At this point the significant digits are in string |s|, and |s[0]!='0'|.
1288
If we put a decimal point at the left of~|s|, the result should
1289
be multiplied by $10^e$.
1290
 
1291
We prefer the output `\.{300.}' to the form `\.{3e2}', and we prefer
1292
`\.{.03}' to `\.{3e-2}'. In general, the output will use an
1293
explicit exponent only if the alternative would take more than
1294
18~characters.
1295
 
1296
@=
1297
if (e>17 || e<(int)strlen(s)-17)
1298
  printf("%c%s%se%d",s[0],(s[1]? ".": ""),s+1,e-1);
1299
else if (e<0) printf(".%0*d%s",-e,0,s);
1300
else if (strlen(s)>=e) printf("%.*s.%s",e,s,s+e);
1301
else printf("%s%0*d.",s,e-(int)strlen(s),0);
1302
 
1303
@*Floating point input conversion. Going the other way, we want to
1304
be able to convert a given decimal number into its floating binary
1305
@^decimal-to-binary conversion@>
1306
@^radix conversion@>
1307
@^multiprecision conversion@>
1308
equivalent. The following syntax is supported:
1309
$$\vbox{\halign{$#$\hfil\cr
1310
\\is\.0\mid\.1\mid\.2\mid\.3\mid\.4\mid
1311
        \.5\mid\.6\mid\.7\mid\.8\mid\.9\cr
1312
\\is\\mid\\\cr
1313
\\is\\..\mid\..\\mid
1314
                      \\..\\cr
1315
\\is\\mid\.+\mid\.-\cr
1316
\\is\.e\\\cr
1317
\\is\\mid\\cr
1318
\\is\\\mid
1319
                    \\\mid\cr
1320
\hskip12em          \.{Inf}\mid\.{NaN}\mid\.{NaN.}\\cr
1321
\\is\\\cr
1322
\\is\\\cr
1323
}}$$
1324
For example, `\.{-3.}' is the floating constant \Hex{c008000000000000}\thinspace;
1325
`\.{1e3}' and `\.{1000}' are both equivalent to \Hex{408f400000000000}\thinspace;
1326
`\.{NaN}' and `\.{+NaN.5}' are both equivalent to \Hex{7ff8000000000000}.
1327
 
1328
The |scan_const| routine looks at a given string and finds the
1329
longest initial substring that matches the syntax of either \
1330
constant> or \. It puts the corresponding value
1331
into the global octabyte variable~|val|; it also puts the position of the first
1332
unscanned character in the global pointer variable |next_char|.
1333
It returns 1 if a floating constant was found, 0~if a decimal constant
1334
was found, $-1$ if nothing was found. A decimal constant that doesn't
1335
fit in an octabyte is computed modulo~$2^{64}$.
1336
@^syntax of floating point constants@>
1337
 
1338
The value of |exceptions| set by |scan_const| is not necessarily correct.
1339
 
1340
@=
1341
static void bignum_double @,@,@[ARGS((bignum*))@];
1342
int scan_const @,@,@[ARGS((char*))@];@+@t}\6{@>
1343
int scan_const(s)
1344
  char *s;
1345
{
1346
  @;
1347
  val.h=val.l=0;
1348
  p=s;
1349
  if (*p=='+' || *p=='-') sign=*p++;@+else sign='+';
1350
  if (strncmp(p,"NaN",3)==0) NaN=true, p+=3;
1351
  else NaN=false;
1352
  if ((isdigit(*p)&&!NaN) || (*p=='.' && isdigit(*(p+1))))
1353
    @;
1354
  if (NaN) @;
1355
  if (strncmp(p,"Inf",3)==0) @;
1356
 no_const_found: next_char=s;@+return -1;
1357
}
1358
 
1359
@ @=
1360
octa val; /* value returned by |scan_const| */
1361
char *next_char; /* pointer returned by |scan_const| */
1362
 
1363
@ @=
1364
register char *p,*q; /* for string manipulations */
1365
register bool NaN; /* are we processing a NaN? */
1366
int sign; /* |'+'| or |'-'| */
1367
 
1368
@ @=
1369
{
1370
  next_char=p;
1371
  val.h=0x600000, exp=0x3fe;
1372
  goto packit;
1373
}
1374
 
1375
@ @=
1376
{
1377
  next_char=p+3;
1378
  goto make_it_infinite;
1379
}
1380
 
1381
@ We saw above that a string of at most 17 digits is enough to characterize
1382
a floating point number, for purposes of output. But a much longer buffer
1383
for digits is needed when we're doing input. For example, consider the
1384
borderline quantity $(1+2^{-53})/2^{1022}$; its decimal expansion, when
1385
written out exactly, is a number with more than 750 significant digits:
1386
\.{2.2250738585...8125e-308}.
1387
If {\it any one\/} of those digits is increased, or if
1388
additional nonzero digits are added as in
1389
\.{2.2250738585...81250000001e-308},
1390
the rounded value is supposed to change from \Hex{0010000000000000}
1391
to \Hex{0010000000000001}.
1392
 
1393
We assume here that the user prefers a perfectly correct answer to
1394
a speedy almost-correct one, so we implement the most general case.
1395
 
1396
@=
1397
{
1398
  for (q=buf0,dec_pt=(char*)0;isdigit(*p);p++) {
1399
    val=oplus(val,shift_left(val,2)); /* multiply by 5 */
1400
    val=incr(shift_left(val,1),*p-'0');
1401
    if (q>buf0 || *p!='0')
1402
       if (q
1403
       else if (*(q-1)=='0') *(q-1)=*p;
1404
  }
1405
  if (NaN) *q++='1';
1406
  if (*p=='.') @;
1407
  next_char=p;
1408
  if (*p=='e' && !NaN) @@;
1409
  else exp=0;
1410
  if (dec_pt) @;
1411
  if (sign=='-') val=ominus(zero_octa,val);
1412
  return 0;
1413
}
1414
 
1415
@ @=
1416
{
1417
  dec_pt=q;
1418
  p++;
1419
  for (zeros=0;isdigit(*p);p++)
1420
    if (*p=='0' && q==buf0) zeros++;
1421
    else if (q
1422
    else if (*(q-1)=='0') *(q-1)=*p;
1423
}
1424
 
1425
@ The buffer needs room for eight digits of padding at the left, followed
1426
by up to $1022+53-307$ significant digits, followed by a ``sticky'' digit
1427
at position |buf_max-1|, and eight more digits of padding.
1428
 
1429
@d buf0 (buf+8)
1430
@d buf_max (buf+777)
1431
 
1432
@=
1433
static char buf[785]="00000000"; /* where we put significant input digits */
1434
 
1435
@ @=
1436
register char* dec_pt; /* position of decimal point in |buf| */
1437
register int exp; /* scanned exponent; later used for raw binary exponent */
1438
register int zeros; /* leading zeros removed after decimal point */
1439
 
1440
@ Here we don't advance |next_char| and force a decimal point until we
1441
know that a syntactically correct exponent exists.
1442
 
1443
The code here will convert extra-large inputs like
1444
`\.{9e+9999999999999999}' into $\infty$ and extra-small inputs into zero.
1445
Strange inputs like `\.{-00.0e9999999}' must also be accommodated.
1446
 
1447
@=
1448
{@+register char exp_sign;
1449
  p++;
1450
  if (*p=='+' || *p=='-') exp_sign=*p++;@+else exp_sign='+';
1451
  if (isdigit(*p)) {
1452
    for (exp=*p++ -'0';isdigit(*p);p++)
1453
      if (exp<1000) exp = 10*exp + *p - '0';
1454
    if (!dec_pt) dec_pt=q, zeros=0;
1455
    if (exp_sign=='-') exp=-exp;
1456
    next_char=p;
1457
  }
1458
}
1459
 
1460
@ @=
1461
{
1462
  @;
1463
  @;
1464
packit:  @;
1465
  return 1;
1466
}
1467
 
1468
@ Now we get ready to compute the binary fraction bits, by putting the
1469
scanned input digits into a multiprecision fixed-point
1470
accumulator |ff| that spans the full necessary range.
1471
After this step, the number that we want to convert to floating binary
1472
will appear in |ff.dat[ff.a]|, |ff.dat[ff.a+1]|, \dots,
1473
|ff.dat[ff.b]|.
1474
The radix-$10^9$ digit in ${\it ff}[36-k]$ is understood to be multiplied
1475
by $10^{9k}$, for $36\ge k\ge-120$.
1476
 
1477
@=
1478
x=buf+341+zeros-dec_pt-exp;
1479
if (q==buf0 || x>=1413) {
1480
 make_it_zero: exp=-99999;@+ goto packit;
1481
}
1482
if (x<0) {
1483
 make_it_infinite: exp=99999;@+ goto packit;
1484
}
1485
ff.a=x/9;
1486
for (p=q;p
1487
q=q-1-(q+341+zeros-dec_pt-exp)%9; /* compute stopping place in |buf| */
1488
for (p=buf0-x%9,k=ff.a;p<=q && k<=156; p+=9, k++)
1489
  @
1490
       into |ff.dat[k]|@>;
1491
ff.b=k-1;
1492
for (x=0;p<=q;p+=9) if (strncmp(p,"000000000",9)!=0) x=1;
1493
ff.dat[156]+=x; /* nonzero digits that fall off the right are sticky */
1494
@^sticky bit@>
1495
while (ff.dat[ff.b]==0) ff.b--;
1496
 
1497
@ @=
1498
{
1499
  for (x=*p-'0',pp=p+1;pp
1500
  ff.dat[k]=x;
1501
}
1502
 
1503
@ @=
1504
register int k,x;
1505
register char *pp;
1506
bignum ff,tt;
1507
 
1508
@ Here's a subroutine that is dual to |bignum_times_ten|. It changes $f$
1509
to~$2f$, assuming that overflow will not occur and that the radix is $10^9$.
1510
 
1511
@=
1512
static void bignum_double(f)
1513
  bignum *f;
1514
{
1515
  register tetra *p,*q; register int x,carry;
1516
  for (p=&f->dat[f->b],q=&f->dat[f->a],carry=0; p>=q; p--) {
1517
    x = *p + *p + carry;
1518
    if (x>=1000000000) carry=1, *p=x-1000000000;
1519
    else carry=0, *p=x;
1520
  }
1521
  *p=carry;
1522
  if (carry) f->a--;
1523
  if (f->dat[f->b]==0 && f->b>f->a) f->b--;
1524
}
1525
 
1526
@ @=
1527
val=zero_octa;
1528
if (ff.a>36) {
1529
  for (exp=0x3fe;ff.a>36;exp--) bignum_double(&ff);
1530
  for (k=54;k;k--) {
1531
    if (ff.dat[36]) {
1532
      if (k>=32) val.h |= 1<<(k-32);@+else val.l |= 1<
1533
      ff.dat[36]=0;
1534
      if (ff.b==36) break; /* break if |ff| now zero */
1535
    }
1536
    bignum_double(&ff);
1537
  }
1538
}@+else {
1539
  tt.a=tt.b=36, tt.dat[36]=2;
1540
  for (exp=0x3fe;bignum_compare(&ff,&tt)>=0;exp++) bignum_double(&tt);
1541
  for (k=54;k;k--) {
1542
    bignum_double(&ff);
1543
    if (bignum_compare(&ff,&tt)>=0) {
1544
      if (k>=32) val.h |= 1<<(k-32);@+else val.l |= 1<
1545
      bignum_dec(&ff,&tt,1000000000);
1546
      if (ff.a==bignum_prec-1) break; /* break if |ff| now zero */
1547
    }
1548
  }
1549
}
1550
if (k==0) val.l |= 1; /* add sticky bit if |ff| nonzero */
1551
 
1552
@ We need to be careful that the input `\.{NaN.999999999999999999999}' doesn't
1553
get rounded up; it is supposed to yield \Hex{7fffffffffffffff}.
1554
 
1555
Although the input `\.{NaN.0}' is illegal, strictly speaking, we silently
1556
convert it to \Hex{7ff0000000000001}---a number that would be
1557
output as `\.{NaN.0000000000000002}'.
1558
 
1559
@=
1560
val=fpack(val,exp,sign,ROUND_NEAR);
1561
if (NaN) {
1562
  if ((val.h&0x7fffffff)==0x40000000) val.h |= 0x7fffffff, val.l=0xffffffff;
1563
  else if ((val.h&0x7fffffff)==0x3ff00000 && !val.l) val.h|=0x40000000,val.l=1;
1564
  else val.h |= 0x40000000;
1565
}
1566
 
1567
@*Floating point remainders. In this section we implement the remainder
1568
of the floating point operations---one of which happens to be the
1569
operation of taking the remainder.
1570
 
1571
The easiest task remaining is to compare two floating point quantities.
1572
Routine |fcomp| returns $-1$~if~$yz$, and
1573
$+2$~if $y$ and~$z$ are unordered.
1574
 
1575
@=
1576
int fcomp @,@,@[ARGS((octa,octa))@];@+@t}\6{@>
1577
int fcomp(y,z)
1578
  octa y,z;
1579
{
1580
  ftype yt,zt;
1581
  int ye,ze;
1582
  char ys,zs;
1583
  octa yf,zf;
1584
  register int x;
1585
  yt=funpack(y,&yf,&ye,&ys);
1586
  zt=funpack(z,&zf,&ze,&zs);
1587
  switch (4*yt+zt) {
1588
 case 4*nan+nan: case 4*zro+nan: case 4*num+nan: case 4*inf+nan:
1589
 case 4*nan+zro: case 4*nan+num: case 4*nan+inf: return 2;
1590
 case 4*zro+zro: return 0;
1591
 case 4*zro+num: case 4*num+zro: case 4*zro+inf: case 4*inf+zro:
1592
 case 4*num+num: case 4*num+inf: case 4*inf+num: case 4*inf+inf:
1593
  if (ys!=zs) x=1;
1594
  else if (y.h>z.h) x=1;
1595
  else if (y.h
1596
  else if (y.l>z.l) x=1;
1597
  else if (y.l
1598
  else return 0;
1599
  break;
1600
  }
1601
  return (ys=='-'? -x: x);
1602
}
1603
 
1604
@ Several \MMIX\ operations act on a single floating point number and
1605
accept an arbitrary rounding mode. For example, consider the
1606
operation of rounding to the nearest floating point integer:
1607
 
1608
@=
1609
octa fintegerize @,@,@[ARGS((octa,int))@];@+@t}\6{@>
1610
octa fintegerize(z,r)
1611
  octa z; /* the operand */
1612
  int r; /* the rounding mode */
1613
{
1614
  ftype zt;
1615
  int ze;
1616
  char zs;
1617
  octa xf,zf;
1618
  zt=funpack(z,&zf,&ze,&zs);
1619
  if (!r) r=cur_round;
1620
  switch (zt) {
1621
 case nan:@+if (!(z.h&0x80000)) {@+exceptions|=I_BIT;@+z.h|=0x80000;@+}
1622
 case inf: case zro: return z;
1623
 case num: @;
1624
  }
1625
}
1626
 
1627
@ @=
1628
if (ze>=1074) return fpack(zf,ze,zs,ROUND_OFF); /* already an integer */
1629
if (ze<=1020) xf.h=0,xf.l=1;
1630
else {@+octa oo;
1631
  xf=shift_right(zf,1074-ze,1);
1632
  oo=shift_left(xf,1074-ze);
1633
  if (oo.l!=zf.l || oo.h!=zf.h) xf.l|=1; /* sticky bit */
1634
@^sticky bit@>
1635
}
1636
switch (r) {
1637
 case ROUND_DOWN:@+ if (zs=='-') xf=incr(xf,3);@+break;
1638
 case ROUND_UP:@+ if (zs!='-') xf=incr(xf,3);
1639
 case ROUND_OFF: break;
1640
 case ROUND_NEAR: xf=incr(xf, xf.l&4? 2: 1);@+break;
1641
}
1642
xf.l&=0xfffffffc;
1643
if (ze>=1022) return fpack(shift_left(xf,1074-ze),ze,zs,ROUND_OFF);
1644
if (xf.l) xf.h=0x3ff00000, xf.l=0;
1645
if (zs=='-') xf.h|=sign_bit;
1646
return xf;
1647
 
1648
@ To convert floating point to fixed point, we use |fixit|.
1649
 
1650
@=
1651
octa fixit @,@,@[ARGS((octa,int))@];@+@t}\6{@>
1652
octa fixit(z,r)
1653
  octa z; /* the operand */
1654
  int r; /* the rounding mode */
1655
{
1656
  ftype zt;
1657
  int ze;
1658
  char zs;
1659
  octa zf,o;
1660
  zt=funpack(z,&zf,&ze,&zs);
1661
  if (!r) r=cur_round;
1662
  switch (zt) {
1663
 case nan: case inf: exceptions|=I_BIT;@+return z;
1664
 case zro: return zero_octa;
1665
 case num:@+if (funpack(fintegerize(z,r),&zf,&ze,&zs)==zro) return zero_octa;
1666
   if (ze<=1076) o=shift_right(zf,1076-ze,1);
1667
   else {
1668
     if (ze>1085 || (ze==1085 && (zf.h>0x400000 || @|
1669
           (zf.h==0x400000 && (zf.l || zs!='-'))))) exceptions|=W_BIT;
1670
     if (ze>=1140) return zero_octa;
1671
     o=shift_left(zf,ze-1076);
1672
    }
1673
  return (zs=='-'? ominus(zero_octa,o): o);
1674
  }
1675
}
1676
 
1677
@ Going the other way, we can specify not only a rounding mode but whether
1678
the given fixed point octabyte is signed or unsigned, and whether the
1679
result should be rounded to short precision.
1680
 
1681
@=
1682
octa floatit @,@,@[ARGS((octa,int,int,int))@];@+@t}\6{@>
1683
octa floatit(z,r,u,p)
1684
  octa z; /* octabyte to float */
1685
  int r; /* rounding mode */
1686
  int u; /* unsigned? */
1687
  int p; /* short precision? */
1688
{
1689
  int e;@+char s;
1690
  register int t;
1691
  exceptions=0;
1692
  if (!z.h && !z.l) return zero_octa;
1693
  if (!r) r=cur_round;
1694
  if (!u && (z.h&sign_bit)) s='-', z=ominus(zero_octa,z);@+ else s='+';
1695
  e=1076;
1696
  while (z.h<0x400000) e--,z=shift_left(z,1);
1697
  while (z.h>=0x800000) {
1698
    e++;
1699
    t=z.l&1;
1700
    z=shift_right(z,1,1);
1701
    z.l|=t;
1702
  }
1703
  if (p) @;
1704
  return fpack(z,e,s,r);
1705
}
1706
 
1707
@ @=
1708
{
1709
  register int ex;@+register tetra t;
1710
  t=sfpack(z,e,s,r);
1711
  ex=exceptions;
1712
  sfunpack(t,&z,&e,&s);
1713
  exceptions=ex;
1714
}
1715
 
1716
@ The square root operation is more interesting.
1717
 
1718
@=
1719
octa froot @,@,@[ARGS((octa,int))@];@+@t}\6{@>
1720
octa froot(z,r)
1721
  octa z; /* the operand */
1722
  int r; /* the rounding mode */
1723
{
1724
  ftype zt;
1725
  int ze;
1726
  char zs;
1727
  octa x,xf,rf,zf;
1728
  register int xe,k;
1729
  if (!r) r=cur_round;
1730
  zt=funpack(z,&zf,&ze,&zs);
1731
  if (zs=='-' && zt!=zro) exceptions|=I_BIT, x=standard_NaN;
1732
  else@+switch (zt) {
1733
 case nan:@+ if (!(z.h&0x80000)) exceptions|=I_BIT, z.h|=0x80000;
1734
  return z;
1735
 case inf: case zro: x=z;@+break;
1736
 case num: @;
1737
  }
1738
  if (zs=='-') x.h|=sign_bit;
1739
  return x;
1740
}
1741
 
1742
@ The square root can be found by an adaptation of the old pencil-and-paper
1743
method. If $n=\lfloor\sqrt s\rfloor$, where $s$ is an integer,
1744
we have $s=n^2+r$ where $0\le r\le2n$;
1745
this invariant can be maintained if we replace $s$ by $4s+(0,1,2,3)$
1746
and $n$ by $2n+(0,1)$. The following code implements this idea with
1747
$2n$ in~|xf| and $r$ in~|rf|. (It could easily be made to run about
1748
twice as fast.)
1749
 
1750
@=
1751
xf.h=0, xf.l=2;
1752
xe=(ze+0x3fe)>>1;
1753
if (ze&1) zf=shift_left(zf,1);
1754
rf.h=0, rf.l=(zf.h>>22)-1;
1755
for (k=53;k;k--) {
1756
  rf=shift_left(rf,2);@+ xf=shift_left(xf,1);
1757
  if (k>=43) rf=incr(rf,(zf.h>>(2*(k-43)))&3);
1758
  else if (k>=27) rf=incr(rf,(zf.l>>(2*(k-27)))&3);
1759
  if ((rf.l>xf.l && rf.h>=xf.h) || rf.h>xf.h) {
1760
    xf.l++;@+rf=ominus(rf,xf);@+xf.l++;
1761
  }
1762
}
1763
if (rf.h || rf.l) xf.l++; /* sticky bit */
1764
return fpack(xf,xe,'+',r);
1765
 
1766
@ And finally, the genuine floating point remainder. Subroutine |fremstep|
1767
either calculates $y\,{\rm rem}\,z$ or reduces $y$ to a smaller number
1768
having the same remainder with respect to~$z$. In the latter case
1769
the |E_BIT| is set in |exceptions|. A third parameter, |delta|,
1770
gives a decrease in exponent that is acceptable for incomplete results;
1771
if |delta| is sufficiently large, say 2500, the correct result will
1772
always be obtained in one step of |fremstep|.
1773
 
1774
@=
1775
octa fremstep @,@,@[ARGS((octa,octa,int))@];@+@t}\6{@>
1776
octa fremstep(y,z,delta)
1777
  octa y,z;
1778
  int delta;
1779
{
1780
  ftype yt,zt;
1781
  int ye,ze;
1782
  char xs,ys,zs;
1783
  octa x,xf,yf,zf;
1784
  register int xe,thresh,odd;
1785
  yt=funpack(y,&yf,&ye,&ys);
1786
  zt=funpack(z,&zf,&ze,&zs);
1787
  switch (4*yt+zt) {
1788
 @t\4@>@;
1789
 case 4*zro+zro: case 4*num+zro: case 4*inf+zro:
1790
 case 4*inf+num: case 4*inf+inf: x=standard_NaN;
1791
  exceptions|=I_BIT;@+break;
1792
 case 4*zro+num: case 4*zro+inf: case 4*num+inf: return y;
1793
 case 4*num+num: @;
1794
 zero_out: x=zero_octa;
1795
  }
1796
  if (ys=='-') x.h|=sign_bit;
1797
  return x;
1798
}
1799
 
1800
@ If there's a huge difference in exponents and the remainder is nonzero,
1801
this computation will take a long time. One could compute
1802
$(2^ny)\,{\rm rem}\,z$ much more quickly for large~$n$ by using $O(\log n)$
1803
multiplications modulo~$z$, but the floating remainder operation isn't
1804
important enough to justify such expensive hardware.
1805
 
1806
Results of floating remainder are always exact, so the rounding mode
1807
is immaterial.
1808
 
1809
@=
1810
odd=0; /* will be 1 if we've subtracted an odd multiple of~$z$ from $y$ */
1811
thresh=ye-delta;
1812
if (thresh
1813
while (ye>=thresh) @
1814
   |goto zero_out| if the remainder is zero,
1815
   |goto try_complement| if appropriate@>;
1816
if (ye>=ze) {
1817
  exceptions|=E_BIT;@+return fpack(yf,ye,ys,ROUND_OFF);
1818
}
1819
if (ye
1820
yf=shift_right(yf,1,1);
1821
try_complement: xf=ominus(zf,yf), xe=ze, xs='+' + '-' - ys;
1822
if (xf.h>yf.h || (xf.h==yf.h && (xf.l>yf.l || (xf.l==yf.l && !odd))))
1823
  xf=yf, xs=ys;
1824
while (xf.h<0x400000) xe--, xf=shift_left(xf,1);
1825
return fpack(xf,xe,xs,ROUND_OFF);
1826
 
1827
@ Here we are careful not to change the sign of |y|, because a remainder
1828
of~0 is supposed to inherit the original sign of~|y|.
1829
 
1830
@=
1831
{
1832
  if (yf.h==zf.h && yf.l==zf.l) goto zero_out;
1833
  if (yf.h
1834
    if (ye==ze) goto try_complement;
1835
    ye--, yf=shift_left(yf,1);
1836
  }
1837
  yf=ominus(yf,zf);
1838
  if (ye==ze) odd=1;
1839
  while (yf.h<0x400000) ye--, yf=shift_left(yf,1);
1840
}
1841
 
1842
@* Index.
1843
 

powered by: WebSVN 2.1.0

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