1 /*
2 * Copyright (c) 2014-2018, NVIDIA CORPORATION. All rights reserved.
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 *
16 */
17
18 #include "dblint64.h"
19
20 typedef union {
21 DBLINT64 wd; /* canonical msw & lsw view of long long values */
22 int hf[2]; /* native msw & lsw signed view of long long values */
23 unsigned uhf[2]; /* native msw & lsw unsigned view of long long values */
24 long long value;
25 unsigned long long uvalue;
26 } LL_SHAPE;
27
28 #undef FL
29 #undef UFL
30 #undef LSW
31 #undef MSW
32 #undef ULSW
33 #undef UMSW
34 #define FL(s) s.value
35 #define UFL(s) s.uvalue
36 /* Little endian view of a long long as two ints */
37 #define LSW(s) s.hf[0]
38 #define MSW(s) s.hf[1]
39 #define ULSW(s) s.uhf[0]
40 #define UMSW(s) s.uhf[1]
41
42 #define VOID void
43
44 #define UTL_I_I64RET(m, l) return (__utl_i_i64ret(m, l))
45 /*extern VOID UTL_I_I64RET();*/
46
47 long long
__utl_i_i64ret(int msw,int lsw)48 __utl_i_i64ret(int msw, int lsw)
49 {
50 union {
51 int i[2];
52 long long z;
53 } uu;
54 uu.i[0] = lsw;
55 uu.i[1] = msw;
56 return uu.z;
57 }
58
59 typedef int INT;
60 typedef unsigned int UINT;
61
62 #define __mth_i_kmul(a, b) (a) * (b)
63
64 INT __mth_i_kcmp(), __mth_i_kucmp(), __mth_i_kcmpz(), __mth_i_kucmpz();
65 VOID __mth_i_krshift();
66 VOID __mth_i_klshift();
67 VOID __mth_i_kurshift();
68 VOID __utl_i_add64(), __utl_i_div64();
69 VOID __utl_i_sub64(DBLINT64 arg1, DBLINT64 arg2, DBLINT64 result);
70 VOID __utl_i_mul64(), __utl_i_udiv64();
71 static VOID neg64(), shf64(), shf128by1();
72
73 /*
74 * Add 2 64-bit integers
75 *
76 * Arguments:
77 * arg1, arg2 64-bit addends
78 * result 64-bit result of arg1 + arg2
79 *
80 * Return value:
81 * none
82 */
__utl_i_add64(arg1,arg2,result)83 VOID __utl_i_add64(arg1, arg2, result) DBLINT64 arg1, arg2, result;
84 {
85 int carry; /* value to be carried from adding the lower
86 * 32 bits */
87 int sign_sum; /* sum of the sign bits of the low-order
88 * words of arg1 and arg2 */
89 sign_sum = ((unsigned)arg1[1] >> 31) + ((unsigned)arg2[1] >> 31);
90 result[1] = arg1[1] + arg2[1];
91 /*
92 * sign_sum = 2 -> carry guaranteed
93 * = 1 -> if sum sign bit off then carry
94 * = 0 -> carry not possible
95 */
96 carry = sign_sum > (int)((unsigned)result[1] >> 31);
97 result[0] = arg1[0] + arg2[0] + carry;
98 return;
99 }
100
101 /** \brief Subtract two 64-bit integers
102 *
103 * \param arg1 minuend
104 * \param arg2 subtrahend
105 * \param result arg1 - arg2
106 */
107 VOID
__utl_i_sub64(DBLINT64 arg1,DBLINT64 arg2,DBLINT64 result)108 __utl_i_sub64(DBLINT64 arg1, DBLINT64 arg2, DBLINT64 result)
109 {
110 int borrow; /* value to be borrowed from adding the lower
111 * 32 bits */
112 int sign_diff; /* difference of the sign bits of the
113 * low-order words of arg1 and arg2 */
114
115 sign_diff = ((unsigned)arg1[1] >> 31) - ((unsigned)arg2[1] >> 31);
116 result[1] = arg1[1] - arg2[1];
117 /*
118 * sign_diff = -1 -> borrow guaranteed
119 * = 0 -> if diff sign bit on (arg2 > arg1) then borrow
120 * 1 -> borrow not possible
121 */
122 borrow = sign_diff < (int)((unsigned)result[1] >> 31);
123 result[0] = (arg1[0] - borrow) - arg2[0];
124 return;
125 }
126
127 /*
128 *
129 * Multiply two 64-bit integers to produce a 64-bit
130 * integer product.
131 */
__utl_i_mul64(arg1,arg2,result)132 VOID __utl_i_mul64(arg1, arg2, result) DBLINT64 arg1, arg2, result;
133 {
134 LL_SHAPE v1, v2, r;
135
136 /* Canonical 64-bit (big endian) form -> little endian */
137 v1.wd[0] = arg1[1];
138 v1.wd[1] = arg1[0];
139
140 v2.wd[0] = arg2[1];
141 v2.wd[1] = arg2[0];
142
143 r.value = __mth_i_kmul(v1.value, v2.value);
144
145 /* Little endian form -> canonical form */
146 result[0] = r.wd[1];
147 result[1] = r.wd[0];
148 }
149
150 /*
151 * Divide two long long ints to produce a long long quotient.
152 */
153 long long
__mth_i_kdiv(long long x,long long y)154 __mth_i_kdiv(long long x, long long y)
155 {
156 int negate; /* flag indicating the result needs to be negated */
157 LL_SHAPE a, b, r;
158
159 FL(a) = x;
160 if (MSW(a) >= 0) {
161 negate = 0;
162 } else {
163 FL(a) = -FL(a);
164 negate = 1;
165 }
166
167 FL(b) = y;
168 if (MSW(b) < 0) {
169 FL(b) = -FL(b);
170 negate = !negate;
171 }
172
173 if (MSW(a) == 0 && MSW(b) == 0) {
174 MSW(r) = 0;
175 *(unsigned *)&LSW(r) = (unsigned)LSW(a) / (unsigned)LSW(b);
176 } else {
177 DBLINT64 arg1, arg2; /* DBLINT64 is big endian!! */
178 DBLINT64 result;
179 arg1[1] = LSW(a);
180 arg1[0] = MSW(a);
181 arg2[1] = LSW(b);
182 arg2[0] = MSW(b);
183 __utl_i_div64(arg1, arg2, result);
184 LSW(r) = result[1];
185 MSW(r) = result[0];
186 }
187
188 if (negate)
189 return -FL(r);
190 else
191 return FL(r);
192 }
193
194 /*
195 * Divide two unsigned long long ints to produce a unsigned long long quotient.
196 */
197 unsigned long long
__mth_i_ukdiv(unsigned long long x,unsigned long long y)198 __mth_i_ukdiv(unsigned long long x, unsigned long long y)
199 {
200 LL_SHAPE a, b, r;
201
202 UFL(a) = x;
203 UFL(b) = y;
204
205 if (UMSW(a) == 0 && UMSW(b) == 0) {
206 UMSW(r) = 0;
207 ULSW(r) = ULSW(a) / ULSW(b);
208 } else {
209 DBLUINT64 arg1, arg2; /* DBLUINT64 is big endian!! */
210 DBLUINT64 result;
211 arg1[1] = ULSW(a);
212 arg1[0] = UMSW(a);
213 arg2[1] = ULSW(b);
214 arg2[0] = UMSW(b);
215 __utl_i_udiv64(arg1, arg2, result);
216 ULSW(r) = result[1];
217 UMSW(r) = result[0];
218 }
219 return UFL(r);
220 }
221
222 /*
223 * Divide two 64-bit integers to produce a 64-bit
224 * integer quotient.
225 */
__utl_i_div64(arg1,arg2,result)226 VOID __utl_i_div64(arg1, arg2, result) DBLINT64 arg1, arg2, result;
227
228 {
229 DBLINT64 den; /* denominator used in calculating the
230 * quotient */
231 int i; /* for loop control variable */
232 int temp_result[4]; /* temporary result used in
233 * calculating the quotient */
234 int negate; /* flag which indicates the result needs to
235 * be negated */
236
237 if ((arg1[0] == 0 && arg1[1] == 0) || (arg2[0] == 0 && arg2[1] == 0)) {
238 result[0] = 0;
239 result[1] = 0;
240 return;
241 }
242 temp_result[0] = 0;
243 temp_result[1] = 0;
244
245 if (arg1[0] < 0) {
246 neg64(arg1, &temp_result[2]);
247 negate = 1;
248 } else {
249 temp_result[2] = arg1[0];
250 temp_result[3] = arg1[1];
251 negate = 0;
252 }
253
254 if (arg2[0] < 0) {
255 neg64(arg2, den);
256 negate = !negate;
257 } else {
258 den[0] = arg2[0];
259 den[1] = arg2[1];
260 }
261
262 for (i = 0; i < 64; ++i) {
263 shf128by1(temp_result, temp_result);
264 if (__mth_i_kucmp(temp_result[1], temp_result[0], den[1], den[0]) >= 0) {
265 __utl_i_sub64(temp_result, den, temp_result);
266 temp_result[3] = temp_result[3] + 1;
267 }
268 }
269
270 if (negate)
271 neg64(&temp_result[2], result);
272 else {
273 result[0] = temp_result[2];
274 result[1] = temp_result[3];
275 }
276 }
277
278 /*
279 * negate a 64-bit integer
280 *
281 * Arguments:
282 * arg 64-bit value to be negated
283 * result - arg
284 *
285 * Return value:
286 * none.
287 */
288
neg64(arg,result)289 static VOID neg64(arg, result) DBLINT64 arg, result;
290
291 {
292 int sign; /* sign of the low-order word of arg prior to
293 * being complemented */
294
295 sign = (unsigned)arg[1] >> 31;
296 result[0] = ~arg[0];
297 result[1] = (~arg[1]) + 1;
298 if (sign == 0 && result[1] >= 0)
299 result[0]++;
300 }
301
302 /** \brief Divide two 64-bit unsigned integers to produce a 64-bit unsigned
303 * integer quotient.
304 */
305 VOID
__utl_i_udiv64(DBLUINT64 arg1,DBLUINT64 arg2,DBLUINT64 result)306 __utl_i_udiv64(DBLUINT64 arg1, DBLUINT64 arg2, DBLUINT64 result)
307 {
308 DBLINT64 den; /* denominator used in calculating the
309 * quotient */
310 int i; /* for loop control variable */
311 int temp_result[4]; /* temporary result used in
312 * calculating the quotient */
313
314 if ((arg1[0] == 0 && arg1[1] == 0) || (arg2[0] == 0 && arg2[1] == 0)) {
315 result[0] = 0;
316 result[1] = 0;
317 return;
318 }
319 temp_result[0] = 0;
320 temp_result[1] = 0;
321 temp_result[2] = arg1[0];
322 temp_result[3] = arg1[1];
323
324 den[0] = arg2[0];
325 den[1] = arg2[1];
326
327 for (i = 0; i < 64; ++i) {
328 shf128by1(temp_result, temp_result);
329 if (__mth_i_kucmp(temp_result[1], temp_result[0], den[1], den[0]) >= 0) {
330 __utl_i_sub64(temp_result, den, temp_result);
331 temp_result[3] = temp_result[3] + 1;
332 }
333 }
334
335 result[0] = temp_result[2];
336 result[1] = temp_result[3];
337 }
338
339 /* this is in assembly now. New routines are krshift,kurshift,klshift */
340 /* leave these in here for now. So old objects/libs get resolved. */
341 /*
342 * shift a 64-bit signed integer
343 *
344 * Arguments:
345 * arg INT [2] operand to be shifted
346 * count int shift count (- => right shift; o.w., left shift).
347 * A count outside of the range -63 to 64 results in a result
348 * of zero; the caller takes into consideration machine
349 * dependicies (such as the shift count is modulo 64).
350 * result
351 *
352 * Return value:
353 * none.
354 */
355 /* new version - inline shf64 */
__mth_i_kicshft(op1,op2,count,direct)356 long long __mth_i_kicshft(op1, op2, count, direct) UINT op1,
357 op2; /* really INT */
358 INT count, direct;
359 {
360 DBLINT64 result;
361 if (count < 0 || count >= 64) {
362 UTL_I_I64RET(0, 0);
363 }
364 if (count == 0) {
365 result[0] = op2;
366 result[1] = op1;
367 UTL_I_I64RET(result[0], result[1]);
368 }
369 if (direct) {/* right shift */
370 if (count < 32) {
371 result[0] = (INT)op2 >> count; /* sign extend */
372 result[1] = (op1 >> count) | (op2 << (32 - count));
373 } else {
374 result[0] = (INT)op2 >> 31; /* sign extend */
375 result[1] = (INT)op2 >> (count - 32);
376 }
377 } else {/* left ushift */
378 if (count < 32) {
379 result[0] = (op2 << count) | (op1 >> (32 - count));
380 result[1] = op1 << count;
381 } else {
382 result[0] = op1 << (count - 32);
383 result[1] = 0;
384 }
385 }
386 UTL_I_I64RET(result[0], result[1]);
387 }
388
389 /* new version - inline ushf64 */
__mth_i_ukicshft(op1,op2,count,direct)390 long long __mth_i_ukicshft(op1, op2, count, direct) UINT op1, op2;
391 INT count, direct;
392 {
393 DBLINT64 result;
394 if (count < 0 || count >= 64) {
395 UTL_I_I64RET(0, 0);
396 }
397 if (count == 0) {
398 result[0] = op2;
399 result[1] = op1;
400 UTL_I_I64RET(result[0], result[1]);
401 }
402 if (direct) {/* right ushift */
403 if (count < 32) {
404 result[0] = op2 >> count;
405 result[1] = (op1 >> count) | (op2 << (32 - count));
406 } else {
407 result[0] = 0;
408 result[1] = op2 >> (count - 32);
409 }
410 } else {/* left ushift */
411 if (count < 32) {
412 result[0] = (op2 << count) | (op1 >> (32 - count));
413 result[1] = op1 << count;
414 } else {
415 result[0] = op1 << (count - 32);
416 result[1] = 0;
417 }
418 }
419 UTL_I_I64RET(result[0], result[1]);
420 }
421
__mth_i_kishft(op1,op2,arg2)422 long long __mth_i_kishft(op1, op2, arg2) INT op1, op2, arg2;
423 {
424 DBLINT64 arg1;
425 DBLINT64 result;
426 arg1[1] = op1;
427 arg1[0] = op2;
428 shf64(arg1, arg2, result);
429 UTL_I_I64RET(result[0], result[1]);
430 }
431
shf64(arg,count,result)432 static VOID shf64(arg, count, result) DBLINT64 arg;
433 int count;
434 DBLINT64 result;
435 {
436 DBLUINT64 u_arg; /* 'copy-in' unsigned value of arg */
437
438 if (count >= 64 || count <= -64) {
439 result[0] = 0;
440 result[1] = 0;
441 return;
442 }
443 if (count == 0) {
444 result[0] = arg[0];
445 result[1] = arg[1];
446 return;
447 }
448 u_arg[0] = arg[0];
449 u_arg[1] = arg[1];
450 if (count >= 0) {
451 if (count < 32) {
452 result[0] = (u_arg[0] << count) | (u_arg[1] >> (32 - count));
453 result[1] = u_arg[1] << count;
454 } else {
455 result[0] = u_arg[1] << (count - 32);
456 result[1] = 0;
457 }
458 } else if (count > -32) {
459 result[0] = arg[0] >> -count; /* sign extend */
460 result[1] = (u_arg[1] >> -count) | (u_arg[0] << (count + 32));
461 } else {
462 result[0] = arg[0] >> 31; /* sign extend */
463 result[1] = arg[0] >> (-count - 32);
464 }
465 }
466
shf128by1(arg,result)467 static VOID shf128by1(arg, result) INT arg[4];
468 INT result[4];
469 {
470 UINT u_arg[4];
471 int i; /* for loop control variable */
472
473 /* to get unsigned, and because arg and result may be same loc */
474 for (i = 0; i < 4; ++i)
475 u_arg[i] = arg[i];
476
477 for (i = 0; i < 3; ++i) {
478 result[i] = (u_arg[i] << 1) | (u_arg[i + 1] >> (32 - 1));
479 } /* end_for */
480 result[3] = (u_arg[3] << 1);
481 }
482
483 /* this is all in assembly now. */
484 /*
485 * compare 64-bit unsigned integer against zero
486 *
487 * Arguments:
488 * arg1 operand 1 of compare
489 *
490 * Return value:
491 * 0 = operand equal to zero
492 * 1 = operand greater than zero (ie. nonzero)
493 */
__mth_i_kucmpz(op1,op2)494 INT __mth_i_kucmpz(op1, op2) register UINT op1, op2;
495 {
496 if ((op1 | op2) == 0)
497 return 0;
498 return 1;
499 }
500
501 /*
502 * compare two 64-bit unsigned integers
503 *
504 * Arguments:
505 * arg1 operand 1 of compare
506 * arg2 operand 2 of compare
507 *
508 * Return value:
509 * -1 = first operand less than second
510 * 0 = first operand equal to second
511 * 1 = first operand greater than second
512 */
__mth_i_kucmp(op1,op2,op3,op4)513 INT __mth_i_kucmp(op1, op2, op3, op4) UINT op1, op2, op3, op4;
514 {
515 if (op2 == op4) {
516 if (op1 == op3)
517 return 0;
518 if (op1 < op3)
519 return -1;
520 return 1;
521 }
522 if (op2 < op4)
523 return -1;
524 return 1;
525 }
526
527 /*
528 * compare 64-bit signed integer against zero
529 */
__mth_i_kcmpz(op1,op2)530 INT __mth_i_kcmpz(op1, op2) INT op1, op2;
531 {
532 if (op2 == 0) {
533 if (op1 == 0)
534 return 0; /* is zero */
535 return 1; /* > zero */
536 }
537 if (op2 < 0)
538 return -1; /* < zero */
539 return 1; /* > zero */
540 }
541
__mth_i_kcmp(op1,op2,op3,op4)542 INT __mth_i_kcmp(op1, op2, op3, op4) INT op1, op2, op3, op4;
543 {
544 if (op2 == op4) {
545 if (op1 == op3)
546 return 0;
547 if ((unsigned)op1 < (unsigned)op3)
548 return -1;
549 return 1;
550 }
551 if (op2 < op4)
552 return -1;
553 return 1;
554 }
555
556 /**********************************************************
557 *
558 * Unpacked floating point routines.
559 * These routines will convert IEEE
560 * floating point to an unpacked internal
561 * format and back again.
562 *
563 *********************************************************/
564
565 /****************************************************************
566 *
567 * The IEEE floating point format is as follows:
568 *
569 * single:
570 * 31 30 23 22 0
571 * +----+-----------+ +----------------------+
572 * W0 |sign| exp + 127 | 1. | mantissa |
573 * +----+-----------+ +----------------------+
574 *
575 * double:
576 * 31 30 20 19 0
577 * +----+-------------+ +--------------------+
578 * W0 |sign| exp + 1023 | 1. | mantissa |
579 * +----+-------------+ +--------------------+
580 * 31 0
581 * +---------------------------------------------+
582 * W1 | mantissa |
583 * +---------------------------------------------+
584 *
585 ***************************************************************/
586
587 /* define bit manipulation macros: */
588
589 #define bitmask(w) ((UINT)(1UL << (w)) - 1)
590 #define extract(i, a, b) (((INT)(i) >> (b)) & bitmask((a) - (b) + 1))
591 #define bit_fld(i, a, b) (((INT)(i)&bitmask((a) - (b) + 1)) << (b))
592 #define bit_insrt(x, i, a, b) x = (((INT)(x) & ~bit_fld(~0L,(a),(b)))
593
594 typedef int IEEE32; /* IEEE single precision float number */
595 typedef int IEEE64[2]; /* IEEE double precision float number */
596
597 #define IEEE32_sign(f) extract(f, 31, 31)
598 #define IEEE32_exp(f) (extract(f, 30, 23) - 127)
599 #define IEEE32_man(f) (bit_fld(1, 23, 23) | extract(f, 22, 0))
600 #define IEEE32_zero(f) (extract(f, 30, 0) == 0)
601 #define IEEE32_infin(f) (extract(f, 30, 0) == bit_fld(255, 30, 23))
602 #define IEEE32_nan(f) (!IEEE32_infin(f) && IEEE32_exp(f) == 128)
603 #define IEEE32_pack(f, sign, exp, man) \
604 f = bit_fld(sign, 31, 31) | bit_fld(exp + 127, 30, 23) | \
605 bit_fld(man[0], 22, 0)
606
607 #define IEEE64_sign(d) extract(d[1], 31, 31)
608 #define IEEE64_exp(d) (extract(d[1], 30, 20) - 1023)
609 #define IEEE64_manhi(d) (bit_fld(1, 20, 20) | extract(d[1], 19, 0))
610 #define IEEE64_manlo(d) (d[0])
611 #define IEEE64_zero(d) (extract(d[1], 30, 0) == 0L && IEEE64_manlo(d) == 0L)
612 #define IEEE64_infin(d) \
613 (extract(d[1], 30, 0) == bit_fld(255, 30, 23) && d[0] == 0)
614 #define IEEE64_nan(d) (!IEEE64_infin(d) && IEEE64_exp(d) == 1024)
615 #define IEEE64_pack(d, sign, exp, man) \
616 d[1] = bit_fld(sign, 31, 31) | bit_fld(exp + 1023, 30, 20) | \
617 bit_fld(man[0], 19, 0); \
618 d[0] = man[1]
619
620 /****************************************************************
621 *
622 * The unpacked floating point format is as follows:
623 *
624 * +-----------+ +-----------+ +-------------+
625 * | val | | sign | | exp |
626 * +-----------+ +-----------+ +-------------+
627 * 31 20 19 0
628 * +------------------+ +--------------------+
629 * W0 | mantissa | . | mantissa |
630 * +------------------+ +--------------------+
631 * 31 0
632 * +---------------------------------------------+
633 * W1 | mantissa |
634 * +---------------------------------------------+
635 * 31 0
636 * +---------------------------------------------+
637 * W2 | mantissa |
638 * +---------------------------------------------+
639 * 31 0
640 * +---------------------------------------------+
641 * W3 | mantissa |
642 * +---------------------------------------------+
643 *
644 ***************************************************************/
645
646 typedef enum { ZERO, NIL, NORMAL, BIG, INFIN, NAN, DIVZ } VAL;
647
648 #define POS 0
649 #define NEG 1
650
651 typedef struct {
652 VAL fval; /* Value */
653 int fsgn; /* Sign */
654 int fexp; /* Exponent */
655 INT fman[4]; /* Mantissa */
656 } UFP;
657
manadd(m1,m2)658 static VOID manadd(m1, m2)
659 /* add mantissas of two unpacked floating point numbers. */
660 INT m1[4]; /* first and result */
661 INT m2[4]; /* second */
662 {
663 INT t1, t2, carry;
664 INT lo, hi;
665 int i;
666
667 carry = 0;
668 for (i = 3; i >= 0; i--) {
669 /* add low halves + carry */
670 t1 = m1[i] & 0x0000FFFFL;
671 t2 = m2[i] & 0x0000FFFFL;
672 lo = t1 + t2 + carry;
673 carry = (lo >> 16) & 0x0000FFFFL;
674 lo &= 0x0000FFFFL;
675
676 /* add high halves + carry */
677 t1 = (m1[i] >> 16) & 0x0000FFFFL;
678 t2 = (m2[i] >> 16) & 0x0000FFFFL;
679 hi = t1 + t2 + carry;
680 carry = (hi >> 16) & 0x0000FFFFL;
681 hi <<= 16;
682
683 /* merge halves */
684 m1[i] = hi | lo;
685 }
686 /* we assume no carry is
687 * generated from the last
688 * add */
689 }
690
manshftl(m,n)691 static VOID manshftl(m, n)
692 /* shift mantissa left */
693 INT m[4]; /* number */
694 int n; /* amount to shift */
695 {
696 register int i;
697 register int j;
698 int mask;
699
700 /* do whole word shifts first */
701 for (i = n; i >= 32; i -= 32) {
702 m[0] = m[1];
703 m[1] = m[2];
704 m[2] = m[3];
705 m[3] = 0L;
706 }
707 if (i > 0) {
708 j = 32 - i;
709 mask = bitmask(i);
710 m[0] = (m[0] << i) | ((m[1] >> j) & mask);
711 m[1] = (m[1] << i) | ((m[2] >> j) & mask);
712 m[2] = (m[2] << i) | ((m[3] >> j) & mask);
713 m[3] = (m[3] << i);
714 }
715 }
716
manshftr(m,n)717 static VOID manshftr(m, n)
718 /* shift mantissa right */
719 INT m[4]; /* number */
720 int n; /* amount to shift */
721 {
722 register int i;
723 register int j;
724 int mask;
725
726 /* do whole word shifts first */
727 for (i = n; i >= 32; i -= 32) {
728 m[3] = m[2];
729 m[2] = m[1];
730 m[1] = m[0];
731 m[0] = 0L;
732 }
733 if (i > 0) {
734 j = 32 - i;
735 mask = bitmask(j);
736 m[3] = ((m[3] >> i) & mask) | (m[2] << j);
737 m[2] = ((m[2] >> i) & mask) | (m[1] << j);
738 m[1] = ((m[1] >> i) & mask) | (m[0] << j);
739 m[0] = (m[0] >> i) & mask;
740 }
741 }
742
manrnd(m,bits)743 static VOID manrnd(m, bits)
744 /* Round the mantissa to the number of bits specified. */
745 INT m[4]; /* Mantissa to round */
746 int bits; /* Number of bits of precision */
747 {
748 int rndwrd, rndbit;
749 int oddwrd, oddbit;
750 INT round[4];
751 static INT one[4] = {0L, 0L, 0L, 1L};
752
753 rndwrd = bits / 32;
754 rndbit = 31 - (bits % 32);
755 oddwrd = (bits - 1) / 32;
756 oddbit = 31 - ((bits - 1) % 32);
757
758 /* check round bit */
759 if (extract(m[rndwrd], rndbit, rndbit) == 1) {
760 round[0] = 0xFFFFFFFFL;
761 round[1] = 0xFFFFFFFFL;
762 round[2] = 0xFFFFFFFFL;
763 round[3] = 0xFFFFFFFFL;
764 manshftr(round, bits + 1);
765 /* Round up by just under
766 * 1/2.
767 */
768 manadd(m, round);
769
770 /* If exactly 1/2 then round
771 * only if it's an odd number.
772 */
773 if (extract(m[rndwrd], rndbit, rndbit) == 1 &&
774 extract(m[oddwrd], oddbit, oddbit) == 1) {
775 manadd(m, one);
776 }
777 }
778 /* Get rid of rounded off bits. */
779 manshftr(m, 128 - bits);
780 manshftl(m, 128 - bits);
781 }
782
ufpnorm(u)783 static VOID ufpnorm(u)
784 /* normalize the unpacked number. */
785 register UFP *u; /* unpacked number */
786 {
787 /* If it's zero then there's no
788 * need to normalize. */
789 if (u->fman[0] == 0 && u->fman[1] == 0 && u->fman[2] == 0 && u->fman[3] == 0)
790 return;
791
792 /* We may have to right shift */
793 while (extract(u->fman[0], 31, 21) != 0) {
794 manshftr(u->fman, 1);
795 u->fexp++;
796 }
797
798 /* repeatedly left shift and decrement
799 * exp until normalized */
800 while (extract(u->fman[0], 20, 20) == 0) {
801 manshftl(u->fman, 1);
802 u->fexp--;
803 }
804 }
805
ufprnd(u,bits)806 static VOID ufprnd(u, bits)
807 /* Round the unpacked floating point to the number of bits of binary
808 * fraction specified.
809 */
810 UFP *u; /* Number to round */
811 int bits; /* Number of bits of precision */
812 {
813 VOID ufpnorm();
814
815 ufpnorm(u);
816 manrnd(u->fman, bits + 12); /* Binary point is 12 bits over */
817 ufpnorm(u); /* Might have denormalized it. */
818 }
819
820 /** \brief Floating point double unpack
821 * \param d double precision number
822 * \param u pointer to area for unpack
823 */
824 static VOID
dtoufp(IEEE64 d,register UFP * u)825 dtoufp(IEEE64 d, register UFP *u)
826 {
827 u->fval = NORMAL;
828 u->fexp = IEEE64_exp(d);
829 u->fsgn = IEEE64_sign(d);
830 u->fman[0] = IEEE64_manhi(d);
831 u->fman[1] = IEEE64_manlo(d);
832 u->fman[2] = 0L;
833 u->fman[3] = 0L;
834
835 /* special case checks */
836 if (IEEE64_nan(d))
837 u->fval = NAN;
838
839 if (IEEE64_infin(d))
840 u->fval = INFIN;
841
842 if (IEEE64_zero(d)) {
843 u->fval = ZERO;
844 u->fexp = 0;
845 u->fman[0] = 0L;
846 u->fman[1] = 0L;
847 }
848 }
849
ftoufp(f,u)850 static VOID ftoufp(f, u)
851 /* -- floating point single unpack */
852 IEEE32 *f; /* single precision number */
853 register UFP *u; /* unpacked result */
854 {
855 u->fval = NORMAL;
856 u->fexp = IEEE32_exp(*f);
857 u->fsgn = IEEE32_sign(*f);
858 u->fman[0] = IEEE32_man(*f);
859 u->fman[1] = 0L;
860 u->fman[2] = 0L;
861 u->fman[3] = 0L;
862 manshftr(u->fman, 3);
863
864 /* special case checks */
865 if (IEEE32_nan(*f))
866 u->fval = NAN;
867
868 if (IEEE32_infin(*f))
869 u->fval = INFIN;
870
871 if (IEEE32_zero(*f)) {
872 u->fval = ZERO;
873 u->fexp = 0;
874 u->fman[0] = 0L;
875 u->fman[1] = 0L;
876 }
877 }
878
i64toufp(i,u)879 static VOID i64toufp(i, u)
880 /* 64-bit integer to unpacked float */
881 DBLINT64 i;
882 UFP *u;
883 {
884 DBLINT64 tmp;
885
886 if (i[0] == 0L && i[1] == 0L) {
887 u->fsgn = POS;
888 u->fval = ZERO;
889 u->fexp = 0;
890 u->fman[0] = 0L;
891 u->fman[1] = 0L;
892 return;
893 }
894 tmp[0] = i[0];
895 tmp[1] = i[1];
896 u->fval = NORMAL;
897 u->fexp = 52;
898 if ((*i & 0x80000000L) != 0) {
899 u->fsgn = NEG;
900 neg64(i, tmp);
901 } else
902 u->fsgn = POS;
903 u->fman[0] = tmp[0];
904 u->fman[1] = tmp[1];
905 u->fman[2] = 0L;
906 u->fman[3] = 0L;
907 }
908
ufptod(u,r)909 static VOID ufptod(u, r)
910 /* floating point double pack */
911 register UFP *u; /* unpacked number */
912 IEEE64 r; /* packed result */
913 {
914 /* Round and normalize the unpacked
915 * number first. */
916 ufprnd(u, 52);
917
918 if (u->fval == ZERO || u->fval == NIL) {/* zero */
919 u->fexp = -1023;
920 u->fman[0] = 0L;
921 u->fman[1] = 0L;
922 }
923
924 if (u->fval == NAN) {/* Not a number */
925 u->fexp = 1024;
926 u->fman[0] = ~0L;
927 u->fman[1] = ~0L;
928 }
929
930 if (u->fval == INFIN || u->fval == BIG || u->fval == DIVZ) {/* infinity */
931 u->fexp = 1024;
932 u->fman[0] = 0L;
933 u->fman[1] = 0L;
934 }
935
936 if (u->fval == NORMAL && u->fexp <= -1023) {/* underflow */
937 u->fval = NIL;
938 u->fexp = -1023;
939 u->fman[0] = 0L;
940 u->fman[1] = 0L;
941 }
942
943 if (u->fval == NORMAL && u->fexp >= 1024) {/* overflow */
944 u->fval = BIG;
945 u->fexp = 1024;
946 u->fman[0] = 0L;
947 u->fman[1] = 0L;
948 }
949
950 IEEE64_pack(r, u->fsgn, u->fexp, u->fman);
951 }
952
ufptof(u,r)953 static VOID ufptof(u, r)
954 /* unpacked floating point single float */
955 register UFP *u; /* unpacked number */
956 IEEE32 *r; /* packed result */
957 {
958 /* Round and normalize the unpacked
959 * number first. */
960 ufprnd(u, 23);
961 manshftl(u->fman, 3);
962
963 if (u->fval == ZERO || u->fval == NIL) {/* zero */
964 u->fexp = -127;
965 u->fman[0] = 0L;
966 u->fman[1] = 0L;
967 }
968
969 if (u->fval == INFIN || u->fval == BIG || u->fval == DIVZ) {/* infinity */
970 u->fexp = 128;
971 u->fman[0] = 0L;
972 u->fman[1] = 0L;
973 }
974
975 if (u->fval == NAN) {/* Not a number */
976 u->fexp = 128;
977 u->fman[0] = ~0L;
978 u->fman[1] = ~0L;
979 }
980
981 if (u->fval == NORMAL && u->fexp <= -127) {/* underflow */
982 u->fval = NIL;
983 u->fexp = -127;
984 u->fman[0] = 0L;
985 u->fman[1] = 0L;
986 }
987
988 if (u->fval == NORMAL && u->fexp >= 128) {/* overflow */
989 u->fval = BIG;
990 u->fexp = 128;
991 u->fman[0] = 0L;
992 u->fman[1] = 0L;
993 }
994
995 IEEE32_pack(*r, u->fsgn, u->fexp, u->fman);
996 }
997
ufptoi64(u,i)998 static VOID ufptoi64(u, i)
999 /* unpacked float to 64-bit integer */
1000 UFP *u;
1001 DBLINT64 i;
1002 {
1003 /* Normalize the unpacked
1004 * number first. */
1005 ufpnorm(u);
1006 if (u->fexp - 52 > 0)
1007 manshftl(u->fman, u->fexp - 52);
1008 else
1009 manshftr(u->fman, 52 - u->fexp);
1010
1011 if (u->fval == ZERO || u->fval == NIL) {
1012 i[0] = 0;
1013 i[1] = 0;
1014 return;
1015 }
1016
1017 if (u->fval == NAN) {/* Not a number */
1018 i[0] = 0;
1019 i[1] = 0;
1020 return;
1021 }
1022
1023 if (u->fval == INFIN || u->fval == BIG || u->fexp > 62 ||
1024 ((u->fman[0] & 0x80000000L) != 0L && u->fman[1] == 0L)) {/* overflow */
1025 u->fval = BIG;
1026 if (u->fsgn == NEG) {
1027 i[0] = 0x80000000L;
1028 i[1] = 0x00000000L;
1029 } else {
1030 i[0] = 0x7FFFFFFFL;
1031 i[1] = 0xFFFFFFFFL;
1032 }
1033 return;
1034 }
1035
1036 i[0] = u->fman[0];
1037 i[1] = u->fman[1];
1038 if (u->fsgn == NEG)
1039 neg64(i, i);
1040 }
1041
__utl_i_dfix64(d,i)1042 VOID __utl_i_dfix64(d, i)
1043 /* double precision to 64-bit integer */
1044 double d; /*IEEE64 format and double are LITTLE_ENDIAN */
1045 DBLINT64 i;
1046 {
1047 UFP u;
1048
1049 dtoufp((int*)&d, &u);
1050 ufptoi64(&u, i);
1051 }
1052
__utl_i_dflt64(i)1053 double __utl_i_dflt64(i)
1054 /* 64 -- 64-bit integer to double */
1055 DBLINT64 i;
1056 {
1057 UFP u;
1058 IEEE64 d;
1059
1060 i64toufp(i, &u);
1061 ufptod(&u, d);
1062 return *((double *)d); /*IEEE64 format and double are LITTLE_ENDIAN */
1063 }
1064
__utl_i_fix64(float ff,DBLINT64 i)1065 VOID __utl_i_fix64(float ff, DBLINT64 i) /* use prototype to pass as float */
1066 /* single float to 64-bit */
1067 {
1068 IEEE32 f;
1069 UFP u;
1070
1071 f = *(IEEE32 *)&ff; /*IEEE32 format and float are LITTLE_ENDIAN*/
1072 ftoufp(&f, &u);
1073 ufptoi64(&u, i);
1074 }
1075
__utl_i_flt64(DBLINT64 i)1076 float __utl_i_flt64(DBLINT64 i) /* use prototype to return as float */
1077 /* 64-bit integer to single precision */
1078 {
1079 UFP u;
1080 IEEE32 f;
1081
1082 i64toufp(i, &u);
1083 ufptof(&u, &f);
1084 return *((float *)&f); /*IEEE32 format and float are LITTLE_ENDIAN */
1085 }
1086