|
|
1.1 root 1: /* flonum_multip.c - multiply two flonums
2: Copyright (C) 1987 Free Software Foundation, Inc.
3:
4: This file is part of Gas, the GNU Assembler.
5:
6: The GNU assembler is distributed in the hope that it will be
7: useful, but WITHOUT ANY WARRANTY. No author or distributor
8: accepts responsibility to anyone for the consequences of using it
9: or for whether it serves any particular purpose or works at all,
10: unless he says so in writing. Refer to the GNU Assembler General
11: Public License for full details.
12:
13: Everyone is granted permission to copy, modify and redistribute
14: the GNU Assembler, but only under the conditions described in the
15: GNU Assembler General Public License. A copy of this license is
16: supposed to have been given to you along with the GNU Assembler
17: so you can know your rights and responsibilities. It should be
18: in a file named COPYING. Among other things, the copyright
19: notice and this notice must be preserved on all copies. */
20:
21: #include "flonum.h"
22:
23: /* plan for a . b => p(roduct)
24:
25:
26: +-------+-------+-/ /-+-------+-------+
27: | a | a | ... | a | a |
28: | A | A-1 | | 1 | 0 |
29: +-------+-------+-/ /-+-------+-------+
30:
31:
32: +-------+-------+-/ /-+-------+-------+
33: | b | b | ... | b | b |
34: | B | B-1 | | 1 | 0 |
35: +-------+-------+-/ /-+-------+-------+
36:
37:
38: +-------+-------+-/ /-+-------+-/ /-+-------+-------+
39: | p | p | ... | p | ... | p | p |
40: | A+B+1| A+B | | N | | 1 | 0 |
41: +-------+-------+-/ /-+-------+-/ /-+-------+-------+
42:
43: /^\
44: (carry) a .b ... | ... a .b a .b
45: A B | 0 1 0 0
46: |
47: ... | ... a .b
48: | 1 0
49: |
50: | ...
51: |
52: |
53: |
54: | ___
55: | \
56: +----- P = > a .b
57: N /__ i j
58:
59: N = 0 ... A+B
60:
61: for all i,j where i+j=N
62: [i,j integers > 0]
63:
64: a[], b[], p[] may not intersect.
65: Zero length factors signify 0 significant bits: treat as 0.0.
66: 0.0 factors do the right thing.
67: Zero length product OK.
68:
69: I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
70: because I felt the ForTran way was more intuitive. The C way would
71: probably yield better code on most C compilers. Dean Elsner.
72: (C style also gives deeper insight [to me] ... oh well ...)
73: */
74:
75: void
76: flonum_multip(
77: const_FLONUM_TYPE *a,
78: FLONUM_TYPE *b,
79: FLONUM_TYPE *product)
80: {
81: int size_of_a; /* 0 origin */
82: int size_of_b; /* 0 origin */
83: int size_of_product; /* 0 origin */
84: int size_of_sum; /* 0 origin */
85: int extra_product_positions;/* 1 origin */
86: unsigned long int work;
87: unsigned long int carry;
88: long int exponent;
89: LITTLENUM_TYPE * q;
90: long int significant; /* TRUE when we emit a non-0 littlenum */
91: /* ForTran accent follows. */
92: int P; /* Scan product low-order -> high. */
93: int N; /* As in sum above. */
94: int A; /* Which [] of a? */
95: int B; /* Which [] of b? */
96:
97: if((a->sign!='-' && a->sign!='+') || (b->sign!='-' && b->sign!='+')) {
98: /* ...
99: Got to fail somehow. Any suggestions? */
100: product->sign=0;
101: return;
102: }
103: product -> sign = (a->sign == b->sign) ? '+' : '-';
104: size_of_a = a -> leader - a -> low;
105: size_of_b = b -> leader - b -> low;
106: exponent = a -> exponent + b -> exponent;
107: size_of_product = product -> high - product -> low;
108: size_of_sum = size_of_a + size_of_b;
109: extra_product_positions = size_of_product - size_of_sum;
110: if (extra_product_positions < 0)
111: {
112: P = extra_product_positions; /* P < 0 */
113: exponent -= extra_product_positions; /* Increases exponent. */
114: }
115: else
116: {
117: P = 0;
118: }
119: carry = 0;
120: significant = 0;
121: for (N = 0;
122: N <= size_of_sum;
123: N++)
124: {
125: work = carry;
126: carry = 0;
127: for (A = 0;
128: A <= N;
129: A ++)
130: {
131: B = N - A;
132: if (A <= size_of_a && B <= size_of_b && B >= 0)
133: {
134: #ifdef TRACE
135: printf("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n", A, a->low[A], B, b->low[B], work);
136: #endif
137: work += a -> low [A] * b -> low [B];
138: carry += work >> LITTLENUM_NUMBER_OF_BITS;
139: work &= LITTLENUM_MASK;
140: #ifdef TRACE
141: printf("work=%08x carry=%04x\n", work, carry);
142: #endif
143: }
144: }
145: significant |= work;
146: if (significant || P<0)
147: {
148: if (P >= 0)
149: {
150: product -> low [P] = work;
151: #ifdef TRACE
152: printf("P=%d. work[p]:=%04x\n", P, work);
153: #endif
154: }
155: P ++;
156: }
157: else
158: {
159: extra_product_positions ++;
160: exponent ++;
161: }
162: }
163: /*
164: * [P]-> position # size_of_sum + 1.
165: * This is where 'carry' should go.
166: */
167: #ifdef TRACE
168: printf("final carry =%04x\n", carry);
169: #endif
170: if (carry)
171: {
172: if (extra_product_positions > 0)
173: {
174: product -> low [P] = carry;
175: }
176: else
177: {
178: /* No room at high order for carry littlenum. */
179: /* Shift right 1 to make room for most significant littlenum. */
180: exponent ++;
181: P --;
182: for (q = product -> low + P;
183: q >= product -> low;
184: q --)
185: {
186: work = * q;
187: * q = carry;
188: carry = work;
189: }
190: }
191: }
192: else
193: {
194: P --;
195: }
196: product -> leader = product -> low + P;
197: product -> exponent = exponent;
198: }
199:
200: /* end: flonum_multip.c */
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.