Single limb floating point $p$-adics

Makefile.in

@@ -182,7 +182,7 @@ HEADER_DIRS	:=	\
 		fq_zech					fq_zech_mat	fq_zech_poly	\
 											n_poly			\
 		fq			fq_vec		fq_mat		fq_poly			\
-		padic					padic_mat	padic_poly		\
+		padic		padic_nmod	padic_mat	padic_poly		\
 		qadic												\
 															\
 		fmpz_extras	fmpz_factor	fmpzi						\

doc/source/index.rst

@@ -297,6 +297,7 @@ p-adic numbers
    padic.rst
    padic_poly.rst
    padic_mat.rst
+   padic_nmod.rst
    qadic.rst
 
 Floating-point support code

doc/source/padic_nmod.rst

@@ -0,0 +1,161 @@
+.. _padic-nmod:
+
+**padic_nmod.h** -- floating point p-adic numbers
+===============================================================================
+
+
+Introduction
+--------------------------------------------------------------------------------
+
+The ``padic_nmod_t`` data type represents elements of `\mathbf{Q}_p` to
+precision `N`, stored in the form `x = p^v u` with `u \in \mathbf{Z}` and
+`v \in \mathbf{Z} / p^N \mathbf{Z}`.
+Arithmetic operations can be carried out with respect to a context
+containing the prime number `p` and various pieces of pre-computed data.
+Multiple precisions `N` are possible.
+
+Independent of the context, we consider a `p`-adic number
+`x = u p^v` to be in canonical form whenever either
+`p \nmid u` or `u = v = 0`.
+
+We briefly describe the interface:
+
+The functions in this module expect arguments of type ``padic_nmod_t``, and
+each variable carries its own precision.  The functions have an interface that
+is similar to the MPFR functions.  In particular, they have the same semantics,
+specified as follows:  Compute the requested operation exactly and then reduce
+the result to the precision of the output variable.
+
+Data structures
+--------------------------------------------------------------------------------
+
+A `p`-adic number of type ``padic_nmod_t`` comprises a mantissa `man`, and a
+valuation `val`.
+
+
+Context
+--------------------------------------------------------------------------------
+
+A context object for `p`-adic arithmetic contains data pertinent to `p`-adic
+computations, but which we choose not to store with each element individually.
+Currently, this includes the prime number `p`, its inverse and the quotient of
+`UWORD_MAX` by `p`, precomputed powers of `p` in the range between `1` and `N`,
+and the integers modulo `p` and `p^N`.
+
+.. function:: int padic_nmod_ctx_init(gr_ctx_t ctx, ulong p, slong n)
+
+    Initialises the context ``ctx`` with the given data.
+
+    Assumes that `p` is a prime.  This is not verified but the subsequent 
+    behaviour is undefined if `p` is a composite number.
+
+.. function:: void padic_nmod_ctx_clear(gr_ctx_t ctx)
+
+
+Memory management
+--------------------------------------------------------------------------------
+
+
+.. function:: void padic_nmod_init(padic_nmod_t res, gr_ctx_t ctx)
+
+    Initialises the `p`-adic number.
+
+.. function:: void _padic_nmod_canonicalise(padic_nmod_t x, gr_ctx_t ctx)
+
+    Brings the `p`-adic number ``x`` into canonical form.
+
+    That is to say, ensures that either `u = v = 0` or 
+    `p \nmid u`.  There is no reduction modulo a power 
+    of `p`.
+
+
+Randomisation
+--------------------------------------------------------------------------------
+
+
+.. function:: int padic_nmod_randtest(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx)
+
+    Sets ``rop`` to a random `p`-adic number.
+
+.. function:: int padic_nmod_randtest_not_zero(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx)
+
+    Sets ``rop`` to a random non-zero `p`-adic number.
+
+
+Assignments and conversions
+--------------------------------------------------------------------------------
+
+All assignment functions set the value of ``res`` from ``x``, 
+reduced to the precision of ``ctx``.
+
+.. function:: int padic_nmod_set(padic_nmod_t res, const padic_nmod_t x, gr_ctx_t ctx)
+
+    Sets ``res`` to the `p`-adic number ``x``.
+
+.. function:: int padic_nmod_set_ui(padic_nmod_t res, ulong x, gr_ctx_t ctx)
+
+    Sets the `p`-adic number ``res`` to the ``ulong`` 
+    integer ``x``.
+
+.. function:: int padic_nmod_zero(padic_nmod_t res, gr_ctx_t ctx)
+
+    Sets the `p`-adic number ``res`` to zero.
+
+.. function:: int padic_nmod_one(padic_nmod_t res, gr_ctx_t ctx)
+
+    Sets the `p`-adic number ``res`` to one, reduced modulo the 
+    precision of ``res``.
+
+
+Comparison
+--------------------------------------------------------------------------------
+
+
+.. function:: truth_t padic_nmod_is_zero(const padic_nmod_t x, gr_ctx_t ctx)
+
+    Returns whether ``x`` is equal to zero.
+
+.. function:: truth_t padic_nmod_is_one(const padic_nmod_t x, gr_ctx_t ctx)
+
+    Returns whether ``x`` is equal to one, that is, whether 
+    `u = 1` and `v = 0`.
+
+
+Arithmetic operations
+--------------------------------------------------------------------------------
+
+
+.. function:: int padic_nmod_add(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx)
+
+    Sets ``res`` to the sum of ``a`` and ``b``.
+
+.. function:: int padic_nmod_div(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx)
+
+    Sets ``res`` to the quotient of ``a`` and ``b``.
+
+.. function:: int padic_nmod_inv(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx)
+
+    Sets ``res`` to the inverse of ``a``.
+
+.. function:: int padic_nmod_mul(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx)
+
+    Sets ``res`` to the product of ``a`` and ``b``.
+
+.. function:: int padic_nmod_neg(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx)
+
+    Sets ``res`` to the additive inverse of ``a``.
+
+.. function:: int padic_nmod_sub(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx)
+
+    Sets ``res`` to the difference of ``a`` and ``b``.
+
+
+Input and output
+--------------------------------------------------------------------------------
+
+
+.. function:: void padic_nmod_println(const padic_nmod_t x, gr_ctx_t ctx)
+
+    Prints the string representation of the `p`-adic number ``x`` 
+    to the stream ``stdout``.
+

