hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

more code and test code

Browse source
This commit is contained in:
Fredrik Johansson 2014-05-08 14:48:36 +02:00
parent 42ba244929
commit 0998aa9821
19 changed files with 1092 additions and 160 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

112
arb/get_mag_lbound.c Normal file
View file

@ -0,0 +1,112 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "arb.h"
/*
arb_get_mag_upper
arb_get_mag_lower
arf_get_mag_upper
arf_get_mag_lower
max(0, |mid|-rad)
*/
static __inline__ void
_arb_get_mag_lower(mag_t z, const arf_t mid, const mag_t rad)
{
if (mag_is_inf(rad) || arf_is_zero(mid))
{
mag_zero(z);
}
else if (mag_is_zero(rad) || !arf_is_finite(mid))
{
abort();
/* arf_get_mag_lower(z, mid); */
}
else
{
long shift, fix;
shift = _fmpz_sub_small(MAG_EXPREF(mid), MAG_EXPREF(rad));
/* mid < rad */
if (shift < 0)
{
mag_zero(z);
}
else if (shift <= 1) /* can be cancellation */
{
arf_t t;
arf_init(t);
arf_set_mag(t, rad);
if (arf_sgn(mid) > 0)
{
arf_sub(t, mid, t, MAG_BITS, ARF_RND_DOWN);
}
else
{
arf_add(t, mid, t, MAG_BITS, ARF_RND_DOWN);
arf_neg(t, t);
}
if (arf_sgn(t) <= 0)
mag_zero(z);
else
abort(); /* arf_get_mag_lower(z, t); */
arf_clear(t);
}
else
{
mp_limb_t m;
ARF_GET_TOP_LIMB(m, mid);
if (shift <= MAG_BITS)
m = m - (MAG_MAN(rad) >> shift) - 1;
else
m = m - 1;
fix = !(m >> (MAG_BITS - 1));
m <<= fix;
_fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(mid), -fix);
}
}
}
void
arb_get_mag_lower(mag_t z, const arb_t x)
{
_arb_get_mag_lower(z, ARB_MIDREF(x), ARB_RADREF(x));
}

140
arf.h
View file

@ -32,17 +32,12 @@
#include <math.h>
#include "flint.h"
#include "fmpr.h"
#include "mag.h"
#ifdef __cplusplus
extern "C" {
#endif
#define LIMB_ONE ((mp_limb_t) 1)
#define LIMB_ONES (-(mp_limb_t) 1)
#define LIMB_TOP (((mp_limb_t) 1) << (FLINT_BITS - 1))
#define MASK_LIMB(n, c) ((n) & (LIMB_ONES << (c)))
#define arf_rnd_t fmpr_rnd_t
#define ARF_RND_DOWN FMPR_RND_DOWN
#define ARF_RND_UP FMPR_RND_UP
@ -87,8 +82,8 @@ arf_rnd_to_mpfr(arf_rnd_t rnd)
#define ARF_RESULT_INEXACT 1
/* Range where we can skip fmpz overflow checks for exponent manipulation. */
#define ARF_MAX_LAGOM_EXP (COEFF_MAX / 4)
#define ARF_MIN_LAGOM_EXP (-ARF_MAX_LAGOM_EXP)
#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
@ -192,6 +187,15 @@ typedef const arf_struct * arf_srcptr;
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 { \
@ -742,74 +746,6 @@ void arf_randtest_not_zero(arf_t x, flint_rand_t state, long bits, long mag_bits
void arf_randtest_special(arf_t x, flint_rand_t state, long bits, long mag_bits);
#define ADD2_FAST_MAX (COEFF_MAX / 4)
#define ADD2_FAST_MIN (-ADD2_FAST_MAX)
static __inline__ void
_fmpz_set_fast(fmpz_t f, const fmpz_t g)
{
if (!COEFF_IS_MPZ(*f) && !COEFF_IS_MPZ(*g))
*f = *g;
else
fmpz_set(f, g);
}
static __inline__ void
_fmpz_add_fast(fmpz_t z, const fmpz_t x, long c)
{
fmpz ze, xe;
ze = *z;
xe = *x;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX))
*z = xe + c;
else
fmpz_add_si(z, x, c);
}
static __inline__ void
_fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
{
fmpz ze, xe, ye;
ze = *z;
xe = *x;
ye = *y;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) &&
(ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX))
{
*z = xe + ye + c;
}
else
{
fmpz_add(z, x, y);
fmpz_add_si(z, z, c);
}
}
static __inline__ void
_fmpz_sub2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
{
fmpz ze, xe, ye;
ze = *z;
xe = *x;
ye = *y;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) &&
(ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX))
{
*z = xe - ye + c;
}
else
{
fmpz_sub(z, x, y);
fmpz_add_si(z, z, c);
}
}
#define MUL_MPFR_MIN_LIMBS 25
#define MUL_MPFR_MAX_LIMBS 10000
@ -1143,6 +1079,58 @@ int arf_sqrt_fmpz(arf_t z, const fmpz_t x, long prec, arf_rnd_t rnd);
int arf_rsqrt(arf_ptr z, arf_srcptr x, long prec, arf_rnd_t rnd);
/* Magnitude bounds */
void arf_get_mag(mag_t y, const arf_t x);
static __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);
}
}
static __inline__ void
mag_init_set_arf(mag_t y, const arf_t x)
{
mag_init(y);
arf_get_mag(y, x);
}
static __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);
}
}
#ifdef __cplusplus
}
#endif

