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