xref: /netbsd-src/external/lgpl3/mpc/dist/src/dot.c (revision 90a8ff2142ed565a73c3c0859f0b1e7d216aeb7b)
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