18
mag/set_arf.c → arf/get_mag.c
View file

@ -23,14 +23,11 @@
******************************************************************************/
#include "mag.h"
#include "arf.h"
void
mag_set_arf(mag_t y, const arf_t x)
arf_get_mag(mag_t y, const arf_t x)
{
mp_srcptr xp;
mp_size_t xn;
if (arf_is_zero(x))
{
mag_zero(y);
@ -41,11 +38,14 @@ mag_set_arf(mag_t y, const arf_t x)
}
else
{
ARF_GET_MPN_READONLY(xp, xn, x);
mp_limb_t t, u;
MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE;
_fmpz_set_fast(MAG_EXPREF(y), ARF_EXPREF(x));
MAG_ADJUST_ONE_TOO_LARGE(y); /* TODO: combine */
ARF_GET_TOP_LIMB(t, x);
t = (t >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE;
u = t >> MAG_BITS;
t = (t >> u) + u;
_fmpz_add_fast(MAG_EXPREF(y), ARF_EXPREF(x), u);
MAG_MAN(y) = t;
}
}

7
doc/source/arb.rst Normal file
View file

@ -0,0 +1,7 @@
.. _arb:
**arb.h** -- arbitrary-precision real floating-point balls
===============================================================================
TBD.

7
doc/source/arf.rst Normal file
View file

@ -0,0 +1,7 @@
.. _arf:
**arf.h** -- arbitrary-precision real floating-point numbers
===============================================================================
TBD.

10
doc/source/index.rst
View file

@ -47,6 +47,16 @@ Module documentation
hypgeom.rst
partitions.rst
Module documentation (Arb 2.x types)
::::::::::::::::::::::::::::::::::::
.. toctree::
:maxdepth: 1
mag.rst
arf.rst
arb.rst
Credits and references
::::::::::::::::::::::::

122
doc/source/mag.rst Normal file
View file

@ -0,0 +1,122 @@
.. _mag:
**mag.h** -- fixed-precision unsigned floating-point numbers for bounds
===============================================================================
The :type:`mag_t` type is an unsigned floating-point type with a
fixed-precision mantissa (30 bits) and unlimited exponent range, suited for
representing and rigorously manipulating magnitude bounds efficiently.
Operations always produce a strict upper or lower bound, but for performance
reasons, no attempt is made to compute the best possible bound
(in general, a result may a few ulps larger/smaller than the optimal value).
The special values zero and positive infinity are supported (but not NaN).
Applications requiring more flexibility (such as correct rounding, or
higher precision) should use the :type:`arf_t` type instead.
The exponent is represented as an :type:`fmpz`. Special "fast" methods are
provided for manipulating :type:`mag_t` objects under the assumption that
no infinities are present and that exponents stay within a certain range.
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: mag_struct
A :type:`mag_struct` holds a mantissa and an exponent.
Special values are encoded by the mantissa being set to zero.
.. type:: mag_t
A :type:`mag_t` is defined as an array of length one of type
:type:`mag_struct`, permitting a :type:`mag_t` to be passed by reference.
Memory management
-------------------------------------------------------------------------------
.. function:: void mag_init(mag_t x)
Initializes the variable *x* for use. Its value is set to zero.
.. function:: void mag_clear(mag_t x)
Clears the variable *x*, freeing or recycling its allocated memory.
.. function:: void mag_init_set(mag_t x, const mag_t y)
.. function:: void mag_swap(mag_t x, mag_t y)
.. function:: void mag_set(mag_t x, const mag_t y)
Special values
-------------------------------------------------------------------------------
.. function:: void mag_zero(mag_t x)
Sets *x* to zero.
.. function:: void mag_inf(mag_t x)
Sets *x* to positive infinity.
.. function:: int mag_is_special(const mag_t x)
Returns nonzero iff *x* is zero or positive infinity.
.. function:: int mag_is_zero(const mag_t x)
Returns nonzero iff *x* is zero.
.. function:: int mag_is_inf(const mag_t x)
Returns nonzero iff *x* is positive infinity.
Arithmetic
-------------------------------------------------------------------------------
.. function:: void mag_mul(mag_t z, const mag_t x, const mag_t y)
Sets `z` to an upper bound for `x + y`.
.. function:: void mag_addmul(mag_t z, const mag_t x, const mag_t y)
Sets `z` to an upper bound for `z + xy`.
.. function:: void mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e)
Sets `z` to an upper bound for `x + 2^e`.
.. function:: void mag_div(mag_t z, const mag_t x, const mag_t y)
Sets `z` to an upper bound for `x / y`.
Fast versions
-------------------------------------------------------------------------------
.. function:: void mag_fast_init_set(mag_t x, const mag_t y)
.. function:: void mag_fast_zero(mag_t x)
.. function:: int mag_fast_is_zero(const mag_t x)
.. function:: void mag_fast_init_set_arf(mag_t y, const arf_t x)
.. function:: void mag_fast_mul(mag_t z, const mag_t x, const mag_t y)
.. function:: void mag_fast_addmul(mag_t z, const mag_t x, const mag_t y)
.. function:: void mag_fast_add_2exp_si(mag_t z, const mag_t x, long e)
Conversions
-------------------------------------------------------------------------------
These functions are intended for debugging purposes: see the :doc:`arf <arf>`
module for other conversion functions.
.. function:: void mag_set_fmpr(mag_t y, const fmpr_t x)
Sets *y* to an upper bound for *x*.
.. function:: void mag_get_fmpr(fmpr_t y, const mag_t x)
Sets *y* to exactly *x*.

