1 /*
2 Copyright (C) 2021 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb.h"
13
14 void
acb_dot_si(acb_t res,const acb_t initial,int subtract,acb_srcptr x,slong xstep,const slong * y,slong ystep,slong len,slong prec)15 acb_dot_si(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec)
16 {
17 arb_ptr t;
18 slong i;
19 slong v;
20 ulong av;
21 unsigned int bc;
22 TMP_INIT;
23
24 /* todo: fast fma and fmma (len=2) code */
25 if (len <= 1)
26 {
27 if (initial == NULL)
28 {
29 if (len <= 0)
30 acb_zero(res);
31 else
32 {
33 acb_mul_si(res, x, y[0], prec);
34 if (subtract)
35 acb_neg(res, res);
36 }
37 return;
38 }
39 else if (len <= 0)
40 {
41 acb_set_round(res, initial, prec);
42 return;
43 }
44 }
45
46 TMP_START;
47 t = TMP_ALLOC(sizeof(arb_struct) * len);
48
49 for (i = 0; i < len; i++)
50 {
51 v = y[i * ystep];
52
53 if (v == 0)
54 {
55 ARF_XSIZE(arb_midref(t + i)) = 0;
56 ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
57 }
58 else
59 {
60 av = FLINT_ABS(v);
61 count_leading_zeros(bc, av);
62
63 ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
64 ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
65 ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0);
66 }
67
68 MAG_EXP(arb_radref(t + i)) = 0;
69 MAG_MAN(arb_radref(t + i)) = 0;
70 }
71
72 arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
73 arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
74
75 TMP_END;
76 }
77
78