hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 08:51:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions
arb/arf.h
Fredrik Johansson 8b2c1704f0 use flint fft mpn_mul
2022-05-18 15:43:00 +02:00

1317 lines
32 KiB
C
Raw Permalink Blame History

/*
Copyright (C) 2014 Fredrik Johansson
2x2 mul code taken from MPFR 2.3.0
(Copyright (C) 1991-2007 Free Software Foundation, Inc.)
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 ARF_H
#define ARF_H
#ifdef ARF_INLINES_C
#define ARF_INLINE
#else
#define ARF_INLINE static __inline__
#endif
#include <stdio.h>
#include <math.h>
#include "flint/flint.h"
#include "fmpr.h"
#include "mag.h"
#ifndef flint_abort
#if __FLINT_RELEASE <= 20502
#define flint_abort abort
#endif
#endif
#ifdef __cplusplus
extern "C" {
#endif
#define arf_rnd_t fmpr_rnd_t
#define ARF_RND_DOWN FMPR_RND_DOWN
#define ARF_RND_UP FMPR_RND_UP
#define ARF_RND_FLOOR FMPR_RND_FLOOR
#define ARF_RND_CEIL FMPR_RND_CEIL
#define ARF_RND_NEAR FMPR_RND_NEAR
ARF_INLINE int
arf_rounds_down(arf_rnd_t rnd, int sgnbit)
{
if (rnd == ARF_RND_DOWN) return 1;
if (rnd == ARF_RND_UP) return 0;
if (rnd == ARF_RND_FLOOR) return !sgnbit;
return sgnbit;
}
ARF_INLINE int
arf_rounds_up(arf_rnd_t rnd, int sgnbit)
{
if (rnd == ARF_RND_DOWN) return 0;
if (rnd == ARF_RND_UP) return 1;
if (rnd == ARF_RND_FLOOR) return sgnbit;
return !sgnbit;
}
ARF_INLINE mpfr_rnd_t
arf_rnd_to_mpfr(arf_rnd_t rnd)
{
if (rnd == ARF_RND_DOWN) return MPFR_RNDZ;
else if (rnd == ARF_RND_UP) return MPFR_RNDA;
else if (rnd == ARF_RND_FLOOR) return MPFR_RNDD;
else if (rnd == ARF_RND_CEIL) return MPFR_RNDU;
else return MPFR_RNDN;
}
/* Allow 'infinite' precision for operations where a result can be computed exactly. */
#define ARF_PREC_EXACT WORD_MAX
#define ARF_PREC_ADD(prec,extra) ((prec) == ARF_PREC_EXACT ? ARF_PREC_EXACT : (prec) + (extra))
#define ARF_RESULT_EXACT 0
#define ARF_RESULT_INEXACT 1
/* Range where we can skip fmpz overflow checks for exponent manipulation. */
#define ARF_MAX_LAGOM_EXP MAG_MAX_LAGOM_EXP
#define ARF_MIN_LAGOM_EXP MAG_MIN_LAGOM_EXP
/* Exponent values used to encode special values. */
#define ARF_EXP_ZERO 0
#define ARF_EXP_NAN COEFF_MIN
#define ARF_EXP_POS_INF (COEFF_MIN+1)
#define ARF_EXP_NEG_INF (COEFF_MIN+2)
/* Direct access to the exponent. */
#define ARF_EXP(x) ((x)->exp)
#define ARF_EXPREF(x) (&(x)->exp)
/* Finite and with lagom big exponents. */
#define ARF_IS_LAGOM(x) (ARF_EXP(x) >= ARF_MIN_LAGOM_EXP && \
ARF_EXP(x) <= ARF_MAX_LAGOM_EXP)
/* More than two limbs (needs pointer). */
#define ARF_HAS_PTR(x) ((x)->size > ((2 << 1) + 1))
/* Raw size field (encodes both limb size and sign). */
#define ARF_XSIZE(x) ((x)->size)
/* Construct size field value from size in limbs and sign bit. */
#define ARF_MAKE_XSIZE(size, sgnbit) ((((mp_size_t) size) << 1) | (sgnbit))
/* The limb size, and the sign bit. */
#define ARF_SIZE(x) (ARF_XSIZE(x) >> 1)
#define ARF_SGNBIT(x) (ARF_XSIZE(x) & 1)
/* Assumes non-special value */
#define ARF_NEG(x) (ARF_XSIZE(x) ^= 1)
/* Note: may also be hardcoded in a few places. */
#define ARF_NOPTR_LIMBS 2
/* Direct access to the limb data. */
#define ARF_NOPTR_D(x) ((x)->d.noptr.d)
#define ARF_PTR_D(x) ((x)->d.ptr.d)
#define ARF_PTR_ALLOC(x) ((x)->d.ptr.alloc)
/* Encoding for special values. */
#define ARF_IS_SPECIAL(x) (ARF_XSIZE(x) == 0)
/* Value is +/- a power of two */
#define ARF_IS_POW2(x) (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP)
/* To set a special value, first call this and then set the exponent. */
#define ARF_MAKE_SPECIAL(x) \
do { \
fmpz_clear(ARF_EXPREF(x)); \
ARF_DEMOTE(x); \
ARF_XSIZE(x) = 0; \
} while (0)
typedef struct
{
mp_limb_t d[ARF_NOPTR_LIMBS];
}
mantissa_noptr_struct;
typedef struct
{
mp_size_t alloc;
mp_ptr d;
}
mantissa_ptr_struct;
typedef union
{
mantissa_noptr_struct noptr;
mantissa_ptr_struct ptr;
}
mantissa_struct;
typedef struct
{
fmpz exp;
mp_size_t size;
mantissa_struct d;
}
arf_struct;
typedef arf_struct arf_t[1];
typedef arf_struct * arf_ptr;
typedef const arf_struct * arf_srcptr;
void _arf_promote(arf_t x, mp_size_t n);
void _arf_demote(arf_t x);
/* Warning: does not set size! -- also doesn't demote exponent. */
#define ARF_DEMOTE(x) \
do { \
if (ARF_HAS_PTR(x)) \
_arf_demote(x); \
} while (0)
/* Get mpn pointer and size (xptr, xn) for read-only use. */
#define ARF_GET_MPN_READONLY(xptr, xn, x) \
do { \
xn = ARF_SIZE(x); \
if (xn <= ARF_NOPTR_LIMBS) \
xptr = ARF_NOPTR_D(x); \
else \
xptr = ARF_PTR_D(x); \
} while (0)
/* Assumes non-special! */
#define ARF_GET_TOP_LIMB(lmb, x) \
do { \
mp_srcptr __xptr; \
mp_size_t __xn; \
ARF_GET_MPN_READONLY(__xptr, __xn, (x)); \
(lmb) = __xptr[__xn - 1]; \
} while (0)
/* Get mpn pointer xptr for writing *exactly* xn limbs to x. */
#define ARF_GET_MPN_WRITE(xptr, xn, x) \
do { \
mp_size_t __xn = (xn); \
if ((__xn) <= ARF_NOPTR_LIMBS) \
{ \
ARF_DEMOTE(x); \
xptr = ARF_NOPTR_D(x); \
} \
else \
{ \
if (!ARF_HAS_PTR(x)) \
{ \
_arf_promote(x, __xn); \
} \
else if (ARF_PTR_ALLOC(x) < (__xn)) \
{ \
ARF_PTR_D(x) = (mp_ptr) \
flint_realloc(ARF_PTR_D(x), \
(xn) * sizeof(mp_limb_t)); \
ARF_PTR_ALLOC(x) = (__xn); \
} \
xptr = ARF_PTR_D(x); \
} \
ARF_XSIZE(x) = ARF_MAKE_XSIZE(__xn, 0); \
} while (0)
ARF_INLINE void
arf_init(arf_t x)
{
fmpz_init(ARF_EXPREF(x));
ARF_XSIZE(x) = 0;
}
void arf_clear(arf_t x);
ARF_INLINE void
arf_zero(arf_t x)
{
ARF_MAKE_SPECIAL(x);
ARF_EXP(x) = ARF_EXP_ZERO;
}
ARF_INLINE void
arf_pos_inf(arf_t x)
{
ARF_MAKE_SPECIAL(x);
ARF_EXP(x) = ARF_EXP_POS_INF;
}
ARF_INLINE void
arf_neg_inf(arf_t x)
{
ARF_MAKE_SPECIAL(x);
ARF_EXP(x) = ARF_EXP_NEG_INF;
}
ARF_INLINE void
arf_nan(arf_t x)
{
ARF_MAKE_SPECIAL(x);
ARF_EXP(x) = ARF_EXP_NAN;
}
ARF_INLINE int
arf_is_special(const arf_t x)
{
return ARF_IS_SPECIAL(x);
}
ARF_INLINE int
arf_is_zero(const arf_t x)
{
return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_ZERO);
}
ARF_INLINE int
arf_is_pos_inf(const arf_t x)
{
return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_POS_INF);
}
ARF_INLINE int
arf_is_neg_inf(const arf_t x)
{
return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_NEG_INF);
}
ARF_INLINE int
arf_is_nan(const arf_t x)
{
return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_NAN);
}
ARF_INLINE int
arf_is_normal(const arf_t x)
{
return !ARF_IS_SPECIAL(x);
}
ARF_INLINE int
arf_is_finite(const arf_t x)
{
return !ARF_IS_SPECIAL(x) || (ARF_EXP(x) == ARF_EXP_ZERO);
}
ARF_INLINE int
arf_is_inf(const arf_t x)
{
return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_POS_INF ||
ARF_EXP(x) == ARF_EXP_NEG_INF);
}
ARF_INLINE void
arf_one(arf_t x)
{
fmpz_clear(ARF_EXPREF(x));
ARF_DEMOTE(x);
ARF_EXP(x) = 1;
ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
ARF_NOPTR_D(x)[0] = LIMB_TOP;
}
ARF_INLINE int
arf_is_one(const arf_t x)
{
return (ARF_EXP(x) == 1) && (ARF_XSIZE(x) == ARF_MAKE_XSIZE(1, 0))
&& ARF_NOPTR_D(x)[0] == LIMB_TOP;
}
ARF_INLINE int
arf_sgn(const arf_t x)
{
if (arf_is_special(x))
{
if (arf_is_zero(x) || arf_is_nan(x))
return 0;
return arf_is_pos_inf(x) ? 1 : -1;
}
else
{
return ARF_SGNBIT(x) ? -1 : 1;
}
}
int arf_cmp(const arf_t x, const arf_t y);
int arf_cmpabs(const arf_t x, const arf_t y);
int arf_cmpabs_ui(const arf_t x, ulong y);
int arf_cmpabs_d(const arf_t x, double y);
int arf_cmp_si(const arf_t x, slong y);
int arf_cmp_ui(const arf_t x, ulong y);
int arf_cmp_d(const arf_t x, double y);
ARF_INLINE void
arf_swap(arf_t y, arf_t x)
{
if (x != y)
{
arf_struct t = *x;
*x = *y;
*y = t;
}
}
ARF_INLINE void
arf_set(arf_t y, const arf_t x)
{
if (x != y)
{
/* Fast path */
if (!COEFF_IS_MPZ(ARF_EXP(x)) && !COEFF_IS_MPZ(ARF_EXP(y)))
ARF_EXP(y) = ARF_EXP(x);
else
fmpz_set(ARF_EXPREF(y), ARF_EXPREF(x));
/* Fast path */
if (!ARF_HAS_PTR(x))
{
ARF_DEMOTE(y);
(y)->d = (x)->d;
}
else
{
mp_ptr yptr;
mp_srcptr xptr;
mp_size_t n;
ARF_GET_MPN_READONLY(xptr, n, x);
ARF_GET_MPN_WRITE(yptr, n, y);
flint_mpn_copyi(yptr, xptr, n);
}
/* Important. */
ARF_XSIZE(y) = ARF_XSIZE(x);
}
}
ARF_INLINE void
arf_neg(arf_t y, const arf_t x)
{
arf_set(y, x);
if (arf_is_special(y))
{
if (arf_is_pos_inf(y))
arf_neg_inf(y);
else if (arf_is_neg_inf(y))
arf_pos_inf(y);
}
else
{
ARF_NEG(y);
}
}
ARF_INLINE void
arf_init_set_ui(arf_t x, ulong v)
{
if (v == 0)
{
ARF_EXP(x) = ARF_EXP_ZERO;
ARF_XSIZE(x) = 0;
}
else
{
unsigned int c;
count_leading_zeros(c, v);
ARF_EXP(x) = FLINT_BITS - c;
ARF_NOPTR_D(x)[0] = v << c;
ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
}
}
/* FLINT_ABS is unsafe for x = WORD_MIN */
#define UI_ABS_SI(x) (((slong)(x) < 0) ? (-(ulong)(x)) : ((ulong)(x)))
ARF_INLINE void
arf_init_set_si(arf_t x, slong v)
{
arf_init_set_ui(x, UI_ABS_SI(v));
if (v < 0)
ARF_NEG(x);
}
ARF_INLINE void
arf_set_ui(arf_t x, ulong v)
{
ARF_DEMOTE(x);
_fmpz_demote(ARF_EXPREF(x));
if (v == 0)
{
ARF_EXP(x) = ARF_EXP_ZERO;
ARF_XSIZE(x) = 0;
}
else
{
unsigned int c;
count_leading_zeros(c, v);
ARF_EXP(x) = FLINT_BITS - c;
ARF_NOPTR_D(x)[0] = v << c;
ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
}
}
ARF_INLINE void
arf_set_si(arf_t x, slong v)
{
arf_set_ui(x, UI_ABS_SI(v));
if (v < 0)
ARF_NEG(x);
}
ARF_INLINE void
arf_init_set_shallow(arf_t z, const arf_t x)
{
*z = *x;
}
ARF_INLINE void
arf_init_neg_shallow(arf_t z, const arf_t x)
{
*z = *x;
arf_neg(z, z);
}
ARF_INLINE void
arf_init_set_mag_shallow(arf_t y, const mag_t x)
{
mp_limb_t t = MAG_MAN(x);
ARF_XSIZE(y) = ARF_MAKE_XSIZE(t != 0, 0);
ARF_EXP(y) = MAG_EXP(x);
ARF_NOPTR_D(y)[0] = t << (FLINT_BITS - MAG_BITS);
}
ARF_INLINE void
arf_init_neg_mag_shallow(arf_t z, const mag_t x)
{
arf_init_set_mag_shallow(z, x);
arf_neg(z, z);
}
ARF_INLINE int
arf_cmpabs_mag(const arf_t x, const mag_t y)
{
arf_t t;
arf_init_set_mag_shallow(t, y); /* no need to free */
return arf_cmpabs(x, t);
}
ARF_INLINE int
arf_mag_cmpabs(const mag_t x, const arf_t y)
{
arf_t t;
arf_init_set_mag_shallow(t, x); /* no need to free */
return arf_cmpabs(t, y);
}
/* Assumes xn > 0, x[xn-1] != 0. */
/* TBD: 1, 2 limb versions */
void arf_set_mpn(arf_t y, mp_srcptr x, mp_size_t xn, int sgnbit);
ARF_INLINE void
arf_set_mpz(arf_t y, const mpz_t x)
{
slong size = x->_mp_size;
if (size == 0)
arf_zero(y);
else
arf_set_mpn(y, x->_mp_d, FLINT_ABS(size), size < 0);
}
ARF_INLINE void
arf_set_fmpz(arf_t y, const fmpz_t x)
{
if (!COEFF_IS_MPZ(*x))
arf_set_si(y, *x);
else
arf_set_mpz(y, COEFF_TO_PTR(*x));
}
int _arf_set_round_ui(arf_t x, ulong v, int sgnbit, slong prec, arf_rnd_t rnd);
int _arf_set_round_uiui(arf_t z, slong * fix, mp_limb_t hi, mp_limb_t lo, int sgnbit, slong prec, arf_rnd_t rnd);
int
_arf_set_round_mpn(arf_t y, slong * exp_shift, mp_srcptr x, mp_size_t xn,
int sgnbit, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_set_round_ui(arf_t x, ulong v, slong prec, arf_rnd_t rnd)
{
return _arf_set_round_ui(x, v, 0, prec, rnd);
}
ARF_INLINE int
arf_set_round_si(arf_t x, slong v, slong prec, arf_rnd_t rnd)
{
return _arf_set_round_ui(x, UI_ABS_SI(v), v < 0, prec, rnd);
}
ARF_INLINE int
arf_set_round_mpz(arf_t y, const mpz_t x, slong prec, arf_rnd_t rnd)
{
int inexact;
slong size = x->_mp_size;
slong fix;
if (size == 0)
{
arf_zero(y);
return 0;
}
inexact = _arf_set_round_mpn(y, &fix, x->_mp_d, FLINT_ABS(size),
(size < 0), prec, rnd);
_fmpz_demote(ARF_EXPREF(y));
ARF_EXP(y) = FLINT_ABS(size) * FLINT_BITS + fix;
return inexact;
}
ARF_INLINE int
arf_set_round_fmpz(arf_t y, const fmpz_t x, slong prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*x))
return arf_set_round_si(y, *x, prec, rnd);
else
return arf_set_round_mpz(y, COEFF_TO_PTR(*x), prec, rnd);
}
int arf_set_round(arf_t y, const arf_t x, slong prec, arf_rnd_t rnd);
int arf_neg_round(arf_t y, const arf_t x, slong prec, arf_rnd_t rnd);
void arf_get_fmpr(fmpr_t y, const arf_t x);
void arf_set_fmpr(arf_t y, const fmpr_t x);
int arf_get_mpfr(mpfr_t x, const arf_t y, mpfr_rnd_t rnd);
void arf_set_mpfr(arf_t x, const mpfr_t y);
int arf_equal(const arf_t x, const arf_t y);
int arf_equal_si(const arf_t x, slong y);
int arf_equal_ui(const arf_t x, ulong y);
int arf_equal_d(const arf_t x, double y);
ARF_INLINE void
arf_min(arf_t z, const arf_t a, const arf_t b)
{
if (arf_cmp(a, b) <= 0)
arf_set(z, a);
else
arf_set(z, b);
}
ARF_INLINE void
arf_max(arf_t z, const arf_t a, const arf_t b)
{
if (arf_cmp(a, b) > 0)
arf_set(z, a);
else
arf_set(z, b);
}
ARF_INLINE void
arf_abs(arf_t y, const arf_t x)
{
if (arf_sgn(x) < 0)
arf_neg(y, x);
else
arf_set(y, x);
}
ARF_INLINE slong
arf_bits(const arf_t x)
{
if (arf_is_special(x))
return 0;
else
{
mp_srcptr xp;
mp_size_t xn;
slong c;
ARF_GET_MPN_READONLY(xp, xn, x);
count_trailing_zeros(c, xp[0]);
return xn * FLINT_BITS - c;
}
}
ARF_INLINE void
arf_bot(fmpz_t e, const arf_t x)
{
if (arf_is_special(x))
fmpz_zero(e);
else
fmpz_sub_si(e, ARF_EXPREF(x), arf_bits(x));
}
int arf_is_int(const arf_t x);
int arf_is_int_2exp_si(const arf_t x, slong e);
int arf_cmp_2exp_si(const arf_t x, slong e);
int arf_cmpabs_2exp_si(const arf_t x, slong e);
ARF_INLINE void
arf_set_si_2exp_si(arf_t x, slong man, slong exp)
{
arf_set_si(x, man);
if (man != 0)
fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
}
ARF_INLINE void
arf_set_ui_2exp_si(arf_t x, ulong man, slong exp)
{
arf_set_ui(x, man);
if (man != 0)
fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
}
ARF_INLINE void
arf_mul_2exp_si(arf_t y, const arf_t x, slong e)
{
arf_set(y, x);
if (!arf_is_special(y))
fmpz_add_si_inline(ARF_EXPREF(y), ARF_EXPREF(y), e);
}
ARF_INLINE void
arf_mul_2exp_fmpz(arf_t y, const arf_t x, const fmpz_t e)
{
arf_set(y, x);
if (!arf_is_special(y))
fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), e);
}
ARF_INLINE int
arf_set_round_fmpz_2exp(arf_t y, const fmpz_t x, const fmpz_t exp, slong prec, arf_rnd_t rnd)
{
if (fmpz_is_zero(x))
{
arf_zero(y);
return 0;
}
else
{
int r = arf_set_round_fmpz(y, x, prec, rnd);
fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), exp);
return r;
}
}
ARF_INLINE void
arf_abs_bound_lt_2exp_fmpz(fmpz_t b, const arf_t x)
{
if (arf_is_special(x))
fmpz_zero(b);
else
fmpz_set(b, ARF_EXPREF(x));
}
ARF_INLINE void
arf_abs_bound_le_2exp_fmpz(fmpz_t b, const arf_t x)
{
if (arf_is_special(x))
fmpz_zero(b);
else if (ARF_IS_POW2(x))
fmpz_sub_ui(b, ARF_EXPREF(x), 1);
else
fmpz_set(b, ARF_EXPREF(x));
}
slong arf_abs_bound_lt_2exp_si(const arf_t x);
void arf_frexp(arf_t man, fmpz_t exp, const arf_t x);
void arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x);
int _arf_get_integer_mpn(mp_ptr y, mp_srcptr x, mp_size_t xn, slong exp);
int _arf_set_mpn_fixed(arf_t z, mp_srcptr xp, mp_size_t xn,
mp_size_t fixn, int negative, slong prec, arf_rnd_t rnd);
int arf_get_fmpz(fmpz_t z, const arf_t x, arf_rnd_t rnd);
slong arf_get_si(const arf_t x, arf_rnd_t rnd);
int arf_get_fmpz_fixed_fmpz(fmpz_t y, const arf_t x, const fmpz_t e);
int arf_get_fmpz_fixed_si(fmpz_t y, const arf_t x, slong e);
ARF_INLINE void
arf_set_fmpz_2exp(arf_t x, const fmpz_t man, const fmpz_t exp)
{
arf_set_fmpz(x, man);
if (!arf_is_zero(x))
fmpz_add_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
}
void arf_floor(arf_t z, const arf_t x);
void arf_ceil(arf_t z, const arf_t x);
void arf_debug(const arf_t x);
char * arf_get_str(const arf_t x, slong d);
void arf_fprint(FILE * file, const arf_t x);
void arf_fprintd(FILE * file, const arf_t y, slong d);
ARF_INLINE void
arf_print(const arf_t x)
{
arf_fprint(stdout, x);
}
ARF_INLINE void
arf_printd(const arf_t y, slong d)
{
arf_fprintd(stdout, y, d);
}
void arf_randtest(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
void arf_randtest_not_zero(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
void arf_randtest_special(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
void arf_urandom(arf_t x, flint_rand_t state, slong bits, arf_rnd_t rnd);
#define MUL_MPFR_MIN_LIMBS 25
#define MUL_MPFR_MAX_LIMBS 10000
#define nn_mul_2x1(r2, r1, r0, a1, a0, b0) \
do { \
mp_limb_t t1; \
umul_ppmm(r1, r0, a0, b0); \
umul_ppmm(r2, t1, a1, b0); \
add_ssaaaa(r2, r1, r2, r1, 0, t1); \
} while (0)
#define nn_mul_2x2(r3, r2, r1, r0, a1, a0, b1, b0) \
do { \
mp_limb_t t1, t2, t3; \
umul_ppmm(r1, r0, a0, b0); \
umul_ppmm(r2, t1, a1, b0); \
add_ssaaaa(r2, r1, r2, r1, 0, t1); \
umul_ppmm(t1, t2, a0, b1); \
umul_ppmm(r3, t3, a1, b1); \
add_ssaaaa(r3, t1, r3, t1, 0, t3); \
add_ssaaaa(r2, r1, r2, r1, t1, t2); \
r3 += r2 < t1; \
} while (0)
#define ARF_USE_FFT_MUL(_xn) ((_xn) > 10000)
/* from flint/fft.h */
void flint_mpn_mul_fft_main(mp_ptr r1, mp_srcptr i1, mp_size_t n1, mp_srcptr i2, mp_size_t n2);
#define ARF_MPN_MUL(_z, _x, _xn, _y, _yn) \
if ((_xn) == (_yn)) \
{ \
if ((_xn) == 1) \
{ \
umul_ppmm((_z)[1], (_z)[0], (_x)[0], (_y)[0]); \
} \
else if ((_xn) == 2) \
{ \
mp_limb_t __arb_x1, __arb_x0, __arb_y1, __arb_y0; \
__arb_x0 = (_x)[0]; \
__arb_x1 = (_x)[1]; \
__arb_y0 = (_y)[0]; \
__arb_y1 = (_y)[1]; \
nn_mul_2x2((_z)[3], (_z)[2], (_z)[1], (_z)[0], __arb_x1, __arb_x0, __arb_y1, __arb_y0); \
} \
else if (ARF_USE_FFT_MUL(_xn)) \
flint_mpn_mul_fft_main((_z), (_x), (_xn), (_y), (_yn)); \
else if ((_x) == (_y)) \
mpn_sqr((_z), (_x), (_xn)); \
else \
mpn_mul_n((_z), (_x), (_y), (_xn)); \
} \
else if ((_xn) > (_yn)) \
{ \
if ((_yn) == 1) \
(_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_x), (_xn), (_y)[0]); \
else if (ARF_USE_FFT_MUL(_yn)) \
flint_mpn_mul_fft_main((_z), (_x), (_xn), (_y), (_yn)); \
else \
mpn_mul((_z), (_x), (_xn), (_y), (_yn)); \
} \
else \
{ \
if ((_xn) == 1) \
(_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_y), (_yn), (_x)[0]); \
else if (ARF_USE_FFT_MUL(_xn)) \
flint_mpn_mul_fft_main((_z), (_y), (_yn), (_x), (_xn)); \
else \
mpn_mul((_z), (_y), (_yn), (_x), (_xn)); \
}
#define ARF_MUL_STACK_ALLOC 40
#define ARF_MUL_TLS_ALLOC 1000
extern TLS_PREFIX mp_ptr __arf_mul_tmp;
extern TLS_PREFIX slong __arf_mul_alloc;
ARB_DLL extern void _arf_mul_tmp_cleanup(void);
#define ARF_MUL_TMP_DECL \
mp_limb_t tmp_stack[ARF_MUL_STACK_ALLOC]; \
#define ARF_MUL_TMP_ALLOC(tmp, alloc) \
if (alloc <= ARF_MUL_STACK_ALLOC) \
{ \
tmp = tmp_stack; \
} \
else if (alloc <= ARF_MUL_TLS_ALLOC) \
{ \
if (__arf_mul_alloc < alloc) \
{ \
if (__arf_mul_alloc == 0) \
{ \
flint_register_cleanup_function(_arf_mul_tmp_cleanup); \
} \
__arf_mul_tmp = flint_realloc(__arf_mul_tmp, sizeof(mp_limb_t) * alloc); \
__arf_mul_alloc = alloc; \
} \
tmp = __arf_mul_tmp; \
} \
else \
{ \
tmp = flint_malloc(sizeof(mp_limb_t) * alloc); \
}
#define ARF_MUL_TMP_FREE(tmp, alloc) \
if (alloc > ARF_MUL_TLS_ALLOC) \
flint_free(tmp);
void arf_mul_special(arf_t z, const arf_t x, const arf_t y);
int arf_mul_via_mpfr(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd);
int arf_mul_rnd_any(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
int arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec);
#define arf_mul(z, x, y, prec, rnd) \
((rnd == FMPR_RND_DOWN) \
? arf_mul_rnd_down(z, x, y, prec) \
: arf_mul_rnd_any(z, x, y, prec, rnd))
ARF_INLINE int
arf_neg_mul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)
{
if (arf_is_special(y))
{
arf_mul(z, x, y, prec, rnd);
arf_neg(z, z);
return 0;
}
else
{
arf_t t;
*t = *y;
ARF_NEG(t);
return arf_mul(z, x, t, prec, rnd);
}
}
ARF_INLINE int
arf_mul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_mul(z, x, t, prec, rnd);
}
ARF_INLINE int
arf_mul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_mul(z, x, t, prec, rnd);
}
int arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_mul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_mul_si(z, x, *y, prec, rnd);
else
return arf_mul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
#define ARF_ADD_STACK_ALLOC 40
#define ARF_ADD_TLS_ALLOC 1000
extern TLS_PREFIX mp_ptr __arf_add_tmp;
extern TLS_PREFIX slong __arf_add_alloc;
ARB_DLL extern void _arf_add_tmp_cleanup(void);
#define ARF_ADD_TMP_DECL \
mp_limb_t tmp_stack[ARF_ADD_STACK_ALLOC]; \
#define ARF_ADD_TMP_ALLOC(tmp, alloc) \
if (alloc <= ARF_ADD_STACK_ALLOC) \
{ \
tmp = tmp_stack; \
} \
else if (alloc <= ARF_ADD_TLS_ALLOC) \
{ \
if (__arf_add_alloc < alloc) \
{ \
if (__arf_add_alloc == 0) \
{ \
flint_register_cleanup_function(_arf_add_tmp_cleanup); \
} \
__arf_add_tmp = flint_realloc(__arf_add_tmp, sizeof(mp_limb_t) * alloc); \
__arf_add_alloc = alloc; \
} \
tmp = __arf_add_tmp; \
} \
else \
{ \
tmp = flint_malloc(sizeof(mp_limb_t) * alloc); \
}
#define ARF_ADD_TMP_FREE(tmp, alloc) \
if (alloc > ARF_ADD_TLS_ALLOC) \
flint_free(tmp);
int _arf_add_mpn(arf_t z, mp_srcptr xp, mp_size_t xn, int xsgnbit,
const fmpz_t xexp, mp_srcptr yp, mp_size_t yn, int ysgnbit,
flint_bitcnt_t shift, slong prec, arf_rnd_t rnd);
int arf_add(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
int arf_add_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd);
int arf_add_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd);
int arf_add_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd);
int arf_add_fmpz_2exp(arf_ptr z, arf_srcptr x, const fmpz_t y, const fmpz_t exp, slong prec, arf_rnd_t rnd);
int arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
int arf_sub_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd);
int arf_sub_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd);
int arf_sub_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd);
int arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_addmul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_addmul(z, x, t, prec, rnd);
}
ARF_INLINE int
arf_addmul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_addmul(z, x, t, prec, rnd);
}
int arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_addmul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_addmul_si(z, x, *y, prec, rnd);
else
return arf_addmul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
int arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_submul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_submul(z, x, t, prec, rnd);
}
ARF_INLINE int
arf_submul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_submul(z, x, t, prec, rnd);
}
int arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_submul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_submul_si(z, x, *y, prec, rnd);
else
return arf_submul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
int arf_fma(arf_ptr res, arf_srcptr x, arf_srcptr y, arf_srcptr z, slong prec, arf_rnd_t rnd);
int arf_sosq(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd);
int arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
ARF_INLINE int
arf_div_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_div(z, x, t, prec, rnd);
}
ARF_INLINE int
arf_ui_div(arf_ptr z, ulong x, arf_srcptr y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, x); /* no need to free */
return arf_div(z, t, y, prec, rnd);
}
ARF_INLINE int
arf_div_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_div(z, x, t, prec, rnd);
}
ARF_INLINE int
arf_si_div(arf_ptr z, slong x, arf_srcptr y, slong prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, x); /* no need to free */
return arf_div(z, t, y, prec, rnd);
}
ARF_INLINE int
arf_div_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
{
arf_t t;
int r;
arf_init(t);
arf_set_fmpz(t, y);
r = arf_div(z, x, t, prec, rnd);
arf_clear(t);
return r;
}
ARF_INLINE int
arf_fmpz_div(arf_ptr z, const fmpz_t x, arf_srcptr y, slong prec, arf_rnd_t rnd)
{
arf_t t;
int r;
arf_init(t);
arf_set_fmpz(t, x);
r = arf_div(z, t, y, prec, rnd);
arf_clear(t);
return r;
}
ARF_INLINE int
arf_fmpz_div_fmpz(arf_ptr z, const fmpz_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)
{
arf_t t, u;
int r;
arf_init(t);
arf_init(u);
arf_set_fmpz(t, x);
arf_set_fmpz(u, y);
r = arf_div(z, t, u, prec, rnd);
arf_clear(t);
arf_clear(u);
return r;
}
int arf_sqrt(arf_ptr z, arf_srcptr x, slong prec, arf_rnd_t rnd);
int arf_sqrt_ui(arf_t z, ulong x, slong prec, arf_rnd_t rnd);
int arf_sqrt_fmpz(arf_t z, const fmpz_t x, slong prec, arf_rnd_t rnd);
int arf_rsqrt(arf_ptr z, arf_srcptr x, slong prec, arf_rnd_t rnd);
int arf_root(arf_t z, const arf_t x, ulong k, slong prec, arf_rnd_t rnd);
/* Magnitude bounds */
void arf_get_mag(mag_t y, const arf_t x);
void arf_get_mag_lower(mag_t y, const arf_t x);
ARF_INLINE void
arf_set_mag(arf_t y, const mag_t x)
{
if (mag_is_zero(x))
{
arf_zero(y);
}
else if (mag_is_inf(x))
{
arf_pos_inf(y);
}
else
{
_fmpz_set_fast(ARF_EXPREF(y), MAG_EXPREF(x));
ARF_DEMOTE(y);
ARF_XSIZE(y) = ARF_MAKE_XSIZE(1, 0);
ARF_NOPTR_D(y)[0] = MAG_MAN(x) << (FLINT_BITS - MAG_BITS);
}
}
ARF_INLINE void
mag_init_set_arf(mag_t y, const arf_t x)
{
mag_init(y);
arf_get_mag(y, x);
}
ARF_INLINE void
mag_fast_init_set_arf(mag_t y, const arf_t x)
{
if (ARF_IS_SPECIAL(x)) /* x == 0 */
{
mag_fast_zero(y);
}
else
{
mp_srcptr xp;
mp_size_t xn;
ARF_GET_MPN_READONLY(xp, xn, x);
MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE;
MAG_EXP(y) = ARF_EXP(x);
MAG_FAST_ADJUST_ONE_TOO_LARGE(y);
}
}
ARF_INLINE void
arf_mag_fast_add_ulp(mag_t z, const mag_t x, const arf_t y, slong prec)
{
mag_fast_add_2exp_si(z, x, ARF_EXP(y) - prec);
}
ARF_INLINE void
arf_mag_add_ulp(mag_t z, const mag_t x, const arf_t y, slong prec)
{
if (ARF_IS_SPECIAL(y))
{
flint_printf("error: ulp error not defined for special value!\n");
flint_abort();
}
else if (MAG_IS_LAGOM(z) && MAG_IS_LAGOM(x) && ARF_IS_LAGOM(y))
{
arf_mag_fast_add_ulp(z, x, y, prec);
}
else
{
fmpz_t e;
fmpz_init(e);
fmpz_sub_ui(e, ARF_EXPREF(y), prec);
mag_add_2exp_fmpz(z, x, e);
fmpz_clear(e);
}
}
ARF_INLINE void
arf_mag_set_ulp(mag_t z, const arf_t y, slong prec)
{
if (ARF_IS_SPECIAL(y))
{
flint_printf("error: ulp error not defined for special value!\n");
flint_abort();
}
else
{
_fmpz_add_fast(MAG_EXPREF(z), ARF_EXPREF(y), 1 - prec);
MAG_MAN(z) = MAG_ONE_HALF;
}
}
void arf_get_fmpq(fmpq_t y, const arf_t x);
ARF_INLINE int
arf_set_fmpq(arf_t y, const fmpq_t x, slong prec, arf_rnd_t rnd)
{
return arf_fmpz_div_fmpz(y, fmpq_numref(x), fmpq_denref(x), prec, rnd);
}
int arf_complex_mul(arf_t e, arf_t f, const arf_t a, const arf_t b,
const arf_t c, const arf_t d,
slong prec, arf_rnd_t rnd);
int arf_complex_mul_fallback(arf_t e, arf_t f,
const arf_t a, const arf_t b,
const arf_t c, const arf_t d,
slong prec, arf_rnd_t rnd);
int arf_complex_sqr(arf_t e, arf_t f, const arf_t a, const arf_t b,
slong prec, arf_rnd_t rnd);
int arf_sum(arf_t s, arf_srcptr terms, slong len, slong prec, arf_rnd_t rnd);
double arf_get_d(const arf_t x, arf_rnd_t rnd);
void arf_set_d(arf_t x, double v);
ARF_INLINE slong
arf_allocated_bytes(const arf_t x)
{
slong size = fmpz_allocated_bytes(ARF_EXPREF(x));
if (ARF_HAS_PTR(x))
size += ARF_PTR_ALLOC(x) * sizeof(mp_limb_t);
return size;
}
int arf_load_str(arf_t res, const char * data);
char * arf_dump_str(const arf_t x);
int arf_load_file(arf_t res, FILE *stream);
int arf_dump_file(FILE* stream, const arf_t x);
#ifdef __cplusplus
}
#endif
#endif
Reference in a new issue View git blame Copy permalink