1/* mpfr_dot -- dot product of two array of numbers
2
3Copyright 2018-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include "mpfr-impl.h"
24
25/* FIXME: handle intermediate overflows and underflows. */
26
27/* res <- a[0]*b[0] + ... + a[n-1]*b[n-1] */
28int
29mpfr_dot (mpfr_ptr res, const mpfr_ptr *a, const mpfr_ptr *b,
30          unsigned long n, mpfr_rnd_t rnd)
31{
32  mpfr_t *c;
33  size_t i;
34  int inex;
35  mpfr_ptr *tab;
36
37  if (MPFR_UNLIKELY (n == 0))
38    {
39      MPFR_SET_ZERO (res);
40      MPFR_SET_POS (res);
41      MPFR_RET (0);
42    }
43
44  c = (mpfr_t *) mpfr_allocate_func (n * sizeof (mpfr_t));
45  tab = (mpfr_ptr *) mpfr_allocate_func (n * sizeof (mpfr_ptr));
46  for (i = 0; i < n; i++)
47    {
48      mpfr_init2 (c[i], mpfr_get_prec (a[i]) + mpfr_get_prec (b[i]));
49      inex = mpfr_mul (c[i], a[i], b[i], MPFR_RNDZ); /* exact except... */
50      MPFR_ASSERTN (inex == 0); /* failure in case of overflow/underflow */
51      tab[i] = c[i];
52    }
53  inex = mpfr_sum (res, tab, n, rnd);
54  for (i = 0; i < n; i++)
55    mpfr_clear (c[i]);
56  mpfr_free_func (c, n * sizeof (mpfr_t));
57  mpfr_free_func (tab, n * sizeof (mpfr_ptr));
58  return inex;
59}
60