1 /*
2 * $Id: ieee.i,v 1.1 2005-09-18 22:06:15 dhmunro Exp $
3 * work out native numeric formats, provide ieee NaN and Inf tests
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 local ieee;
12 /* DOCUMENT ieee.i
13 native_align, native_fix, native_flt, native_flim, native_dlim
14 describe the binary formats of the native primitive numeric types
15 as_chars(x) -- gets/sets bits of x as char array
16 ieee_test(x) -- tests for ieee754 special values
17 ieee_set(x,what) -- sets ieee754 special values
18 */
19
20 local native_align;
21 /* DOCUMENT native_align =
22 [char_align, short_align, int_align, long_align,
23 float_align, double_align, pointer_align, struct_align]
24
25 struct_align may be -1 for ppc/ibm interpretation
26 -2 for ppc/gcc interpretation
27 */
28
29 local native_fix;
30 /* DOCUMENT native_fix = [short_format, int_format, long_format]
31 format = [sizeof(type), order]
32 order = 1 for big-endian (most significant byte first)
33 -1 for little-endian (least significant byte first)
34 */
35
36 local native_flt;
37 /* DOCUMENT native_flt = [float_format, double_format]
38 format = [sizeof(type), order, sgn_addr, exp_addr, exp_size,
39 man_addr, man_size, man_norm, exp_bias, denorm]
40 order = 1 for big-endian (most significant byte first)
41 -1 for little-endian (least significant byte first)
42 (interpreted code here does not yet handle VAX middle-endian)
43 addr, size are bit address (from msb), number of bits
44 man_norm is 0 if mantissa 1-bit implied, 1 if 1-bit explicit
45 exp_bias is the exponent field for 1.0
46 denorm is 0 if there are no denormals
47 1 if there are functioning denormals
48 2 if denormals exist but are treated as NaNs
49 */
50
51 local native_flim, native_dlim;
52 /* DOCUMENT native_flim = [big, tiny, eps]
53 big is the biggest number
54 tiny is the smallest non-denormal number
55 eps is smallest number such that 1+eps != 1
56 */
57
as_chars(x,xnew)58 func as_chars(x, xnew)
59 /* DOCUMENT as_chars(x)
60 as_chars, x, xnew
61 return the bits of X as an array of char
62 return value has leading dimension of sizeof(x(1)), otherwise
63 same dimensions as X
64 in second form, sets bits of X to char array XNEW
65 SEE ALSO: ieee_test, ieee_set
66 */
67 {
68 size = sizeof(x)/numberof(x);
69 if (size < 2) return x;
70 dims = dimsof(x);
71 dims = grow([dims(1)+1],dims);
72 dims(2) = size;
73 xc = [];
74 reshape, xc, &x, char, dims;
75 if (!am_subroutine()) x = xc;
76 if (!is_void(xnew)) xc = xnew;
77 reshape, xc;
78 if (am_subroutine()) return;
79 return x;
80 }
81
82 /* note: struct and func names are bogus --
83 * intended to be overwritten later
84 */
85
86 struct _numfmt2 { char x; char y; }
87 native_align = [sizeof(_numfmt2)-sizeof(char)];
88 struct _numfmt2 { char x; short y; }
89 grow, native_align, [sizeof(_numfmt2)-sizeof(short)];
90 struct _numfmt2 { char x; int y; }
91 grow, native_align, [sizeof(_numfmt2)-sizeof(int)];
92 struct _numfmt2 { char x; long y; }
93 grow, native_align, [sizeof(_numfmt2)-sizeof(long)];
94 struct _numfmt2 { char x; float y; }
95 grow, native_align, [sizeof(_numfmt2)-sizeof(float)];
96 struct _numfmt2 { char x; double y; }
97 grow, native_align, [sizeof(_numfmt2)-sizeof(double)];
98 struct _numfmt2 { char x; pointer y; }
99 grow, native_align, [sizeof(_numfmt2)-sizeof(pointer)];
100 struct _numfmt1 { char c; }
101 struct _numfmt2 { char x; _numfmt1 y; }
102 grow, native_align, [sizeof(_numfmt2)-sizeof(_numfmt1)];
103
104 struct _numfmt1 { double d; char c; }
105 struct _numfmt2 { char x; _numfmt1 y; }
106 native_fix = [sizeof(_numfmt1), sizeof(_numfmt2)];
107 struct _numfmt1 { double d(1); char c; }
108 grow, native_fix, [sizeof(_numfmt1)];
109 if (native_fix(1)!=sizeof(double)+native_align(6) &&
110 native_fix(1)==2*sizeof(double)) {
111 if (native_fix(2) == 2*sizeof(double)+native_align(6))
112 native_align(8) = -1;
113 else if (native_fix(2) == 3*sizeof(double) &&
114 native_fix(3) == sizeof(double)+native_align(6))
115 native_align(8) = -2;
116 }
117
_numfmt1(type)118 func _numfmt1(type) /* computes integer (fix point) order */
119 {
120 x = type();
121 size = sizeof(x);
122 if (size<2 || structof(x+0)!=long) return 0;
123 xc = [];
124 reshape, xc, &x, char, [1,size];
125 for (i=1 ; i<=size ; i++) xc(i) = i;
126 reshape, xc;
127 for (i=0 ; i<size ; i++) {
128 y = (x>>(8*i)) & 0xff;
129 if (y == 1) break;
130 }
131 if (i >= size) error, "this can't be happening?";
132 if (2*i < size) {
133 word = order = i+1;
134 if (word == 1) order = -1;
135 } else {
136 word = order = size-i;
137 if (word != 1) order = -word;
138 }
139 if (word != 1) error, "middle endian integers (previously) unknown!";
140 return order;
141 }
142
143 native_fix = [[sizeof(short),_numfmt1(short)],
144 [sizeof(int),_numfmt1(int)],
145 [sizeof(long),_numfmt1(long)]];
146
_numfmt1(aa,bb,msb)147 func _numfmt1(aa, bb, msb) /* returns bit address where aa!=bb */
148 {
149 if (is_void(msb)) msb = 1;
150 size = sizeof(aa);
151 if (size<2) {
152 a = aa;
153 b = bb;
154 } else {
155 reshape, a, &aa, char, [1,size];
156 reshape, b, &bb, char, [1,size];
157 }
158 i = where(a != b);
159 diff = (a ~ b) & 0xff;
160 reshape, a;
161 reshape, b;
162 if (!numberof(i)) return -1;
163 i = i(msb);
164 diff = diff(i);
165 addr = msb? 8*(i-1) : 8*(size-i);
166 for (mask=0x80 ; !(mask&diff) ; mask>>=1) addr++;
167 return addr;
168 }
169
170 func _numfmt2(type, &epsilon) /* analyzes floating point type */
171 {
172 x = type();
173 size = sizeof(x);
174 if (structof(x+0)==long || structof(x+0.)!=double) return 0;
175
176 /* 1.0, 1.5, 1.5+0.5/256., ... 1.5+0.5/256.^n
177 * each term differs from previous by one bit, and the
178 * bit moves by one byte toward lower significance each time
179 */
180 one = type(1.0);
181 p256 = array(type, size);
182 for (i=1,s=0.5 ; i<=size ; i++,s*=1./256.) p256(i) = s;
183 p256c = one+p256(cum);
184 addr = array(0, size);
185 for (i=1 ; i<=size ; i++)
186 addr(i) = _numfmt1(p256c(i), p256c(i+1));
187 len = where(addr>=0)(0);
188 step = addr(dif:1:len);
189 if (anyof(step(dif))) error, "can't handle VAX float format yet";
190 if (step(1)==8) order = 1;
191 else if (step(1)==-8) order = -1;
192 else error, "unknown floating point format";
193
194 if (order > 0) msb = 1;
195 else msb = 0;
196
197 /* sign address is where 1.0 differs from -1.0 */
198 sgn_addr = _numfmt1(type(1.0), type(-1.0), msb);
199
200 /* exponent address is where 2.0^96 differs from 0.5^96 */
201 exp_addr = _numfmt1(type(65536.^6), type((1.0/65536.)^6), msb);
202
203 /* mantissa address is where 1.0 differs from 1.5 */
204 man_addr = _numfmt1(type(1.0), type(1.5), msb);
205
206 /* mantissa size p256c(len+1) is last different 1.5+0.5/256.^n
207 * scan down in steps of 0.5 to find epsilon */
208 base = p256c(len+1);
209 epsilon = p256(len);
210 man_size = addr(1)-addr(len)+1;
211 for (i=1 ; i<=8 ; epsilon*=type(0.5),man_size++,i++)
212 if (_numfmt1(base, base+type(0.5)*epsilon, msb) < 0) break;
213
214 /* mantissa normalization: check whether bit before man_addr
215 * is always set in series 1.5*2.0^n, -1.75*1.5*2.0^n, n=0 to 15 */
216 i = (man_addr-1)%8;
217 mask = 0x80 >> i;
218 i = (man_addr-1)/8;
219 if (order < 0) i = size-i;
220 else i++;
221 norm = type(1.5)*type(2.0)^indgen(0:15);
222 grow, norm, type(-1.75)*norm;
223 norm = as_chars(norm) & mask;
224 man_norm = allof(norm);
225 if (man_norm) {
226 man_addr--;
227 man_size++;
228 }
229
230 /* exponent size is where 1.0 and 2.0 differ, scanned from lsb */
231 test = as_chars(type(1.0)) ~ as_chars(type(2.0));
232 if (order < 0) test = test(::-1);
233 a = where(test)(0);
234 test = test(a);
235 for (i=0 ; i<8 ; i++) if ((1<<i)&test) break;
236 a = 8*(a-1) + 7-i;
237 exp_size = a - exp_addr + 1;
238
239 /* exponent bias is exponent for 1.0 */
240 one = as_chars(type(1.0));
241 if (order < 0) one = one(::-1);
242 i1 = a/8 + 1;
243 i0 = exp_addr/8 + 1;
244 one(i0) &= (1<<(8-exp_addr%8))-1;
245 s = 7-(a%8);
246 bias = long(one(i1)) >> s;
247 for (i=i1-1,s=8-s ; i>=i0 ; i--,s+=8) { bias |= long(one(i))<<s; }
248
249 /* try to decide if this is ieee754 type
250 * idea is to look at whether denormals exist */
251 denorm = 0;
252 if (!man_norm) {
253 small = [type(0.)];
254 c = array(char, size);
255 c(i1) = 1<<(7-(a%8));
256 if (order < 0) c = c(::-1);
257 as_chars, small, c;
258 denorm = (type(0.25)*small(1)) != 0.;
259 if (!denorm && (type(0.5)*small(1)) == 0.) denorm = 2;
260 }
261
262 return [size, order, sgn_addr, exp_addr, exp_size,
263 man_addr, man_size, man_norm, bias, denorm];
264 }
265
266 native_flt = [_numfmt2(float, native_flim), _numfmt2(double, native_dlim)];
267
_numfmt2(type)268 func _numfmt2(type)
269 {
270 if (type==float) fmt = native_flt(,1);
271 else if (type==double) fmt = native_flt(,2);
272
273 sgn_addr = fmt(3);
274 exp_addr = fmt(4);
275 exp_size = fmt(5);
276 man_addr = fmt(6);
277 man_size = fmt(7);
278 man_norm = fmt(8);
279
280 c = array(char, fmt(1), 2);
281 if (fmt(2) < 0) c = c(::-1,..);
282
283 a = exp_size+exp_addr-1;
284 if (fmt(10)) a--;
285 i1 = a/8 + 1;
286 i0 = exp_addr/8 + 1;
287 mask = (1<<(8-(exp_addr%8)))-1;
288 if (i1==i0) {
289 mask &= ~((1<<(7-(a%8)))-1);
290 } else {
291 c(i1,1) |= ~((1<<(7-(a%8)))-1);
292 }
293 c(i0,1) |= mask;
294 for (i=i0+1 ; i<i1 ; i++) c(i,1) |= 0xff;
295
296 a = man_size+man_addr-1;
297 i1 = a/8 + 1;
298 i0 = man_addr/8 + 1;
299 mask = (1<<(8-(man_addr%8)))-1;
300 if (i1==i0) {
301 mask &= ~((1<<(7-(a%8)))-1);
302 } else {
303 c(i1,1) |= ~((1<<(7-(a%8)))-1);
304 }
305 c(i0,1) |= mask;
306 for (i=i0+1 ; i<i1 ; i++) c(i,1) |= 0xff;
307
308 if (man_norm) {
309 i = man_addr/8 + 1;
310 mask = 1 << (7-man_addr%8);
311 c(i,2) |= mask;
312 }
313 if (fmt(10) || !man_norm) {
314 a = exp_size+exp_addr-1;
315 i = a/8 + 1;
316 mask = 1 << (7-a%8);
317 c(i,2) |= mask;
318 }
319
320 if (fmt(2) < 0) c = c(::-1,..);
321 x = array(type, 2);
322 as_chars, x, c;
323 return x;
324 }
325
326 native_flim = grow(_numfmt2(float), native_flim);
327 native_dlim = grow(_numfmt2(double), native_dlim);
328
329 _numfmt1 = _numfmt2 = [];
330
ieee_test(x)331 func ieee_test(x)
332 /* DOCUMENT ieee_test(x)
333 return values:
334 0 if this is an ordinary number
335 -1 if this is -Inf
336 1 if this is Inf
337 2 if this is qNaN
338 3 if this is sNaN
339 4 if this is a denormal
340 5 if this is a denormal which will be treated as NaN
341
342 Warning-- apparently there is no universal standard for what
343 constitutes signalling versus quiet NaN
344 on MIPS and HPPA architectures, qNaN and sNaN are reversed
345
346 SEE ALSO: ieee_set, as_chars
347 */
348 {
349 if (structof(x)==float) fmt = native_flt(,1);
350 else if (structof(x)==double) fmt = native_flt(,2);
351 else error, "can only test float or double arrays";
352 if (!fmt(10)) return array(0n, dimsof(x));
353
354 x = as_chars(x);
355
356 if (fmt(2) < 0) x = x(::-1,..);
357 sgn_addr = fmt(3);
358 exp_addr = fmt(4);
359 exp_size = fmt(5);
360 man_addr = fmt(6);
361 man_size = fmt(7);
362
363 i0 = sgn_addr/8+1;
364 i1 = 7-sgn_addr%8;
365 p = int(x(i0,..)>>i1) & 1n; /* sign bit */
366 x(i0,..) &= ~(1<<i1);
367 z = !x(sum,..);
368
369 a = exp_size+exp_addr-1;
370 i1 = a/8 + 1;
371 i0 = exp_addr/8 + 1;
372 s = 7-(a%8);
373 e = long(x(i1,..)) >> s;
374 for (i=i1-1,s=8-s ; i>=i0 ; i--,s+=8) { e |= long(x(i,..))<<s; }
375 mask = (1<<exp_size) - 1;
376 e &= mask;
377
378 d = !e & !z; /* denormals */
379 i = (e==mask); /* inf and nan */
380
381 i0 = man_addr/8+1;
382 i1 = 7-man_addr%8;
383 q = int(x(i0,..)>>i1) & 1n; /* first mantissa bit */
384
385 a = man_size+man_addr-1;
386 i2 = a/8+1;
387 x = x(i0:i2,..);
388 x(1,..) &= (1<<(1+i1))-1;
389 x(0,..) &= ~((1<<(7-a%8))-1);
390 n = i & (x(sum,..)!=0);
391 q &= n;
392 i &= !n;
393
394 if (dimsof(i)(1)==0) {
395 scalar = 1;
396 x = [i];
397 } else {
398 x = i;
399 }
400 list = where(p & i);
401 if (numberof(list)) x(list) = -1;
402 list = where(q);
403 if (numberof(list)) x(list) = 2;
404 list = where(n & !q);
405 if (numberof(list)) x(list) = 3;
406 list = where(d);
407 if (numberof(list)) x(list) = (fmt(10)==2? 5 : 4);
408 if (scalar) x = x(1);
409
410 return x;
411 }
412
ieee_set(x,what)413 func ieee_set(x, what)
414 /* DOCUMENT ieee_set, x, what
415 set X to ieee754 special value WHAT
416 X must be an array of float or double values
417 (note that X cannot be a scalar double value)
418 WHAT = 0 means leave unchanged
419 WHAT = 1 means set to Inf
420 WHAT = 2 means set to qNaN
421 WHAT = 3 means set to sNaN
422 WHAT = 4 means set to 0.0
423 negate WHAT to set the sign bit of X as well
424
425 WHAT may be an array conformable with X, in order to set only
426 some values of X
427
428 this routine is a no-op if this machine is not known to
429 support these ieee754 special values
430
431 Warning-- apparently there is no universal standard for what
432 constitutes signalling versus quiet NaN
433 on MIPS and HPPA architectures, qNaN and sNaN are reversed
434
435 SEE ALSO: ieee_test, as_chars
436 */
437 {
438 if (structof(x)==float) fmt = native_flt(,1);
439 else if (structof(x)==double) fmt = native_flt(,2);
440 else error, "can only test float or double arrays";
441 if (!fmt(10) || noneof(what)) return;
442 sgn_addr = fmt(3);
443 exp_addr = fmt(4);
444 exp_size = fmt(5);
445 man_addr = fmt(6);
446 man_size = fmt(7);
447
448 dims = dimsof(x);
449 what += array(0, dims);
450 s = what < 0;
451 what = abs(what);
452 c = as_chars(x);
453 if (fmt(2) < 0) c = c(::-1,..);
454
455 list = where((what>=1) & (what<=5));
456 if (numberof(list)) c(,list) = 0;
457
458 list = where((what>=1) & (what<=3));
459 if (numberof(list)) {
460 a = exp_size+exp_addr-1;
461 i1 = a/8 + 1;
462 i0 = exp_addr/8 + 1;
463 mask = (1<<(8-(exp_addr%8)))-1;
464 if (i1==i0) {
465 mask &= ~((1<<(7-(a%8)))-1);
466 } else {
467 c(i1,list) |= ~((1<<(7-(a%8)))-1);
468 }
469 c(i0,list) |= mask;
470 for (i=i0+1 ; i<i1 ; i++) c(i,list) |= 0xff;
471
472 list = where(what==2);
473 if (numberof(list)) {
474 i = man_addr/8 + 1;
475 mask = 1 << (7-man_addr%8);
476 c(i,list) |= mask;
477 }
478
479 list = where(what==3);
480 if (numberof(list)) {
481 a = man_size+man_addr-1;
482 i1 = a/8 + 1;
483 i0 = (man_addr+1)/8 + 1;
484 mask = (1<<(8-((man_addr+1)%8)))-1;
485 if (i1==i0) {
486 mask &= ~((1<<(7-(a%8)))-1);
487 } else {
488 c(i1,list) |= ~((1<<(7-(a%8)))-1);
489 }
490 c(i0,list) |= mask;
491 for (i=i0+1 ; i<i1 ; i++) c(i,list) |= 0xff;
492 }
493 }
494
495 list = where(s);
496 if (numberof(list)) {
497 i = sgn_addr/8 + 1;
498 mask = 1 << (7-sgn_addr%8);
499 c(i,list) |= mask;
500 }
501
502 if (fmt(2) < 0) c = c(::-1,..);
503 as_chars, x, c;
504 }
505