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 "acb_hypgeom.h"
13
main()14 int main()
15 {
16 slong iter;
17 flint_rand_t state;
18
19 flint_printf("li_series....");
20 fflush(stdout);
21
22 flint_randinit(state);
23
24 for (iter = 0; iter < 200 * arb_test_multiplier(); iter++)
25 {
26 slong m, n1, n2, n3, bits1, bits2, bits3;
27 acb_poly_t S, A, B, C, T, U;
28 int offset;
29
30 bits1 = 2 + n_randint(state, 200);
31 bits2 = 2 + n_randint(state, 200);
32 bits3 = 2 + n_randint(state, 200);
33
34 m = 1 + n_randint(state, 10);
35 n1 = 1 + n_randint(state, 10);
36 n2 = 1 + n_randint(state, 10);
37 n3 = FLINT_MIN(n1, n2);
38 offset = n_randint(state, 2);
39
40 acb_poly_init(S);
41 acb_poly_init(A);
42 acb_poly_init(B);
43 acb_poly_init(C);
44 acb_poly_init(T);
45 acb_poly_init(U);
46
47 acb_poly_randtest(S, state, m, bits1, 3);
48 acb_poly_randtest(A, state, m, bits1, 3);
49 acb_poly_randtest(B, state, m, bits1, 3);
50
51 acb_hypgeom_li_series(A, S, offset, n1, bits2);
52 acb_hypgeom_li_series(B, S, offset, n2, bits3);
53
54 acb_poly_set(C, A);
55 acb_poly_truncate(C, n3);
56 acb_poly_truncate(B, n3);
57
58 /* [li(h(x))]' log(h(x)) = h'(x) */
59 acb_poly_derivative(T, A, bits2);
60 acb_poly_log_series(U, S, n3, bits2);
61 acb_poly_mullow(T, T, U, FLINT_MAX(0, n3 - 1), bits2);
62 acb_poly_derivative(U, S, bits2);
63 acb_poly_truncate(U, FLINT_MAX(0, n3 - 1));
64
65 if (!acb_poly_overlaps(B, C) || !acb_poly_overlaps(T, U))
66 {
67 flint_printf("FAIL\n\n");
68 flint_printf("S = "); acb_poly_printd(S, 15); flint_printf("\n\n");
69 flint_printf("A = "); acb_poly_printd(A, 15); flint_printf("\n\n");
70 flint_printf("B = "); acb_poly_printd(B, 15); flint_printf("\n\n");
71 flint_printf("T = "); acb_poly_printd(T, 15); flint_printf("\n\n");
72 flint_printf("U = "); acb_poly_printd(U, 15); flint_printf("\n\n");
73 flint_abort();
74 }
75
76 acb_hypgeom_li_series(S, S, offset, n1, bits2);
77
78 if (!acb_poly_overlaps(A, S))
79 {
80 flint_printf("FAIL (aliasing)\n\n");
81 flint_abort();
82 }
83
84 acb_poly_clear(S);
85 acb_poly_clear(A);
86 acb_poly_clear(B);
87 acb_poly_clear(C);
88 acb_poly_clear(T);
89 acb_poly_clear(U);
90 }
91
92 flint_randclear(state);
93 flint_cleanup();
94 flint_printf("PASS\n");
95 return EXIT_SUCCESS;
96 }
97