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 void
acb_mul_naive(acb_t z,const acb_t x,const acb_t y,slong prec)15 acb_mul_naive(acb_t z, const acb_t x, const acb_t y, slong prec)
16 {
17 #define a acb_realref(x)
18 #define b acb_imagref(x)
19 #define c acb_realref(y)
20 #define d acb_imagref(y)
21 #define e acb_realref(z)
22 #define f acb_imagref(z)
23
24 if (arb_is_zero(b))
25 {
26 arb_mul(f, d, a, prec);
27 arb_mul(e, c, a, prec);
28 }
29 else if (arb_is_zero(d))
30 {
31 arb_mul(f, b, c, prec);
32 arb_mul(e, a, c, prec);
33 }
34 else if (arb_is_zero(a))
35 {
36 arb_mul(e, c, b, prec);
37 arb_mul(f, d, b, prec);
38 acb_mul_onei(z, z);
39 }
40 else if (arb_is_zero(c))
41 {
42 arb_mul(e, a, d, prec);
43 arb_mul(f, b, d, prec);
44 acb_mul_onei(z, z);
45 }
46 /* squaring = a^2-b^2, 2ab */
47 else if (x == y)
48 {
49 /* aliasing */
50 if (z == x)
51 {
52 arb_t t;
53 arb_init(t);
54
55 arb_mul(t, a, b, prec);
56 arb_mul_2exp_si(t, t, 1);
57 arb_mul(e, a, a, prec);
58 arb_mul(f, b, b, prec);
59 arb_sub(e, e, f, prec);
60 arb_swap(f, t);
61
62 arb_clear(t);
63 }
64 else
65 {
66 arb_mul(e, a, a, prec);
67 arb_mul(f, b, b, prec);
68 arb_sub(e, e, f, prec);
69 arb_mul(f, a, b, prec);
70 arb_mul_2exp_si(f, f, 1);
71 }
72 }
73 else
74 {
75 /* aliasing */
76 if (z == x)
77 {
78 arb_t t, u;
79
80 arb_init(t);
81 arb_init(u);
82
83 arb_mul(t, a, c, prec);
84 arb_mul(u, a, d, prec);
85
86 arb_mul(e, b, d, prec);
87 arb_sub(e, t, e, prec);
88
89 arb_mul(f, b, c, prec);
90 arb_add(f, u, f, prec);
91
92 arb_clear(t);
93 arb_clear(u);
94 }
95 else if (z == y)
96 {
97 arb_t t, u;
98
99 arb_init(t);
100 arb_init(u);
101
102 arb_mul(t, a, c, prec);
103 arb_mul(u, b, c, prec);
104
105 arb_mul(e, b, d, prec);
106 arb_sub(e, t, e, prec);
107
108 arb_mul(f, a, d, prec);
109 arb_add(f, u, f, prec);
110
111 arb_clear(t);
112 arb_clear(u);
113 }
114 else
115 {
116 arb_t t;
117 arb_init(t);
118
119 arb_mul(e, a, c, prec);
120 arb_mul(t, b, d, prec);
121 arb_sub(e, e, t, prec);
122
123 arb_mul(f, a, d, prec);
124 arb_mul(t, b, c, prec);
125 arb_add(f, f, t, prec);
126
127 arb_clear(t);
128 }
129 }
130
131 #undef a
132 #undef b
133 #undef c
134 #undef d
135 #undef e
136 #undef f
137 }
138