Skip to content

Commit

Permalink
Merge pull request #2004 from flintlib/nfloat9
Browse files Browse the repository at this point in the history 
Complex nfloats
  • Loading branch information
@fredrik-johansson
fredrik-johansson authored May 27, 2024
2 parents 3d6eaab + a6f6fb8 commit 80e9b24
Show file tree
Hide file tree
Showing 15 changed files with 2,473 additions and 18 deletions.
75 changes: 72 additions & 3 deletions doc/source/nfloat.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,9 @@ Types, macros and constants
NFLOAT_MAX_LIMBS
The number of limbs `n` permitted as precision. The current
limits are are `1 \le n \le 33` on a 64-bit machine and
`1 \le n \le 66` on a 32-bit machine, permitting precision
up to 2112 bits. The upper limit exists so that elements and
limits are are `1 \le n \le 66` on a 64-bit machine and
`1 \le n \le 132` on a 32-bit machine, permitting precision
up to 4224 bits. The upper limit exists so that elements and
temporary buffers are safe to allocate on the stack and so that
simple operations like swapping are not too expensive.
Expand All @@ -84,6 +84,7 @@ Types, macros and constants
nfloat512_struct
nfloat1024_struct
nfloat2048_struct
nfloat4096_struct
nfloat64_t
nfloat128_t
nfloat192_t
Expand All @@ -92,6 +93,7 @@ Types, macros and constants
nfloat512_t
nfloat1024_t
nfloat2048_t
nfloat4096_t

For convenience we define types of the correct structure size for
some common levels of bit precision. An ``nfloatX_t`` is defined as
Expand Down Expand Up @@ -254,6 +256,7 @@ These methods are interchangeable with their ``gr`` counterparts.
int nfloat_mul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx)
int nfloat_submul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx)
int nfloat_addmul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx)
int nfloat_sqr(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx)

.. function:: int nfloat_mul_2exp_si(nfloat_ptr res, nfloat_srcptr x, slong y, gr_ctx_t ctx)

Expand Down Expand Up @@ -329,3 +332,69 @@ Internal functions
int _nfloat_sub_4(nfloat_ptr res, nn_srcptr x, slong xexp, int xsgnbit, nn_srcptr y, slong delta, gr_ctx_t ctx)
int _nfloat_add_n(nfloat_ptr res, nn_srcptr xd, slong xexp, int xsgnbit, nn_srcptr yd, slong delta, slong nlimbs, gr_ctx_t ctx)
int _nfloat_sub_n(nfloat_ptr res, nn_srcptr xd, slong xexp, int xsgnbit, nn_srcptr yd, slong delta, slong nlimbs, gr_ctx_t ctx)

Complex numbers
-------------------------------------------------------------------------------

Complex floating-point numbers have the obvious representation as
real pairs.

.. type:: nfloat_complex_ptr
nfloat_complex_srcptr

.. function:: int nfloat_complex_ctx_init(gr_ctx_t ctx, slong prec, int flags)

.. macro:: NFLOAT_COMPLEX_CTX_DATA_NLIMBS(ctx)
NFLOAT_COMPLEX_RE(ptr, ctx)
NFLOAT_COMPLEX_IM(ptr, ctx)
NFLOAT_COMPLEX_IS_SPECIAL(x, ctx)
NFLOAT_COMPLEX_IS_ZERO(x, ctx)

