Updates to the fixed module

doc/source/fixed.rst

@@ -18,6 +18,62 @@ This module is mainly optimized for 64-bit systems. With 32-bit limbs, some
 generated straight-line and register implementations are disabled and
 evaluation goes through generic and fallback code paths.
 
+Arithmetic
+-------------------------------------------------------------------------------
+
+Newton-based division and square root
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. function:: void fixed_inv_newton_basecase(nn_ptr q, nn_srcptr a, slong an, slong n)
+              void fixed_inv_newton(nn_ptr q, nn_srcptr a, slong an, slong n)
+
+    Given `(a, an)` with `a_{an-1} \ne 0` representing a fixed-point number
+    `a \in [1/B, 1)` with `an` fraction limbs, sets `(q, n+2)` to an
+    approximation of `1/a \in (1, B]` with `n` fraction limbs and two
+    integral limbs (the highest limb may be zero). The absolute error is
+    bounded by `4 B^{-n} / a`. The *newton* suffix flags that the result
+    is not ulp-accurate. The basecase divides by ``mpn_tdiv_qr``; the
+    main function runs a Newton iteration on middle products, ported
+    from :func:`radix_inv_approx`.
+
+.. function:: void fixed_div_newton_invmul(nn_ptr q, nn_srcptr b, slong bn, nn_srcptr a, slong an, slong n)
+              void fixed_div_newton(nn_ptr q, nn_srcptr b, slong bn, nn_srcptr a, slong an, slong n)
+
+    Given a numerator `(b, bn)` with `bn \ge 1` fraction limbs representing
+    `b \in [0, 1)` and a denominator `(a, an)` with `a_{an-1} \ne 0`
+    representing `a \in [1/B, 1)`, sets `(q, n+2)` to an approximation of
+    `b/a` with `n` fraction limbs and two integral limbs. The absolute
+    error is bounded by `4 B^{-n} / a`. The *invmul* variant multiplies
+    the numerator by :func:`fixed_inv_newton`; the main function performs
+    a Karp-Markstein iteration, ported from :func:`radix_div_approx`.
+
+.. function:: void fixed_rsqrt_ui_newton_basecase(nn_ptr res, ulong a, slong n)
+              void fixed_rsqrt_ui_newton(nn_ptr res, ulong a, slong n)
+
+    Sets `(res, n)` to the fraction limbs of an approximation of
+    `1/\sqrt{a}`, requiring `2 \le a < B`. The error is bounded by
+    `2 B^{-n}`. Ported from :func:`radix_rsqrt_1_approx`.
+
+.. function:: void fixed_rsqrt_newton_basecase(nn_ptr q, nn_srcptr a, slong an, slong n)
+              void fixed_rsqrt_newton(nn_ptr q, nn_srcptr a, slong an, slong n)
+
+    Given `(a, an)` representing `a \in [B^{-2}, 1)` with `an` fraction
+    limbs (at least one of the two highest limbs must be nonzero), sets
+    `(q, n+2)` to an approximation of `1/\sqrt{a} \in (1, B]` with `n`
+    fraction limbs and two integral limbs. The absolute error is bounded
+    by `4 B^{-n} / \sqrt{a}`. Ported from :func:`radix_rsqrt_approx`;
+    the basecase combines ``mpn_sqrtrem`` and ``mpn_tdiv_qr``.
+
+.. function:: void fixed_sqrt_newton_rsqrtmul(nn_ptr q, nn_srcptr a, slong an, slong n)
+              void fixed_sqrt_newton(nn_ptr q, nn_srcptr a, slong an, slong n)
+
+    Input as for :func:`fixed_rsqrt_newton`; sets `(q, n+2)` to an
+    approximation of `\sqrt{a} \in [1/B, 1)` (the computed value can
+    round up to 1) with absolute error bounded by `4 B^{-n} / \sqrt{a}`.
+    Note that the error is proportional to `1/\sqrt{a}` rather than to
+    the output. The main function performs a Karp-Markstein iteration,
+    ported from :func:`radix_sqrt_approx`.
+
 Elementary functions
 -------------------------------------------------------------------------------
 
@@ -97,6 +153,23 @@ On 32-bit systems, it is assumes that `n \ge 2`.
     all work happens at the output precision; callers wanting sub-ulp
     accuracy should pad the precision by one limb themselves.
 
