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