doc/source/ulong_extras.rst

@@ -1344,6 +1344,19 @@ Factorisation
     `p` we make repeated use of :func:`n_divrem2_precomp` until division
     by `p` is no longer possible.
 
+.. function:: int n_remove2_prime_inv(ulong * n, ulong p, ulong inv1, ulong inv2)
+
+    Removes the highest possible power of `p` from `n`, replacing `n` with the
+    quotient. The return value is the highest power of `p` that divided `n`.
+    Assumes `n` is not `0`.
+
+    For `p = 2` trailing zeroes are counted. For other primes `p` we require
+    `inv1` to be set to a precomputed inverse of `p` computed with
+    :func:`n_binvert`, and `inv2` to be the quotient of ``UWORD_MAX`` by `p`.
+    Then repeated divisions by `p` are done until it is no longer possible.
+
+    Assumes that `p` is prime.
+
 .. function:: void n_factor_insert(n_factor_t * factors, ulong p, ulong exp)
 
     Inserts the given prime power factor ``p^exp`` into ``factors``.

src/padic/randtest.c

@@ -13,7 +13,8 @@
 
 #define PADIC_RANDTEST_TRIES  10
 
-void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
+void
+padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
 {
     const slong N = padic_prec(rop);
 
@@ -23,7 +24,7 @@ void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
 
     if (N > 0)
     {
-        min = - ((N + 9) / 10);
+        min = -((N + 9) / 10);
         max = N;
     }
     else if (N < 0)
@@ -46,15 +47,17 @@ void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
         fmpz_clear(pow);
 }
 
-void padic_randtest_not_zero(padic_t rop, flint_rand_t state,
-                             const padic_ctx_t ctx)
+void
+padic_randtest_not_zero(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
 {
-    slong i;
+    slong i = 1;
 
-    padic_randtest(rop, state, ctx);
-
-    for (i = 1; !padic_is_zero(rop) && i < PADIC_RANDTEST_TRIES; i++)
+    do
+    {
         padic_randtest(rop, state, ctx);
+        i++;
+    }
+    while (padic_is_zero(rop) && i < PADIC_RANDTEST_TRIES);
 
     if (padic_is_zero(rop))
     {
@@ -63,8 +66,8 @@ void padic_randtest_not_zero(padic_t rop, flint_rand_t state,
     }
 }
 
-void padic_randtest_int(padic_t rop, flint_rand_t state,
-                        const padic_ctx_t ctx)
+void
+padic_randtest_int(padic_t rop, flint_rand_t state, const padic_ctx_t ctx)
 {
     const slong N = padic_prec(rop);