+.. function:: void fixed_exp_reduced(nn_ptr y, nn_srcptr t, slong wn, flint_bitcnt_t r, int alg)
+
+    Sets `(y, wn + 1)` (``wn`` fraction limbs and a units limb) to an
+    approximation of `\exp(t)` for a reduced argument `(t, wn)` with
+    `t < 2^{-r}`, `r \ge 16`, independent of any particular argument
+    reduction (algorithms 1 and 2 additionally require
+    `r \ge 32`).  The error is at most ``FIXED_EXP_REDUCED_MAX_ERR``
+    ulps.  *alg* selects the internal method: 0 the tuned automatic
+    choice, 1 the direct rectangular-splitting series, 2 the sinh
+    series plus a square root, 3 one bit-burst step (the leading
+    slice of *t*, on limb boundaries, evaluated by binary splitting
+    and the remainder by the sinh series at the doubled rate), 4 the
+    full bit-burst algorithm with slice lengths doubling, which is
+    asymptotically quasi-optimal for very large ``wn``.  The
+    automatic thresholds can be recalibrated with
+    ``tune/tune-exp-reduced``.
+
 .. function:: void fixed_exp_bitwise_rs(nn_ptr res, nn_srcptr x, slong n, int r)
 
     Sets `(res, n + 1)` to an approximation of `\exp((x, n))` for any
@@ -327,6 +400,32 @@ On 32-bit systems, it is assumes that `n \ge 2`.
     four small multiplications, and `\cos t'` cancels in every ratio
     just as `|W|` does.
 
+.. function:: void fixed_log1p_2mexp_ui_bs(nn_ptr res, ulong i, slong n)
+              void fixed_atan_2mexp_ui_bs(nn_ptr res, ulong i, slong n)
+
+    Sets `(res, n)` to a one-sided fixed-point approximation of
+    `\log(1 + 2^{-i})` resp. `\operatorname{atan}(2^{-i})`, `i \ge 1`:
+    at most the true value, short by no more than a couple of ulps.
+    These build the entries of the cached reduction tables (which call
+    them with one guard limb) by binary splitting in mpn arithmetic:
+    the split products are carried as truncated mantissas with
+    limb-radix exponents, an iterative basecase accumulates blocks of
+    terms with shifts and single-limb multiplications, and every
+    truncation rounds numerators down and denominators up so that the
+    final truncating division stays one-sided.  The logarithm sums
+    `2\operatorname{atanh}(1/(2^{i+1}+1))`, carrying one `q^2` factor
+    per TERM in the materialized denominator so that ranges compose
+    without any `q`-power at the merges; the `q^2` multiplications
+    all happen at the leaves, each a single ``mpn_mul_1`` whenever
+    `q^2` fits in a limb.  At high precision the direct `\log(1+x)`
+    series takes over, its larger term count offset by purely dyadic
+    merges.
+
+    The smallest indices of the tables (`i \le 6` for the logarithms,
+    `i \le 3` for the angles) are instead combined from
+    `\log 2, \log 3, \ldots` resp. Gauss-machin style atans, computed
+    by generic binary splitting.
+
 
 The reduction parameter selected by `r = 0` is tuned in two tiers.
 
@@ -381,3 +480,4 @@ size is called once before timing so that table precomputation stays
 out of the measurement, and the timing loop cycles over an array of
 random inputs so that the branchy reductions pay their real
 misprediction costs.
+

src/fixed.h

@@ -13,6 +13,7 @@
 #define FIXED_H
 
 #include "flint.h"
+#include "arb_types.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -67,6 +68,16 @@ void fixed_exp_rs(nn_ptr res, nn_srcptr x, slong n);
    the output precision, so the error bound grows linearly with r;
    callers wanting sub-ulp accuracy should pad the precision by one
    limb themselves. */
+/* exp(t) of a reduced argument t < 2^-r, r >= 32, into
+   (y, wn + 1): wn fraction limbs and a units limb.  alg: 0 = tuned
+   automatic choice, 1 = direct rectangular-splitting series,
+   2 = sinh series + square root, 3 = one bit-burst step + sinh,
+   4 = full bit-burst.  Error at most FIXED_EXP_REDUCED_MAX_ERR
+   ulps. */
+#define FIXED_EXP_REDUCED_MAX_ERR 96
+void fixed_exp_reduced(nn_ptr y, nn_srcptr t, slong wn,
+    flint_bitcnt_t r, int alg);
+
 void fixed_exp_bitwise_rs(nn_ptr res, nn_srcptr x, slong n, int r);
 
 /* fully specialized per-size implementations (default dispatch;
@@ -237,6 +248,7 @@ void _fixed_exp_rs_fallback(nn_ptr res, nn_srcptr x, slong n);
    this thread. */
 void _fixed_exp_logs_ensure(slong nv, slong rc);
 nn_srcptr _fixed_exp_logs_entry(slong i, slong n);
