1 /*
2 Copyright (C) 2016 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpq_poly.h"
13
14 void
_fmpq_poly_sin_cos_series_basecase_can(fmpz * S,fmpz_t Sden,fmpz * C,fmpz_t Cden,const fmpz * A,const fmpz_t Aden,slong Alen,slong n,int can)15 _fmpq_poly_sin_cos_series_basecase_can(fmpz * S, fmpz_t Sden,
16 fmpz * C, fmpz_t Cden, const fmpz * A, const fmpz_t Aden, slong Alen, slong n, int can)
17 {
18 fmpz_t t, u, v;
19 slong j, k;
20
21 Alen = FLINT_MIN(Alen, n);
22
23 if (Alen == 1 || n == 1)
24 {
25 fmpz_zero(S);
26 fmpz_one(C);
27 _fmpz_vec_zero(S + 1, n - 1);
28 _fmpz_vec_zero(C + 1, n - 1);
29 fmpz_one(Sden);
30 fmpz_one(Cden);
31 return;
32 }
33
34 /* support aliasing */
35 if (A == S || A == C)
36 {
37 fmpz * tmp = _fmpz_vec_init(Alen + 1);
38 _fmpz_vec_set(tmp, A, Alen);
39 fmpz_set(tmp + Alen, Aden);
40 _fmpq_poly_sin_cos_series_basecase_can(S, Sden, C, Cden,
41 tmp, tmp + Alen, Alen, n, can);
42 _fmpz_vec_clear(tmp, Alen + 1);
43 return;
44 }
45
46 fmpz_init(t);
47 fmpz_init(u);
48 fmpz_init(v);
49
50 fmpz_fac_ui(t, n - 1);
51 fmpz_pow_ui(v, Aden, n - 1);
52 fmpz_mul(Sden, t, v);
53 fmpz_set(C, Sden);
54 fmpz_set(Cden, Sden);
55 fmpz_zero(S);
56
57 for (k = 1; k < n; k++)
58 {
59 fmpz_zero(t);
60 fmpz_zero(u);
61
62 for (j = 1; j < FLINT_MIN(Alen, k + 1); j++)
63 {
64 fmpz_mul_ui(v, A + j, j);
65 fmpz_submul(t, v, S + k - j);
66 fmpz_addmul(u, v, C + k - j);
67 }
68
69 fmpz_mul_ui(v, Aden, k);
70 fmpz_divexact(C + k, t, v);
71 fmpz_divexact(S + k, u, v);
72 }
73
74 if (can & 1) _fmpq_poly_canonicalise(S, Sden, n);
75 if (can & 2) _fmpq_poly_canonicalise(C, Cden, n);
76
77 fmpz_clear(t);
78 fmpz_clear(u);
79 fmpz_clear(v);
80 }
81
82 void
_fmpq_poly_sin_cos_series_basecase(fmpz * S,fmpz_t Sden,fmpz * C,fmpz_t Cden,const fmpz * A,const fmpz_t Aden,slong Alen,slong n)83 _fmpq_poly_sin_cos_series_basecase(fmpz * S, fmpz_t Sden,
84 fmpz * C, fmpz_t Cden, const fmpz * A, const fmpz_t Aden, slong Alen, slong n)
85 {
86 _fmpq_poly_sin_cos_series_basecase_can(S, Sden, C, Cden, A, Aden, Alen, n, 3);
87 }
88
89 void
_fmpq_poly_sin_cos_series_tangent(fmpz * S,fmpz_t Sden,fmpz * C,fmpz_t Cden,const fmpz * A,const fmpz_t Aden,slong Alen,slong n)90 _fmpq_poly_sin_cos_series_tangent(fmpz * S, fmpz_t Sden,
91 fmpz * C, fmpz_t Cden, const fmpz * A, const fmpz_t Aden,
92 slong Alen, slong n)
93 {
94 fmpz * t;
95 fmpz * u;
96 fmpz_t tden;
97 fmpz_t uden;
98
99 Alen = FLINT_MIN(Alen, n);
100
101 t = _fmpz_vec_init(n);
102 u = _fmpz_vec_init(n);
103 fmpz_init(tden);
104 fmpz_init(uden);
105
106 /* sin(x) = 2*tan(x/2)/(1+tan(x/2)^2) */
107 /* cos(x) = (1-tan(x/2)^2)/(1+tan(x/2)^2) */
108
109 /* t = tan(x/2), u = 1+tan(x/2)^2 */
110 fmpz_mul_ui(uden, Aden, 2);
111 _fmpq_poly_tan_series(t, tden, A, uden, Alen, n);
112 _fmpq_poly_mullow(u, uden, t, tden, n, t, tden, n, n);
113 fmpz_set(u, uden);
114 _fmpq_poly_canonicalise(u, uden, n);
115
116 /* C = 1/(1+tan(x/2))^2 */
117 _fmpq_poly_inv_series(C, Cden, u, uden, n, n);
118
119 /* S = sin(x)/2 */
120 _fmpq_poly_mullow(S, Sden, t, tden, n, C, Cden, n, n);
121 _fmpq_poly_canonicalise(S, Sden, n);
122
123 /* u = sin(x)/2 * tan(x/2) */
124 /* C = C - u */
125 _fmpq_poly_mullow(u, uden, S, Sden, n, t, tden, n, n);
126 _fmpq_poly_canonicalise(u, uden, n);
127 _fmpq_poly_sub(C, Cden, C, Cden, n, u, uden, n);
128
129 /* S = sin(x) */
130 _fmpq_poly_scalar_mul_ui(S, Sden, S, Sden, n, 2);
131
132 _fmpz_vec_clear(t, n);
133 _fmpz_vec_clear(u, n);
134 fmpz_clear(tden);
135 fmpz_clear(uden);
136 }
137
138 void
_fmpq_poly_sin_cos_series(fmpz * S,fmpz_t Sden,fmpz * C,fmpz_t Cden,const fmpz * A,const fmpz_t Aden,slong Alen,slong n)139 _fmpq_poly_sin_cos_series(fmpz * S, fmpz_t Sden,
140 fmpz * C, fmpz_t Cden, const fmpz * A, const fmpz_t Aden,
141 slong Alen, slong n)
142 {
143 if (Alen < 20 || n < 20)
144 _fmpq_poly_sin_cos_series_basecase(S, Sden, C, Cden, A, Aden, Alen, n);
145 else
146 _fmpq_poly_sin_cos_series_tangent(S, Sden, C, Cden, A, Aden, Alen, n);
147 }
148
149 void
fmpq_poly_sin_cos_series(fmpq_poly_t res1,fmpq_poly_t res2,const fmpq_poly_t poly,slong n)150 fmpq_poly_sin_cos_series(fmpq_poly_t res1, fmpq_poly_t res2, const fmpq_poly_t poly, slong n)
151 {
152 if (n == 0)
153 {
154 fmpq_poly_zero(res1);
155 fmpq_poly_zero(res2);
156 return;
157 }
158
159 if (fmpq_poly_is_zero(poly) || n == 1)
160 {
161 fmpq_poly_zero(res1);
162 fmpq_poly_one(res2);
163 return;
164 }
165
166 if (!fmpz_is_zero(poly->coeffs))
167 {
168 flint_printf("Exception (fmpq_poly_sin_cos_series). Constant term != 0.\n");
169 flint_abort();
170 }
171
172 fmpq_poly_fit_length(res1, n);
173 fmpq_poly_fit_length(res2, n);
174 _fmpq_poly_sin_cos_series(res1->coeffs, res1->den,
175 res2->coeffs, res2->den, poly->coeffs, poly->den, poly->length, n);
176 _fmpq_poly_set_length(res1, n);
177 _fmpq_poly_normalise(res1);
178 _fmpq_poly_set_length(res2, n);
179 _fmpq_poly_normalise(res2);
180 }
181
182