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 #include "arb_poly.h"
14 #include "arb_fmpz_poly.h" /* for minpoly */
15
16 void
_arb_cos_pi_fmpq_algebraic(arb_t c,ulong p,ulong q,slong prec)17 _arb_cos_pi_fmpq_algebraic(arb_t c, ulong p, ulong q, slong prec)
18 {
19 /* handle simple angles using exact formulas */
20 if (q <= 6)
21 {
22 if (p == 0)
23 {
24 arb_one(c);
25 }
26 else if (q == 2) /* p/q must be 1/2 */
27 {
28 arb_zero(c);
29 }
30 else if (q == 3) /* p/q must be 1/3 */
31 {
32 arb_set_ui(c, 1);
33 arb_mul_2exp_si(c, c, -1);
34 }
35 else if (q == 4) /* p/q must be 1/4 */
36 {
37 arb_sqrt_ui(c, 2, prec);
38 arb_mul_2exp_si(c, c, -1);
39 }
40 else if (q == 5) /* p/q must be 1/5 or 2/5 */
41 {
42 arb_sqrt_ui(c, 5, prec + 3);
43 arb_add_si(c, c, (p == 1) ? 1 : -1, prec);
44 arb_mul_2exp_si(c, c, -2);
45 }
46 else if (q == 6) /* p/q must be 1/6 */
47 {
48 arb_sqrt_ui(c, 3, prec);
49 arb_mul_2exp_si(c, c, -1);
50 }
51 }
52 /* reduce even denominator */
53 else if (q % 2 == 0)
54 {
55 slong extra = 2 * FLINT_BIT_COUNT(q) + 2;
56
57 if (4 * p <= q)
58 {
59 _arb_cos_pi_fmpq_algebraic(c, p, q / 2, prec + extra);
60 arb_add_ui(c, c, 1, prec + extra);
61 }
62 else
63 {
64 _arb_cos_pi_fmpq_algebraic(c, q / 2 - p, q / 2, prec + extra);
65 arb_sub_ui(c, c, 1, prec + extra);
66 arb_neg(c, c);
67 }
68
69 arb_mul_2exp_si(c, c, -1);
70 arb_sqrt(c, c, prec);
71 }
72 else
73 {
74 /* compute root of the minimal polynomial */
75 slong start_prec, eval_extra_prec;
76 fmpz_poly_t poly;
77 arb_poly_t fpoly;
78 arf_t interval_bound;
79 arb_t interval;
80
81 arf_init(interval_bound);
82 arb_init(interval);
83 fmpz_poly_init(poly);
84 arb_poly_init(fpoly);
85
86 if (p % 2 == 0)
87 arb_fmpz_poly_cos_minpoly(poly, q);
88 else
89 arb_fmpz_poly_cos_minpoly(poly, 2 * q);
90
91 eval_extra_prec = fmpz_poly_max_bits(poly) * 2; /* heuristic */
92 eval_extra_prec = FLINT_ABS(eval_extra_prec);
93 arb_poly_set_fmpz_poly(fpoly, poly, ARF_PREC_EXACT);
94
95 /* todo: smallify for accuracy */
96 start_prec = 100 + eval_extra_prec;
97 arb_const_pi(c, start_prec);
98 arb_mul_ui(c, c, p, start_prec);
99 arb_div_ui(c, c, q, start_prec);
100 arb_cos(c, c, start_prec);
101 arb_mul_2exp_si(c, c, 1); /* poly is for 2*cos */
102
103 if (100 + eval_extra_prec - 10 < prec)
104 {
105 arb_set(interval, c);
106 mag_mul_2exp_si(arb_radref(interval), arb_radref(interval), 1);
107 _arb_poly_newton_convergence_factor(interval_bound,
108 fpoly->coeffs, fpoly->length, interval, start_prec);
109 _arb_poly_newton_refine_root(c, fpoly->coeffs, fpoly->length,
110 c, interval, interval_bound, eval_extra_prec, prec);
111 }
112
113 arb_mul_2exp_si(c, c, -1);
114
115 fmpz_poly_clear(poly);
116 arb_poly_clear(fpoly);
117 arf_clear(interval_bound);
118 arb_clear(interval);
119 }
120 }
121
122 void
_arb_sin_pi_fmpq_algebraic(arb_t s,ulong p,ulong q,slong prec)123 _arb_sin_pi_fmpq_algebraic(arb_t s, ulong p, ulong q, slong prec)
124 {
125 if (q % 2 == 0)
126 {
127 p = q / 2 - p;
128
129 while ((p % 2 == 0) && (q % 2 == 0))
130 {
131 p /= 2;
132 q /= 2;
133 }
134
135 _arb_cos_pi_fmpq_algebraic(s, p, q, prec);
136 }
137 else
138 {
139 _arb_cos_pi_fmpq_algebraic(s, q - 2 * p, 2 * q, prec);
140 }
141 }
142
143 void
_arb_sin_cos_pi_fmpq_algebraic(arb_t s,arb_t c,ulong p,ulong q,slong prec)144 _arb_sin_cos_pi_fmpq_algebraic(arb_t s, arb_t c, ulong p, ulong q, slong prec)
145 {
146 slong wp;
147
148 if (q <= 6)
149 {
150 if (p == 0)
151 {
152 arb_one(c);
153 arb_zero(s);
154 return;
155 }
156 else if (q == 2) /* p/q must be 1/2 */
157 {
158 arb_zero(c);
159 arb_one(s);
160 return;
161 }
162 else if (q == 4) /* p/q must be 1/4 */
163 {
164 arb_sqrt_ui(c, 2, prec);
165 arb_mul_2exp_si(c, c, -1);
166 arb_set(s, c);
167 return;
168 }
169 }
170
171 wp = prec + 3;
172
173 /* prefer the formula with less cancellation */
174 if (p <= q / 4)
175 {
176 _arb_sin_pi_fmpq_algebraic(s, p, q, wp);
177 arb_mul(c, s, s, wp);
178 arb_sub_ui(c, c, 1, wp);
179 arb_neg(c, c);
180 arb_sqrt(c, c, prec);
181 }
182 else
183 {
184 _arb_cos_pi_fmpq_algebraic(c, p, q, wp);
185 arb_mul(s, c, c, wp);
186 arb_sub_ui(s, s, 1, wp);
187 arb_neg(s, s);
188 arb_sqrt(s, s, prec);
189 }
190 }
191
192