1/* flonum_multip.c - multiply two flonums
2   Copyright (C) 1987 Free Software Foundation, Inc.
3
4This file is part of Gas, the GNU Assembler.
5
6The GNU assembler is distributed in the hope that it will be
7useful, but WITHOUT ANY WARRANTY.  No author or distributor
8accepts responsibility to anyone for the consequences of using it
9or for whether it serves any particular purpose or works at all,
10unless he says so in writing.  Refer to the GNU Assembler General
11Public License for full details.
12
13Everyone is granted permission to copy, modify and redistribute
14the GNU Assembler, but only under the conditions described in the
15GNU Assembler General Public License.  A copy of this license is
16supposed to have been given to you along with the GNU Assembler
17so you can know your rights and responsibilities.  It should be
18in a file named COPYING.  Among other things, the copyright
19notice 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
64a[], b[], p[] may not intersect.
65Zero length factors signify 0 significant bits: treat as 0.0.
660.0 factors do the right thing.
67Zero length product OK.
68
69I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
70because I felt the ForTran way was more intuitive. The C way would
71probably yield better code on most C compilers. Dean Elsner.
72(C style also gives deeper insight [to me] ... oh well ...)
73*/
74
75void
76flonum_multip(
77const_FLONUM_TYPE *a,
78FLONUM_TYPE *b,
79FLONUM_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  uint32_t		work;
87  uint32_t		carry;
88  int32_t		exponent;
89  LITTLENUM_TYPE *	q;
90  int32_t		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
135printf("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
141printf("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
152printf("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
168printf("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 */
201