1 /*
2 Copyright (C) 2012 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 #define a acb_realref(x)
15 #define b acb_imagref(x)
16 #define c acb_realref(y)
17 #define d acb_imagref(y)
18 #define e acb_realref(z)
19 #define f acb_imagref(z)
20 #define ar arb_radref(a)
21 #define br arb_radref(b)
22 #define cr arb_radref(c)
23 #define dr arb_radref(d)
24
25 static void
_acb_sqr_fast(acb_t z,const acb_t x,slong prec)26 _acb_sqr_fast(acb_t z, const acb_t x, slong prec)
27 {
28 int inexact;
29
30 mag_t am, bm, er, fr;
31
32 mag_fast_init_set_arf(am, arb_midref(a));
33 mag_fast_init_set_arf(bm, arb_midref(b));
34
35 mag_init(er);
36 mag_init(fr);
37
38 mag_fast_addmul(er, am, ar);
39 mag_fast_addmul(er, bm, br);
40 mag_fast_mul_2exp_si(er, er, 1);
41 mag_fast_addmul(er, ar, ar);
42 mag_fast_addmul(er, br, br);
43
44 mag_fast_addmul(fr, bm, ar);
45 mag_fast_addmul(fr, am, br);
46 mag_fast_addmul(fr, ar, br);
47 mag_fast_mul_2exp_si(fr, fr, 1);
48
49 inexact = arf_complex_sqr(arb_midref(e), arb_midref(f),
50 arb_midref(a), arb_midref(b), prec, ARB_RND);
51
52 if (inexact & 1)
53 arf_mag_add_ulp(arb_radref(e), er, arb_midref(e), prec);
54 else
55 mag_set(arb_radref(e), er);
56
57 if (inexact & 2)
58 arf_mag_add_ulp(arb_radref(f), fr, arb_midref(f), prec);
59 else
60 mag_set(arb_radref(f), fr);
61 }
62
63 static void
_acb_sqr_slow(acb_t z,const acb_t x,slong prec)64 _acb_sqr_slow(acb_t z, const acb_t x, slong prec)
65 {
66 int inexact;
67
68 mag_t am, bm, er, fr;
69
70 mag_init_set_arf(am, arb_midref(a));
71 mag_init_set_arf(bm, arb_midref(b));
72
73 mag_init(er);
74 mag_init(fr);
75
76 mag_addmul(er, am, ar);
77 mag_addmul(er, bm, br);
78 mag_mul_2exp_si(er, er, 1);
79 mag_addmul(er, ar, ar);
80 mag_addmul(er, br, br);
81
82 mag_addmul(fr, bm, ar);
83 mag_addmul(fr, am, br);
84 mag_addmul(fr, ar, br);
85 mag_mul_2exp_si(fr, fr, 1);
86
87 inexact = arf_complex_sqr(arb_midref(e), arb_midref(f),
88 arb_midref(a), arb_midref(b), prec, ARB_RND);
89
90 if (inexact & 1)
91 arf_mag_add_ulp(arb_radref(e), er, arb_midref(e), prec);
92 else
93 mag_swap(arb_radref(e), er);
94
95 if (inexact & 2)
96 arf_mag_add_ulp(arb_radref(f), fr, arb_midref(f), prec);
97 else
98 mag_swap(arb_radref(f), fr);
99
100 mag_clear(am);
101 mag_clear(bm);
102
103 mag_clear(er);
104 mag_clear(fr);
105 }
106
107 static void
_acb_mul_fast(acb_t z,const acb_t x,const acb_t y,slong prec)108 _acb_mul_fast(acb_t z, const acb_t x, const acb_t y, slong prec)
109 {
110 int inexact;
111
112 mag_t am, bm, cm, dm, er, fr;
113
114 mag_fast_init_set_arf(am, arb_midref(a));
115 mag_fast_init_set_arf(bm, arb_midref(b));
116 mag_fast_init_set_arf(cm, arb_midref(c));
117 mag_fast_init_set_arf(dm, arb_midref(d));
118
119 mag_init(er);
120 mag_init(fr);
121
122 mag_fast_addmul(er, am, cr);
123 mag_fast_addmul(er, bm, dr);
124 mag_fast_addmul(er, cm, ar);
125 mag_fast_addmul(er, dm, br);
126 mag_fast_addmul(er, ar, cr);
127 mag_fast_addmul(er, br, dr);
128
129 mag_fast_addmul(fr, am, dr);
130 mag_fast_addmul(fr, bm, cr);
131 mag_fast_addmul(fr, cm, br);
132 mag_fast_addmul(fr, dm, ar);
133 mag_fast_addmul(fr, br, cr);
134 mag_fast_addmul(fr, ar, dr);
135
136 inexact = arf_complex_mul(arb_midref(e), arb_midref(f),
137 arb_midref(a), arb_midref(b),
138 arb_midref(c), arb_midref(d), prec, ARB_RND);
139
140 if (inexact & 1)
141 arf_mag_add_ulp(arb_radref(e), er, arb_midref(e), prec);
142 else
143 mag_set(arb_radref(e), er);
144
145 if (inexact & 2)
146 arf_mag_add_ulp(arb_radref(f), fr, arb_midref(f), prec);
147 else
148 mag_set(arb_radref(f), fr);
149 }
150
151 static void
_acb_mul_slow(acb_t z,const acb_t x,const acb_t y,slong prec)152 _acb_mul_slow(acb_t z, const acb_t x, const acb_t y, slong prec)
153 {
154 int inexact;
155
156 mag_t am, bm, cm, dm, er, fr;
157
158 mag_init_set_arf(am, arb_midref(a));
159 mag_init_set_arf(bm, arb_midref(b));
160 mag_init_set_arf(cm, arb_midref(c));
161 mag_init_set_arf(dm, arb_midref(d));
162
163 mag_init(er);
164 mag_init(fr);
165
166 mag_addmul(er, am, cr);
167 mag_addmul(er, bm, dr);
168 mag_addmul(er, cm, ar);
169 mag_addmul(er, dm, br);
170 mag_addmul(er, ar, cr);
171 mag_addmul(er, br, dr);
172
173 mag_addmul(fr, am, dr);
174 mag_addmul(fr, bm, cr);
175 mag_addmul(fr, cm, br);
176 mag_addmul(fr, dm, ar);
177 mag_addmul(fr, br, cr);
178 mag_addmul(fr, ar, dr);
179
180 inexact = arf_complex_mul(arb_midref(e), arb_midref(f),
181 arb_midref(a), arb_midref(b),
182 arb_midref(c), arb_midref(d), prec, ARB_RND);
183
184 if (inexact & 1)
185 arf_mag_add_ulp(arb_radref(e), er, arb_midref(e), prec);
186 else
187 mag_swap(arb_radref(e), er);
188
189 if (inexact & 2)
190 arf_mag_add_ulp(arb_radref(f), fr, arb_midref(f), prec);
191 else
192 mag_swap(arb_radref(f), fr);
193
194 mag_clear(am);
195 mag_clear(bm);
196 mag_clear(cm);
197 mag_clear(dm);
198
199 mag_clear(er);
200 mag_clear(fr);
201 }
202
203 void
acb_mul(acb_t z,const acb_t x,const acb_t y,slong prec)204 acb_mul(acb_t z, const acb_t x, const acb_t y, slong prec)
205 {
206 if (arb_is_zero(b))
207 {
208 arb_mul(f, d, a, prec);
209 arb_mul(e, c, a, prec);
210 }
211 else if (arb_is_zero(d))
212 {
213 arb_mul(f, b, c, prec);
214 arb_mul(e, a, c, prec);
215 }
216 else if (arb_is_zero(a))
217 {
218 arb_mul(e, c, b, prec);
219 arb_mul(f, d, b, prec);
220 acb_mul_onei(z, z);
221 }
222 else if (arb_is_zero(c))
223 {
224 arb_mul(e, a, d, prec);
225 arb_mul(f, b, d, prec);
226 acb_mul_onei(z, z);
227 }
228 /* squaring = a^2-b^2, 2ab */
229 else if (x == y)
230 {
231 if (ARB_IS_LAGOM(a) && ARB_IS_LAGOM(b))
232 _acb_sqr_fast(z, x, prec);
233 else
234 _acb_sqr_slow(z, x, prec);
235 }
236 else
237 {
238 if (ARB_IS_LAGOM(a) && ARB_IS_LAGOM(b) &&
239 ARB_IS_LAGOM(c) && ARB_IS_LAGOM(d))
240 _acb_mul_fast(z, x, y, prec);
241 else
242 _acb_mul_slow(z, x, y, prec);
243 }
244 }
245
246 #undef a
247 #undef b
248 #undef c
249 #undef d
250 #undef e
251 #undef f
252 #undef ar
253 #undef br
254 #undef cr
255 #undef dr
256
257