[9] | 1 | /*
|
---|
| 2 | (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
|
---|
| 3 | See the copyright notice in the ACK home directory, in the file "Copyright".
|
---|
| 4 | */
|
---|
| 5 |
|
---|
| 6 | /* $Header: /cvsup/minix/src/lib/ack/float/mul_ext.fc,v 1.1 2005/10/10 15:27:43 beng Exp $ */
|
---|
| 7 |
|
---|
| 8 | /*
|
---|
| 9 | ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
|
---|
| 10 | */
|
---|
| 11 |
|
---|
| 12 | # include "FP_bias.h"
|
---|
| 13 | # include "FP_trap.h"
|
---|
| 14 | # include "FP_types.h"
|
---|
| 15 | # include "FP_shift.h"
|
---|
| 16 |
|
---|
| 17 | void
|
---|
| 18 | mul_ext(e1,e2)
|
---|
| 19 | EXTEND *e1,*e2;
|
---|
| 20 | {
|
---|
| 21 | register int i,j; /* loop control */
|
---|
| 22 | unsigned short mp[4]; /* multiplier */
|
---|
| 23 | unsigned short mc[4]; /* multipcand */
|
---|
| 24 | unsigned short result[8]; /* result */
|
---|
| 25 | register unsigned short *pres;
|
---|
| 26 |
|
---|
| 27 | /* first save the sign (XOR) */
|
---|
| 28 | e1->sign ^= e2->sign;
|
---|
| 29 |
|
---|
| 30 | /* compute new exponent */
|
---|
| 31 | e1->exp += e2->exp + 1;
|
---|
| 32 | /* 128 bit multiply of mantissas */
|
---|
| 33 |
|
---|
| 34 | /* assign unknown long formats */
|
---|
| 35 | /* to known unsigned word formats */
|
---|
| 36 | mp[0] = e1->m1 >> 16;
|
---|
| 37 | mp[1] = (unsigned short) e1->m1;
|
---|
| 38 | mp[2] = e1->m2 >> 16;
|
---|
| 39 | mp[3] = (unsigned short) e1->m2;
|
---|
| 40 | mc[0] = e2->m1 >> 16;
|
---|
| 41 | mc[1] = (unsigned short) e2->m1;
|
---|
| 42 | mc[2] = e2->m2 >> 16;
|
---|
| 43 | mc[3] = (unsigned short) e2->m2;
|
---|
| 44 | for (i = 8; i--;) {
|
---|
| 45 | result[i] = 0;
|
---|
| 46 | }
|
---|
| 47 | /*
|
---|
| 48 | * fill registers with their components
|
---|
| 49 | */
|
---|
| 50 | for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
|
---|
| 51 | unsigned short k = 0;
|
---|
| 52 | unsigned long mpi = mp[i];
|
---|
| 53 | for(j=4;j--;) {
|
---|
| 54 | unsigned long tmp = (unsigned long)pres[j] + k;
|
---|
| 55 | if (mc[j]) tmp += mpi * mc[j];
|
---|
| 56 | pres[j] = tmp;
|
---|
| 57 | k = tmp >> 16;
|
---|
| 58 | }
|
---|
| 59 | pres[-1] = k;
|
---|
| 60 | }
|
---|
| 61 | if (! (result[0] & 0x8000)) {
|
---|
| 62 | e1->exp--;
|
---|
| 63 | for (i = 0; i <= 3; i++) {
|
---|
| 64 | result[i] <<= 1;
|
---|
| 65 | if (result[i+1]&0x8000) result[i] |= 1;
|
---|
| 66 | }
|
---|
| 67 | result[4] <<= 1;
|
---|
| 68 | }
|
---|
| 69 |
|
---|
| 70 | /*
|
---|
| 71 | * combine the registers to a total
|
---|
| 72 | */
|
---|
| 73 | e1->m1 = ((unsigned long)(result[0]) << 16) + result[1];
|
---|
| 74 | e1->m2 = ((unsigned long)(result[2]) << 16) + result[3];
|
---|
| 75 | if (result[4] & 0x8000) {
|
---|
| 76 | if (++e1->m2 == 0)
|
---|
| 77 | if (++e1->m1 == 0) {
|
---|
| 78 | e1->m1 = NORMBIT;
|
---|
| 79 | e1->exp++;
|
---|
| 80 | }
|
---|
| 81 | }
|
---|
| 82 |
|
---|
| 83 | /* check for overflow */
|
---|
| 84 | if (e1->exp >= EXT_MAX) {
|
---|
| 85 | trap(EFOVFL);
|
---|
| 86 | /* if caught */
|
---|
| 87 | /* return signed infinity */
|
---|
| 88 | e1->exp = EXT_MAX;
|
---|
| 89 | infinity: e1->m1 = e1->m2 =0L;
|
---|
| 90 | return;
|
---|
| 91 | }
|
---|
| 92 | /* check for underflow */
|
---|
| 93 | if (e1->exp < EXT_MIN) {
|
---|
| 94 | trap(EFUNFL);
|
---|
| 95 | e1->exp = EXT_MIN;
|
---|
| 96 | goto infinity;
|
---|
| 97 | }
|
---|
| 98 | }
|
---|