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 "arb.h"
13
14 void
arb_dot_ui(arb_t res,const arb_t initial,int subtract,arb_srcptr x,slong xstep,const ulong * y,slong ystep,slong len,slong prec)15 arb_dot_ui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
16 {
17 arb_ptr t;
18 slong i;
19 ulong v;
20 unsigned int bc;
21 TMP_INIT;
22
23 /* todo: fast fma and fmma (len=2) code */
24 if (len <= 1)
25 {
26 if (initial == NULL)
27 {
28 if (len <= 0)
29 arb_zero(res);
30 else
31 {
32 arb_mul_ui(res, x, y[0], prec);
33 if (subtract)
34 arb_neg(res, res);
35 }
36 return;
37 }
38 else if (len <= 0)
39 {
40 arb_set_round(res, initial, prec);
41 return;
42 }
43 }
44
45 TMP_START;
46 t = TMP_ALLOC(sizeof(arb_struct) * len);
47
48 for (i = 0; i < len; i++)
49 {
50 v = y[i * ystep];
51
52 if (v == 0)
53 {
54 ARF_XSIZE(arb_midref(t + i)) = 0;
55 ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
56 }
57 else
58 {
59 count_leading_zeros(bc, v);
60
61 ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
62 ARF_NOPTR_D(arb_midref(t + i))[0] = v << bc;
63 ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, 0);
64 }
65
66 MAG_EXP(arb_radref(t + i)) = 0;
67 MAG_MAN(arb_radref(t + i)) = 0;
68 }
69
70 arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
71
72 TMP_END;
73 }
74
75