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 static void
arf_shallow_set_siui(arf_t res,ulong vhi,ulong vlo)15 arf_shallow_set_siui(arf_t res, ulong vhi, ulong vlo)
16 {
17 int negative;
18 unsigned int bc;
19
20 negative = ((slong) vhi) < 0;
21
22 if (negative)
23 {
24 vhi = -vhi - (vlo != 0);
25 vlo = -vlo;
26 }
27
28 if (vhi == 0)
29 {
30 if (vlo == 0)
31 {
32 ARF_XSIZE(res) = 0;
33 ARF_EXP(res) = ARF_EXP_ZERO;
34 }
35 else
36 {
37 count_leading_zeros(bc, vlo);
38 ARF_EXP(res) = FLINT_BITS - bc;
39 ARF_NOPTR_D(res)[0] = vlo << bc;
40 ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
41 }
42 }
43 else if (vlo == 0)
44 {
45 count_leading_zeros(bc, vhi);
46 ARF_EXP(res) = 2 * FLINT_BITS - bc;
47 ARF_NOPTR_D(res)[0] = vhi << bc;
48 ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
49 }
50 else
51 {
52 count_leading_zeros(bc, vhi);
53 ARF_EXP(res) = 2 * FLINT_BITS - bc;
54 ARF_NOPTR_D(res)[0] = vlo << bc;
55 if (bc == 0)
56 ARF_NOPTR_D(res)[1] = vhi;
57 else
58 ARF_NOPTR_D(res)[1] = (vhi << bc) | (vlo >> (FLINT_BITS - bc));
59 ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, negative);
60 }
61 }
62
63 void
arb_dot_siui(arb_t res,const arb_t initial,int subtract,arb_srcptr x,slong xstep,const ulong * y,slong ystep,slong len,slong prec)64 arb_dot_siui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
65 {
66 arb_ptr t;
67 slong i;
68 ulong vhi, vlo;
69 TMP_INIT;
70
71 /* todo: fast fma and fmma (len=2) code */
72 if (len <= 1)
73 {
74 if (initial == NULL)
75 {
76 if (len <= 0)
77 arb_zero(res);
78 else
79 {
80 arf_t t;
81 arf_shallow_set_siui(t, y[1], y[0]);
82 arb_mul_arf(res, x, t, prec);
83 if (subtract)
84 arb_neg(res, res);
85 }
86 return;
87 }
88 else if (len <= 0)
89 {
90 arb_set_round(res, initial, prec);
91 return;
92 }
93 }
94
95 TMP_START;
96 t = TMP_ALLOC(sizeof(arb_struct) * len);
97
98 for (i = 0; i < len; i++)
99 {
100 vlo = y[2 * i * ystep];
101 vhi = y[2 * i * ystep + 1];
102
103 arf_shallow_set_siui(arb_midref(t + i), vhi, vlo);
104
105 MAG_EXP(arb_radref(t + i)) = 0;
106 MAG_MAN(arb_radref(t + i)) = 0;
107 }
108
109 arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
110
111 TMP_END;
112 }
113
114