1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
7Software Foundation; either version 3, or (at your option) any later
8version.
9
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13for more details.
14
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22<http://www.gnu.org/licenses/>.  */
23
24#include "bid_internal.h"
25
26#if DECIMAL_CALL_BY_REFERENCE
27void
28bid64dq_mul (UINT64 * pres, UINT64 * px, UINT128 * py
29	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
30	     _EXC_INFO_PARAM) {
31  UINT64 x = *px;
32#if !DECIMAL_GLOBAL_ROUNDING
33  unsigned int rnd_mode = *prnd_mode;
34#endif
35#else
36UINT64
37bid64dq_mul (UINT64 x, UINT128 y
38	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
39	     _EXC_INFO_PARAM) {
40#endif
41  UINT64 res = 0xbaddbaddbaddbaddull;
42  UINT128 x1;
43
44#if DECIMAL_CALL_BY_REFERENCE
45  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
46  bid64qq_mul (&res, &x1, py
47	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
48	       _EXC_INFO_ARG);
49#else
50  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
51  res = bid64qq_mul (x1, y
52		     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
53		     _EXC_INFO_ARG);
54#endif
55  BID_RETURN (res);
56}
57
58
59#if DECIMAL_CALL_BY_REFERENCE
60void
61bid64qd_mul (UINT64 * pres, UINT128 * px, UINT64 * py
62	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
63	     _EXC_INFO_PARAM) {
64  UINT64 y = *py;
65#if !DECIMAL_GLOBAL_ROUNDING
66  unsigned int rnd_mode = *prnd_mode;
67#endif
68#else
69UINT64
70bid64qd_mul (UINT128 x, UINT64 y
71	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
72	     _EXC_INFO_PARAM) {
73#endif
74  UINT64 res = 0xbaddbaddbaddbaddull;
75  UINT128 y1;
76
77#if DECIMAL_CALL_BY_REFERENCE
78  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
79  bid64qq_mul (&res, px, &y1
80	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
81	       _EXC_INFO_ARG);
82#else
83  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
84  res = bid64qq_mul (x, y1
85		     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
86		     _EXC_INFO_ARG);
87#endif
88  BID_RETURN (res);
89}
90
91
92#if DECIMAL_CALL_BY_REFERENCE
93void
94bid64qq_mul (UINT64 * pres, UINT128 * px, UINT128 * py
95	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
96	     _EXC_INFO_PARAM) {
97  UINT128 x = *px, y = *py;
98#if !DECIMAL_GLOBAL_ROUNDING
99  unsigned int rnd_mode = *prnd_mode;
100#endif
101#else
102UINT64
103bid64qq_mul (UINT128 x, UINT128 y
104	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
105	     _EXC_INFO_PARAM) {
106#endif
107
108  UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
109  };
110  UINT64 res = 0xbaddbaddbaddbaddull;
111  UINT64 x_sign, y_sign, p_sign;
112  UINT64 x_exp, y_exp, p_exp;
113  int true_p_exp;
114  UINT128 C1, C2;
115
116  BID_SWAP128 (z);
117  // skip cases where at least one operand is NaN or infinity
118  if (!(((x.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
119	((y.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
120	((x.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF) ||
121	((y.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF))) {
122    // x, y are 0 or f but not inf or NaN => unpack the arguments and check
123    // for non-canonical values
124
125    x_sign = x.w[HIGH_128W] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
126    C1.w[1] = x.w[HIGH_128W] & MASK_COEFF;
127    C1.w[0] = x.w[LOW_128W];
128    // check for non-canonical values - treated as zero
129    if ((x.w[HIGH_128W] & 0x6000000000000000ull) ==
130	0x6000000000000000ull) {
131      // G0_G1=11 => non-canonical
132      x_exp = (x.w[HIGH_128W] << 2) & MASK_EXP;	// biased and shifted left 49 bits
133      C1.w[1] = 0;	// significand high
134      C1.w[0] = 0;	// significand low
135    } else {	// G0_G1 != 11
136      x_exp = x.w[HIGH_128W] & MASK_EXP;	// biased and shifted left 49 bits
137      if (C1.w[1] > 0x0001ed09bead87c0ull ||
138	  (C1.w[1] == 0x0001ed09bead87c0ull &&
139	   C1.w[0] > 0x378d8e63ffffffffull)) {
140	// x is non-canonical if coefficient is larger than 10^34 -1
141	C1.w[1] = 0;
142	C1.w[0] = 0;
143      } else {	// canonical
144	;
145      }
146    }
147    y_sign = y.w[HIGH_128W] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
148    C2.w[1] = y.w[HIGH_128W] & MASK_COEFF;
149    C2.w[0] = y.w[LOW_128W];
150    // check for non-canonical values - treated as zero
151    if ((y.w[HIGH_128W] & 0x6000000000000000ull) ==
152	0x6000000000000000ull) {
153      // G0_G1=11 => non-canonical
154      y_exp = (y.w[HIGH_128W] << 2) & MASK_EXP;	// biased and shifted left 49 bits
155      C2.w[1] = 0;	// significand high
156      C2.w[0] = 0;	// significand low
157    } else {	// G0_G1 != 11
158      y_exp = y.w[HIGH_128W] & MASK_EXP;	// biased and shifted left 49 bits
159      if (C2.w[1] > 0x0001ed09bead87c0ull ||
160	  (C2.w[1] == 0x0001ed09bead87c0ull &&
161	   C2.w[0] > 0x378d8e63ffffffffull)) {
162	// y is non-canonical if coefficient is larger than 10^34 -1
163	C2.w[1] = 0;
164	C2.w[0] = 0;
165      } else {	// canonical
166	;
167      }
168    }
169    p_sign = x_sign ^ y_sign;	// sign of the product
170
171    true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
172    // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0
173    if (true_p_exp < -398)
174      p_exp = 0;	// cannot be less than EXP_MIN
175    else if (true_p_exp > 369)
176      p_exp = (UINT64) (369 + 398) << 53;	// cannot be more than EXP_MAX
177    else
178      p_exp = (UINT64) (true_p_exp + 398) << 53;
179
180    if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
181	(C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
182      // x = 0 or y = 0
183      // the result is 0
184      res = p_sign | p_exp;	// preferred exponent in [EXP_MIN, EXP_MAX]
185      BID_RETURN (res)
186    }	// else continue
187  }
188  // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
189#if DECIMAL_CALL_BY_REFERENCE
190  bid64qqq_fma (&res, &y, &x, &z
191		_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
192		_EXC_INFO_ARG);
193#else
194  res = bid64qqq_fma (y, x, z
195		      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
196		      _EXC_INFO_ARG);
197#endif
198  BID_RETURN (res);
199}
200
201
202#if DECIMAL_CALL_BY_REFERENCE
203void
204bid128dd_mul (UINT128 * pres, UINT64 * px, UINT64 * py
205	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
206	      _EXC_INFO_PARAM) {
207  UINT64 x = *px, y = *py;
208#if !DECIMAL_GLOBAL_ROUNDING
209  unsigned int rnd_mode = *prnd_mode;
210#endif
211#else
212UINT128
213bid128dd_mul (UINT64 x, UINT64 y
214	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
215	      _EXC_INFO_PARAM) {
216#endif
217  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
218  };
219  UINT128 x1, y1;
220
221#if DECIMAL_CALL_BY_REFERENCE
222  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
223  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
224  bid128_mul (&res, &x1, &y1
225	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
226	      _EXC_INFO_ARG);
227#else
228  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
229  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
230  res = bid128_mul (x1, y1
231		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
232		    _EXC_INFO_ARG);
233#endif
234  BID_RETURN (res);
235}
236
237
238#if DECIMAL_CALL_BY_REFERENCE
239void
240bid128dq_mul (UINT128 * pres, UINT64 * px, UINT128 * py
241	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
242	      _EXC_INFO_PARAM) {
243  UINT64 x = *px;
244#if !DECIMAL_GLOBAL_ROUNDING
245  unsigned int rnd_mode = *prnd_mode;
246#endif
247#else
248UINT128
249bid128dq_mul (UINT64 x, UINT128 y
250	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
251	      _EXC_INFO_PARAM) {
252#endif
253  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
254  };
255  UINT128 x1;
256
257#if DECIMAL_CALL_BY_REFERENCE
258  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
259  bid128_mul (&res, &x1, py
260	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
261	      _EXC_INFO_ARG);
262#else
263  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
264  res = bid128_mul (x1, y
265		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
266		    _EXC_INFO_ARG);
267#endif
268  BID_RETURN (res);
269}
270
271
272#if DECIMAL_CALL_BY_REFERENCE
273void
274bid128qd_mul (UINT128 * pres, UINT128 * px, UINT64 * py
275	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
276	      _EXC_INFO_PARAM) {
277  UINT64 y = *py;
278#if !DECIMAL_GLOBAL_ROUNDING
279  unsigned int rnd_mode = *prnd_mode;
280#endif
281#else
282UINT128
283bid128qd_mul (UINT128 x, UINT64 y
284	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
285	      _EXC_INFO_PARAM) {
286#endif
287  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
288  };
289  UINT128 y1;
290
291#if DECIMAL_CALL_BY_REFERENCE
292  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
293  bid128_mul (&res, px, &y1
294	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
295	      _EXC_INFO_ARG);
296#else
297  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
298  res = bid128_mul (x, y1
299		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
300		    _EXC_INFO_ARG);
301#endif
302  BID_RETURN (res);
303}
304
305
306// bid128_mul stands for bid128qq_mul
307#if DECIMAL_CALL_BY_REFERENCE
308void
309bid128_mul (UINT128 * pres, UINT128 * px,
310	    UINT128 *
311	    py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
312	    _EXC_INFO_PARAM) {
313  UINT128 x = *px, y = *py;
314
315#if !DECIMAL_GLOBAL_ROUNDING
316  unsigned int rnd_mode = *prnd_mode;
317
318#endif
319#else
320UINT128
321bid128_mul (UINT128 x,
322	    UINT128 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
323	    _EXC_INFO_PARAM) {
324
325#endif
326  UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
327  };
328  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
329  };
330  UINT64 x_sign, y_sign, p_sign;
331  UINT64 x_exp, y_exp, p_exp;
332  int true_p_exp;
333  UINT128 C1, C2;
334
335  BID_SWAP128 (x);
336  BID_SWAP128 (y);
337  // skip cases where at least one operand is NaN or infinity
338  if (!(((x.w[1] & MASK_NAN) == MASK_NAN) ||
339	((y.w[1] & MASK_NAN) == MASK_NAN) ||
340	((x.w[1] & MASK_ANY_INF) == MASK_INF) ||
341	((y.w[1] & MASK_ANY_INF) == MASK_INF))) {
342    // x, y are 0 or f but not inf or NaN => unpack the arguments and check
343    // for non-canonical values
344
345    x_sign = x.w[1] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
346    C1.w[1] = x.w[1] & MASK_COEFF;
347    C1.w[0] = x.w[0];
348    // check for non-canonical values - treated as zero
349    if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
350      // G0_G1=11 => non-canonical
351      x_exp = (x.w[1] << 2) & MASK_EXP;	// biased and shifted left 49 bits
352      C1.w[1] = 0;	// significand high
353      C1.w[0] = 0;	// significand low
354    } else {	// G0_G1 != 11
355      x_exp = x.w[1] & MASK_EXP;	// biased and shifted left 49 bits
356      if (C1.w[1] > 0x0001ed09bead87c0ull ||
357	  (C1.w[1] == 0x0001ed09bead87c0ull &&
358	   C1.w[0] > 0x378d8e63ffffffffull)) {
359	// x is non-canonical if coefficient is larger than 10^34 -1
360	C1.w[1] = 0;
361	C1.w[0] = 0;
362      } else {	// canonical
363	;
364      }
365    }
366    y_sign = y.w[1] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
367    C2.w[1] = y.w[1] & MASK_COEFF;
368    C2.w[0] = y.w[0];
369    // check for non-canonical values - treated as zero
370    if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
371      // G0_G1=11 => non-canonical
372      y_exp = (y.w[1] << 2) & MASK_EXP;	// biased and shifted left 49 bits
373      C2.w[1] = 0;	// significand high
374      C2.w[0] = 0;	// significand low
375    } else {	// G0_G1 != 11
376      y_exp = y.w[1] & MASK_EXP;	// biased and shifted left 49 bits
377      if (C2.w[1] > 0x0001ed09bead87c0ull ||
378	  (C2.w[1] == 0x0001ed09bead87c0ull &&
379	   C2.w[0] > 0x378d8e63ffffffffull)) {
380	// y is non-canonical if coefficient is larger than 10^34 -1
381	C2.w[1] = 0;
382	C2.w[0] = 0;
383      } else {	// canonical
384	;
385      }
386    }
387    p_sign = x_sign ^ y_sign;	// sign of the product
388
389    true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
390    // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0
391    if (true_p_exp < -6176)
392      p_exp = 0;	// cannot be less than EXP_MIN
393    else if (true_p_exp > 6111)
394      p_exp = (UINT64) (6111 + 6176) << 49;	// cannot be more than EXP_MAX
395    else
396      p_exp = (UINT64) (true_p_exp + 6176) << 49;
397
398    if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
399	(C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
400      // x = 0 or y = 0
401      // the result is 0
402      res.w[1] = p_sign | p_exp;	// preferred exponent in [EXP_MIN, EXP_MAX]
403      res.w[0] = 0x0;
404      BID_SWAP128 (res);
405      BID_RETURN (res)
406    }	// else continue
407  }
408
409  BID_SWAP128 (x);
410  BID_SWAP128 (y);
411  BID_SWAP128 (z);
412  // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
413#if DECIMAL_CALL_BY_REFERENCE
414  bid128_fma (&res, &y, &x, &z
415	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
416	      _EXC_INFO_ARG);
417#else
418  res = bid128_fma (y, x, z
419		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
420		    _EXC_INFO_ARG);
421#endif
422  BID_RETURN (res);
423}
424