Exponential and logarithm

The $p$-adic exponential and logarithm are defined by their power series

which converge for $\operatorname{ord}_p(x) > 1/(p - 1)$ and for $\operatorname{ord}_p(x - 1) > 0$, respectively. For $x \in \mathbf{Q}_p$ this means $\operatorname{ord}_p(x) \ge 1$ ($\ge 2$ if $p = 2$) for the exponential and $x \equiv 1 \pmod p$ for the logarithm. On these domains they are inverse isomorphisms between $p\mathbf{Z}_p$ (or $4\mathbf{Z}_2$) and $1 + p\mathbf{Z}_p$ (or $1 + 4\mathbf{Z}_2$). The functions return 0 and leave the output unchanged when the series diverges, and otherwise return 1 and set the output to the value reduced modulo $p^N$. Since the result of the exponential is a unit and that of the logarithm has valuation at least that of $x - 1$, the precision $N$ of the output is also, up to the valuation, its number of correct digits.

Algorithms#

Both series are truncated after the first term whose valuation, and that of every later term, is at least $N$ (_padic_exp_bound and _padic_log_bound). The tail bounds use $\operatorname{ord}_p(i!) \le (i - 1)/(p - 1)$ and $\operatorname{ord}_p(i) \le \log_p i$; for an argument of valuation $v$ this gives about $N(p - 1)/((p - 1)v - 1)$ terms for the exponential and about $N/v$ for the logarithm.

• Rectangular splitting evaluates the truncated series with about $2\sqrt n$ multiplications of full-size integers modulo $p^{N + k}$, the $k$ extra digits absorbing the powers of $p$ in the denominators $i!$ or $i$. It is used by padic_exp for $N < 1024$ and by padic_log for $N < 512/b$, where $b$ is the bit length of $p$.
• The balanced algorithms split the argument into pieces $x = x_0 + x_1 + \cdots$ with $x_j$ having valuation at least $2^j$ and $2^j$ digits, and use $\exp(x) = \prod \exp(x_j)$, respectively $\log(x) = \sum \log(1 - y_j)$ for a matching factorisation $x = \prod (1 - y_j)$. Each short series is summed exactly by binary splitting. This is the $p$-adic analogue of Brent's bit-burst method and costs $O(M(N \log p) \log^2 N)$, quasi-linear in $N$.
• Satoh's method (padic_log_satoh) uses $\log(x) = p^{-k} \log(x^{p^k})$: since $x^{p^k} - 1$ has valuation at least $k + 1$, the series for the right-hand side has about $N/(k + 1)$ terms, and with $k = \lfloor \sqrt N \rfloor$ the power (about $k \log_2 p$ squarings) and the shorter series together cost $O(\sqrt N \log p)$ multiplications [SST2003]. It is never selected automatically, and for $N < 16$ it falls back to rectangular splitting.

Example#

With $p = 5$ and precision 20, the logarithm inverts the exponential at $15 = 3 \cdot 5$, by each of the three algorithms. Neither series converges at the unit 3.

#include <flint/flint.h>
#include <flint/padic.h>

int main(void)
{
    fmpz_t p;
    padic_ctx_t ctx;
    padic_t x, y, z;
    int ok;

    fmpz_init_set_ui(p, 5);
    padic_ctx_init(ctx, p, 0, 20, PADIC_TERSE);
    padic_init2(x, 20);
    padic_init2(y, 20);
    padic_init2(z, 20);

    padic_set_ui(x, 15, ctx);            /* x = 3 * 5, valuation 1 */
    ok = padic_exp(y, x, ctx);
    flint_printf("exp(15) converges: %d\ny = ", ok);
    padic_print(y, ctx); flint_printf("\n");

    ok = padic_log(z, y, ctx);
    flint_printf("log(y) converges: %d, log(y) == 15: %d\n",
                 ok, padic_equal(z, x));

    padic_log_balanced(z, y, ctx);
    flint_printf("balanced agrees: %d\n", padic_equal(z, x));
    padic_log_satoh(z, y, ctx);
    flint_printf("Satoh agrees: %d\n", padic_equal(z, x));

    padic_set_ui(x, 3, ctx);             /* a unit, not 1 mod 5 */
    flint_printf("exp(3) converges: %d, log(3) converges: %d\n",
                 padic_exp(y, x, ctx), padic_log(z, x, ctx));

    padic_clear(x); padic_clear(y); padic_clear(z);
    padic_ctx_clear(ctx);
    fmpz_clear(p);
    return 0;
}

exp(15) converges: 1
y = 87481552191941
log(y) converges: 1, log(y) == 15: 1
balanced agrees: 1
Satoh agrees: 1
exp(3) converges: 0, log(3) converges: 0

The value printed is the representative of $\exp(15)$ modulo $5^{20}$. Applying the logarithm to it recovers 15 exactly, because $\log(\exp(x)) = x$ on the domain of convergence and both functions are computed to the full precision of the output.