add basic support for Gaussian integers

flintlib · Jun 21, 2022 · bce3b40 · bce3b40
1 parent 7cd542a
commit bce3b40
Show file tree

Hide file tree

Showing 11 changed files with 819 additions and 1 deletion.
diff --git a/Makefile.in b/Makefile.in
@@ -14,7 +14,7 @@ AT=@
 
 BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
    acb_dft acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \
-   arb_hypgeom bernoulli hypgeom fmpz_extras bool_mat partitions dlog \
+   arb_hypgeom bernoulli hypgeom fmpz_extras fmpzi bool_mat partitions dlog \
    double_interval arb_fmpz_poly arb_fpwrap \
    $(EXTRA_BUILD_DIRS)
 

diff --git a/doc/source/fmpzi.rst b/doc/source/fmpzi.rst
@@ -0,0 +1,67 @@
+.. _fmpzi:
+
+**fmpzi.h** -- Gaussian integers
+===============================================================================
+
+This module allows working with elements of the ring `\mathbb{Z}[i]`.
+At present, only a minimal interface is provided.
+
+Types, macros and constants
+-------------------------------------------------------------------------------
+
+.. type:: fmpzi_struct
+
+.. type:: fmpzi_t
+
+    Contains a pairs of integers representing the real and imaginary
+    parts. An *fmpzi_t* is defined as an array
+    of length one of type *fmpzi_struct*, permitting an *fmpzi_t* to
+    be passed by reference.
+
+.. macro:: fmpzi_realref(x)
+
+    Macro giving a pointer to the real part of *x*.
+
+.. macro:: fmpzi_imagref(x)
+
+    Macro giving a pointer to the imaginary part of *x*.
+
+Basic manipulation
+-------------------------------------------------------------------------------
+
+.. function:: void fmpzi_init(fmpzi_t x)
+
+.. function:: void fmpzi_clear(fmpzi_t x)
+
+.. function:: int fmpzi_equal(const fmpzi_t x, const fmpzi_t y)
+
+.. function:: void fmpzi_zero(fmpzi_t x)
+
+.. function:: void fmpzi_one(fmpzi_t x)
+
+.. function:: void fmpzi_set(fmpzi_t res, const fmpzi_t x)
+
+.. function:: void fmpzi_set_si_si(fmpzi_t res, slong a, slong b)
+
+.. function:: void fmpzi_swap(fmpzi_t x, fmpzi_t y)
+
+.. function:: void fmpzi_print(const fmpzi_t x)
+
+.. function:: void fmpzi_randtest(fmpzi_t res, flint_rand_t state, mp_bitcnt_t bits)
+
+Arithmetic
+-------------------------------------------------------------------------------
+
+.. function:: void fmpzi_conj(fmpzi_t res, const fmpzi_t x)
+
+.. function:: void fmpzi_neg(fmpzi_t res, const fmpzi_t x)
+
+.. function:: void fmpzi_add(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+
+.. function:: void fmpzi_sub(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+
+.. function:: void fmpzi_sqr(fmpzi_t res, const fmpzi_t x)
+
+.. function:: void fmpzi_mul(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+
+.. function:: void fmpzi_pow_ui(fmpzi_t res, const fmpzi_t x, ulong exp)
diff --git a/doc/source/index.rst b/doc/source/index.rst
@@ -217,6 +217,7 @@ Mainly for internal use.
 .. toctree::
    :maxdepth: 1
 
+   fmpzi.rst
    double_interval.rst
    fmpz_extras.rst
    bool_mat.rst

diff --git a/fmpzi.h b/fmpzi.h
@@ -0,0 +1,150 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef FMPZI_H
+#define FMPZI_H
+
+#ifdef FMPZI_INLINES_C
+#define FMPZI_INLINE
+#else
+#define FMPZI_INLINE static __inline__
+#endif
+
+#include "flint/fmpz.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct
+{
+    fmpz a;
+    fmpz b;
+}
+fmpzi_struct;
+
+typedef fmpzi_struct fmpzi_t[1];
+
+#define fmpzi_realref(x) (&((x)->a))
+#define fmpzi_imagref(x) (&((x)->b))
+
+FMPZI_INLINE void
+fmpzi_init(fmpzi_t x)
+{
+    fmpz_init(fmpzi_realref(x));
+    fmpz_init(fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_clear(fmpzi_t x)
+{
+    fmpz_clear(fmpzi_realref(x));
+    fmpz_clear(fmpzi_imagref(x));
+}
+
+FMPZI_INLINE int
+fmpzi_equal(const fmpzi_t x, const fmpzi_t y)
+{
+    return fmpz_equal(fmpzi_realref(x), fmpzi_realref(y)) &&
+           fmpz_equal(fmpzi_imagref(x), fmpzi_imagref(y));
+}
+
+FMPZI_INLINE void
+fmpzi_zero(fmpzi_t x)
+{
+    fmpz_zero(fmpzi_realref(x));
+    fmpz_zero(fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_one(fmpzi_t x)
+{
+    fmpz_one(fmpzi_realref(x));
+    fmpz_zero(fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_set(fmpzi_t res, const fmpzi_t x)
+{
+    fmpz_set(fmpzi_realref(res), fmpzi_realref(x));
+    fmpz_set(fmpzi_imagref(res), fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_conj(fmpzi_t res, const fmpzi_t x)
+{
+    fmpz_set(fmpzi_realref(res), fmpzi_realref(x));
+    fmpz_neg(fmpzi_imagref(res), fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_swap(fmpzi_t x, fmpzi_t y)
+{
+    fmpzi_struct t;
+    t = *x;
+    *x = *y;
+    *y = t;
+}
+
+FMPZI_INLINE void
+fmpzi_print(const fmpzi_t x)
+{
+    fmpz_print(fmpzi_realref(x));
+    if (fmpz_sgn(fmpzi_imagref(x)) >= 0)
+        flint_printf("+");
+    fmpz_print(fmpzi_imagref(x));
+    flint_printf("*I");
+}
+
+FMPZI_INLINE void
+fmpzi_set_si_si(fmpzi_t res, slong a, slong b)
+{
+    fmpz_set_si(fmpzi_realref(res), a);
+    fmpz_set_si(fmpzi_imagref(res), b);
+}
+
+FMPZI_INLINE void
+fmpzi_randtest(fmpzi_t res, flint_rand_t state, mp_bitcnt_t bits)
+{
+    fmpz_randtest(fmpzi_realref(res), state, bits);
+    fmpz_randtest(fmpzi_imagref(res), state, bits);
+}
+
+FMPZI_INLINE void
+fmpzi_neg(fmpzi_t res, const fmpzi_t x)
+{
+    fmpz_neg(fmpzi_realref(res), fmpzi_realref(x));
+    fmpz_neg(fmpzi_imagref(res), fmpzi_imagref(x));
+}
+
+FMPZI_INLINE void
+fmpzi_add(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+{
+    fmpz_add(fmpzi_realref(res), fmpzi_realref(x), fmpzi_realref(y));
+    fmpz_add(fmpzi_imagref(res), fmpzi_imagref(x), fmpzi_imagref(y));
+}
+
+FMPZI_INLINE void
+fmpzi_sub(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+{
+    fmpz_sub(fmpzi_realref(res), fmpzi_realref(x), fmpzi_realref(y));
+    fmpz_sub(fmpzi_imagref(res), fmpzi_imagref(x), fmpzi_imagref(y));
+}
+
+void fmpzi_sqr(fmpzi_t res, const fmpzi_t x);
+void fmpzi_mul(fmpzi_t res, const fmpzi_t x, const fmpzi_t y);
+void fmpzi_pow_ui(fmpzi_t res, const fmpzi_t x, ulong exp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/fmpzi/inlines.c b/fmpzi/inlines.c
@@ -0,0 +1,14 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#define FMPZI_INLINES_C
+#include "fmpzi.h"
+
diff --git a/fmpzi/mul.c b/fmpzi/mul.c
@@ -0,0 +1,116 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "fmpzi.h"
+
+void
+fmpzi_mul(fmpzi_t res, const fmpzi_t x, const fmpzi_t y)
+{
+    slong asize, bsize, csize, dsize;
+    fmpzi_struct * rp;
+    fmpzi_t tmp;
+    fmpz * t;
+    fmpz * u;
+    int xsmall, ysmall;
+
+    const fmpz *a = fmpzi_realref(x);
+    const fmpz *b = fmpzi_imagref(x);
+    const fmpz *c = fmpzi_realref(y);
+    const fmpz *d = fmpzi_imagref(y);
+
+    const fmpz ca = *a;
+    const fmpz cb = *b;
+    const fmpz cc = *c;
+    const fmpz cd = *d;
+
+    if (x == y)
+    {
+        fmpzi_sqr(res, x);
+        return;
+    }
+
+    xsmall = !COEFF_IS_MPZ(ca) && !COEFF_IS_MPZ(cb);
+    ysmall = !COEFF_IS_MPZ(cc) && !COEFF_IS_MPZ(cd);
+
+    if (xsmall && ysmall)
+    {
+        ulong thi, tlo, uhi, ulo, ahi, alo, bhi, blo;
+
+        smul_ppmm(thi, tlo, ca, cc);
+        smul_ppmm(uhi, ulo, cb, cd);
+        sub_ddmmss(ahi, alo, thi, tlo, uhi, ulo);
+
+        smul_ppmm(thi, tlo, ca, cd);
+        smul_ppmm(uhi, ulo, cb, cc);
+        add_ssaaaa(bhi, blo, thi, tlo, uhi, ulo);
+
+        fmpz_set_signed_uiui(fmpzi_realref(res), ahi, alo);
+        fmpz_set_signed_uiui(fmpzi_imagref(res), bhi, blo);
+        return;
+    }
+
+    /* todo: detect pure real and imaginary operands */
+
+    if (res == x || res == y)
+    {
+        fmpzi_init(tmp);
+        rp = tmp;
+    }
+    else
+    {
+        rp = res;
+    }
+
+    t = fmpzi_realref(rp);
+    u = fmpzi_imagref(rp);
+
+    if (!xsmall && !ysmall)
+    {
+        asize = fmpz_size(a);
+
+        if (asize >= 13)
+        {
+            bsize = fmpz_size(b);
+            csize = fmpz_size(c);
+            dsize = fmpz_size(d);
+
+            if (csize >= 13 && FLINT_ABS(asize - bsize) <= 2 && FLINT_ABS(csize - dsize) <= 2)
+            {
+                fmpz_t v;
+
+                fmpz_init(v);
+                fmpz_add(t, a, b);
+                fmpz_add(v, c, d);
+                fmpz_mul(u, t, v);
+                fmpz_mul(t, a, c);
+                fmpz_mul(v, b, d);
+                fmpz_sub(u, u, t);
+                fmpz_sub(u, u, v);
+                fmpz_sub(t, t, v);
+                fmpz_clear(v);
+
+                goto cleanup;
+            }
+        }
+    }
+
+    fmpz_mul(t, a, c);
+    fmpz_submul(t, b, d);
+    fmpz_mul(u, a, d);
+    fmpz_addmul(u, b, c);
+
+cleanup:
+    if (res == x || res == y)
+    {
+        fmpzi_swap(res, tmp);
+        fmpzi_clear(tmp);
+    }
+}