1 /*
2 Copyright (C) 2013 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 int
use_algebraic(const fmpz_t v,const fmpz_t w,slong prec)15 use_algebraic(const fmpz_t v, const fmpz_t w, slong prec)
16 {
17 fmpz q = *w;
18 int r;
19
20 if (COEFF_IS_MPZ(q))
21 return 0;
22
23 if (q <= 6)
24 return 1;
25
26 count_trailing_zeros(r, q);
27 q >>= r;
28
29 if (r >= 4 && prec < (r - 3) * 300)
30 return 0;
31
32 if (q > 1000)
33 return 0;
34
35 if (prec < 1500 + 150 * q)
36 return 0;
37
38 return 1;
39 }
40
41 void
_arb_sin_cos_pi_fmpq_oct(arb_t s,arb_t c,const fmpz_t v,const fmpz_t w,slong prec)42 _arb_sin_cos_pi_fmpq_oct(arb_t s, arb_t c,
43 const fmpz_t v, const fmpz_t w, slong prec)
44 {
45 if (use_algebraic(v, w, prec))
46 {
47 _arb_sin_cos_pi_fmpq_algebraic(s, c, *v, *w, prec);
48 }
49 else
50 {
51 arb_const_pi(s, prec);
52 arb_mul_fmpz(s, s, v, prec);
53 arb_div_fmpz(s, s, w, prec);
54 arb_sin_cos(s, c, s, prec);
55 }
56 }
57
58 void
_arb_sin_pi_fmpq_oct(arb_t s,const fmpz_t v,const fmpz_t w,slong prec)59 _arb_sin_pi_fmpq_oct(arb_t s, const fmpz_t v, const fmpz_t w, slong prec)
60 {
61 if (use_algebraic(v, w, prec))
62 {
63 _arb_sin_pi_fmpq_algebraic(s, *v, *w, prec);
64 }
65 else
66 {
67 arb_const_pi(s, prec);
68 arb_mul_fmpz(s, s, v, prec);
69 arb_div_fmpz(s, s, w, prec);
70 arb_sin(s, s, prec);
71 }
72 }
73
74 void
_arb_cos_pi_fmpq_oct(arb_t c,const fmpz_t v,const fmpz_t w,slong prec)75 _arb_cos_pi_fmpq_oct(arb_t c, const fmpz_t v, const fmpz_t w, slong prec)
76 {
77 if (use_algebraic(v, w, prec))
78 {
79 _arb_cos_pi_fmpq_algebraic(c, *v, *w, prec);
80 }
81 else
82 {
83 arb_const_pi(c, prec);
84 arb_mul_fmpz(c, c, v, prec);
85 arb_div_fmpz(c, c, w, prec);
86 arb_cos(c, c, prec);
87 }
88 }
89
90 static unsigned int
reduce_octant(fmpz_t v,fmpz_t w,const fmpq_t x)91 reduce_octant(fmpz_t v, fmpz_t w, const fmpq_t x)
92 {
93 const fmpz * p = fmpq_numref(x);
94 const fmpz * q = fmpq_denref(x);
95 unsigned int octant;
96 flint_bitcnt_t vval, wval;
97
98 if (*p > COEFF_MIN / 8 &&
99 *p < COEFF_MAX / 8 &&
100 *q > 0 &&
101 *q < COEFF_MAX / 4)
102 {
103 slong pp, qq, ww, vv, tt;
104
105 pp = *p;
106 qq = *q;
107
108 tt = pp * 4;
109 ww = tt / qq;
110 vv = tt - qq * ww;
111 /* compute correct (floor) quotient and remainder */
112 if (vv < 0)
113 {
114 ww--;
115 vv += qq;
116 }
117
118 octant = ((ulong) ww) % 8;
119 ww = qq * 4;
120 if (octant % 2 != 0)
121 vv = qq - vv;
122
123 if (vv != 0)
124 {
125 count_trailing_zeros(vval, vv);
126 count_trailing_zeros(wval, ww);
127 vval = FLINT_MIN(vval, wval);
128 vv >>= vval;
129 ww >>= vval;
130 }
131
132 fmpz_set_si(v, vv);
133 fmpz_set_si(w, ww);
134 }
135 else
136 {
137 fmpz_mul_2exp(w, p, 2);
138 fmpz_fdiv_qr(w, v, w, q);
139 octant = fmpz_fdiv_ui(w, 8);
140 fmpz_mul_2exp(w, q, 2);
141
142 if (octant % 2 != 0)
143 fmpz_sub(v, q, v);
144
145 vval = fmpz_val2(v);
146 wval = fmpz_val2(w);
147 vval = FLINT_MIN(vval, wval);
148
149 if (vval != 0)
150 {
151 fmpz_tdiv_q_2exp(v, v, vval);
152 fmpz_tdiv_q_2exp(w, w, vval);
153 }
154 }
155
156 return octant;
157 }
158
159 void
arb_sin_cos_pi_fmpq(arb_t s,arb_t c,const fmpq_t x,slong prec)160 arb_sin_cos_pi_fmpq(arb_t s, arb_t c, const fmpq_t x, slong prec)
161 {
162 fmpz_t v, w;
163 unsigned int octant;
164
165 fmpz_init(v);
166 fmpz_init(w);
167
168 octant = reduce_octant(v, w, x);
169
170 if ((octant + 1) % 4 < 2)
171 _arb_sin_cos_pi_fmpq_oct(s, c, v, w, prec);
172 else
173 _arb_sin_cos_pi_fmpq_oct(c, s, v, w, prec);
174
175 if ((octant + 6) % 8 < 4) arb_neg(c, c);
176 if (octant >= 4) arb_neg(s, s);
177
178 fmpz_clear(v);
179 fmpz_clear(w);
180 }
181
182 void
arb_sin_pi_fmpq(arb_t s,const fmpq_t x,slong prec)183 arb_sin_pi_fmpq(arb_t s, const fmpq_t x, slong prec)
184 {
185 fmpz_t v, w;
186 unsigned int octant;
187
188 fmpz_init(v);
189 fmpz_init(w);
190
191 octant = reduce_octant(v, w, x);
192
193 if ((octant + 1) % 4 < 2)
194 _arb_sin_pi_fmpq_oct(s, v, w, prec);
195 else
196 _arb_cos_pi_fmpq_oct(s, v, w, prec);
197
198 if (octant >= 4)
199 arb_neg(s, s);
200
201 fmpz_clear(v);
202 fmpz_clear(w);
203 }
204
205 void
arb_cos_pi_fmpq(arb_t c,const fmpq_t x,slong prec)206 arb_cos_pi_fmpq(arb_t c, const fmpq_t x, slong prec)
207 {
208 fmpz_t v, w;
209 unsigned int octant;
210
211 fmpz_init(v);
212 fmpz_init(w);
213
214 octant = reduce_octant(v, w, x);
215
216 if ((octant + 1) % 4 < 2)
217 _arb_cos_pi_fmpq_oct(c, v, w, prec);
218 else
219 _arb_sin_pi_fmpq_oct(c, v, w, prec);
220
221 if ((octant + 6) % 8 < 4)
222 arb_neg(c, c);
223
224 fmpz_clear(v);
225 fmpz_clear(w);
226 }
227
228