+void _fixed_exp_logs_clear(void);
 slong _fixed_exp_logs_max_index(void);
 
 /* Internal: number of slots the used array of _fixed_bitwise_reduce
@@ -261,7 +273,49 @@ slong _fixed_bitwise_reduce(nn_ptr t, slong wn, int r, slong istart,
    fixed_tan_bitwise_rs and fixed_atan_bitwise_rs.  Storage and
    accessors work exactly as for the logarithm table above. */
 void _fixed_atans_ensure(slong nv, slong rc);
+
+/* Internal: n-limb one-sided fixed-point approximations of the table
+   values by mpn binary splitting (at most the true value, short by a
+   couple of ulps); i >= 1. */
+void fixed_atan_2mexp_ui_bs(nn_ptr res, ulong i, slong n);
+void fixed_log1p_2mexp_ui_bs(nn_ptr res, ulong i, slong n);
+
+/* Approximate (not ulp-accurate) fixed-point inversion,
+   division and square roots by Newton / Karp-Markstein iteration on
+   middle products; ports of the radix_*_approx functions.  Writing
+   B = 2^64:
+
+   fixed_inv_newton: given (a, an), a_{an-1} != 0, representing
+   a in [1/B, 1) with an fraction limbs, sets (q, n+2) to 1/a in
+   (1, B] with n fraction limbs and two integral limbs (the top limb
+   may be zero); |error| <= 4 B^-n / a.
+
+   fixed_div_newton: numerator (b, bn) in [0, 1), denominator
+   (a, an) in [1/B, 1) with a_{an-1} != 0; sets (q, n+2) to b/a with
+   n fraction limbs and two integral limbs; |error| <= 4 B^-n / a.
+
+   fixed_rsqrt_ui_newton: 2 <= a < B; sets (res, n) to the fraction
+   limbs of 1/sqrt(a); |error| <= 2 B^-n.
+
+   fixed_rsqrt_newton: (a, an) in [B^-2, 1), one of the two top limbs
+   nonzero; sets (q, n+2) to 1/sqrt(a) in (1, B] with n fraction
+   limbs and two integral limbs; |error| <= 4 B^-n / sqrt(a).
+
+   fixed_sqrt_newton: input as for fixed_rsqrt_newton; sets (q, n+2)
+   to sqrt(a) in [1/B, 1) (the value can round up to 1);
+   |error| <= 4 B^-n / sqrt(a). */
+void fixed_inv_newton_basecase(nn_ptr q, nn_srcptr a, slong an, slong n);
+void fixed_inv_newton(nn_ptr q, nn_srcptr a, slong an, slong n);
+void fixed_div_newton_invmul(nn_ptr q, nn_srcptr b, slong bn, nn_srcptr a, slong an, slong n);
+void fixed_div_newton(nn_ptr q, nn_srcptr b, slong bn, nn_srcptr a, slong an, slong n);
+void fixed_rsqrt_ui_newton_basecase(nn_ptr res, ulong a, slong n);
+void fixed_rsqrt_ui_newton(nn_ptr res, ulong a, slong n);
+void fixed_rsqrt_newton_basecase(nn_ptr q, nn_srcptr a, slong an, slong n);
+void fixed_rsqrt_newton(nn_ptr q, nn_srcptr a, slong an, slong n);
+void fixed_sqrt_newton_rsqrtmul(nn_ptr q, nn_srcptr a, slong an, slong n);
+void fixed_sqrt_newton(nn_ptr q, nn_srcptr a, slong an, slong n);
 nn_srcptr _fixed_atans_entry(slong i, slong n);
+void _fixed_atans_clear(void);
 slong _fixed_atans_max_index(void);
 void _fixed_sin_cos_rs_fallback(nn_ptr ysin, nn_ptr ycos, nn_srcptr x,
     slong n, int alternating);