172
mag.h
View file

@ -28,37 +28,102 @@
#include <math.h>
#include "flint.h"
#include "arf.h"
#include "fmpz.h"
#include "fmpz_extras.h"
#ifdef __cplusplus
extern "C" {
#endif
/*
The mag_t type is an unsigned floating-point type with a fixed-precision
mantissa (30 bits) and unlimited exponent range, suited for representing
magnitude bounds efficiently.
#define LIMB_ONE ((mp_limb_t) 1)
#define LIMB_ONES (-(mp_limb_t) 1)
#define LIMB_TOP (((mp_limb_t) 1) << (FLINT_BITS - 1))
#define MASK_LIMB(n, c) ((n) & (LIMB_ONES << (c)))
Operations always produce a strict upper/lower bound, but for performance
reasons, no attempt is made to compute the best possible bound
(in general, a result may a few ulps larger/smaller than the optimal value).
#define MAG_MAX_LAGOM_EXP (COEFF_MAX / 4)
#define MAG_MIN_LAGOM_EXP (-MAG_MAX_LAGOM_EXP)
The special values zero and positive infinity are supported (but not NaN).
#define ADD2_FAST_MAX (COEFF_MAX / 4)
#define ADD2_FAST_MIN (-ADD2_FAST_MAX)
Applications requiring more flexibility (such as correct rounding, or
higher precision) should use the arf_t type instead.
*/
/* TODO: rename these and move them to fmpz_extras */
/* TODO: arf.h should depend on this, not the other way around */
static __inline__ void
_fmpz_set_fast(fmpz_t f, const fmpz_t g)
{
if (!COEFF_IS_MPZ(*f) && !COEFF_IS_MPZ(*g))
*f = *g;
else
fmpz_set(f, g);
}
static __inline__ void
_fmpz_add_fast(fmpz_t z, const fmpz_t x, long c)
{
fmpz ze, xe;
ze = *z;
xe = *x;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX))
*z = xe + c;
else
fmpz_add_si(z, x, c);
}
static __inline__ void
_fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
{
fmpz ze, xe, ye;
ze = *z;
xe = *x;
ye = *y;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) &&
(ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX))
{
*z = xe + ye + c;
}
else
{
fmpz_add(z, x, y);
fmpz_add_si(z, z, c);
}
}
static __inline__ void
_fmpz_sub2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
{
fmpz ze, xe, ye;
ze = *z;
xe = *x;
ye = *y;
if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) &&
(ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX))
{
*z = xe - ye + c;
}
else
{
fmpz_sub(z, x, y);
fmpz_add_si(z, z, c);
}
}
#define MAG_EXP_POS_INF (COEFF_MIN+1)
/* Finite and with lagom big exponents. */
#define MAG_IS_LAGOM(x) (MAG_EXP(x) >= MAG_MIN_LAGOM_EXP && \
MAG_EXP(x) <= MAG_MAX_LAGOM_EXP)
#define MAG_EXPREF(x) (&(x)->exp)
#define MAG_EXP(x) ((x)->exp)
#define MAG_MAN(x) ((x)->man)
/* Finite and with lagom big exponents. */
#define MAG_IS_LAGOM(x) (MAG_EXP(x) >= ARF_MIN_LAGOM_EXP && \
MAG_EXP(x) <= ARF_MAX_LAGOM_EXP)
#define MAG_BITS 30
#define MAG_ONE_HALF (1UL << (MAG_BITS - 1))
@ -155,7 +220,7 @@ mag_swap(mag_t x, mag_t y)
static __inline__ void
mag_set(mag_t x, const mag_t y)
{
_fmpz_set_fast(ARF_EXPREF(x), ARF_EXPREF(y));
_fmpz_set_fast(MAG_EXPREF(x), MAG_EXPREF(y));
x->man = y->man;
}
@ -166,13 +231,6 @@ mag_zero(mag_t x)
MAG_MAN(x) = 0;
}
static __inline__ void
mag_inf(mag_t x)
{
fmpz_clear(MAG_EXPREF(x));
MAG_EXP(x) = ARF_EXP_POS_INF;
}
static __inline__ int
mag_is_special(const mag_t x)
{
@ -185,6 +243,13 @@ mag_is_zero(const mag_t x)
return (MAG_MAN(x) == 0) && (MAG_EXP(x) == 0);
}
static __inline__ void
mag_inf(mag_t x)
{
fmpz_clear(MAG_EXPREF(x));
MAG_EXP(x) = MAG_EXP_POS_INF;
}
static __inline__ int
mag_is_inf(const mag_t x)
{
@ -193,23 +258,14 @@ mag_is_inf(const mag_t x)
/* general versions */
void mag_set_fmpr(mag_t x, const fmpr_t y);
void mag_set_arf(mag_t y, const arf_t x);
static __inline__ void
mag_init_set_arf(mag_t y, const arf_t x)
{
mag_init(y);
mag_set_arf(y, x);
}
void mag_mul(mag_t z, const mag_t x, const mag_t y);
void mag_addmul(mag_t z, const mag_t x, const mag_t y);
void mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e);
void mag_div(mag_t z, const mag_t x, const mag_t y);
/* Fast versions (no infs/nans, small exponents). Note that this
applies to outputs too! */
@ -233,27 +289,6 @@ mag_fast_is_zero(const mag_t x)
return MAG_MAN(x) == 0;
}
static __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);
}
}
static __inline__ void
mag_fast_mul(mag_t z, const mag_t x, const mag_t y)
{
@ -293,8 +328,7 @@ mag_fast_addmul(mag_t z, const mag_t x, const mag_t y)
if (shift >= MAG_BITS)
MAG_MAN(z)++;
else
MAG_MAN(z) = MAG_MAN(z) + (MAG_FIXMUL(MAG_MAN(x),
MAG_MAN(y)) >> shift) + LIMB_ONE;
MAG_MAN(z) = MAG_MAN(z) + (MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) >> shift) + 1;
}
else
{
@ -302,11 +336,11 @@ mag_fast_addmul(mag_t z, const mag_t x, const mag_t y)
MAG_EXP(z) = e;
if (shift >= MAG_BITS)
MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y))
+ (2 * LIMB_ONE);
MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) + 2;
else
MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y))
+ (MAG_MAN(z) >> shift) + (2 * LIMB_ONE);
MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) + (MAG_MAN(z) >> shift) + 2;
MAG_FAST_ADJUST_ONE_TOO_SMALL(z);
}
MAG_FAST_ADJUST_ONE_TOO_LARGE(z);
@ -327,7 +361,7 @@ mag_fast_add_2exp_si(mag_t z, const mag_t x, long e)
long shift;
shift = MAG_EXP(x) - e;
if (shift >= 0)
if (shift > 0)
{
MAG_EXP(z) = MAG_EXP(x);
@ -340,18 +374,21 @@ mag_fast_add_2exp_si(mag_t z, const mag_t x, long e)
{
shift = -shift;
MAG_EXP(z) = e;
MAG_EXP(z) = e + 1;
if (shift >= MAG_BITS)
MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + LIMB_ONE;
MAG_MAN(z) = MAG_ONE_HALF + LIMB_ONE;
else
MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + (MAG_MAN(x) >> shift);
MAG_MAN(z) = MAG_ONE_HALF + (MAG_MAN(x) >> (shift + 1)) + LIMB_ONE;
}
MAG_FAST_ADJUST_ONE_TOO_LARGE(z);
}
}
#include "fmpr.h"
void mag_set_fmpr(mag_t x, const fmpr_t y);
static __inline__ void
mag_get_fmpr(fmpr_t x, const mag_t r)
@ -371,7 +408,6 @@ mag_get_fmpr(fmpr_t x, const mag_t r)
}
}
#ifdef __cplusplus
}
#endif

