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