.. function:: void nfloat_complex_init(nfloat_complex_ptr res, gr_ctx_t ctx)
void nfloat_complex_clear(nfloat_complex_ptr res, gr_ctx_t ctx)
int nfloat_complex_zero(nfloat_complex_ptr res, gr_ctx_t ctx)
int nfloat_complex_get_acf(acf_t res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_set_acf(nfloat_complex_ptr res, const acf_t x, gr_ctx_t ctx)
int nfloat_complex_get_acb(acb_t res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_set_acb(nfloat_complex_ptr res, const acb_t x, gr_ctx_t ctx)
int nfloat_complex_write(gr_stream_t out, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_randtest(nfloat_complex_ptr res, flint_rand_t state, gr_ctx_t ctx)
void nfloat_complex_swap(nfloat_complex_ptr x, nfloat_complex_ptr y, gr_ctx_t ctx)
int nfloat_complex_set(nfloat_complex_ptr res, nfloat_complex_ptr x, gr_ctx_t ctx)
int nfloat_complex_one(nfloat_complex_ptr res, gr_ctx_t ctx)
int nfloat_complex_neg_one(nfloat_complex_ptr res, gr_ctx_t ctx)
truth_t nfloat_complex_is_zero(nfloat_complex_srcptr x, gr_ctx_t ctx)
truth_t nfloat_complex_is_one(nfloat_complex_srcptr x, gr_ctx_t ctx)
truth_t nfloat_complex_is_neg_one(nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_i(nfloat_complex_ptr res, gr_ctx_t ctx)
int nfloat_complex_pi(nfloat_complex_ptr res, gr_ctx_t ctx)
int nfloat_complex_conj(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_re(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_im(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
truth_t nfloat_complex_equal(nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx)
int nfloat_complex_set_si(nfloat_complex_ptr res, slong x, gr_ctx_t ctx)
int nfloat_complex_set_ui(nfloat_complex_ptr res, ulong x, gr_ctx_t ctx)
int nfloat_complex_set_fmpz(nfloat_complex_ptr res, const fmpz_t x, gr_ctx_t ctx)
int nfloat_complex_set_fmpq(nfloat_complex_ptr res, const fmpq_t x, gr_ctx_t ctx)
int nfloat_complex_set_d(nfloat_complex_ptr res, double x, gr_ctx_t ctx)
int nfloat_complex_neg(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_add(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx)
int nfloat_complex_sub(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx)
int _nfloat_complex_sqr_naive(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx)
int _nfloat_complex_sqr_standard(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx)
int _nfloat_complex_sqr_karatsuba(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx)
int _nfloat_complex_sqr(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx)
int nfloat_complex_sqr(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int _nfloat_complex_mul_naive(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, nfloat_srcptr c, nfloat_srcptr d, gr_ctx_t ctx)
int _nfloat_complex_mul_standard(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, nfloat_srcptr c, nfloat_srcptr d, gr_ctx_t ctx)
int _nfloat_complex_mul_karatsuba(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, nfloat_srcptr c, nfloat_srcptr d, gr_ctx_t ctx)
int nfloat_complex_mul(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx)
int nfloat_complex_inv(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx)
int nfloat_complex_div(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx)
void _nfloat_complex_vec_init(nfloat_complex_ptr res, slong len, gr_ctx_t ctx)
void _nfloat_complex_vec_clear(nfloat_complex_ptr res, slong len, gr_ctx_t ctx)
int _nfloat_complex_vec_zero(nfloat_complex_ptr res, slong len, gr_ctx_t ctx)
int _nfloat_complex_vec_set(nfloat_complex_ptr res, nfloat_complex_srcptr x, slong len, gr_ctx_t ctx)
int _nfloat_complex_vec_add(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx)
int _nfloat_complex_vec_sub(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx)

2 changes: 1 addition & 1 deletion src/gr.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -681,7 +681,7 @@ typedef enum
GR_CTX_COMPLEX_EXTENDED_CA,
GR_CTX_RR_ARB, GR_CTX_CC_ACB,
GR_CTX_REAL_FLOAT_ARF, GR_CTX_COMPLEX_FLOAT_ACF,
GR_CTX_NFLOAT,
GR_CTX_NFLOAT, GR_CTX_NFLOAT_COMPLEX,
GR_CTX_FMPZ_POLY, GR_CTX_FMPQ_POLY, GR_CTX_GR_POLY,
GR_CTX_FMPZ_MPOLY, GR_CTX_GR_MPOLY,
GR_CTX_FMPZ_MPOLY_Q,
Expand Down
48 changes: 48 additions & 0 deletions src/gr/acb.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "gr_generic.h"
#include "gr_vec.h"
#include "gr_poly.h"
#include "nfloat.h"

typedef struct
{
Expand Down Expand Up @@ -281,6 +282,18 @@ _gr_acb_set_other(acb_t res, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
return GR_DOMAIN;
}

case GR_CTX_NFLOAT_COMPLEX:
if (NFLOAT_CTX_HAS_INF_NAN(x_ctx)) /* todo */
{
return GR_UNABLE;
}
else
{
nfloat_complex_get_acb(res, x, x_ctx);
acb_set_round(res, res, ACB_CTX_PREC(ctx));
return GR_SUCCESS;
}

case GR_CTX_RR_ARB:
arb_set_round(acb_realref(res), x, ACB_CTX_PREC(ctx));
arb_zero(acb_imagref(res));
Expand Down Expand Up @@ -958,6 +971,39 @@ _gr_acb_arg(acb_t res, const acb_t x, const gr_ctx_t ctx)
return GR_SUCCESS;
}

int
_gr_acb_cmp(int * res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
{
if (arb_is_zero(acb_imagref(x)) && arb_is_zero(acb_imagref(y)) &&
((arb_is_exact(acb_realref(x)) && arb_is_exact(acb_realref(y))) || !arb_overlaps(acb_realref(x), acb_realref(y))))
{
*res = arf_cmp(arb_midref(acb_realref(x)), arb_midref(acb_realref(y)));
return GR_SUCCESS;
}
else
{
*res = 0;
return GR_UNABLE;
}
}

int
_gr_acb_cmpabs(int * res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
{
acb_t t, u;

*t = *x;
*u = *y;

if (arf_sgn(arb_midref(acb_realref(t))) < 0)
ARF_NEG(arb_midref(acb_realref(t)));

if (arf_sgn(arb_midref(acb_realref(u))) < 0)
ARF_NEG(arb_midref(acb_realref(u)));

return _gr_acb_cmp(res, t, u, ctx);
}

int
_gr_acb_pi(acb_t res, const gr_ctx_t ctx)
{
Expand Down Expand Up @@ -2172,6 +2218,8 @@ gr_method_tab_input _acb_methods_input[] =
{GR_METHOD_SGN, (gr_funcptr) _gr_acb_sgn},
{GR_METHOD_CSGN, (gr_funcptr) _gr_acb_csgn},
{GR_METHOD_ARG, (gr_funcptr) _gr_acb_arg},
{GR_METHOD_CMP, (gr_funcptr) _gr_acb_cmp},
{GR_METHOD_CMPABS, (gr_funcptr) _gr_acb_cmpabs},
{GR_METHOD_PI, (gr_funcptr) _gr_acb_pi},
{GR_METHOD_EXP, (gr_funcptr) _gr_acb_exp},
{GR_METHOD_EXPM1, (gr_funcptr) _gr_acb_expm1},
Expand Down
11 changes: 8 additions & 3 deletions src/gr/test_ring.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1369,12 +1369,16 @@ gr_test_zero_one(gr_ctx_t R, flint_rand_t state, int test_flags)
}

status |= gr_randtest(a, state, R);
status |= gr_one(a, R);
status |= gr_neg(a, a, R);
status |= gr_neg_one(a, R);
equal = gr_is_neg_one(a, R);
if (status == GR_SUCCESS && equal == T_FALSE)
status = GR_TEST_FAIL;

status |= gr_neg(a, a, R);
equal = gr_is_one(a, R);
if (status == GR_SUCCESS && equal == T_FALSE)
status = GR_TEST_FAIL;

if ((test_flags & GR_TEST_ALWAYS_ABLE) && (status & GR_UNABLE))
status = GR_TEST_FAIL;

Expand Down Expand Up @@ -3939,7 +3943,8 @@ gr_test_floating_point(gr_ctx_t R, slong iters, int test_flags)
gr_test_iter(R, state, "add: aliasing", gr_test_add_aliasing, iters, test_flags);
gr_test_iter(R, state, "sub: equal neg add", gr_test_sub_equal_neg_add, iters, test_flags);
gr_test_iter(R, state, "sub: aliasing", gr_test_sub_aliasing, iters, test_flags);
gr_test_iter(R, state, "mul: commutative", gr_test_mul_commutative, iters, test_flags);
/* can fail for complex */
/* gr_test_iter(R, state, "mul: commutative", gr_test_mul_commutative, iters, test_flags); */
gr_test_iter(R, state, "mul: aliasing", gr_test_mul_aliasing, iters, test_flags);
gr_test_iter(R, state, "div: aliasing", gr_test_div_aliasing, iters, test_flags);
gr_test_iter(R, state, "pow: aliasing", gr_test_pow_aliasing, iters, test_flags);
Expand Down
14 changes: 14 additions & 0 deletions src/mpn_extras.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,20 @@ char * _flint_mpn_get_str(mp_srcptr x, mp_size_t n);
(r0) = __r0; (r1) = __r1; (r2) = __r2; \
} while (0)

#define FLINT_MPN_SQR_2X2(r3, r2, r1, r0, a1, a0) \
do { \
mp_limb_t __t1, __t2, __t3; \
mp_limb_t __r3, __r2, __r1, __r0; \
mp_limb_t __a1 = (a1), __a0 = (a0); \
umul_ppmm(__t2, __t1, __a0, __a1); \
add_sssaaaaaa(__t3, __t2, __t1, 0, __t2, __t1, 0, __t2, __t1); \
umul_ppmm(__r1, __r0, __a0, __a0); \
umul_ppmm(__r3, __r2, __a1, __a1); \
add_sssaaaaaa(__r3, __r2, __r1, __r3, __r2, __r1, __t3, __t2, __t1); \
(r0) = __r0; (r1) = __r1; (r2) = __r2; (r3) = __r3; \
} while (0)


/* {s0,s1,s2} = u[0]v[n-1] + u[1]v[n-2] + ... */
/* Assumes n >= 2 */
#define NN_DOTREV_S3_1X1(s2, s1, s0, u, v, n) \
Expand Down
Loading

0 comments on commit 80e9b24

Please sign in to comment.