include/flint/nfloat.h

0001 /*
0002     Copyright (C) 2024 Fredrik Johansson
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef NFLOAT_H
0013 #define NFLOAT_H
0014
0015 #ifdef NFLOAT_INLINES_C
0016 #define NFLOAT_INLINE
0017 #else
0018 #define NFLOAT_INLINE static inline
0019 #endif
0020
0021 #include "mpn_extras.h"
0022 #include "gr_types.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 /* The upper limit is so that temporary buffers are safe to allocate on
0029    the stack and so that simple operations like swapping are not too
0030    expensive compared to a pointer-and-size representation. */
0031 #if FLINT_BITS == 64
0032 #define NFLOAT_MIN_LIMBS 1
0033 #define NFLOAT_MAX_LIMBS 66
0034 #else
0035 #define NFLOAT_MIN_LIMBS 1
0036 #define NFLOAT_MAX_LIMBS 132
0037 #endif
0038
0039 /* Number of header limbs used to encode sign + exponent. We use a
0040    whole limb for the sign bit to avoid the overhead of bit fiddling. */
0041 #define NFLOAT_HEADER_LIMBS 2
0042 #define NFLOAT_EXP(x) (((slong *) x)[0])
0043 #define NFLOAT_SGNBIT(x) (((nn_ptr) x)[1])
0044 #define NFLOAT_D(x) (((nn_ptr) x) + NFLOAT_HEADER_LIMBS)
0045 #define NFLOAT_DATA(x) (((nn_ptr) x))
0046
0047 /* Limbs needed to hold a temporary element of any precision. */
0048 #define NFLOAT_MAX_ALLOC (NFLOAT_HEADER_LIMBS + NFLOAT_MAX_LIMBS)
0049
0050 /* The exponent limits should not be too close to the word boundary,
0051    so that we can add two exponents plus some adjustmenets
0052    and safely intercept overflow/underflow. */
0053 #define NFLOAT_MIN_EXP (WORD_MIN / 4)
0054 #define NFLOAT_MAX_EXP (WORD_MAX / 4)
0055
0056 /* Special values (including zero). */
0057 /* TODO: is it better to encode special values in the exponent or in
0058          the sign field? */
0059 /* TODO: do we want signed zero (optionally)? */
0060 #define NFLOAT_EXP_ZERO      WORD_MIN
0061 #define NFLOAT_EXP_POS_INF   (WORD_MIN+1)
0062 #define NFLOAT_EXP_NEG_INF   (WORD_MIN+2)
0063 #define NFLOAT_EXP_NAN       (WORD_MIN+3)
0064 #define NFLOAT_IS_SPECIAL(x) (NFLOAT_EXP(x) < NFLOAT_MIN_EXP)
0065 #define NFLOAT_IS_ZERO(x)    (NFLOAT_EXP(x) == NFLOAT_EXP_ZERO)
0066 #define NFLOAT_IS_POS_INF(x) (NFLOAT_EXP(x) == NFLOAT_EXP_POS_INF)
0067 #define NFLOAT_IS_NEG_INF(x) (NFLOAT_EXP(x) == NFLOAT_EXP_NEG_INF)
0068 #define NFLOAT_IS_INF(x)     (NFLOAT_IS_POS_INF(x) || NFLOAT_IS_NEG_INF(x))
0069 #define NFLOAT_IS_NAN(x)     (NFLOAT_EXP(x) == NFLOAT_EXP_NAN)
0070
0071 /* Context flags to indicate whether we want to allow special values
0072    or return GR_DOMAIN / GR_UNABLE. */
0073 #define NFLOAT_ALLOW_UNDERFLOW 1
0074 #define NFLOAT_ALLOW_INF       2
0075 #define NFLOAT_ALLOW_NAN       4
0076 /* Context flags for directed rounding (experimental) */
0077 #define NFLOAT_RND_FLOOR       8
0078 #define NFLOAT_RND_CEIL        16
0079
0080 typedef struct
0081 {
0082     slong nlimbs;
0083     int flags;
0084     int rnd;    /* Allow rounding modes? Currently unused. */
0085 }
0086 _nfloat_ctx_struct;
0087
0088 typedef void * nfloat_ptr;
0089 typedef const void * nfloat_srcptr;
0090
0091 #define NFLOAT_CTX(ctx) ((_nfloat_ctx_struct *)(ctx))
0092 #define NFLOAT_CTX_NLIMBS(ctx) (NFLOAT_CTX(ctx)->nlimbs)
0093 #define NFLOAT_CTX_PREC(ctx) ((NFLOAT_CTX(ctx)->nlimbs) * FLINT_BITS)
0094 #define NFLOAT_CTX_FLAGS(ctx) (NFLOAT_CTX(ctx)->flags)
0095 #define NFLOAT_CTX_RND(ctx) (NFLOAT_CTX(ctx)->rnd)
0096 #define NFLOAT_CTX_DATA_NLIMBS(ctx) (NFLOAT_CTX_NLIMBS(ctx) + NFLOAT_HEADER_LIMBS)
0097 #define NFLOAT_CTX_HAS_INF_NAN(ctx) ((NFLOAT_CTX_FLAGS(ctx) & (NFLOAT_ALLOW_INF | NFLOAT_ALLOW_NAN)) != 0)
0098 #define NFLOAT_CTX_HAS_DIRECTED_ROUNDING(ctx) ((NFLOAT_CTX_FLAGS(ctx) & (NFLOAT_RND_FLOOR | NFLOAT_RND_CEIL)) != 0)
0099
0100 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[64 / FLINT_BITS]; } nfloat64_struct;
0101 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[128 / FLINT_BITS]; } nfloat128_struct;
0102 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[192 / FLINT_BITS]; } nfloat192_struct;
0103 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[256 / FLINT_BITS]; } nfloat256_struct;
0104 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[384 / FLINT_BITS]; } nfloat384_struct;
0105 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[512 / FLINT_BITS]; } nfloat512_struct;
0106 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[1024 / FLINT_BITS]; } nfloat1024_struct;
0107 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[2048 / FLINT_BITS]; } nfloat2048_struct;
0108 typedef struct { ulong head[NFLOAT_HEADER_LIMBS]; ulong d[4096 / FLINT_BITS]; } nfloat4096_struct;
0109
0110 typedef nfloat64_struct nfloat64_t[1];
0111 typedef nfloat128_struct nfloat128_t[1];
0112 typedef nfloat192_struct nfloat192_t[1];
0113 typedef nfloat256_struct nfloat256_t[1];
0114 typedef nfloat384_struct nfloat384_t[1];
0115 typedef nfloat512_struct nfloat512_t[1];
0116 typedef nfloat1024_struct nfloat1024_t[1];
0117 typedef nfloat2048_struct nfloat2048_t[1];
0118 typedef nfloat4096_struct nfloat4096_t[1];
0119
0120 #define LIMB_MSB_IS_SET(n) ((slong) (n) < 0)
0121
0122 int nfloat_ctx_init(gr_ctx_t ctx, slong prec, int flags);
0123 int nfloat_ctx_write(gr_stream_t out, gr_ctx_t ctx);
0124
0125 NFLOAT_INLINE void
0126 nfloat_init(nfloat_ptr res, gr_ctx_t FLINT_UNUSED(ctx))
0127 {
0128     NFLOAT_EXP(res) = NFLOAT_EXP_ZERO;
0129     NFLOAT_SGNBIT(res) = 0;
0130 }
0131
0132 NFLOAT_INLINE void
0133 nfloat_clear(nfloat_ptr FLINT_UNUSED(res), gr_ctx_t FLINT_UNUSED(ctx))
0134 {
0135     FLINT_ASSERT(NFLOAT_IS_SPECIAL(res) || LIMB_MSB_IS_SET(NFLOAT_D(res)[NFLOAT_CTX_NLIMBS(ctx)-1]));
0136 }
0137
0138 void nfloat_swap(nfloat_ptr x, nfloat_ptr y, gr_ctx_t ctx);
0139 int nfloat_set(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0140
0141 truth_t nfloat_equal(nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0142
0143 NFLOAT_INLINE
0144 int nfloat_ctx_set_real_prec(gr_ctx_t FLINT_UNUSED(ctx), slong FLINT_UNUSED(prec))
0145 {
0146     return GR_UNABLE;
0147 }
0148
0149 NFLOAT_INLINE
0150 int nfloat_ctx_get_real_prec(slong * res, gr_ctx_t ctx)
0151 {
0152     *res = NFLOAT_CTX_PREC(ctx);
0153     return GR_SUCCESS;
0154 }
0155
0156 NFLOAT_INLINE int
0157 nfloat_zero(nfloat_ptr res, gr_ctx_t FLINT_UNUSED(ctx))
0158 {
0159     NFLOAT_EXP(res) = NFLOAT_EXP_ZERO;
0160     NFLOAT_SGNBIT(res) = 0;
0161     return GR_SUCCESS;
0162 }
0163
0164 int nfloat_pos_inf(nfloat_ptr res, gr_ctx_t ctx);
0165 int nfloat_neg_inf(nfloat_ptr res, gr_ctx_t ctx);
0166 int nfloat_nan(nfloat_ptr res, gr_ctx_t ctx);
0167
0168 int _nfloat_underflow(nfloat_ptr res, int sgnbit, gr_ctx_t ctx);
0169 int _nfloat_overflow(nfloat_ptr res, int sgnbit, gr_ctx_t ctx);
0170
0171 #define NFLOAT_HANDLE_UNDERFLOW(res, ctx) \
0172     do { \
0173         if (FLINT_UNLIKELY(NFLOAT_EXP(res) < NFLOAT_MIN_EXP)) \
0174             return _nfloat_underflow(res, NFLOAT_SGNBIT(res), ctx); \
0175     } while (0)
0176
0177 #define NFLOAT_HANDLE_OVERFLOW(res, ctx) \
0178     do { \
0179         if (FLINT_UNLIKELY(NFLOAT_EXP(res) < NFLOAT_MIN_EXP)) \
0180             return _nfloat_underflow(res, NFLOAT_SGNBIT(res), ctx); \
0181     } while (0)
0182
0183 #define NFLOAT_HANDLE_UNDERFLOW_OVERFLOW(res, ctx) \
0184     do { \
0185         NFLOAT_HANDLE_UNDERFLOW(res, ctx); \
0186         NFLOAT_HANDLE_OVERFLOW(res, ctx); \
0187     } while (0)
0188
0189 int nfloat_one(nfloat_ptr res, gr_ctx_t ctx);
0190 int nfloat_neg_one(nfloat_ptr res, gr_ctx_t ctx);
0191
0192 NFLOAT_INLINE truth_t
0193 nfloat_is_zero(nfloat_srcptr x, gr_ctx_t FLINT_UNUSED(ctx))
0194 {
0195     if (NFLOAT_IS_NAN(x))
0196         return T_UNKNOWN;
0197
0198     return NFLOAT_IS_ZERO(x) ? T_TRUE : T_FALSE;
0199 }
0200
0201 truth_t nfloat_is_one(nfloat_srcptr x, gr_ctx_t ctx);
0202 truth_t nfloat_is_neg_one(nfloat_srcptr x, gr_ctx_t ctx);
0203
0204 int nfloat_set_ui(nfloat_ptr res, ulong x, gr_ctx_t ctx);
0205 int nfloat_set_si(nfloat_ptr res, slong x, gr_ctx_t ctx);
0206
0207 /* Here exp is understood such that {x, xn} is a fraction in [0, 1). */
0208 int _nfloat_set_mpn_2exp(nfloat_ptr res, nn_srcptr x, slong xn, slong exp, int xsgnbit, gr_ctx_t ctx);
0209
0210 NFLOAT_INLINE int
0211 nfloat_set_mpn_2exp(nfloat_ptr res, nn_srcptr x, slong xn, slong exp, int xsgnbit, gr_ctx_t ctx)
0212 {
0213     while (xn != 0 && x[xn - 1] == 0)
0214     {
0215         xn--;
0216         exp -= FLINT_BITS;
0217     }
0218
0219     if (xn == 0)
0220         return nfloat_zero(res, ctx);
0221
0222     return _nfloat_set_mpn_2exp(res, x, xn, exp, xsgnbit, ctx);
0223 }
0224
0225 NFLOAT_INLINE int
0226 nfloat_1_set_2_2exp(nfloat_ptr res, ulong x1, ulong x0, slong exp, int xsgnbit, gr_ctx_t ctx)
0227 {
0228     slong norm;
0229
0230     if (x1 == 0)
0231     {
0232         if (x0 == 0)
0233             return nfloat_zero(res, ctx);
0234
0235         norm = flint_clz(x0);
0236         NFLOAT_EXP(res) = exp - FLINT_BITS - norm;
0237         NFLOAT_SGNBIT(res) = xsgnbit;
0238         NFLOAT_D(res)[0] = x0 << norm;
0239     }
0240     else if (LIMB_MSB_IS_SET(x1))
0241     {
0242         NFLOAT_EXP(res) = exp;
0243         NFLOAT_SGNBIT(res) = xsgnbit;
0244         NFLOAT_D(res)[0] = x1;
0245     }
0246     else
0247     {
0248         norm = flint_clz(x1);
0249         NFLOAT_EXP(res) = exp - norm;
0250         NFLOAT_SGNBIT(res) = xsgnbit;
0251         NFLOAT_D(res)[0] = (x1 << norm) | (x0 >> (FLINT_BITS - norm));
0252     }
0253
0254     NFLOAT_HANDLE_UNDERFLOW_OVERFLOW(res, ctx);
0255     return GR_SUCCESS;
0256 }
0257
0258 NFLOAT_INLINE int
0259 nfloat_1_set_3_2exp(nfloat_ptr res, ulong x2, ulong x1, ulong x0, slong exp, int xsgnbit, gr_ctx_t ctx)
0260 {
0261     if (x2 == 0)
0262         return nfloat_1_set_2_2exp(res, x1, x0, exp - FLINT_BITS, xsgnbit, ctx);
0263     else
0264         return nfloat_1_set_2_2exp(res, x2, x1, exp, xsgnbit, ctx);
0265 }
0266
0267 NFLOAT_INLINE int
0268 nfloat_2_set_3_2exp(nfloat_ptr res, ulong x2, ulong x1, ulong x0, slong exp, int xsgnbit, gr_ctx_t ctx)
0269 {
0270     slong norm;
0271
0272     if (x2 == 0)
0273     {
0274         if (x1 == 0)
0275         {
0276             if (x0 == 0)
0277                 return nfloat_zero(res, ctx);
0278
0279             norm = flint_clz(x0);
0280             exp = exp - 2 * FLINT_BITS - norm;
0281             x1 = 0;
0282             x2 = x0 << norm;
0283         }
0284         else if (LIMB_MSB_IS_SET(x1))
0285         {
0286             exp = exp - FLINT_BITS;
0287             x2 = x1;
0288             x1 = x0;
0289         }
0290         else
0291         {
0292             norm = flint_clz(x1);
0293             exp = exp - FLINT_BITS - norm;
0294             x2 = (x1 << norm) | (x0 >> (FLINT_BITS - norm));
0295             x1 = (x0 << norm);
0296         }
0297     }
0298     else if (!LIMB_MSB_IS_SET(x2))
0299     {
0300         norm = flint_clz(x2);
0301         exp = exp - norm;
0302         x2 = (x2 << norm) | (x1 >> (FLINT_BITS - norm));
0303         x1 = (x1 << norm) | (x0 >> (FLINT_BITS - norm));
0304     }
0305
0306     NFLOAT_EXP(res) = exp;
0307     NFLOAT_SGNBIT(res) = xsgnbit;
0308     NFLOAT_D(res)[0] = x1;
0309     NFLOAT_D(res)[1] = x2;
0310     NFLOAT_HANDLE_UNDERFLOW_OVERFLOW(res, ctx);
0311     return GR_SUCCESS;
0312 }
0313
0314 NFLOAT_INLINE int
0315 nfloat_2_set_4_2exp(nfloat_ptr res, ulong x3, ulong x2, ulong x1, ulong x0, slong exp, int xsgnbit, gr_ctx_t ctx)
0316 {
0317     if (x3 == 0)
0318         return nfloat_2_set_3_2exp(res, x2, x1, x0, exp - FLINT_BITS, xsgnbit, ctx);
0319     else
0320         return nfloat_2_set_3_2exp(res, x3, x2, x1, exp, xsgnbit, ctx);
0321 }
0322
0323
0324
0325 int nfloat_set_fmpz(nfloat_ptr res, const fmpz_t x, gr_ctx_t ctx);
0326
0327 #ifdef ARF_H
0328 int nfloat_set_arf(nfloat_ptr res, const arf_t x, gr_ctx_t ctx);
0329 int nfloat_get_arf(arf_t res, nfloat_srcptr x, gr_ctx_t ctx);
0330 #endif
0331
0332 int nfloat_get_d_2exp_si(double * mant, slong * expo, nfloat_srcptr x, gr_ctx_t ctx);
0333 int nfloat_get_fmpz(fmpz_t res, nfloat_srcptr x, gr_ctx_t ctx);
0334
0335 int nfloat_set_fmpq(nfloat_ptr res, const fmpq_t v, gr_ctx_t ctx);
0336 int nfloat_set_d(nfloat_ptr res, double x, gr_ctx_t ctx);
0337 int nfloat_set_str(nfloat_ptr res, const char * x, gr_ctx_t ctx);
0338 int nfloat_set_other(nfloat_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx);
0339
0340 int nfloat_write(gr_stream_t out, nfloat_srcptr x, gr_ctx_t ctx);
0341 int nfloat_randtest(nfloat_ptr res, flint_rand_t state, gr_ctx_t ctx);
0342
0343 int nfloat_neg(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0344 int nfloat_abs(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0345 int _nfloat_cmp(nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0346 int _nfloat_cmpabs(nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0347 int nfloat_cmp(int * res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0348 int nfloat_cmpabs(int * res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0349 int nfloat_sgn(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0350 int nfloat_im(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0351
0352 int _nfloat_add_1(nfloat_ptr res, ulong x0, slong xexp, int xsgnbit, ulong y0, slong delta, gr_ctx_t ctx);
0353 int _nfloat_sub_1(nfloat_ptr res, ulong x0, slong xexp, int xsgnbit, ulong y0, slong delta, gr_ctx_t ctx);
0354 int _nfloat_add_2(nfloat_ptr res, nn_srcptr xd, slong xexp, int xsgnbit, nn_srcptr yd, slong delta, gr_ctx_t ctx);
0355 int _nfloat_sub_2(nfloat_ptr res, nn_srcptr xd, slong xexp, int xsgnbit, nn_srcptr yd, slong delta, gr_ctx_t ctx);
0356 int _nfloat_add_3(nfloat_ptr res, nn_srcptr x, slong xexp, int xsgnbit, nn_srcptr y, slong delta, gr_ctx_t ctx);
0357 int _nfloat_sub_3(nfloat_ptr res, nn_srcptr x, slong xexp, int xsgnbit, nn_srcptr y, slong delta, gr_ctx_t ctx);
0358 int _nfloat_add_4(nfloat_ptr res, nn_srcptr x, slong xexp, int xsgnbit, nn_srcptr y, slong delta, gr_ctx_t ctx);
0359 int _nfloat_sub_4(nfloat_ptr res, nn_srcptr x, slong xexp, int xsgnbit, nn_srcptr y, slong delta, gr_ctx_t ctx);
0360 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);
0361 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);
0362
0363 int nfloat_add(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0364 int nfloat_sub(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0365 int nfloat_mul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0366 int nfloat_addmul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0367 int nfloat_submul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0368 int nfloat_sqr(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0369
0370 int nfloat_mul_2exp_si(nfloat_ptr res, nfloat_srcptr x, slong y, gr_ctx_t ctx);
0371
0372 int nfloat_inv(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0373 int nfloat_div(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0374 int nfloat_div_ui(nfloat_ptr res, nfloat_srcptr x, ulong y, gr_ctx_t ctx);
0375 int nfloat_div_si(nfloat_ptr res, nfloat_srcptr x, slong y, gr_ctx_t ctx);
0376
0377 int nfloat_sqrt(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0378 int nfloat_rsqrt(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0379
0380 int nfloat_floor(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0381 int nfloat_ceil(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0382 int nfloat_trunc(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0383 int nfloat_nint(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0384
0385 int nfloat_pi(nfloat_ptr res, gr_ctx_t ctx);
0386 int nfloat_pow(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, gr_ctx_t ctx);
0387 int nfloat_exp(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0388 int nfloat_expm1(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0389 int nfloat_log(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0390 int nfloat_log1p(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0391 int nfloat_sin(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0392 int nfloat_cos(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0393 int nfloat_tan(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0394 int nfloat_sinh(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0395 int nfloat_cosh(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0396 int nfloat_tanh(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0397 int nfloat_atan(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0398 int nfloat_gamma(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0399 int nfloat_zeta(nfloat_ptr res, nfloat_srcptr x, gr_ctx_t ctx);
0400
0401 void _nfloat_vec_init(nfloat_ptr res, slong len, gr_ctx_t ctx);
0402 void _nfloat_vec_clear(nfloat_ptr res, slong len, gr_ctx_t ctx);
0403 int _nfloat_vec_set(nfloat_ptr res, nfloat_srcptr x, slong len, gr_ctx_t ctx);
0404 int _nfloat_vec_zero(nfloat_ptr res, slong len, gr_ctx_t ctx);
0405
0406 int _nfloat_vec_add(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx);
0407 int _nfloat_vec_sub(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx);
0408 int _nfloat_vec_mul(nfloat_ptr res, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx);
0409 int _nfloat_vec_mul_scalar(nfloat_ptr res, nfloat_srcptr x, slong len, nfloat_srcptr y, gr_ctx_t ctx);
0410 int _nfloat_vec_addmul_scalar(nfloat_ptr res, nfloat_srcptr x, slong len, nfloat_srcptr y, gr_ctx_t ctx);
0411 int _nfloat_vec_submul_scalar(nfloat_ptr res, nfloat_srcptr x, slong len, nfloat_srcptr y, gr_ctx_t ctx);
0412
0413 int _nfloat_vec_dot(nfloat_ptr res, nfloat_srcptr initial, int subtract, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx);
0414 int _nfloat_vec_dot_rev(nfloat_ptr res, nfloat_srcptr initial, int subtract, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx);
0415
0416 int nfloat_mat_mul_fixed(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, slong max_extra_prec, gr_ctx_t ctx);
0417 int nfloat_mat_mul_block(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, slong min_block_size, gr_ctx_t ctx);
0418 int nfloat_mat_mul(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0419
0420 int nfloat_mat_nonsingular_solve_tril(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);
0421 int nfloat_mat_nonsingular_solve_triu(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);
0422 int nfloat_mat_lu(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx);
0423 int nfloat_mat_lq(gr_mat_t L, gr_mat_t Q, const gr_mat_t A, gr_ctx_t ctx);
0424
0425 /* Complex numbers */
0426 /* Note: we use the same context data for real and complex rings
0427    (only which_ring and sizeof_elem differ). This allows us to call
0428    nfloat methods on the real and imaginary parts without creating
0429    a temporary nfloat context object, as long as the nfloat methods
0430    don't call any generic gr methods internally.
0431 */
0432
0433 typedef nfloat_ptr nfloat_complex_ptr;
0434 typedef nfloat_srcptr nfloat_complex_srcptr;
0435
0436 int nfloat_complex_ctx_init(gr_ctx_t ctx, slong prec, int flags);
0437
0438 #define NFLOAT_COMPLEX_CTX_DATA_NLIMBS(ctx) (2 * NFLOAT_CTX_DATA_NLIMBS(ctx))
0439
0440 #define NFLOAT_COMPLEX_RE(ptr, ctx) (ptr)
0441 #define NFLOAT_COMPLEX_IM(ptr, ctx) ((nn_ptr) (ptr) + NFLOAT_CTX_DATA_NLIMBS(ctx))
0442
0443 #define NFLOAT_COMPLEX_IS_SPECIAL(x, ctx) (NFLOAT_IS_SPECIAL(NFLOAT_COMPLEX_RE(x, ctx)) || NFLOAT_IS_SPECIAL(NFLOAT_COMPLEX_IM(x, ctx)))
0444 #define NFLOAT_COMPLEX_IS_ZERO(x, ctx) (NFLOAT_IS_ZERO(NFLOAT_COMPLEX_RE(x, ctx)) && NFLOAT_IS_ZERO(NFLOAT_COMPLEX_IM(x, ctx)))
0445
0446
0447 NFLOAT_INLINE void
0448 nfloat_complex_init(nfloat_complex_ptr res, gr_ctx_t ctx)
0449 {
0450     nfloat_init(NFLOAT_COMPLEX_RE(res, ctx), ctx);
0451     nfloat_init(NFLOAT_COMPLEX_IM(res, ctx), ctx);
0452 }
0453
0454 NFLOAT_INLINE void
0455 nfloat_complex_clear(nfloat_complex_ptr FLINT_UNUSED(res), gr_ctx_t FLINT_UNUSED(ctx))
0456 {
0457     FLINT_ASSERT(NFLOAT_IS_SPECIAL(NFLOAT_COMPLEX_RE(res, ctx)) || LIMB_MSB_IS_SET(NFLOAT_D(NFLOAT_COMPLEX_RE(res, ctx))[NFLOAT_CTX_NLIMBS(ctx)-1]));
0458     FLINT_ASSERT(NFLOAT_IS_SPECIAL(NFLOAT_COMPLEX_IM(res, ctx)) || LIMB_MSB_IS_SET(NFLOAT_D(NFLOAT_COMPLEX_IM(res, ctx))[NFLOAT_CTX_NLIMBS(ctx)-1]));
0459 }
0460
0461 NFLOAT_INLINE int
0462 nfloat_complex_zero(nfloat_complex_ptr res, gr_ctx_t ctx)
0463 {
0464     nfloat_zero(NFLOAT_COMPLEX_RE(res, ctx), ctx);
0465     nfloat_zero(NFLOAT_COMPLEX_IM(res, ctx), ctx);
0466     return GR_SUCCESS;
0467 }
0468
0469 #ifdef ACF_H
0470 int nfloat_complex_get_acf(acf_t res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0471 int nfloat_complex_set_acf(nfloat_complex_ptr res, const acf_t x, gr_ctx_t ctx);
0472 #endif
0473
0474 #ifdef ACB_H
0475 int nfloat_complex_get_acb(acb_t res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0476 int nfloat_complex_set_acb(nfloat_complex_ptr res, const acb_t x, gr_ctx_t ctx);
0477 #endif
0478
0479 int nfloat_complex_write(gr_stream_t out, nfloat_complex_srcptr x, gr_ctx_t ctx);
0480 int nfloat_complex_randtest(nfloat_complex_ptr res, flint_rand_t state, gr_ctx_t ctx);
0481
0482 void nfloat_complex_swap(nfloat_complex_ptr x, nfloat_complex_ptr y, gr_ctx_t ctx);
0483 int nfloat_complex_set(nfloat_complex_ptr res, nfloat_complex_ptr x, gr_ctx_t ctx);
0484 int nfloat_complex_one(nfloat_complex_ptr res, gr_ctx_t ctx);
0485 int nfloat_complex_neg_one(nfloat_complex_ptr res, gr_ctx_t ctx);
0486 truth_t nfloat_complex_is_zero(nfloat_complex_srcptr x, gr_ctx_t ctx);
0487 truth_t nfloat_complex_is_one(nfloat_complex_srcptr x, gr_ctx_t ctx);
0488 truth_t nfloat_complex_is_neg_one(nfloat_complex_srcptr x, gr_ctx_t ctx);
0489 int nfloat_complex_i(nfloat_complex_ptr res, gr_ctx_t ctx);
0490 int nfloat_complex_pi(nfloat_complex_ptr res, gr_ctx_t ctx);
0491 int nfloat_complex_conj(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0492 int nfloat_complex_re(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0493 int nfloat_complex_im(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0494 truth_t nfloat_complex_equal(nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0495 int nfloat_complex_set_si(nfloat_complex_ptr res, slong x, gr_ctx_t ctx);
0496 int nfloat_complex_set_ui(nfloat_complex_ptr res, ulong x, gr_ctx_t ctx);
0497 int nfloat_complex_set_fmpz(nfloat_complex_ptr res, const fmpz_t x, gr_ctx_t ctx);
0498 int nfloat_complex_set_fmpq(nfloat_complex_ptr res, const fmpq_t x, gr_ctx_t ctx);
0499 int nfloat_complex_set_d(nfloat_complex_ptr res, double x, gr_ctx_t ctx);
0500 int nfloat_complex_set_other(nfloat_complex_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx);
0501 int nfloat_complex_neg(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0502 int nfloat_complex_add(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0503 int nfloat_complex_sub(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0504 int _nfloat_complex_sqr_naive(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx);
0505 int _nfloat_complex_sqr_standard(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx);
0506 int _nfloat_complex_sqr_karatsuba(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx);
0507 int _nfloat_complex_sqr(nfloat_ptr res1, nfloat_ptr res2, nfloat_srcptr a, nfloat_srcptr b, gr_ctx_t ctx);
0508 int nfloat_complex_sqr(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0509 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);
0510 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);
0511 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);
0512 int nfloat_complex_mul(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0513 int nfloat_complex_inv(nfloat_complex_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0514 int nfloat_complex_div(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0515 int nfloat_complex_mul_2exp_si(nfloat_complex_ptr res, nfloat_complex_srcptr x, slong y, gr_ctx_t ctx);
0516 int nfloat_complex_cmp(int * res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0517 int nfloat_complex_cmpabs(int * res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, gr_ctx_t ctx);
0518 int nfloat_complex_abs(nfloat_ptr res, nfloat_complex_srcptr x, gr_ctx_t ctx);
0519
0520 void _nfloat_complex_vec_init(nfloat_complex_ptr res, slong len, gr_ctx_t ctx);
0521 void _nfloat_complex_vec_clear(nfloat_complex_ptr res, slong len, gr_ctx_t ctx);
0522 int _nfloat_complex_vec_zero(nfloat_complex_ptr res, slong len, gr_ctx_t ctx);
0523 int _nfloat_complex_vec_set(nfloat_complex_ptr res, nfloat_complex_srcptr x, slong len, gr_ctx_t ctx);
0524 int _nfloat_complex_vec_add(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx);
0525 int _nfloat_complex_vec_sub(nfloat_complex_ptr res, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx);
0526
0527 int _nfloat_complex_vec_dot(nfloat_complex_ptr res, nfloat_complex_srcptr initial, int subtract, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx);
0528 int _nfloat_complex_vec_dot_rev(nfloat_complex_ptr res, nfloat_complex_srcptr initial, int subtract, nfloat_complex_srcptr x, nfloat_complex_srcptr y, slong len, gr_ctx_t ctx);
0529
0530 int nfloat_complex_mat_mul_fixed(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, slong max_extra_prec, gr_ctx_t ctx);
0531 int nfloat_complex_mat_mul_block(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, slong min_block_size, gr_ctx_t ctx);
0532 int nfloat_complex_mat_mul_reorder(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0533 int nfloat_complex_mat_mul(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0534
0535 int nfloat_complex_mat_nonsingular_solve_tril(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);
0536 int nfloat_complex_mat_nonsingular_solve_triu(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);
0537 int nfloat_complex_mat_lu(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx);
0538
0539 /* Fixed-point arithmetic */
0540
0541 void _nfixed_print(nn_srcptr x, slong nlimbs, slong exp);
0542 void _nfixed_vec_add(nn_ptr res, nn_srcptr a, nn_srcptr b, slong len, slong nlimbs);
0543 void _nfixed_vec_sub(nn_ptr res, nn_srcptr a, nn_srcptr b, slong len, slong nlimbs);
0544
0545 void _nfixed_dot_2(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0546 void _nfixed_dot_3(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0547 void _nfixed_dot_4(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0548 void _nfixed_dot_5(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0549 void _nfixed_dot_6(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0550 void _nfixed_dot_7(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0551 void _nfixed_dot_8(nn_ptr res, nn_srcptr x, slong xstride, nn_srcptr y, slong ystride, slong len);
0552
0553 void _nfixed_mat_mul_classical_precise(nn_ptr C, nn_srcptr A, nn_srcptr B, slong m, slong n, slong p, slong nlimbs);
0554 void _nfixed_mat_mul_classical(nn_ptr C, nn_srcptr A, nn_srcptr B, slong m, slong n, slong p, slong nlimbs);
0555 void _nfixed_mat_mul_waksman(nn_ptr C, nn_srcptr A, nn_srcptr B, slong m, slong n, slong p, slong nlimbs);
0556 void _nfixed_mat_mul_strassen(nn_ptr C, nn_srcptr A, nn_srcptr B, slong m, slong n, slong p, slong cutoff, slong nlimbs);
0557 void _nfixed_mat_mul(nn_ptr C, nn_srcptr A, nn_srcptr B, slong m, slong n, slong p, slong nlimbs);
0558
0559 void _nfixed_mat_mul_bound_classical(double * bound, double * error, slong m, slong n, slong p, double A, double B, slong nlimbs);
0560 void _nfixed_mat_mul_bound_waksman(double * bound, double * error, slong m, slong n, slong p, double A, double B, slong nlimbs);
0561 void _nfixed_mat_mul_bound_strassen(double * bound, double * error, slong m, slong n, slong p, double A, double B, slong cutoff, slong nlimbs);
0562 void _nfixed_mat_mul_bound(double * bound, double * error, slong m, slong n, slong p, double A, double B, slong nlimbs);
0563 void _nfixed_complex_mat_mul_bound(double * bound, double * error, slong m, slong n, slong p, double A, double B, double C, double D, slong nlimbs);
0564
0565
0566 #ifdef __cplusplus
0567 }
0568 #endif
0569
0570 #endif