@@ -289,6 +343,110 @@ void _fixed_atan_rs_fallback(nn_ptr res, nn_srcptr x, slong n,
    it in tan_bitwise_rs.c. */
 #endif
 
+
+/* -------------------------------------------------------------- */
+/* Internal declarations (formerly fixed/impl.h), placed here so
+   that the test, tune and profile subdirectories build without
+   extra include paths. */
+
+/* Library-internal view of the cached reduction tables.
+
+   These thread-local objects are shared across the translation units
+   of the module (the dispatch files, the reduction, and the
+   specialized per-size implementations) but are deliberately NOT
+   declared in fixed.h: Windows DLLs cannot export thread-local data,
+   so external consumers -- the test suite -- go through the
+   _fixed_exp_logs_entry / _fixed_atans_entry accessors instead.
+
+   Each entry occupies _fixed_{exp_logs,atans}_n limbs: the value
+   limbs with one guard limb below them.  Consumers wanting the top n
+   limbs of entry i read tab + i * stride + (stride - n). */
+
+#define FIXED_STATIC_TAB_INLINE static inline
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern FLINT_TLS_PREFIX nn_ptr _fixed_exp_logs;
+extern FLINT_TLS_PREFIX slong _fixed_exp_logs_n;
+extern FLINT_TLS_PREFIX slong _fixed_exp_logs_r;
+
+extern FLINT_TLS_PREFIX nn_ptr _fixed_atans;
+extern FLINT_TLS_PREFIX slong _fixed_atans_n;
+extern FLINT_TLS_PREFIX slong _fixed_atans_r;
+
+/* Static prefixes of the two tables covering all reductions with
+   r <= FIXED_STATIC_TAB_R whose per-entry reads fit in
+   FIXED_STATIC_TAB_N limbs.  Unlike the dynamic tables, whose
+   bottom limb is a guard used for in-place generation, every stored
+   limb here is a value limb (the entries are the top limbs of a
+   dynamic table built one limb deeper).  The accessors below hand
+   out the static data when it suffices -- avoiding the
+   precomputation, the TLS lookup, and the wide entry stride of a
+   high-precision dynamic table -- and fall back to ensuring the
+   dynamic one. */
+
+#define FIXED_STATIC_TAB_N 12
+#define FIXED_STATIC_TAB_R 32
+
+/* internal: mpn binary splitting for sum_{k=1}^N x^k / (k! 2^(rk))
+   (exp_sum_bs.c); T needs (N (r + 128))/64 + 4 limbs, Q needs
+   (N bits(N+1))/64 + 3 */
+slong _fixed_exp_bs_num_terms(flint_bitcnt_t r, slong prec);
+void _fixed_exp_burst_factor(nn_ptr F, slong * fn, slong * fexp,
+    nn_srcptr T, slong tn, nn_srcptr Q, slong qn,
+    flint_bitcnt_t Qexp, slong cap);
+void _fixed_exp_sum_bs_powtab(nn_ptr T, slong * tn, nn_ptr Q,
+    slong * qn, flint_bitcnt_t * Qexp, nn_srcptr xp, slong xn,
+    flint_bitcnt_t r, slong N);
+
+/* internal: exact-floor entry helpers (tab_exact.c) */
+int _fixed_tab_store_floor(nn_ptr e, const arb_t x, slong nc, slong prec);
+void _fixed_tab_entry_exact(nn_ptr e, int which, ulong i, slong nc);
+
+/* fast-path entries sit at most a few guard ulps below the truth;
+   a guard limb this close to wrapping means the deficit may have
+   borrowed into the value limbs, so the entry is recomputed exactly */
+#define FIXED_TAB_GUARD_SLACK UWORD(1024)
+
+FLINT_DLL extern const ulong _fixed_exp_logs_static[(FIXED_STATIC_TAB_R + 1) * FIXED_STATIC_TAB_N];
+FLINT_DLL extern const ulong _fixed_atans_static[(FIXED_STATIC_TAB_R + 1) * FIXED_STATIC_TAB_N];
+
+/* NOTE: the two inline _tab accessors below dereference the
+   thread-local table data and are therefore usable only from
+   translation units compiled INTO the library.  Code linking
+   against a Windows DLL (tests, tuning and profiling programs)
+   must use the exported _fixed_*_entry / _fixed_*_ensure /
+   _fixed_*_clear functions instead: thread-local data cannot be
+   DLL-exported, so direct references fail to link there. */
+
+FIXED_STATIC_TAB_INLINE nn_srcptr
+_fixed_exp_logs_tab(slong nv, slong rc, slong * nc)
+{
+    if (rc <= FIXED_STATIC_TAB_R && nv + 1 <= FIXED_STATIC_TAB_N)
+    {
+        *nc = FIXED_STATIC_TAB_N;
+        return _fixed_exp_logs_static;
+    }
+    _fixed_exp_logs_ensure(nv, rc);
+    *nc = _fixed_exp_logs_n;
+    return _fixed_exp_logs;
+}
+
+FIXED_STATIC_TAB_INLINE nn_srcptr
+_fixed_atans_tab(slong nv, slong rc, slong * nc)
+{
+    if (rc <= FIXED_STATIC_TAB_R && nv + 1 <= FIXED_STATIC_TAB_N)
+    {
+        *nc = FIXED_STATIC_TAB_N;
+        return _fixed_atans_static;
+    }
+    _fixed_atans_ensure(nv, rc);
+    *nc = _fixed_atans_n;
+    return _fixed_atans;
+}
+
 #ifdef __cplusplus
 }
 #endif

