1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22 #pragma ident "%Z%%M% %I% %E% SMI"
23
24 /*
25 * Copyright (c) 1988 by Sun Microsystems, Inc.
26 */
27
28 #include "_Qquad.h"
29 #include "_Qglobals.h"
30
31
32 void
_fp_mul(px,py,pz)33 _fp_mul(px, py, pz)
34 unpacked *px, *py, *pz;
35
36 {
37 unpacked *pt;
38 unsigned acc[4]; /* Product accumulator. */
39 unsigned i,j,y,*x,s,r,c;
40
41 if ((int) px->fpclass < (int) py->fpclass) {
42 pt = px;
43 px = py;
44 py = pt;
45 }
46 /* Now class(x) >= class(y). */
47
48 *pz = *px;
49 pz->sign = px->sign ^ py->sign;
50
51 switch (px->fpclass) {
52 case fp_quiet:
53 case fp_signaling:
54 case fp_zero:
55 return;
56 case fp_infinity:
57 if (py->fpclass == fp_zero) {
58 fpu_error_nan(pz);
59 pz->fpclass = fp_quiet;
60 }
61 return;
62 case fp_normal:
63 if (py->fpclass == fp_zero) {
64 pz->fpclass = fp_zero;
65 return;
66 }
67 }
68
69 /* Now x and y are both normal or subnormal. */
70
71 x = px->significand; /* save typing */
72
73 s=r=acc[0]=acc[1]=acc[2]=acc[3]=0; /* intialize acc to zero */
74
75 y = py->significand[3]; /* py->significand[3] * x */
76 if(y!=0) {
77 j=1;
78 do {
79 s |= r; /* shift acc right one bit */
80 r = acc[3]&1;
81 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
82 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
83 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
84 acc[0] = (acc[0]>>1);
85 if(j&y) { /* bit i of y != 0, add x to acc */
86 c = 0;
87 c = fpu_add3wc(&acc[3],acc[3],x[3],c);
88 c = fpu_add3wc(&acc[2],acc[2],x[2],c);
89 c = fpu_add3wc(&acc[1],acc[1],x[1],c);
90 c = fpu_add3wc(&acc[0],acc[0],x[0],c);
91 }
92 j += j;
93 } while (j!=0);
94 }
95
96 y = py->significand[2]; /* py->significand[2] * x */
97 if(y!=0) {
98 j=1;
99 do {
100 s |= r; /* shift acc right one bit */
101 r = acc[3]&1;
102 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
103 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
104 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
105 acc[0] = (acc[0]>>1);
106 if(j&y) { /* bit i of y != 0, add x to acc */
107 c = 0;
108 c = fpu_add3wc(&acc[3],acc[3],x[3],c);
109 c = fpu_add3wc(&acc[2],acc[2],x[2],c);
110 c = fpu_add3wc(&acc[1],acc[1],x[1],c);
111 c = fpu_add3wc(&acc[0],acc[0],x[0],c);
112 }
113 j += j;
114 } while (j!=0);
115 } else {
116 s |= r|(acc[3]&0x7fffffff);
117 r = (acc[3]&0x80000000)>>31;
118 acc[3]=acc[2];acc[2]=acc[1];acc[1]=acc[0];acc[0]=0;
119 }
120
121 y = py->significand[1]; /* py->significand[1] * x */
122 if(y!=0) {
123 j=1;
124 do {
125 s |= r; /* shift acc right one bit */
126 r = acc[3]&1;
127 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
128 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
129 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
130 acc[0] = (acc[0]>>1);
131 if(j&y) { /* bit i of y != 0, add x to acc */
132 c = 0;
133 c = fpu_add3wc(&acc[3],acc[3],x[3],c);
134 c = fpu_add3wc(&acc[2],acc[2],x[2],c);
135 c = fpu_add3wc(&acc[1],acc[1],x[1],c);
136 c = fpu_add3wc(&acc[0],acc[0],x[0],c);
137 }
138 j += j;
139 } while (j!=0);
140 } else {
141 s |= r|(acc[3]&0x7fffffff);
142 r = (acc[3]&0x80000000)>>31;
143 acc[3]=acc[2];acc[2]=acc[1];acc[1]=acc[0];acc[0]=0;
144 }
145
146 /* py->significand[0] * x */
147 y = py->significand[0]; /* y is of form 0x0001???? */
148 j=1;
149 do {
150 s |= r; /* shift acc right one bit */
151 r = acc[3]&1;
152 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
153 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
154 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
155 acc[0] = (acc[0]>>1);
156 if(j&y) { /* bit i of y != 0, add x to acc */
157 c = 0;
158 c = fpu_add3wc(&acc[3],acc[3],x[3],c);
159 c = fpu_add3wc(&acc[2],acc[2],x[2],c);
160 c = fpu_add3wc(&acc[1],acc[1],x[1],c);
161 c = fpu_add3wc(&acc[0],acc[0],x[0],c);
162 }
163 j += j;
164 } while (j<=y);
165
166 if(acc[0]>=0x20000) { /* right shift one bit to normalize */
167 pz->exponent = px->exponent + py->exponent + 1;
168 pz->sticky = s|r;
169 pz->rounded = acc[3]&1;
170 pz->significand[3]=((acc[2]&1)<<31)|(acc[3]>>1);
171 pz->significand[2]=((acc[1]&1)<<31)|(acc[2]>>1);
172 pz->significand[1]=((acc[0]&1)<<31)|(acc[1]>>1);
173 pz->significand[0]=(acc[0]>>1);
174 } else {
175 pz->exponent = px->exponent + py->exponent;
176 pz->sticky = s;
177 pz->rounded = r;
178 pz->significand[3]=acc[3];
179 pz->significand[2]=acc[2];
180 pz->significand[1]=acc[1];
181 pz->significand[0]=acc[0];
182 }
183 }
184