Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-11-07 10:15:08 

0001 /*
0002     Copyright (C) 2014 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 ARF_H
0013 #define ARF_H
0014 
0015 #ifdef ARF_INLINES_C
0016 #define ARF_INLINE
0017 #else
0018 #define ARF_INLINE static inline
0019 #endif
0020 
0021 #include <gmp.h>
0022 #include "longlong.h"
0023 #include "mag.h"
0024 #include "arf_types.h"
0025 
0026 #ifdef __cplusplus
0027 extern "C" {
0028 #endif
0029 
0030 typedef enum
0031 {
0032     ARF_RND_DOWN  = 0,
0033     ARF_RND_UP    = 1,
0034     ARF_RND_FLOOR = 2,
0035     ARF_RND_CEIL  = 3,
0036     ARF_RND_NEAR  = 4
0037 }
0038 arf_rnd_t;
0039 
0040 ARF_INLINE int
0041 arf_rounds_down(arf_rnd_t rnd, int sgnbit)
0042 {
0043     if (rnd == ARF_RND_DOWN) return 1;
0044     if (rnd == ARF_RND_UP) return 0;
0045     if (rnd == ARF_RND_FLOOR) return !sgnbit;
0046     return sgnbit;
0047 }
0048 
0049 ARF_INLINE int
0050 arf_rounds_up(arf_rnd_t rnd, int sgnbit)
0051 {
0052     if (rnd == ARF_RND_DOWN) return 0;
0053     if (rnd == ARF_RND_UP) return 1;
0054     if (rnd == ARF_RND_FLOOR) return sgnbit;
0055     return !sgnbit;
0056 }
0057 
0058 #ifdef __MPFR_H
0059 ARF_INLINE mpfr_rnd_t
0060 arf_rnd_to_mpfr(arf_rnd_t rnd)
0061 {
0062     if (rnd == ARF_RND_DOWN) return MPFR_RNDZ;
0063     else if (rnd == ARF_RND_UP) return MPFR_RNDA;
0064     else if (rnd == ARF_RND_FLOOR) return MPFR_RNDD;
0065     else if (rnd == ARF_RND_CEIL) return MPFR_RNDU;
0066     else return MPFR_RNDN;
0067 }
0068 #endif
0069 
0070 /* Allow 'infinite' precision for operations where a result can be computed exactly. */
0071 #define ARF_PREC_EXACT WORD_MAX
0072 
0073 #define ARF_PREC_ADD(prec,extra) ((prec) == ARF_PREC_EXACT ? ARF_PREC_EXACT : (prec) + (extra))
0074 
0075 #define ARF_RESULT_EXACT 0
0076 #define ARF_RESULT_INEXACT 1
0077 
0078 /* Range where we can skip fmpz overflow checks for exponent manipulation. */
0079 #define ARF_MAX_LAGOM_EXP MAG_MAX_LAGOM_EXP
0080 #define ARF_MIN_LAGOM_EXP MAG_MIN_LAGOM_EXP
0081 
0082 /* Exponent values used to encode special values. */
0083 #define ARF_EXP_ZERO 0
0084 #define ARF_EXP_NAN COEFF_MIN
0085 #define ARF_EXP_POS_INF (COEFF_MIN+1)
0086 #define ARF_EXP_NEG_INF (COEFF_MIN+2)
0087 
0088 /* Direct access to the exponent. */
0089 #define ARF_EXP(x) ((x)->exp)
0090 #define ARF_EXPREF(x) (&(x)->exp)
0091 
0092 /* Finite and with lagom big exponents. */
0093 #define ARF_IS_LAGOM(x) (ARF_EXP(x) >= ARF_MIN_LAGOM_EXP && \
0094                          ARF_EXP(x) <= ARF_MAX_LAGOM_EXP)
0095 
0096 /* More than two limbs (needs pointer). */
0097 #define ARF_HAS_PTR(x) ((x)->size > ((2 << 1) + 1))
0098 
0099 /* Raw size field (encodes both limb size and sign). */
0100 #define ARF_XSIZE(x) ((x)->size)
0101 
0102 /* Construct size field value from size in limbs and sign bit. */
0103 #define ARF_MAKE_XSIZE(size, sgnbit) ((((slong) size) << 1) | (sgnbit))
0104 
0105 /* The limb size, and the sign bit. */
0106 #define ARF_SIZE(x) (ARF_XSIZE(x) >> 1)
0107 #define ARF_SGNBIT(x) (ARF_XSIZE(x) & 1)
0108 
0109 /* Assumes non-special value */
0110 #define ARF_NEG(x) (ARF_XSIZE(x) ^= 1)
0111 
0112 /* Direct access to the limb data. */
0113 #define ARF_NOPTR_D(x)   ((x)->d.noptr.d)
0114 #define ARF_PTR_D(x)     ((x)->d.ptr.d)
0115 #define ARF_PTR_ALLOC(x) ((x)->d.ptr.alloc)
0116 
0117 /* Encoding for special values. */
0118 #define ARF_IS_SPECIAL(x) (ARF_XSIZE(x) == 0)
0119 
0120 /* Value is +/- a power of two */
0121 #define ARF_IS_POW2(x) (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP)
0122 
0123 /* To set a special value, first call this and then set the exponent. */
0124 #define ARF_MAKE_SPECIAL(x)         \
0125     do {                            \
0126         fmpz_clear(ARF_EXPREF(x));  \
0127         ARF_DEMOTE(x);              \
0128         ARF_XSIZE(x) = 0;           \
0129     } while (0)
0130 
0131 
0132 void _arf_promote(arf_t x, slong n);
0133 
0134 void _arf_demote(arf_t x);
0135 
0136 
0137 /* Warning: does not set size! -- also doesn't demote exponent. */
0138 #define ARF_DEMOTE(x)                 \
0139     do {                              \
0140         if (ARF_HAS_PTR(x))           \
0141             _arf_demote(x);           \
0142     } while (0)
0143 
0144 /* Get mpn pointer and size (xptr, xn) for read-only use. */
0145 #define ARF_GET_MPN_READONLY(xptr, xn, x)   \
0146     do {                                    \
0147         xn = ARF_SIZE(x);                   \
0148         if (xn <= ARF_NOPTR_LIMBS)          \
0149             xptr = ARF_NOPTR_D(x);          \
0150         else                                \
0151             xptr = ARF_PTR_D(x);            \
0152     } while (0)
0153 
0154 /* Assumes non-special! */
0155 #define ARF_GET_TOP_LIMB(lmb, x)                     \
0156     do {                                             \
0157         nn_srcptr __xptr;                            \
0158         slong __xn;                              \
0159         ARF_GET_MPN_READONLY(__xptr, __xn, (x));     \
0160         (lmb) = __xptr[__xn - 1];                    \
0161     } while (0)
0162 
0163 /* Get mpn pointer xptr for writing *exactly* xn limbs to x. */
0164 #define ARF_GET_MPN_WRITE(xptr, xn, x)                      \
0165     do {                                                    \
0166         slong __xn = (xn);                              \
0167         if ((__xn) <= ARF_NOPTR_LIMBS)                      \
0168         {                                                   \
0169             ARF_DEMOTE(x);                                  \
0170             xptr = ARF_NOPTR_D(x);                          \
0171         }                                                   \
0172         else                                                \
0173         {                                                   \
0174             if (!ARF_HAS_PTR(x))                            \
0175             {                                               \
0176                 _arf_promote(x, __xn);                      \
0177             }                                               \
0178             else if (ARF_PTR_ALLOC(x) < (__xn))             \
0179             {                                               \
0180                 ARF_PTR_D(x) = (nn_ptr)                     \
0181                         flint_realloc(ARF_PTR_D(x),         \
0182                         (xn) * sizeof(ulong));          \
0183                 ARF_PTR_ALLOC(x) = (__xn);                  \
0184             }                                               \
0185             xptr = ARF_PTR_D(x);                            \
0186         }                                                   \
0187         ARF_XSIZE(x) = ARF_MAKE_XSIZE(__xn, 0);             \
0188     } while (0)
0189 
0190 ARF_INLINE void
0191 arf_init(arf_t x)
0192 {
0193     fmpz_init(ARF_EXPREF(x));
0194     ARF_XSIZE(x) = 0;
0195 }
0196 
0197 void arf_clear(arf_t x);
0198 
0199 arf_ptr _arf_vec_init(slong n);
0200 void _arf_vec_clear(arf_ptr v, slong n);
0201 
0202 ARF_INLINE void
0203 arf_zero(arf_t x)
0204 {
0205     ARF_MAKE_SPECIAL(x);
0206     ARF_EXP(x) = ARF_EXP_ZERO;
0207 }
0208 
0209 ARF_INLINE void
0210 arf_pos_inf(arf_t x)
0211 {
0212     ARF_MAKE_SPECIAL(x);
0213     ARF_EXP(x) = ARF_EXP_POS_INF;
0214 }
0215 
0216 ARF_INLINE void
0217 arf_neg_inf(arf_t x)
0218 {
0219     ARF_MAKE_SPECIAL(x);
0220     ARF_EXP(x) = ARF_EXP_NEG_INF;
0221 }
0222 
0223 ARF_INLINE void
0224 arf_nan(arf_t x)
0225 {
0226     ARF_MAKE_SPECIAL(x);
0227     ARF_EXP(x) = ARF_EXP_NAN;
0228 }
0229 
0230 ARF_INLINE int
0231 arf_is_special(const arf_t x)
0232 {
0233     return ARF_IS_SPECIAL(x);
0234 }
0235 
0236 ARF_INLINE int
0237 arf_is_zero(const arf_t x)
0238 {
0239     return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_ZERO);
0240 }
0241 
0242 ARF_INLINE int
0243 arf_is_pos_inf(const arf_t x)
0244 {
0245     return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_POS_INF);
0246 }
0247 
0248 ARF_INLINE int
0249 arf_is_neg_inf(const arf_t x)
0250 {
0251     return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_NEG_INF);
0252 }
0253 
0254 ARF_INLINE int
0255 arf_is_nan(const arf_t x)
0256 {
0257     return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_NAN);
0258 }
0259 
0260 ARF_INLINE int
0261 arf_is_normal(const arf_t x)
0262 {
0263     return !ARF_IS_SPECIAL(x);
0264 }
0265 
0266 ARF_INLINE int
0267 arf_is_finite(const arf_t x)
0268 {
0269     return !ARF_IS_SPECIAL(x) || (ARF_EXP(x) == ARF_EXP_ZERO);
0270 }
0271 
0272 ARF_INLINE int
0273 arf_is_inf(const arf_t x)
0274 {
0275     return ARF_IS_SPECIAL(x) && (ARF_EXP(x) == ARF_EXP_POS_INF ||
0276                                  ARF_EXP(x) == ARF_EXP_NEG_INF);
0277 }
0278 
0279 ARF_INLINE void
0280 arf_one(arf_t x)
0281 {
0282     fmpz_clear(ARF_EXPREF(x));
0283     ARF_DEMOTE(x);
0284     ARF_EXP(x) = 1;
0285     ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
0286     ARF_NOPTR_D(x)[0] = LIMB_TOP;
0287 }
0288 
0289 ARF_INLINE int
0290 arf_is_one(const arf_t x)
0291 {
0292     return (ARF_EXP(x) == 1) && (ARF_XSIZE(x) == ARF_MAKE_XSIZE(1, 0))
0293                              && ARF_NOPTR_D(x)[0] == LIMB_TOP;
0294 }
0295 
0296 ARF_INLINE int
0297 arf_sgn(const arf_t x)
0298 {
0299     if (arf_is_special(x))
0300     {
0301         if (arf_is_zero(x) || arf_is_nan(x))
0302             return 0;
0303         return arf_is_pos_inf(x) ? 1 : -1;
0304     }
0305     else
0306     {
0307         return ARF_SGNBIT(x) ? -1 : 1;
0308     }
0309 }
0310 
0311 int arf_cmp(const arf_t x, const arf_t y);
0312 
0313 int arf_cmpabs(const arf_t x, const arf_t y);
0314 
0315 int arf_cmpabs_ui(const arf_t x, ulong y);
0316 
0317 int arf_cmpabs_d(const arf_t x, double y);
0318 
0319 int arf_cmp_si(const arf_t x, slong y);
0320 
0321 int arf_cmp_ui(const arf_t x, ulong y);
0322 
0323 int arf_cmp_d(const arf_t x, double y);
0324 
0325 ARF_INLINE void
0326 arf_swap(arf_t y, arf_t x)
0327 {
0328     FLINT_SWAP(arf_struct, *x, *y);
0329 }
0330 
0331 void arf_set(arf_t y, const arf_t x);
0332 
0333 ARF_INLINE void
0334 arf_neg(arf_t y, const arf_t x)
0335 {
0336     arf_set(y, x);
0337 
0338     if (arf_is_special(y))
0339     {
0340         if (arf_is_pos_inf(y))
0341             arf_neg_inf(y);
0342         else if (arf_is_neg_inf(y))
0343             arf_pos_inf(y);
0344     }
0345     else
0346     {
0347         ARF_NEG(y);
0348     }
0349 }
0350 
0351 ARF_INLINE void
0352 arf_init_set_ui(arf_t x, ulong v)
0353 {
0354     if (v == 0)
0355     {
0356         ARF_EXP(x) = ARF_EXP_ZERO;
0357         ARF_XSIZE(x) = 0;
0358     }
0359     else
0360     {
0361         unsigned int c;
0362         c = flint_clz(v);
0363         ARF_EXP(x) = FLINT_BITS - c;
0364         ARF_NOPTR_D(x)[0] = v << c;
0365         ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
0366     }
0367 }
0368 
0369 ARF_INLINE void
0370 arf_init_set_si(arf_t x, slong v)
0371 {
0372     arf_init_set_ui(x, FLINT_UABS(v));
0373     if (v < 0)
0374         ARF_NEG(x);
0375 }
0376 
0377 ARF_INLINE void
0378 arf_set_ui(arf_t x, ulong v)
0379 {
0380     ARF_DEMOTE(x);
0381     _fmpz_demote(ARF_EXPREF(x));
0382 
0383     if (v == 0)
0384     {
0385         ARF_EXP(x) = ARF_EXP_ZERO;
0386         ARF_XSIZE(x) = 0;
0387     }
0388     else
0389     {
0390         unsigned int c;
0391         c = flint_clz(v);
0392         ARF_EXP(x) = FLINT_BITS - c;
0393         ARF_NOPTR_D(x)[0] = v << c;
0394         ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, 0);
0395     }
0396 }
0397 
0398 ARF_INLINE void
0399 arf_set_si(arf_t x, slong v)
0400 {
0401     arf_set_ui(x, FLINT_UABS(v));
0402     if (v < 0)
0403         ARF_NEG(x);
0404 }
0405 
0406 ARF_INLINE void
0407 arf_init_set_shallow(arf_t z, const arf_t x)
0408 {
0409     *z = *x;
0410 }
0411 
0412 ARF_INLINE void
0413 arf_init_neg_shallow(arf_t z, const arf_t x)
0414 {
0415     *z = *x;
0416     arf_neg(z, z);
0417 }
0418 
0419 ARF_INLINE void
0420 arf_init_set_mag_shallow(arf_t y, const mag_t x)
0421 {
0422     ulong t = MAG_MAN(x);
0423     ARF_XSIZE(y) = ARF_MAKE_XSIZE(t != 0, 0);
0424     ARF_EXP(y) = MAG_EXP(x);
0425     ARF_NOPTR_D(y)[0] = t << (FLINT_BITS - MAG_BITS);
0426 }
0427 
0428 ARF_INLINE void
0429 arf_init_neg_mag_shallow(arf_t z, const mag_t x)
0430 {
0431     arf_init_set_mag_shallow(z, x);
0432     arf_neg(z, z);
0433 }
0434 
0435 ARF_INLINE int
0436 arf_cmpabs_mag(const arf_t x, const mag_t y)
0437 {
0438     arf_t t;
0439     arf_init_set_mag_shallow(t, y);  /* no need to free */
0440     return arf_cmpabs(x, t);
0441 }
0442 
0443 ARF_INLINE int
0444 arf_mag_cmpabs(const mag_t x, const arf_t y)
0445 {
0446     arf_t t;
0447     arf_init_set_mag_shallow(t, x);  /* no need to free */
0448     return arf_cmpabs(t, y);
0449 }
0450 
0451 /* Assumes xn > 0, x[xn-1] != 0. */
0452 /* TBD: 1, 2 limb versions */
0453 void arf_set_mpn(arf_t y, nn_srcptr x, slong xn, int sgnbit);
0454 
0455 ARF_INLINE void
0456 arf_set_mpz(arf_t y, const mpz_t x)
0457 {
0458     slong size = x->_mp_size;
0459 
0460     if (size == 0)
0461         arf_zero(y);
0462     else
0463         arf_set_mpn(y, x->_mp_d, FLINT_ABS(size), size < 0);
0464 }
0465 
0466 ARF_INLINE void
0467 arf_set_fmpz(arf_t y, const fmpz_t x)
0468 {
0469     if (!COEFF_IS_MPZ(*x))
0470         arf_set_si(y, *x);
0471     else
0472         arf_set_mpz(y, COEFF_TO_PTR(*x));
0473 }
0474 
0475 int _arf_set_round_ui(arf_t x, ulong v, int sgnbit, slong prec, arf_rnd_t rnd);
0476 
0477 int _arf_set_round_uiui(arf_t z, slong * fix, ulong hi, ulong lo, int sgnbit, slong prec, arf_rnd_t rnd);
0478 
0479 int
0480 _arf_set_round_mpn(arf_t y, slong * exp_shift, nn_srcptr x, slong xn,
0481     int sgnbit, slong prec, arf_rnd_t rnd);
0482 
0483 ARF_INLINE int
0484 arf_set_round_ui(arf_t x, ulong v, slong prec, arf_rnd_t rnd)
0485 {
0486     return _arf_set_round_ui(x, v, 0, prec, rnd);
0487 }
0488 
0489 ARF_INLINE int
0490 arf_set_round_si(arf_t x, slong v, slong prec, arf_rnd_t rnd)
0491 {
0492     return _arf_set_round_ui(x, FLINT_UABS(v), v < 0, prec, rnd);
0493 }
0494 
0495 ARF_INLINE int
0496 arf_set_round_mpz(arf_t y, const mpz_t x, slong prec, arf_rnd_t rnd)
0497 {
0498     int inexact;
0499     slong size = x->_mp_size;
0500     slong fix;
0501 
0502     if (size == 0)
0503     {
0504         arf_zero(y);
0505         return 0;
0506     }
0507 
0508     inexact = _arf_set_round_mpn(y, &fix, x->_mp_d, FLINT_ABS(size),
0509         (size < 0), prec, rnd);
0510     _fmpz_demote(ARF_EXPREF(y));
0511     ARF_EXP(y) = FLINT_ABS(size) * FLINT_BITS + fix;
0512     return inexact;
0513 }
0514 
0515 ARF_INLINE int
0516 arf_set_round_fmpz(arf_t y, const fmpz_t x, slong prec, arf_rnd_t rnd)
0517 {
0518     if (!COEFF_IS_MPZ(*x))
0519         return arf_set_round_si(y, *x, prec, rnd);
0520     else
0521         return arf_set_round_mpz(y, COEFF_TO_PTR(*x), prec, rnd);
0522 }
0523 
0524 int arf_set_round(arf_t y, const arf_t x, slong prec, arf_rnd_t rnd);
0525 
0526 int arf_neg_round(arf_t y, const arf_t x, slong prec, arf_rnd_t rnd);
0527 
0528 #ifdef __MPFR_H
0529 int arf_get_mpfr(mpfr_t x, const arf_t y, mpfr_rnd_t rnd);
0530 void arf_set_mpfr(arf_t x, const mpfr_t y);
0531 int _arf_call_mpfr_func(arf_ptr r1, arf_ptr r2, int (*func)(void), arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0532 #endif
0533 
0534 int arf_equal(const arf_t x, const arf_t y);
0535 int arf_equal_si(const arf_t x, slong y);
0536 int arf_equal_ui(const arf_t x, ulong y);
0537 int arf_equal_d(const arf_t x, double y);
0538 
0539 ARF_INLINE void
0540 arf_min(arf_t z, const arf_t a, const arf_t b)
0541 {
0542     if (arf_cmp(a, b) <= 0)
0543         arf_set(z, a);
0544     else
0545         arf_set(z, b);
0546 }
0547 
0548 ARF_INLINE void
0549 arf_max(arf_t z, const arf_t a, const arf_t b)
0550 {
0551     if (arf_cmp(a, b) > 0)
0552         arf_set(z, a);
0553     else
0554         arf_set(z, b);
0555 }
0556 
0557 ARF_INLINE void
0558 arf_abs(arf_t y, const arf_t x)
0559 {
0560     if (arf_sgn(x) < 0)
0561         arf_neg(y, x);
0562     else
0563         arf_set(y, x);
0564 }
0565 
0566 ARF_INLINE slong
0567 arf_bits(const arf_t x)
0568 {
0569     if (arf_is_special(x))
0570         return 0;
0571     else
0572     {
0573         nn_srcptr xp;
0574         slong xn;
0575         slong c;
0576 
0577         ARF_GET_MPN_READONLY(xp, xn, x);
0578         c = flint_ctz(xp[0]);
0579         return xn * FLINT_BITS - c;
0580     }
0581 }
0582 
0583 ARF_INLINE void
0584 arf_bot(fmpz_t e, const arf_t x)
0585 {
0586     if (arf_is_special(x))
0587         fmpz_zero(e);
0588     else
0589         fmpz_sub_si(e, ARF_EXPREF(x), arf_bits(x));
0590 }
0591 
0592 int arf_is_int(const arf_t x);
0593 
0594 int arf_is_int_2exp_si(const arf_t x, slong e);
0595 
0596 int arf_cmp_2exp_si(const arf_t x, slong e);
0597 
0598 int arf_cmpabs_2exp_si(const arf_t x, slong e);
0599 
0600 ARF_INLINE void
0601 arf_set_si_2exp_si(arf_t x, slong man, slong exp)
0602 {
0603     arf_set_si(x, man);
0604     if (man != 0)
0605         fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
0606 }
0607 
0608 ARF_INLINE void
0609 arf_set_ui_2exp_si(arf_t x, ulong man, slong exp)
0610 {
0611     arf_set_ui(x, man);
0612     if (man != 0)
0613         fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
0614 }
0615 
0616 ARF_INLINE void
0617 arf_mul_2exp_si(arf_t y, const arf_t x, slong e)
0618 {
0619     arf_set(y, x);
0620     if (!arf_is_special(y))
0621         fmpz_add_si_inline(ARF_EXPREF(y), ARF_EXPREF(y), e);
0622 }
0623 
0624 ARF_INLINE void
0625 arf_mul_2exp_fmpz(arf_t y, const arf_t x, const fmpz_t e)
0626 {
0627     arf_set(y, x);
0628     if (!arf_is_special(y))
0629         fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), e);
0630 }
0631 
0632 ARF_INLINE int
0633 arf_set_round_fmpz_2exp(arf_t y, const fmpz_t x, const fmpz_t exp, slong prec, arf_rnd_t rnd)
0634 {
0635     if (fmpz_is_zero(x))
0636     {
0637         arf_zero(y);
0638         return 0;
0639     }
0640     else
0641     {
0642         int r = arf_set_round_fmpz(y, x, prec, rnd);
0643         fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), exp);
0644         return r;
0645     }
0646 }
0647 
0648 ARF_INLINE void
0649 arf_abs_bound_lt_2exp_fmpz(fmpz_t b, const arf_t x)
0650 {
0651     if (arf_is_special(x))
0652         fmpz_zero(b);
0653     else
0654         fmpz_set(b, ARF_EXPREF(x));
0655 }
0656 
0657 ARF_INLINE void
0658 arf_abs_bound_le_2exp_fmpz(fmpz_t b, const arf_t x)
0659 {
0660     if (arf_is_special(x))
0661         fmpz_zero(b);
0662     else if (ARF_IS_POW2(x))
0663         fmpz_sub_ui(b, ARF_EXPREF(x), 1);
0664     else
0665         fmpz_set(b, ARF_EXPREF(x));
0666 }
0667 
0668 slong arf_abs_bound_lt_2exp_si(const arf_t x);
0669 
0670 void arf_frexp(arf_t man, fmpz_t exp, const arf_t x);
0671 
0672 void arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x);
0673 
0674 int _arf_get_integer_mpn(nn_ptr y, nn_srcptr x, slong xn, slong exp);
0675 
0676 int _arf_set_mpn_fixed(arf_t z, nn_srcptr xp, slong xn,
0677         slong fixn, int negative, slong prec, arf_rnd_t rnd);
0678 
0679 int arf_get_fmpz(fmpz_t z, const arf_t x, arf_rnd_t rnd);
0680 
0681 slong arf_get_si(const arf_t x, arf_rnd_t rnd);
0682 
0683 int arf_get_fmpz_fixed_fmpz(fmpz_t y, const arf_t x, const fmpz_t e);
0684 
0685 int arf_get_fmpz_fixed_si(fmpz_t y, const arf_t x, slong e);
0686 
0687 ARF_INLINE void
0688 arf_set_fmpz_2exp(arf_t x, const fmpz_t man, const fmpz_t exp)
0689 {
0690     arf_set_fmpz(x, man);
0691     if (!arf_is_zero(x))
0692         fmpz_add_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
0693 }
0694 
0695 void arf_floor(arf_t z, const arf_t x);
0696 void arf_ceil(arf_t z, const arf_t x);
0697 void arf_nint(arf_t z, const arf_t x);
0698 
0699 void arf_debug(const arf_t x);
0700 
0701 char * arf_get_str(const arf_t x, slong d);
0702 
0703 #ifdef FLINT_HAVE_FILE
0704 void arf_fprint(FILE * file, const arf_t x);
0705 void arf_fprintd(FILE * file, const arf_t y, slong d);
0706 #endif
0707 
0708 void arf_print(const arf_t x);
0709 void arf_printd(const arf_t y, slong d);
0710 
0711 void arf_randtest(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
0712 
0713 void arf_randtest_not_zero(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
0714 
0715 void arf_randtest_special(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
0716 
0717 void arf_urandom(arf_t x, flint_rand_t state, slong bits, arf_rnd_t rnd);
0718 
0719 #define MUL_MPFR_MIN_LIMBS 25
0720 
0721 #if FLINT_HAVE_FFT_SMALL
0722 #define MUL_MPFR_MAX_LIMBS 800
0723 #else
0724 #define MUL_MPFR_MAX_LIMBS 10000
0725 #endif
0726 
0727 #define ARF_MUL_STACK_ALLOC 40
0728 #define ARF_MUL_TLS_ALLOC 1000
0729 
0730 extern FLINT_TLS_PREFIX nn_ptr __arf_mul_tmp;
0731 extern FLINT_TLS_PREFIX slong __arf_mul_alloc;
0732 
0733 extern void _arf_mul_tmp_cleanup(void);
0734 
0735 #define ARF_MUL_TMP_DECL \
0736     ulong tmp_stack[ARF_MUL_STACK_ALLOC]; \
0737 
0738 #define ARF_MUL_TMP_ALLOC(tmp, alloc) \
0739     if (alloc <= ARF_MUL_STACK_ALLOC) \
0740     { \
0741         tmp = tmp_stack; \
0742     } \
0743     else if (alloc <= ARF_MUL_TLS_ALLOC) \
0744     { \
0745         if (__arf_mul_alloc < alloc) \
0746         { \
0747             if (__arf_mul_alloc == 0) \
0748             { \
0749                 flint_register_cleanup_function(_arf_mul_tmp_cleanup); \
0750             } \
0751             __arf_mul_tmp = flint_realloc(__arf_mul_tmp, sizeof(ulong) * alloc); \
0752             __arf_mul_alloc = alloc; \
0753         } \
0754         tmp = __arf_mul_tmp; \
0755     } \
0756     else \
0757     { \
0758         tmp = flint_malloc(sizeof(ulong) * alloc); \
0759     }
0760 
0761 #define ARF_MUL_TMP_FREE(tmp, alloc) \
0762     if (alloc > ARF_MUL_TLS_ALLOC) \
0763         flint_free(tmp);
0764 
0765 void arf_mul_special(arf_t z, const arf_t x, const arf_t y);
0766 
0767 int arf_mul_via_mpfr(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd);
0768 
0769 int arf_mul_rnd_any(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0770 
0771 int arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec);
0772 
0773 #define arf_mul(z, x, y, prec, rnd)              \
0774     ((rnd == ARF_RND_DOWN)                      \
0775         ? arf_mul_rnd_down(z, x, y, prec)        \
0776         : arf_mul_rnd_any(z, x, y, prec, rnd))
0777 
0778 ARF_INLINE int
0779 arf_neg_mul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)
0780 {
0781     if (arf_is_special(y))
0782     {
0783         arf_mul(z, x, y, prec, rnd);
0784         arf_neg(z, z);
0785         return 0;
0786     }
0787     else
0788     {
0789         arf_t t;
0790         *t = *y;
0791         ARF_NEG(t);
0792         return arf_mul(z, x, t, prec, rnd);
0793     }
0794 }
0795 
0796 ARF_INLINE int
0797 arf_mul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
0798 {
0799     arf_t t;
0800     arf_init_set_ui(t, y); /* no need to free */
0801     return arf_mul(z, x, t, prec, rnd);
0802 }
0803 
0804 ARF_INLINE int
0805 arf_mul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
0806 {
0807     arf_t t;
0808     arf_init_set_si(t, y); /* no need to free */
0809     return arf_mul(z, x, t, prec, rnd);
0810 }
0811 
0812 int arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
0813 
0814 ARF_INLINE int
0815 arf_mul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
0816 {
0817     if (!COEFF_IS_MPZ(*y))
0818         return arf_mul_si(z, x, *y, prec, rnd);
0819     else
0820         return arf_mul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
0821 }
0822 
0823 #define ARF_ADD_STACK_ALLOC 40
0824 #define ARF_ADD_TLS_ALLOC 1000
0825 
0826 extern FLINT_TLS_PREFIX nn_ptr __arf_add_tmp;
0827 extern FLINT_TLS_PREFIX slong __arf_add_alloc;
0828 
0829 extern void _arf_add_tmp_cleanup(void);
0830 
0831 #define ARF_ADD_TMP_DECL \
0832     ulong tmp_stack[ARF_ADD_STACK_ALLOC]; \
0833 
0834 #define ARF_ADD_TMP_ALLOC(tmp, alloc) \
0835     if (alloc <= ARF_ADD_STACK_ALLOC) \
0836     { \
0837         tmp = tmp_stack; \
0838     } \
0839     else if (alloc <= ARF_ADD_TLS_ALLOC) \
0840     { \
0841         if (__arf_add_alloc < alloc) \
0842         { \
0843             if (__arf_add_alloc == 0) \
0844             { \
0845                 flint_register_cleanup_function(_arf_add_tmp_cleanup); \
0846             } \
0847             __arf_add_tmp = flint_realloc(__arf_add_tmp, sizeof(ulong) * alloc); \
0848             __arf_add_alloc = alloc; \
0849         } \
0850         tmp = __arf_add_tmp; \
0851     } \
0852     else \
0853     { \
0854         tmp = flint_malloc(sizeof(ulong) * alloc); \
0855     }
0856 
0857 #define ARF_ADD_TMP_FREE(tmp, alloc) \
0858     if (alloc > ARF_ADD_TLS_ALLOC) \
0859         flint_free(tmp);
0860 
0861 int _arf_add_mpn(arf_t z, nn_srcptr xp, slong xn, int xsgnbit,
0862     const fmpz_t xexp, nn_srcptr yp, slong yn, int ysgnbit,
0863     flint_bitcnt_t shift, slong prec, arf_rnd_t rnd);
0864 
0865 int arf_add(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0866 int arf_add_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd);
0867 int arf_add_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd);
0868 int arf_add_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd);
0869 
0870 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);
0871 
0872 int arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0873 int arf_sub_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd);
0874 int arf_sub_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd);
0875 int arf_sub_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd);
0876 
0877 int arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0878 
0879 ARF_INLINE int
0880 arf_addmul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
0881 {
0882     arf_t t;
0883     arf_init_set_ui(t, y); /* no need to free */
0884     return arf_addmul(z, x, t, prec, rnd);
0885 }
0886 
0887 ARF_INLINE int
0888 arf_addmul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
0889 {
0890     arf_t t;
0891     arf_init_set_si(t, y); /* no need to free */
0892     return arf_addmul(z, x, t, prec, rnd);
0893 }
0894 
0895 int arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
0896 
0897 ARF_INLINE int
0898 arf_addmul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
0899 {
0900     if (!COEFF_IS_MPZ(*y))
0901         return arf_addmul_si(z, x, *y, prec, rnd);
0902     else
0903         return arf_addmul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
0904 }
0905 
0906 int arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0907 
0908 ARF_INLINE int
0909 arf_submul_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
0910 {
0911     arf_t t;
0912     arf_init_set_ui(t, y); /* no need to free */
0913     return arf_submul(z, x, t, prec, rnd);
0914 }
0915 
0916 ARF_INLINE int
0917 arf_submul_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
0918 {
0919     arf_t t;
0920     arf_init_set_si(t, y); /* no need to free */
0921     return arf_submul(z, x, t, prec, rnd);
0922 }
0923 
0924 int arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, slong prec, arf_rnd_t rnd);
0925 
0926 ARF_INLINE int
0927 arf_submul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
0928 {
0929     if (!COEFF_IS_MPZ(*y))
0930         return arf_submul_si(z, x, *y, prec, rnd);
0931     else
0932         return arf_submul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
0933 }
0934 
0935 int arf_fma(arf_ptr res, arf_srcptr x, arf_srcptr y, arf_srcptr z, slong prec, arf_rnd_t rnd);
0936 
0937 int arf_sosq(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd);
0938 
0939 int arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd);
0940 
0941 ARF_INLINE int
0942 arf_div_ui(arf_ptr z, arf_srcptr x, ulong y, slong prec, arf_rnd_t rnd)
0943 {
0944     arf_t t;
0945     arf_init_set_ui(t, y); /* no need to free */
0946     return arf_div(z, x, t, prec, rnd);
0947 }
0948 
0949 ARF_INLINE int
0950 arf_ui_div(arf_ptr z, ulong x, arf_srcptr y, slong prec, arf_rnd_t rnd)
0951 {
0952     arf_t t;
0953     arf_init_set_ui(t, x); /* no need to free */
0954     return arf_div(z, t, y, prec, rnd);
0955 }
0956 
0957 ARF_INLINE int
0958 arf_div_si(arf_ptr z, arf_srcptr x, slong y, slong prec, arf_rnd_t rnd)
0959 {
0960     arf_t t;
0961     arf_init_set_si(t, y); /* no need to free */
0962     return arf_div(z, x, t, prec, rnd);
0963 }
0964 
0965 ARF_INLINE int
0966 arf_si_div(arf_ptr z, slong x, arf_srcptr y, slong prec, arf_rnd_t rnd)
0967 {
0968     arf_t t;
0969     arf_init_set_si(t, x); /* no need to free */
0970     return arf_div(z, t, y, prec, rnd);
0971 }
0972 
0973 ARF_INLINE int
0974 arf_div_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t rnd)
0975 {
0976     arf_t t;
0977     int r;
0978     arf_init(t);
0979     arf_set_fmpz(t, y);
0980     r = arf_div(z, x, t, prec, rnd);
0981     arf_clear(t);
0982     return r;
0983 }
0984 
0985 ARF_INLINE int
0986 arf_fmpz_div(arf_ptr z, const fmpz_t x, arf_srcptr y, slong prec, arf_rnd_t rnd)
0987 {
0988     arf_t t;
0989     int r;
0990     arf_init(t);
0991     arf_set_fmpz(t, x);
0992     r = arf_div(z, t, y, prec, rnd);
0993     arf_clear(t);
0994     return r;
0995 }
0996 
0997 ARF_INLINE int
0998 arf_fmpz_div_fmpz(arf_ptr z, const fmpz_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)
0999 {
1000     arf_t t, u;
1001     int r;
1002     arf_init(t);
1003     arf_init(u);
1004     arf_set_fmpz(t, x);
1005     arf_set_fmpz(u, y);
1006     r = arf_div(z, t, u, prec, rnd);
1007     arf_clear(t);
1008     arf_clear(u);
1009     return r;
1010 }
1011 
1012 int arf_sqrt(arf_ptr z, arf_srcptr x, slong prec, arf_rnd_t rnd);
1013 
1014 int arf_sqrt_ui(arf_t z, ulong x, slong prec, arf_rnd_t rnd);
1015 
1016 int arf_sqrt_fmpz(arf_t z, const fmpz_t x, slong prec, arf_rnd_t rnd);
1017 
1018 int arf_rsqrt(arf_ptr z, arf_srcptr x, slong prec, arf_rnd_t rnd);
1019 
1020 int arf_root(arf_ptr z, arf_srcptr x, ulong k, slong prec, arf_rnd_t rnd);
1021 
1022 /* Magnitude bounds */
1023 
1024 void arf_get_mag(mag_t y, const arf_t x);
1025 
1026 void arf_get_mag_lower(mag_t y, const arf_t x);
1027 
1028 ARF_INLINE void
1029 arf_set_mag(arf_t y, const mag_t x)
1030 {
1031     if (mag_is_zero(x))
1032     {
1033         arf_zero(y);
1034     }
1035     else if (mag_is_inf(x))
1036     {
1037         arf_pos_inf(y);
1038     }
1039     else
1040     {
1041         _fmpz_set_fast(ARF_EXPREF(y), MAG_EXPREF(x));
1042         ARF_DEMOTE(y);
1043         ARF_XSIZE(y) = ARF_MAKE_XSIZE(1, 0);
1044         ARF_NOPTR_D(y)[0] = MAG_MAN(x) << (FLINT_BITS - MAG_BITS);
1045     }
1046 }
1047 
1048 ARF_INLINE void
1049 mag_init_set_arf(mag_t y, const arf_t x)
1050 {
1051     mag_init(y);
1052     arf_get_mag(y, x);
1053 }
1054 
1055 ARF_INLINE void
1056 mag_fast_init_set_arf(mag_t y, const arf_t x)
1057 {
1058     if (ARF_IS_SPECIAL(x))   /* x == 0 */
1059     {
1060         mag_fast_zero(y);
1061     }
1062     else
1063     {
1064         nn_srcptr xp;
1065         slong xn;
1066 
1067         ARF_GET_MPN_READONLY(xp, xn, x);
1068 
1069         MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE;
1070         MAG_EXP(y) = ARF_EXP(x);
1071 
1072         MAG_FAST_ADJUST_ONE_TOO_LARGE(y);
1073     }
1074 }
1075 
1076 ARF_INLINE void
1077 arf_mag_fast_add_ulp(mag_t z, const mag_t x, const arf_t y, slong prec)
1078 {
1079     mag_fast_add_2exp_si(z, x, ARF_EXP(y) - prec);
1080 }
1081 
1082 ARF_INLINE void
1083 arf_mag_add_ulp(mag_t z, const mag_t x, const arf_t y, slong prec)
1084 {
1085     if (ARF_IS_SPECIAL(y))
1086     {
1087         flint_throw(FLINT_ERROR, "error: ulp error not defined for special value!\n");
1088     }
1089     else if (MAG_IS_LAGOM(z) && MAG_IS_LAGOM(x) && ARF_IS_LAGOM(y))
1090     {
1091         arf_mag_fast_add_ulp(z, x, y, prec);
1092     }
1093     else
1094     {
1095         fmpz_t e;
1096         fmpz_init(e);
1097         fmpz_sub_ui(e, ARF_EXPREF(y), prec);
1098         mag_add_2exp_fmpz(z, x, e);
1099         fmpz_clear(e);
1100     }
1101 }
1102 
1103 ARF_INLINE void
1104 arf_mag_set_ulp(mag_t z, const arf_t y, slong prec)
1105 {
1106     if (ARF_IS_SPECIAL(y))
1107     {
1108         flint_throw(FLINT_ERROR, "error: ulp error not defined for special value!\n");
1109     }
1110     else
1111     {
1112         _fmpz_add_fast(MAG_EXPREF(z), ARF_EXPREF(y), 1 - prec);
1113         MAG_MAN(z) = MAG_ONE_HALF;
1114     }
1115 }
1116 
1117 void arf_get_fmpq(fmpq_t y, const arf_t x);
1118 
1119 ARF_INLINE int
1120 arf_set_fmpq(arf_t y, const fmpq_t x, slong prec, arf_rnd_t rnd)
1121 {
1122     return arf_fmpz_div_fmpz(y, fmpq_numref(x), fmpq_denref(x), prec, rnd);
1123 }
1124 
1125 int arf_complex_mul(arf_t e, arf_t f, const arf_t a, const arf_t b,
1126                                       const arf_t c, const arf_t d,
1127                                       slong prec, arf_rnd_t rnd);
1128 
1129 int arf_complex_mul_fallback(arf_t e, arf_t f,
1130         const arf_t a, const arf_t b,
1131         const arf_t c, const arf_t d,
1132         slong prec, arf_rnd_t rnd);
1133 
1134 int arf_complex_sqr(arf_t e, arf_t f, const arf_t a, const arf_t b,
1135                                       slong prec, arf_rnd_t rnd);
1136 
1137 int arf_sum(arf_t s, arf_srcptr terms, slong len, slong prec, arf_rnd_t rnd);
1138 
1139 double arf_get_d(const arf_t x, arf_rnd_t rnd);
1140 void arf_set_d(arf_t x, double v);
1141 
1142 ARF_INLINE slong
1143 arf_allocated_bytes(const arf_t x)
1144 {
1145     slong size = fmpz_allocated_bytes(ARF_EXPREF(x));
1146 
1147     if (ARF_HAS_PTR(x))
1148         size += ARF_PTR_ALLOC(x) * sizeof(ulong);
1149 
1150     return size;
1151 }
1152 
1153 int arf_load_str(arf_t res, const char * data);
1154 char * arf_dump_str(const arf_t x);
1155 
1156 #ifdef FLINT_HAVE_FILE
1157 int arf_load_file(arf_t res, FILE *stream);
1158 int arf_dump_file(FILE* stream, const arf_t x);
1159 #endif
1160 
1161 int arf_dot(arf_t res, const arf_t initial, int subtract,
1162     arf_srcptr x, slong xstep, arf_srcptr y, slong ystep, slong len, slong prec, arf_rnd_t rnd);
1163 
1164 void arf_approx_dot(arf_t res, const arf_t initial, int subtract,
1165     arf_srcptr x, slong xstep, arf_srcptr y, slong ystep, slong len, slong prec, arf_rnd_t rnd);
1166 
1167 #ifdef __cplusplus
1168 }
1169 #endif
1170 
1171 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!