1 /*
2 Copyright (C) 2011 Fredrik Johansson
3 Copyright (C) 2019 Daniel Schultz
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include "fmpq.h"
14
15 #define _UI_MAT22_RMUL_ELEM(m11, m12, m21, m22, q) \
16 do { \
17 mp_limb_t __t1 = m12 + q*m11; \
18 mp_limb_t __t2 = m22 + q*m21; \
19 m12 = m11; \
20 m22 = m21; \
21 m11 = __t1; \
22 m21 = __t2; \
23 } while (0)
24
fmpq_dedekind_sum(fmpq_t sum,const fmpz_t h,const fmpz_t k)25 void fmpq_dedekind_sum(fmpq_t sum, const fmpz_t h, const fmpz_t k)
26 {
27 if (fmpz_cmp_ui(k, 2) <= 0 || fmpz_is_zero(h))
28 {
29 fmpq_zero(sum);
30 return;
31 }
32 else if (fmpz_fits_si(k))
33 {
34 /*
35 since the alternating sum of quotients
36 s1 - s2 + - ... + sn - 3, n odd
37 s1 - s2 + - ... - sn, n even
38 is < k in absolute value, we have a version for machine k
39 */
40 ulong a = fmpz_get_ui(k);
41 ulong b = fmpz_fdiv_ui(h, a);
42 ulong m11 = 1, m12 = 0, m21 = 0, m22 = 1;
43 ulong sum_hi, sum_lo, t = 0, q, r;
44
45 while (b != 0)
46 {
47 udiv_qrnnd(q, r, 0, a, b);
48 a = b;
49 b = r;
50 t += q;
51 _UI_MAT22_RMUL_ELEM(m11, m12, m21, m22, q);
52
53 if (b == 0)
54 {
55 /* break with det -1 */
56 t -= 3;
57 smul_ppmm(sum_hi, sum_lo, t, m11);
58 add_ssaaaa(sum_hi, sum_lo, sum_hi, sum_lo, 0, m21 + m12);
59 goto set_sum;
60 }
61
62 udiv_qrnnd(q, r, 0, a, b);
63 a = b;
64 b = r;
65 t -= q;
66 _UI_MAT22_RMUL_ELEM(m11, m12, m21, m22, q);
67 }
68
69 /* break with det +1 */
70 smul_ppmm(sum_hi, sum_lo, t, m11);
71 t = m21 - m12;
72 add_ssaaaa(sum_hi, sum_lo, sum_hi, sum_lo, FLINT_SIGN_EXT(t), t);
73
74 set_sum:
75
76 fmpz_set_signed_uiui(fmpq_numref(sum), sum_hi, sum_lo);
77 fmpz_set_ui(fmpq_denref(sum), m11);
78 }
79 else
80 {
81 _fmpq_cfrac_list_t s;
82 _fmpz_mat22_t M;
83 _fmpq_ball_t x;
84
85 _fmpq_cfrac_list_init(s);
86 s->length = -1;
87 s->want_alt_sum = 1;
88
89 _fmpz_mat22_init(M);
90 _fmpz_mat22_one(M);
91
92 _fmpq_ball_init(x);
93 x->exact = 1;
94 fmpz_set(x->left_num, k);
95 fmpz_fdiv_r(x->left_den, h, k);
96
97 if (!fmpz_is_zero(x->left_den))
98 {
99 _fmpq_ball_get_cfrac(s, M, 1, x);
100
101 /* exactly one extra iteration needed if get_cfrac is working */
102 FLINT_ASSERT(fmpz_divisible(x->left_num, x->left_den));
103
104 do {
105 fmpz_fdiv_qr(x->right_num, x->left_num, x->left_num, x->left_den);
106 _fmpz_mat22_rmul_elem(M, x->right_num);
107 _fmpq_cfrac_list_push_back(s, x->right_num);
108 fmpz_swap(x->left_num, x->left_den);
109 } while (!fmpz_is_zero(x->left_den));
110 }
111
112 FLINT_ASSERT(!fmpz_is_zero(M->_11));
113
114 FLINT_ASSERT(s->want_alt_sum == M->det);
115 if (M->det == 1)
116 {
117 fmpz_sub(fmpq_numref(sum), M->_21, M->_12);
118 }
119 else
120 {
121 FLINT_ASSERT(M->det == -1);
122 fmpz_sub_ui(s->alt_sum, s->alt_sum, 3);
123 fmpz_add(fmpq_numref(sum), M->_21, M->_12);
124 }
125
126 fmpz_swap(fmpq_denref(sum), M->_11);
127 fmpz_addmul(fmpq_numref(sum), s->alt_sum, fmpq_denref(sum));
128
129 _fmpq_ball_clear(x);
130 _fmpq_cfrac_list_clear(s);
131 _fmpz_mat22_clear(M);
132 }
133
134 fmpz_mul_ui(fmpq_denref(sum), fmpq_denref(sum), 12);
135 fmpq_canonicalise(sum); /* extra gcd seems to be unavoidable */
136 }
137