8
mag/add_2exp_fmpz.c
View file

@ -46,7 +46,7 @@ mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e)
shift = _fmpz_sub_small(MAG_EXPREF(x), e);
if (shift >= 0)
if (shift > 0)
{
_fmpz_set_fast(MAG_EXPREF(z), MAG_EXPREF(x));
@ -59,12 +59,12 @@ mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e)
{
shift = -shift;
_fmpz_set_fast(MAG_EXPREF(z), e);
_fmpz_add_fast(MAG_EXPREF(z), e, 1);
if (shift >= MAG_BITS)
MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + LIMB_ONE;
MAG_MAN(z) = MAG_ONE_HALF + LIMB_ONE;
else
MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + (MAG_MAN(x) >> shift);
MAG_MAN(z) = MAG_ONE_HALF + (MAG_MAN(x) >> (shift + 1)) + LIMB_ONE;
}
MAG_ADJUST_ONE_TOO_LARGE(z);

4
mag/addmul.c
View file

@ -38,7 +38,7 @@ mag_addmul(mag_t z, const mag_t x, const mag_t y)
}
else if (mag_is_zero(x) || mag_is_zero(y))
{
mag_zero(z);
return;
}
else
{
@ -70,6 +70,8 @@ mag_addmul(mag_t z, const mag_t x, const mag_t y)
else
MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y))
+ (MAG_MAN(z) >> shift) + (2 * LIMB_ONE);
MAG_ADJUST_ONE_TOO_SMALL(z);
}
MAG_ADJUST_ONE_TOO_LARGE(z);