src/fixed/atan_bitwise_rs.c

@@ -12,7 +12,15 @@
 #include "flint.h"
 #include "mpn_extras.h"
 #include "fixed.h"
-#include "impl.h"
+
+/* quotient length from which the argument-reduction division runs
+   through fixed_div_newton instead of mpn_tdiv_qr; set from
+   in-context measurements (tight-loop benchmarks favor Newton much
+   earlier, but embedded in the evaluation it loses until the FFT
+   range, see dev/notes) */
+#ifndef FIXED_RED_DIV_NEWTON_CUTOFF
+#define FIXED_RED_DIV_NEWTON_CUTOFF 60
+#endif
 
 /* atan(t) for t in [0, 1) by greedy vectoring (the trigonometric
    analog of the fixed_log1p_bitwise_rs reduction): the vector
@@ -53,7 +61,7 @@
    the series error (<= 15 ulp two-sided) are constants.  Total
    below FIXED_ATAN_BITWISE_RS_MAX_ERR(n, r) = 4 r + 64. */
 
-#define ATANS_A(i) (_fixed_atans + (i) * _fixed_atans_n)
+#define ATANS_A(i) (tab + (i) * nc)
 
 /* smallest n at which each reduction parameter becomes optimal
    (generated by src/fixed/tune/tune-bitwise-r.c with nmax deep enough
@@ -134,8 +142,8 @@ fixed_atan_bitwise_rs(nn_ptr res, nn_srcptr x, slong n, int r)
 
     wn = n + 1;
 
-    _fixed_atans_ensure(n, r);
-    nc = _fixed_atans_n;
+    nn_srcptr tab;
+    tab = _fixed_atans_tab(n, r, &nc);
 
     TMP_START;
     X = TMP_ALLOC((wn + n + wn + wn + wn + n + n + 2 * n + 2)
@@ -270,7 +278,25 @@ fixed_atan_bitwise_rs(nn_ptr res, nn_srcptr x, slong n, int r)
     {
         flint_mpn_zero(nd, n);
         flint_mpn_copyi(nd + n, Y, n - qr);
-        mpn_tdiv_qr(t, nd, 0, nd, ds, X, wn);
+        if (ds - wn + 1 < FIXED_RED_DIV_NEWTON_CUTOFF)
+        {
+            mpn_tdiv_qr(t, nd, 0, nd, ds, X, wn);
+        }
+        else
+        {
+            /* Same quotient by Karp-Markstein division: numerator
+               and divisor read as fractions (the numerator's n low
+               limbs are zero, so pass only its significant top
+               part), quotient computed with k + 1 fraction limbs, k
+               the quotient length.  Truncating the two guard limbs
+               leaves t within 1 + eps ulps of the floored quotient
+               (eps = 4 B^(-1) / a' < 2^-62), which the series
+               argument budget absorbs like the tdiv floor. */
+            slong k = ds - wn + 1;
+            nn_ptr q = TMP_ALLOC((k + 4) * sizeof(ulong));
+            fixed_div_newton(q, nd + n, ds - n, X, wn, k + 1);
+            flint_mpn_copyi(t, q + 2, k);
+        }
     }
 
     /* res = acc + atan(t') (the series does not support aliased