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 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
Indexes created Sun Sep 13 00:25:22 UTC 2026