Line data Source code
1 : #pragma once
2 :
3 : /// @file
4 : /// @author hacatu
5 : /// @version 0.2.0
6 : /// @section LICENSE
7 : /// This Source Code Form is subject to the terms of the Mozilla Public
8 : /// License, v. 2.0. If a copy of the MPL was not distributed with this
9 : /// file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 : /// @section DESCRIPTION
11 : /// Dirichlet Hyperbola based functions for computing sums of multiplicative functions at a single value quickly
12 : ///
13 : /// For a function f defined over the natural numbers, we say f is multiplicative if f(ab) = f(a)f(b) whenever a and b are coprime
14 : /// (have gcd 1). There are many functions in number theory that satisfy this condition, with one of the most well known being
15 : /// euler's phi function.
16 : ///
17 : /// A brief list of common multiplicative functions is:
18 : /// (see here https://en.wikipedia.org/wiki/Multiplicative_function)
19 : /// The dirichlet convolution identity (explained later): I(n) = 1 if n == 1; 0 otherwise
20 : /// Sometimes written epsilon(n), (very confusingly) u(n), or delta(n, 0) (since it is just a kroneker delta function of course)
21 : /// The constant function: u(n) = 1
22 : /// Sometimes written 1(n)
23 : /// The identity function: N(n) = n
24 : /// Sometimes written Id(n) or id(n)
25 : /// The power functions: N_k(n) = n**k
26 : /// Sometimes written Id_k(n)
27 : /// The divisor count function: d(n) = #{k | n} (the number of natural numbers including 1 and n which divide n)
28 : /// The generalized divisor functions: d_k(n) = # k-tuples ks of natural numbers such that ks multiplies out to n
29 : /// The divisor power sum functions: sigma_a(n) = sum(k | n, k**a) (note that a can be any real number)
30 : /// The mobius function: mu(n) = 0 if n is not squarefree; (-1)**Omega(n) otherwise, where Omega(n) is the number of prime factors of n with multiplicity
31 : /// Euler's totient function: phi(n) = #{k = 1 ... n : (k coprime n)}
32 : /// Any constant raised to the power of Omega(n) or omega(n) (the number of prime factors of n counted with and without multiplicity respectively)
33 : /// The number of non-isomorphic abelian groups of order n: a(n)
34 : /// The Ramanujan tau function: tau(n)
35 : /// The unitary divisor power sums: sigma_a^*(n) = sum(d | n where (d coprime n/d), d**a)
36 : /// All dirichlet characters, including the legendre symbol (n/p) for fixed p and gcd(n, m) for fixed m
37 : ///
38 : /// Some of these are fully multiplicative, meaning not only do they satisfy f(ab) = f(a)f(b) for a, b coprime, but they satisfy it for all a, b.
39 : /// In particular, I, u, N, and N_k are fully multiplicative.
40 : ///
41 : /// We often want to find the sum of a multiplicative function, often denoted by a capital version, for instance phi and Phi, f and F (for a general function, etc).
42 : /// We also often want to be able to analyze multiplicative functions in terms of simpler functions.
43 : ///
44 : /// Multiplicative functions are fully determined by their values at prime powers, which often leads to efficient sieve based approaches to calculating them.
45 : /// This header provides some functions for doing this, the nut_euler_sieve family of functions, although these are mainly intended as low level
46 : /// subroutines for the upcoming Dirichlet Hyperbola algorithm.
47 : ///
48 : /// However, if we only want to compute the sum of a multiplicative function up to some value n, sieve based methods will basically be O(nlogn),
49 : /// which seems very bad considering how structured multiplicative functions are.
50 : ///
51 : /// Indeed, we can often do better, in particular when we want to find the sum H(n) of some multiplicative function h(n) defined as h = f <*> g
52 : /// where <*> denotes Dirichlet convolution. What is Dirichlet convolution? It is an operation for combining functions, similar to ordinary
53 : /// multiplication, division, addition, subtraction, composition, and so on. It is defined as (f <*> g)(n) = sum(k | n, f(k)g(n/k)).
54 : /// This is mostly useful because when f and g are multiplicative f <*> g will be as well, and it lets us build up more complicated functions in terms
55 : /// of simpler ones.
56 : /// (see here https://en.wikipedia.org/wiki/Dirichlet_convolution)
57 : ///
58 : /// Obviously, not all operations on multiplicative functions will produce multiplicative functions. For example, if f is multiplicative, 2f is not,
59 : /// so multiplicative functions are not closed under scalar multiplication, and thus they are not closed under addition.
60 : /// They are closed under multiplication and dirichlet convolution though.
61 : ///
62 : /// The Dirichlet Hyperbola algorithm lets us rewrite H(x) = sum(n <= x, h(n)) in terms of F and G as well as f and g, but crucially it lets us evaluate them
63 : /// at fewer values.
64 : ///
65 : /// In particular, sum(n <= x, (f <*> g)(n)) = sum(n <= x, ab = n, f(a)g(b)) = sum(ab <= x, f(a)g(b))
66 : /// = sum(a <= A, b <= x/a, f(a)g(b)) + sum(b <= B, a <= x/b, f(a)g(b)) - sum(a <= A, b <= B, f(a)g(b))
67 : /// = sum(a <= A, f(a)G(x/a)) + sum(b <= B, F(x/b)g(b)) - F(A)G(B)
68 : /// where AB = x.
69 : /// (see here https://gbroxey.github.io/blog/2023/04/30/mult-sum-1.html)
70 : /// (see here https://angyansheng.github.io/blog/dirichlet-hyperbola-method)
71 : /// (you may also find this instructive https://en.wikipedia.org/wiki/Marginal_distribution)
72 : /// To accomplish this rearrangement, we first expanded (f <*> g) using the definition of dirichlet convolution.
73 : /// Then, we interpreted the sum both as an a-indexed sum and as a b-indexed sum.
74 : /// We could find the sum either way, but considering both sets us up to cut the work down quadratically.
75 : /// These two ways of interpreting the sum are EXACTLY EQUIVALENT to the trick of transforming an integral of f(x) in terms of x into
76 : /// an integral of f_inverse(y) in terms of y. For integrals, it is common to just pick whichever interpretation is easier, but for our sum,
77 : /// we benefit from realizing that if we sum f(a)g(b) for (a <= A, b <= x/a) and then separately for (b <= B, a <= x/b), then we double count the
78 : /// rectangular region (a <= A, b <= B).
79 : /// But more importantly, if we pick A = B = sqrt(x), then we've turned one sum over (f <*> g) for n from 1 to x, which we could accomplish in
80 : /// about O(xlogx), into two sums over f*G and F*g for n from 1 to sqrt(x), which we can sometimes accomplish in O(sqrt(x)).
81 : ///
82 : /// There is a big catch though, we now need to know values of F and G as well as f and g.
83 : /// Wonderfully though, we only need to know F(x/n) and G(x/n) for integral values of n, and floor(x/n) has at most 2sqrt(x) distinct values.
84 : /// That, arguably, is the real magic, since it lets us store only O(sqrt(x)) values of F and G.
85 : /// When we want to compute only H(x), this is fine, but when we want to repeatedly compute the table of values of H(x/n) for integral values of n,
86 : /// it is often better to make A larger than B, eg A = x**(2/3) and B = x**(1/3).
87 :
88 : #include <stdbool.h>
89 : #include <assert.h>
90 :
91 : #include <nut/modular_math.h>
92 :
93 : /// Wrapper to hold values of some multiplicative function.
94 : /// Stores all values f(n) for n up to (and including) y, then
95 : /// stores values f(x/n) for n less than x/y
96 : typedef struct{
97 : int64_t x;
98 : int64_t y, yinv;
99 : int64_t *buf;
100 : } nut_Diri;
101 :
102 : /// Compute the sum of the divisor count function d from 1 to max.
103 : /// @param [in] max: inclusive upper bound of range to compute sum for
104 : /// @param [in] m: modulus to reduce result by, or 0 to skip reducing
105 : /// @return the sum d(1) + ... + d(max)
106 : NUT_ATTR_CONST
107 : uint64_t nut_dirichlet_D(uint64_t max, uint64_t m);
108 :
109 : /// Compute the sum of the divisor sum function sigma from 1 to max.
110 : /// @param [in] max: inclusive upper bound of range to compute sum for
111 : /// @param [in] m: modulus to reduce result by, or 0 to skip reducing
112 : /// @return the sum sigma(1) + ... + sigma(max)
113 : NUT_ATTR_CONST
114 : uint64_t nut_dirichlet_Sigma(uint64_t max, uint64_t m);
115 :
116 : /// Given a table of values of a multiplicative function f, compute (f <*> u)(x) for all x from 1 to n
117 : /// See { @link nut_euler_sieve_conv_u} for more info}
118 : /// @param [in] n: inclusive upper bound of range to compute f <*> u over
119 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
120 : /// @param [in] f_vals: table of values for f
121 : /// @param [out] f_conv_u_vals: table to store values of f <*> u in
122 : /// @return true on success, false on allocation failure
123 : NUT_ATTR_NONNULL(3, 4)
124 : NUT_ATTR_ACCESS(read_only, 3, 1)
125 : NUT_ATTR_ACCESS(read_write, 4, 1)
126 : bool nut_euler_sieve_conv_u(int64_t n, int64_t m, const int64_t f_vals[static n+1], int64_t f_conv_u_vals[restrict static n+1]);
127 :
128 : /// Given a table of values of a multiplicative function f, compute (f <*> N)(x) for all x from 1 to n
129 : /// See { @link nut_euler_sieve_conv_u} for more info}
130 : /// @param [in] n: inclusive upper bound of range to compute f <*> N over
131 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
132 : /// @param [in] f_vals: table of values for f
133 : /// @param [out] f_conv_N_vals: table to store values of f <*> N in
134 : /// @return true on success, false on allocation failure
135 : NUT_ATTR_NONNULL(3, 4)
136 : NUT_ATTR_ACCESS(read_only, 3, 1)
137 : NUT_ATTR_ACCESS(read_write, 4, 1)
138 : bool nut_euler_sieve_conv_N(int64_t n, int64_t m, const int64_t f_vals[static n+1], int64_t f_conv_N_vals[restrict static n+1]);
139 :
140 : /// Given tables of values of multiplicative functions f and g, compute (f <*> g)(x) for all x from 1 to n
141 : /// This uses Euler's sieve, a variant of the sieve of Eratosthenes that only marks off each multiple once.
142 : /// This generally has worse performance than the sieve of Eratosthenes, but some preliminary tests showed
143 : /// that Eratosthenes is only about 14% faster in release mode than Euler, so currently only Euler's sieve
144 : /// is implemented.
145 : /// @param [in] n: inclusive upper bound of range to compute f <*> g over
146 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
147 : /// @param [in] f_vals: table of values for f
148 : /// @param [in] g_vals: table of values for f
149 : /// @param [out] f_conv_g_vals: table to store values of f <*> g in
150 : /// @return true on success, false on allocation failure
151 : NUT_ATTR_NONNULL(3, 4, 5)
152 : NUT_ATTR_ACCESS(read_only, 3, 1)
153 : NUT_ATTR_ACCESS(read_only, 4, 1)
154 : NUT_ATTR_ACCESS(read_write, 5, 1)
155 : bool nut_euler_sieve_conv(int64_t n, int64_t m, const int64_t f_vals[static n+1], const int64_t g_vals[static n+1], int64_t f_conv_g_vals[restrict static n+1]);
156 :
157 : /// Allocate internal buffers for a diri table
158 : /// self->buf will have f(0) through f(y) at indicies 0 through y,
159 : /// and then f(x/1) through f(x/(yinv - 1)) at indicies y + 1 through y + yinv - 1.
160 : /// The first two indicies, f(0) and f(1), are basically dummies though: f(0) should
161 : /// never be relied on, and f(1) should always be 1
162 : /// @param [in] x: inclusive upper bound of the domain of interest
163 : /// @param [in] y: inclusive upper bound of the dense portion of the domain of interest.
164 : /// Will be increased to sqrt(x) if needed, so 0 can be used to explicitly signal you want that behavior
165 : /// @return true on success, false on allocation failure
166 : NUT_ATTR_NONNULL(1)
167 : NUT_ATTR_ACCESS(write_only, 1)
168 : bool nut_Diri_init(nut_Diri *self, int64_t x, int64_t y);
169 :
170 : /// Copy the values from one diri table to another, which must be initialized
171 : NUT_ATTR_NONNULL(1, 2)
172 : NUT_ATTR_ACCESS(read_write, 1)
173 : NUT_ATTR_ACCESS(read_only, 2)
174 : void nut_Diri_copy(nut_Diri *restrict dest, const nut_Diri *restrict src);
175 :
176 : /// Deallocate internal buffers for a diri table
177 : NUT_ATTR_NONNULL(1)
178 : NUT_ATTR_ACCESS(read_write, 1)
179 : void nut_Diri_destroy(nut_Diri *self);
180 :
181 : NUT_ATTR_PURE
182 : NUT_ATTR_NONNULL(1)
183 : NUT_ATTR_ACCESS(read_only, 1)
184 55230523 : static inline int64_t nut_Diri_get_dense(const nut_Diri *self, int64_t k){
185 55230523 : assert(k >= 0 && k <= self->y);
186 55230523 : return self->buf[k];
187 : }
188 :
189 : NUT_ATTR_NONNULL(1)
190 : NUT_ATTR_ACCESS(read_write, 1)
191 2320856 : static inline void nut_Diri_set_dense(nut_Diri *self, int64_t k, int64_t v){
192 2320856 : assert(k >= 0 && k <= self->y);
193 2320856 : self->buf[k] = v;
194 2320856 : }
195 :
196 : NUT_ATTR_PURE
197 : NUT_ATTR_NONNULL(1)
198 : NUT_ATTR_ACCESS(read_only, 1)
199 92066 : static inline int64_t nut_Diri_get_sparse(const nut_Diri *self, int64_t k){
200 92066 : assert(k > 0 && k <= self->yinv);
201 92066 : return self->buf[self->y + k];
202 : }
203 :
204 : NUT_ATTR_NONNULL(1)
205 : NUT_ATTR_ACCESS(read_write, 1)
206 16961 : static inline void nut_Diri_set_sparse(nut_Diri *self, int64_t k, int64_t v){
207 16961 : assert(k > 0 && k <= self->yinv);
208 16961 : self->buf[self->y + k] = v;
209 16961 : }
210 :
211 : /// Compute the value table for the dirichlet convolution identity I(n) = \{1 if n == 0, 0 otherwise\}
212 : /// Just memset's the dense part, then sets index 1 and y + 1 through y + yinv - 1 to 1 (remember the sparse indicies are sums)
213 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
214 : NUT_ATTR_NONNULL(1)
215 : NUT_ATTR_ACCESS(read_write, 1)
216 : void nut_Diri_compute_I(nut_Diri *self);
217 :
218 : /// Compute the value table for the unit function u(n) = 1
219 : /// Fills the dense part of the table with 1s, and computes the sparse entries with table[y + k] = U(x/k) = x/k
220 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
221 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
222 : NUT_ATTR_NONNULL(1)
223 : NUT_ATTR_ACCESS(read_write, 1)
224 : void nut_Diri_compute_u(nut_Diri *self, int64_t m);
225 :
226 : /// Compute the value table for the identity function N(n) = n
227 : /// Fills the dense part of the table with increasing numbers, and computes the sparse entries with table[y + k] = sum_N(v = x/k) = v*(v+1)/2
228 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
229 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
230 : NUT_ATTR_NONNULL(1)
231 : NUT_ATTR_ACCESS(read_write, 1)
232 : void nut_Diri_compute_N(nut_Diri *self, int64_t m);
233 :
234 : /// Compute the value table for the mobius function mu(n) (whose sum is called the Mertens function)
235 : /// Requires both an initialized nut_Diri from {@link nut_Diri_init} and a packed table of mobius values from {@link nut_sieve_mobius}.
236 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
237 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
238 : /// @param [in] mobius: packed table of mobius values, from {@link nut_sieve_mobius} (with upper bound self->y).
239 : NUT_ATTR_NONNULL(1, 3)
240 : NUT_ATTR_ACCESS(read_write, 1)
241 : NUT_ATTR_ACCESS(read_only, 3)
242 : void nut_Diri_compute_mertens(nut_Diri *restrict self, int64_t m, const uint8_t mobius[restrict static self->y/4 + 1]);
243 :
244 : /// Compute the value table for the kth generalized divisor function dk(n)
245 : /// dk(n) = da(n) <*> db(n) whenever a + b = k.
246 : /// This function uses binary exponentiation to compute dk in only log(k) convolutions.
247 : /// However, if all dk's are needed, calling { @link nut_Diri_compute_conv_u}
248 : /// is more efficient
249 : /// @param [in, out] self: the table to store the result in, initialized by { @link nut_Diri_init}.
250 : /// self->y, self->x, and self->yinv must be set and consistent with the inputs
251 : /// @param [in] k: k, ie which generalized divisor function to compute. dk is u <*> u <*> ... <*> u with k u's
252 : /// @param [in] m: modulus to reduce the result by, or 0 to skip reducing
253 : /// @param [in, out] f_tbl, g_tbl: temporary tables for scratch work, initialized by { @link nut_Diri_init}, fields y, x, and yinv must be set
254 : /// Will still have scratch data stored on return
255 : NUT_ATTR_NONNULL(1, 4, 5)
256 : NUT_ATTR_ACCESS(read_write, 1, 4, 5)
257 : bool nut_Diri_compute_dk(nut_Diri *restrict self, uint64_t k, int64_t m, nut_Diri *restrict f_tbl, nut_Diri *restrict g_tbl);
258 :
259 : /// Compute the value table for the kth power function N^k(n)
260 : /// N^k(n) = n^k, so we can find the sums using https://en.wikipedia.org/wiki/Faulhaber%27s_formula
261 : /// (we actually do this by taking a lower triangular matrix of pascal's triangle and inverting it)
262 : /// This is one of the core components needed to compute the Dirichlet sum table for f where f(p) is a polynomial,
263 : /// the other one being https://en.wikipedia.org/wiki/Jordan%27s_totient_function
264 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
265 : /// @param [in] k: the power of the power function to compute
266 : /// @param [in] m: modulus to reduce the result by, or 0 to skip reducing
267 : NUT_ATTR_NONNULL(1)
268 : NUT_ATTR_ACCESS(read_write, 1)
269 : bool nut_Diri_compute_Nk(nut_Diri *restrict self, uint64_t k, int64_t m);
270 :
271 : /// Compute the value table for the Jacobi symbol jacobi(n/p)
272 : /// Currently, p MUST be and odd prime
273 : /// This is 1 if n is a quadratic residue (perfect square) mod p, -1 if it is not, or 0 if n is a multiple of p
274 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
275 : /// @param [in] p: prime modulus for jacobi symbol
276 : NUT_ATTR_NONNULL(1)
277 : NUT_ATTR_ACCESS(read_write, 1)
278 : bool nut_Diri_compute_J(nut_Diri *restrict self, uint64_t p);
279 :
280 : /// Compute the value table for a given Dirichlet Character Chi, whose total is 0
281 : /// A dirichlet character is some n-periodic function Chi(x) which is
282 : /// - completely multiplicative (Chi(ab) = Chi(a)Chi(b))
283 : /// - 0 if x is not coprime to n
284 : /// And here, we also require that the sum of Chi(x) over its period is 0, so it is trivial to sum
285 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
286 : /// @param [in] n: period of the dirichlet character
287 : /// @param [in] m: modulus to reduce the result by, or 0 to skip reducing. Note that the input values must already be reduced
288 : /// @param [in] period_tbl: table of values of Chi(x) for 0 <= x < n. These values should be reduced mod m and sum to 0 mod m, or if m is 0, some of these values should be 0
289 : bool nut_Diri_compute_Chi_balanced(nut_Diri *restrict self, uint64_t n, int64_t m, const int64_t period_tbl[restrict static n]);
290 :
291 : /// Compute the value table for the prime counting function, aka prime pi
292 : /// The indicator function for primes is effectively multiplicative,
293 : /// and the prime counting function is its sum. But also, Lucy Hedgehog's prime counting
294 : /// algorithm essentially uses Dirichlet tables already, hence its inclusion in this part of the library
295 : /// @param [in, out] self: the table to store the result in, and take the bounds from. Must be initialized
296 : NUT_ATTR_NONNULL(1)
297 : NUT_ATTR_ACCESS(read_write, 1)
298 : bool nut_Diri_compute_pi(nut_Diri *restrict self);
299 :
300 : /// Compute the value table for h = f <*> u, the dirichlet convolution of f and u (the unit function u(n) = 1), given the value table for f
301 : /// See { @link nut_Diri_compute_conv } for details, this is just that function but with several specializations due to u being very simple.
302 : /// @param [in, out] self: the table to store the result in, initialized by { @link nut_Diri_init}.
303 : /// self->y, self->x, and self->yinv must be set and consistent with the inputs
304 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
305 : /// @param [in] f_tbl: table for the first operand
306 : NUT_ATTR_NONNULL(1, 3)
307 : NUT_ATTR_ACCESS(read_write, 1)
308 : NUT_ATTR_ACCESS(read_only, 3)
309 : bool nut_Diri_compute_conv_u(nut_Diri *restrict self, int64_t m, const nut_Diri *restrict f_tbl);
310 :
311 : /// Compute the value table for h = f <*> N, the dirichlet convolution of f and N (the identity function N(n) = n), given the value table for f
312 : /// See { @link nut_Diri_compute_conv } for details, this is just that function but with several specializations due to N being very simple.
313 : /// @param [in, out] self: the table to store the result in, initialized by { @link nut_Diri_init}.
314 : /// self->y, self->x, and self->yinv must be set and consistent with the inputs
315 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
316 : /// @param [in] f_tbl: table for the first operand
317 : NUT_ATTR_NONNULL(1, 3)
318 : NUT_ATTR_ACCESS(read_write, 1)
319 : NUT_ATTR_ACCESS(read_only, 3)
320 : bool nut_Diri_compute_conv_N(nut_Diri *restrict self, int64_t m, const nut_Diri *restrict f_tbl);
321 :
322 : /// Compute the value table for h = f <*> g, the dirichlet convolution of f and g, given value tables for the operands
323 : /// self must have been initialized using { @link nut_Diri_init}, and in particular the lengths and cutoffs for self, f_tbl, and g_tbl
324 : /// must be set and match
325 : /// This uses a sieve to compute h(1) through h(y), then uses dirichlet's hyperbola formula to compute H(x/1) through H(x/(yinv-1))
326 : /// If one of the operands is a simple standard multiplicative function like the unit function u or the mobius function mu, then
327 : /// try to use a specialized function from this library and use that instead.
328 : /// @param [in, out] self: the table to store the result in, initialized by { @link nut_Diri_init}.
329 : /// self->y, self->x, and self->yinv must be set and consistent with the inputs
330 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
331 : /// @param [in] f_tbl: table for the first operand (dirichlet convolution is commutative, so order doesn't matter)
332 : /// @param [in] g_tbl: table for the second operand
333 : NUT_ATTR_NONNULL(1, 3, 4)
334 : NUT_ATTR_ACCESS(read_write, 1)
335 : NUT_ATTR_ACCESS(read_only, 3)
336 : NUT_ATTR_ACCESS(read_only, 4)
337 : bool nut_Diri_compute_conv(nut_Diri *restrict self, int64_t m, const nut_Diri *f_tbl, const nut_Diri *g_tbl);
338 :
339 :
340 : /// Compute the value table for h such that f = g <*> h, aka h = f </> g where the division is in terms of dirichlet convolution
341 : /// self must have been initialized using { @link nut_Diri_init }, and in particular the lengths and cutoffs for self, f_tbl, and g_tbl
342 : /// must be set and match
343 : /// This is the inverse of <*>, in the sense that if h = f </> g, then h <*> g = f. In fact, this function is essentially
344 : /// implemented by solving this for H(v) and applying the dirichlet hyperbola method (see the code comments for details).
345 : /// Unlike the `compute_conv_*` family of functions, this needs to allocate space for scratch work, which can cause it to fail if
346 : /// there is not enough memory. About 16*(self->y + 1) bytes will be allocated.
347 : /// Currently there are not specialized functions for dividing by simple standard functions.
348 : /// @param [in, out] self: the table to store the result in, initialized by { @link nut_Diri_init }.
349 : /// self->y, self->x, and self->yinv must be set ant consistent with the inputs
350 : /// @param [in] m: modulus to reduce results by, or 0 to skip reducing
351 : /// @param [in] f_tbl: table for the first operand (the "numerator" aka "dividend" in the "division")
352 : /// @param [in] g_tbl: the table for the second operand (the "denominator" aka "divisor" in the "division")
353 : NUT_ATTR_NONNULL(1, 3, 4)
354 : NUT_ATTR_ACCESS(read_write, 1, 3, 4)
355 : bool nut_Diri_convdiv(nut_Diri *restrict self, int64_t m, const nut_Diri *restrict f_tbl, const nut_Diri *restrict g_tbl);
356 :
|