58
mag/div.c Normal file
View file

@ -0,0 +1,58 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
void
mag_div(mag_t z, const mag_t x, const mag_t y)
{
if (mag_is_special(x) || mag_is_special(y))
{
if (mag_is_zero(y) || mag_is_inf(x))
mag_inf(z);
else
mag_zero(z);
}
else
{
mp_limb_t q;
long fix;
#if FLINT_BITS == 64
q = (MAG_MAN(x) << MAG_BITS) / MAG_MAN(y) + 1;
#else
mp_limb_t hi, lo, r;
lo = MAG_MAN(x) << MAG_BITS;
hi = MAG_MAN(x) >> (FLINT_BITS - MAG_BITS);
udiv_qrnnd(q, r, hi, lo, MAG_MAN(y));
q += 1;
#endif
fix = q >> MAG_BITS;
q = (q >> fix) + fix;
_fmpz_sub2_fast(MAG_EXPREF(z), MAG_EXPREF(x), MAG_EXPREF(y), fix);
}
}

2
mag/mul.c
View file

@ -30,7 +30,7 @@ mag_mul(mag_t z, const mag_t x, const mag_t y)
{
if (mag_is_special(x) || mag_is_special(y))
{
if (mag_is_inf(x) || mag_is_special(y))
if (mag_is_inf(x) || mag_is_inf(y))
mag_inf(z);
else
mag_zero(z);

3
mag/set_fmpr.c
View file

@ -24,6 +24,7 @@
******************************************************************************/
#include "mag.h"
#include "arf.h"
void
mag_set_fmpr(mag_t x, const fmpr_t y)
@ -31,7 +32,7 @@ mag_set_fmpr(mag_t x, const fmpr_t y)
arf_t t;
arf_init(t);
arf_set_fmpr(t, y);
mag_set_arf(x, t);
arf_get_mag(x, t);
arf_clear(t);
}

100
mag/test/t-add_2exp_fmpz.c Normal file
View file

@ -0,0 +1,100 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("add_2exp_fmpz....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, zb;
fmpz_t e;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
fmpz_init(e);
mag_init(xb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 100);
fmpr_abs(x, x);
fmpz_randtest(e, state, 100);
fmpr_one(y);
fmpr_mul_2exp_fmpz(y, y, e);
fmpr_add(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_set_fmpr(xb, x);
mag_add_2exp_fmpz(zb, xb, e);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_printd(x, 15); printf("\n\n");
printf("y = "); fmpr_printd(y, 15); printf("\n\n");
printf("z = "); fmpr_printd(z, 15); printf("\n\n");
printf("w = "); fmpr_printd(w, 15); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(zb);
fmpz_clear(e);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

102
mag/test/t-addmul.c Normal file
View file

@ -0,0 +1,102 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("addmul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, yb, zb;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
mag_init(xb);
mag_init(yb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 100);
fmpr_randtest(y, state, 2 + n_randint(state, 200), 100);
fmpr_randtest(z, state, 2 + n_randint(state, 200), 100);
fmpr_abs(x, x);
fmpr_abs(y, y);
fmpr_abs(z, z);
mag_set_fmpr(xb, x);
mag_set_fmpr(yb, y);
mag_set_fmpr(zb, z);
fmpr_addmul(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_addmul(zb, xb, yb);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(yb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_printd(x, 15); printf("\n\n");
printf("y = "); fmpr_printd(y, 15); printf("\n\n");
printf("z = "); fmpr_printd(z, 15); printf("\n\n");
printf("w = "); fmpr_printd(w, 15); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(yb);
mag_clear(zb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

95
mag/test/t-fast_add_2exp_si.c Normal file
View file

@ -0,0 +1,95 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("fast_add_2exp_si....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, zb;
long e;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
mag_init(xb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 15);
fmpr_abs(x, x);
e = n_randint(state, 10000) - n_randint(state, 10000);
fmpr_set_ui_2exp_si(y, 1, e);
fmpr_add(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_set_fmpr(xb, x);
mag_fast_add_2exp_si(zb, xb, e);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_printd(x, 15); printf("\n\n");
printf("y = "); fmpr_printd(y, 15); printf("\n\n");
printf("z = "); fmpr_printd(z, 15); printf("\n\n");
printf("w = "); fmpr_printd(w, 15); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(zb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

102
mag/test/t-fast_addmul.c Normal file
View file

@ -0,0 +1,102 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("fast_addmul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, yb, zb;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
mag_init(xb);
mag_init(yb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 15);
fmpr_randtest(y, state, 2 + n_randint(state, 200), 15);
fmpr_randtest(z, state, 2 + n_randint(state, 200), 15);
fmpr_abs(x, x);
fmpr_abs(y, y);
fmpr_abs(z, z);
mag_set_fmpr(xb, x);
mag_set_fmpr(yb, y);
mag_set_fmpr(zb, z);
fmpr_addmul(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_fast_addmul(zb, xb, yb);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(yb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_printd(x, 15); printf("\n\n");
printf("y = "); fmpr_printd(y, 15); printf("\n\n");
printf("z = "); fmpr_printd(z, 15); printf("\n\n");
printf("w = "); fmpr_printd(w, 15); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(yb);
mag_clear(zb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

95
mag/test/t-fast_mul.c Normal file
View file

@ -0,0 +1,95 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("fast_mul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, yb, zb;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
mag_init(xb);
mag_init(yb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 15);
fmpr_randtest(y, state, 2 + n_randint(state, 200), 15);
fmpr_mul(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_set_fmpr(xb, x);
mag_set_fmpr(yb, y);
mag_fast_mul(zb, xb, yb);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(yb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_printd(x, 15); printf("\n\n");
printf("y = "); fmpr_printd(y, 15); printf("\n\n");
printf("z = "); fmpr_printd(z, 15); printf("\n\n");
printf("w = "); fmpr_printd(w, 15); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(yb);
mag_clear(zb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

95
mag/test/t-mul.c Normal file
View file

@ -0,0 +1,95 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("mul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2, w;
mag_t xb, yb, zb;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
fmpr_init(w);
mag_init(xb);
mag_init(yb);
mag_init(zb);
fmpr_randtest(x, state, 2 + n_randint(state, 200), 100);
fmpr_randtest(y, state, 2 + n_randint(state, 200), 100);
fmpr_mul(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
mag_set_fmpr(xb, x);
mag_set_fmpr(yb, y);
mag_mul(zb, xb, yb);
mag_get_fmpr(w, zb);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(yb)
MAG_CHECK_BITS(zb)
if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_print(x); printf("\n\n");
printf("y = "); fmpr_print(y); printf("\n\n");
printf("z = "); fmpr_print(z); printf("\n\n");
printf("w = "); fmpr_print(w); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
fmpr_clear(w);
mag_clear(xb);
mag_clear(yb);
mag_clear(zb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}