1 /* mpc_dot -- Dot product of two arrays of complex numbers.
2
3 Copyright (C) 2018, 2020 INRIA
4
5 This file is part of GNU MPC.
6
7 GNU MPC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with this program. If not, see http://www.gnu.org/licenses/ .
19 */
20
21 #include <stdio.h> /* for MPC_ASSERT */
22 #include "mpc-impl.h"
23
24 /* res <- x[0]*y[0] + ... + x[n-1]*y[n-1] */
25 int
mpc_dot(mpc_ptr res,const mpc_ptr * x,const mpc_ptr * y,unsigned long n,mpc_rnd_t rnd)26 mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
27 unsigned long n, mpc_rnd_t rnd)
28 {
29 int inex_re, inex_im;
30 mpfr_ptr *t;
31 mpfr_t *z;
32 unsigned long i;
33 mpfr_t re_res;
34
35 z = (mpfr_t *) malloc (2 * n * sizeof (mpfr_t));
36 /* warning: when n=0, malloc() might return NULL (e.g., gcc119) */
37 MPC_ASSERT(n == 0 || z != NULL);
38 t = (mpfr_ptr *) malloc (2 * n * sizeof(mpfr_ptr));
39 MPC_ASSERT(n == 0 || t != NULL);
40 for (i = 0; i < 2 * n; i++)
41 t[i] = z[i];
42 /* we first store in z[i] the value of Re(x[i])*Re(y[i])
43 and in z[n+i] that of -Im(x[i])*Im(y[i]) */
44 for (i = 0; i < n; i++)
45 {
46 mpfr_prec_t prec_x_re = mpfr_get_prec (mpc_realref (x[i]));
47 mpfr_prec_t prec_x_im = mpfr_get_prec (mpc_imagref (x[i]));
48 mpfr_prec_t prec_y_re = mpfr_get_prec (mpc_realref (y[i]));
49 mpfr_prec_t prec_y_im = mpfr_get_prec (mpc_imagref (y[i]));
50 mpfr_prec_t prec_y_max = MPC_MAX (prec_y_re, prec_y_im);
51 /* we allocate z[i] with prec_x_re + prec_y_max bits
52 so that the second loop below does not reallocate */
53 mpfr_init2 (z[i], prec_x_re + prec_y_max);
54 mpfr_set_prec (z[i], prec_x_re + prec_y_re);
55 mpfr_mul (z[i], mpc_realref (x[i]), mpc_realref (y[i]), MPFR_RNDZ);
56 /* idem for z[n+i]: we allocate with prec_x_im + prec_y_max bits */
57 mpfr_init2 (z[n+i], prec_x_im + prec_y_max);
58 mpfr_set_prec (z[n+i], prec_x_im + prec_y_im);
59 mpfr_mul (z[n+i], mpc_imagref (x[i]), mpc_imagref (y[i]), MPFR_RNDZ);
60 mpfr_neg (z[n+i], z[n+i], MPFR_RNDZ);
61 }
62 /* copy the real part in a temporary variable, since it might be in the
63 input array */
64 mpfr_init2 (re_res, mpfr_get_prec (mpc_realref (res)));
65 inex_re = mpfr_sum (re_res, t, 2 * n, MPC_RND_RE (rnd));
66 /* we then store in z[i] the value of Re(x[i])*Im(y[i])
67 and in z[n+i] that of Im(x[i])*Re(y[i]) */
68 for (i = 0; i < n; i++)
69 {
70 mpfr_prec_t prec_x_re = mpfr_get_prec (mpc_realref (x[i]));
71 mpfr_prec_t prec_x_im = mpfr_get_prec (mpc_imagref (x[i]));
72 mpfr_prec_t prec_y_re = mpfr_get_prec (mpc_realref (y[i]));
73 mpfr_prec_t prec_y_im = mpfr_get_prec (mpc_imagref (y[i]));
74 mpfr_set_prec (z[i], prec_x_re + prec_y_im);
75 mpfr_mul (z[i], mpc_realref (x[i]), mpc_imagref (y[i]), MPFR_RNDZ);
76 mpfr_set_prec (z[n+i], prec_x_im + prec_y_re);
77 mpfr_mul (z[n+i], mpc_imagref (x[i]), mpc_realref (y[i]), MPFR_RNDZ);
78 }
79 inex_im = mpfr_sum (mpc_imagref (res), t, 2 * n, MPC_RND_IM (rnd));
80 mpfr_swap (mpc_realref (res), re_res);
81 mpfr_clear (re_res);
82 for (i = 0; i < 2 * n; i++)
83 mpfr_clear (z[i]);
84 free (t);
85 free (z);
86
87 return MPC_INEX(inex_re, inex_im);
88 }
89