1*a9fa9459Szrj /* flonum_mult.c - multiply two flonums
2*a9fa9459Szrj Copyright (C) 1987-2016 Free Software Foundation, Inc.
3*a9fa9459Szrj
4*a9fa9459Szrj This file is part of GAS, the GNU Assembler.
5*a9fa9459Szrj
6*a9fa9459Szrj GAS is free software; you can redistribute it and/or modify
7*a9fa9459Szrj it under the terms of the GNU General Public License as published by
8*a9fa9459Szrj the Free Software Foundation; either version 3, or (at your option)
9*a9fa9459Szrj any later version.
10*a9fa9459Szrj
11*a9fa9459Szrj GAS is distributed in the hope that it will be useful, but WITHOUT
12*a9fa9459Szrj ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13*a9fa9459Szrj or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
14*a9fa9459Szrj License for more details.
15*a9fa9459Szrj
16*a9fa9459Szrj You should have received a copy of the GNU General Public License
17*a9fa9459Szrj along with GAS; see the file COPYING. If not, write to the Free
18*a9fa9459Szrj Software Foundation, 51 Franklin Street - Fifth Floor, Boston, MA
19*a9fa9459Szrj 02110-1301, USA. */
20*a9fa9459Szrj
21*a9fa9459Szrj #include "ansidecl.h"
22*a9fa9459Szrj #include "flonum.h"
23*a9fa9459Szrj
24*a9fa9459Szrj /* plan for a . b => p(roduct)
25*a9fa9459Szrj
26*a9fa9459Szrj +-------+-------+-/ /-+-------+-------+
27*a9fa9459Szrj | a | a | ... | a | a |
28*a9fa9459Szrj | A | A-1 | | 1 | 0 |
29*a9fa9459Szrj +-------+-------+-/ /-+-------+-------+
30*a9fa9459Szrj
31*a9fa9459Szrj +-------+-------+-/ /-+-------+-------+
32*a9fa9459Szrj | b | b | ... | b | b |
33*a9fa9459Szrj | B | B-1 | | 1 | 0 |
34*a9fa9459Szrj +-------+-------+-/ /-+-------+-------+
35*a9fa9459Szrj
36*a9fa9459Szrj +-------+-------+-/ /-+-------+-/ /-+-------+-------+
37*a9fa9459Szrj | p | p | ... | p | ... | p | p |
38*a9fa9459Szrj | A+B+1| A+B | | N | | 1 | 0 |
39*a9fa9459Szrj +-------+-------+-/ /-+-------+-/ /-+-------+-------+
40*a9fa9459Szrj
41*a9fa9459Szrj /^\
42*a9fa9459Szrj (carry) a .b ... | ... a .b a .b
43*a9fa9459Szrj A B | 0 1 0 0
44*a9fa9459Szrj |
45*a9fa9459Szrj ... | ... a .b
46*a9fa9459Szrj | 1 0
47*a9fa9459Szrj |
48*a9fa9459Szrj | ...
49*a9fa9459Szrj |
50*a9fa9459Szrj |
51*a9fa9459Szrj |
52*a9fa9459Szrj | ___
53*a9fa9459Szrj | \
54*a9fa9459Szrj +----- P = > a .b
55*a9fa9459Szrj N /__ i j
56*a9fa9459Szrj
57*a9fa9459Szrj N = 0 ... A+B
58*a9fa9459Szrj
59*a9fa9459Szrj for all i,j where i+j=N
60*a9fa9459Szrj [i,j integers > 0]
61*a9fa9459Szrj
62*a9fa9459Szrj a[], b[], p[] may not intersect.
63*a9fa9459Szrj Zero length factors signify 0 significant bits: treat as 0.0.
64*a9fa9459Szrj 0.0 factors do the right thing.
65*a9fa9459Szrj Zero length product OK.
66*a9fa9459Szrj
67*a9fa9459Szrj I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
68*a9fa9459Szrj because I felt the ForTran way was more intuitive. The C way would
69*a9fa9459Szrj probably yield better code on most C compilers. Dean Elsner.
70*a9fa9459Szrj (C style also gives deeper insight [to me] ... oh well ...) */
71*a9fa9459Szrj
72*a9fa9459Szrj void
flonum_multip(const FLONUM_TYPE * a,const FLONUM_TYPE * b,FLONUM_TYPE * product)73*a9fa9459Szrj flonum_multip (const FLONUM_TYPE *a, const FLONUM_TYPE *b,
74*a9fa9459Szrj FLONUM_TYPE *product)
75*a9fa9459Szrj {
76*a9fa9459Szrj int size_of_a; /* 0 origin */
77*a9fa9459Szrj int size_of_b; /* 0 origin */
78*a9fa9459Szrj int size_of_product; /* 0 origin */
79*a9fa9459Szrj int size_of_sum; /* 0 origin */
80*a9fa9459Szrj int extra_product_positions; /* 1 origin */
81*a9fa9459Szrj unsigned long work;
82*a9fa9459Szrj unsigned long carry;
83*a9fa9459Szrj long exponent;
84*a9fa9459Szrj LITTLENUM_TYPE *q;
85*a9fa9459Szrj long significant; /* TRUE when we emit a non-0 littlenum */
86*a9fa9459Szrj /* ForTran accent follows. */
87*a9fa9459Szrj int P; /* Scan product low-order -> high. */
88*a9fa9459Szrj int N; /* As in sum above. */
89*a9fa9459Szrj int A; /* Which [] of a? */
90*a9fa9459Szrj int B; /* Which [] of b? */
91*a9fa9459Szrj
92*a9fa9459Szrj if ((a->sign != '-' && a->sign != '+')
93*a9fa9459Szrj || (b->sign != '-' && b->sign != '+'))
94*a9fa9459Szrj {
95*a9fa9459Szrj /* Got to fail somehow. Any suggestions? */
96*a9fa9459Szrj product->sign = 0;
97*a9fa9459Szrj return;
98*a9fa9459Szrj }
99*a9fa9459Szrj product->sign = (a->sign == b->sign) ? '+' : '-';
100*a9fa9459Szrj size_of_a = a->leader - a->low;
101*a9fa9459Szrj size_of_b = b->leader - b->low;
102*a9fa9459Szrj exponent = a->exponent + b->exponent;
103*a9fa9459Szrj size_of_product = product->high - product->low;
104*a9fa9459Szrj size_of_sum = size_of_a + size_of_b;
105*a9fa9459Szrj extra_product_positions = size_of_product - size_of_sum;
106*a9fa9459Szrj if (extra_product_positions < 0)
107*a9fa9459Szrj {
108*a9fa9459Szrj P = extra_product_positions; /* P < 0 */
109*a9fa9459Szrj exponent -= extra_product_positions; /* Increases exponent. */
110*a9fa9459Szrj }
111*a9fa9459Szrj else
112*a9fa9459Szrj {
113*a9fa9459Szrj P = 0;
114*a9fa9459Szrj }
115*a9fa9459Szrj carry = 0;
116*a9fa9459Szrj significant = 0;
117*a9fa9459Szrj for (N = 0; N <= size_of_sum; N++)
118*a9fa9459Szrj {
119*a9fa9459Szrj work = carry;
120*a9fa9459Szrj carry = 0;
121*a9fa9459Szrj for (A = 0; A <= N; A++)
122*a9fa9459Szrj {
123*a9fa9459Szrj B = N - A;
124*a9fa9459Szrj if (A <= size_of_a && B <= size_of_b && B >= 0)
125*a9fa9459Szrj {
126*a9fa9459Szrj #ifdef TRACE
127*a9fa9459Szrj printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n",
128*a9fa9459Szrj A, a->low[A], B, b->low[B], work);
129*a9fa9459Szrj #endif
130*a9fa9459Szrj /* Watch out for sign extension! Without the casts, on
131*a9fa9459Szrj the DEC Alpha, the multiplication result is *signed*
132*a9fa9459Szrj int, which gets sign-extended to convert to the
133*a9fa9459Szrj unsigned long! */
134*a9fa9459Szrj work += (unsigned long) a->low[A] * (unsigned long) b->low[B];
135*a9fa9459Szrj carry += work >> LITTLENUM_NUMBER_OF_BITS;
136*a9fa9459Szrj work &= LITTLENUM_MASK;
137*a9fa9459Szrj #ifdef TRACE
138*a9fa9459Szrj printf ("work=%08x carry=%04x\n", work, carry);
139*a9fa9459Szrj #endif
140*a9fa9459Szrj }
141*a9fa9459Szrj }
142*a9fa9459Szrj significant |= work;
143*a9fa9459Szrj if (significant || P < 0)
144*a9fa9459Szrj {
145*a9fa9459Szrj if (P >= 0)
146*a9fa9459Szrj {
147*a9fa9459Szrj product->low[P] = work;
148*a9fa9459Szrj #ifdef TRACE
149*a9fa9459Szrj printf ("P=%d. work[p]:=%04x\n", P, work);
150*a9fa9459Szrj #endif
151*a9fa9459Szrj }
152*a9fa9459Szrj P++;
153*a9fa9459Szrj }
154*a9fa9459Szrj else
155*a9fa9459Szrj {
156*a9fa9459Szrj extra_product_positions++;
157*a9fa9459Szrj exponent++;
158*a9fa9459Szrj }
159*a9fa9459Szrj }
160*a9fa9459Szrj /* [P]-> position # size_of_sum + 1.
161*a9fa9459Szrj This is where 'carry' should go. */
162*a9fa9459Szrj #ifdef TRACE
163*a9fa9459Szrj printf ("final carry =%04x\n", carry);
164*a9fa9459Szrj #endif
165*a9fa9459Szrj if (carry)
166*a9fa9459Szrj {
167*a9fa9459Szrj if (extra_product_positions > 0)
168*a9fa9459Szrj product->low[P] = carry;
169*a9fa9459Szrj else
170*a9fa9459Szrj {
171*a9fa9459Szrj /* No room at high order for carry littlenum. */
172*a9fa9459Szrj /* Shift right 1 to make room for most significant littlenum. */
173*a9fa9459Szrj exponent++;
174*a9fa9459Szrj P--;
175*a9fa9459Szrj for (q = product->low + P; q >= product->low; q--)
176*a9fa9459Szrj {
177*a9fa9459Szrj work = *q;
178*a9fa9459Szrj *q = carry;
179*a9fa9459Szrj carry = work;
180*a9fa9459Szrj }
181*a9fa9459Szrj }
182*a9fa9459Szrj }
183*a9fa9459Szrj else
184*a9fa9459Szrj P--;
185*a9fa9459Szrj product->leader = product->low + P;
186*a9fa9459Szrj product->exponent = exponent;
187*a9fa9459Szrj }
188