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 multiplication and addmul/submul methods

Browse source
This commit is contained in:
Fredrik Johansson 2014-05-05 10:52:51 +02:00
parent a47e1580c6
commit 4d66d4968a
15 changed files with 2004 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

85
arf.h
View file

@ -915,6 +915,33 @@ int arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec);
? arf_mul_rnd_down(z, x, y, prec) \
: arf_mul_rnd_any(z, x, y, prec, rnd))
static __inline__ int
arf_mul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_mul(z, x, t, prec, rnd);
}
static __inline__ int
arf_mul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_mul(z, x, t, prec, rnd);
}
int arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd);
static __inline__ int
arf_mul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_mul_si(z, x, *y, prec, rnd);
else
return arf_mul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
#define ARF_ADD_STACK_ALLOC 40
#define ARF_ADD_TLS_ALLOC 1000
@ -967,6 +994,64 @@ int arf_sub_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd);
int arf_sub_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd);
int arf_sub_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd);
int arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
static __inline__ int
arf_addmul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_addmul(z, x, t, prec, rnd);
}
static __inline__ int
arf_addmul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_addmul(z, x, t, prec, rnd);
}
int arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd);
static __inline__ int
arf_addmul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_addmul_si(z, x, *y, prec, rnd);
else
return arf_addmul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
int arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
static __inline__ int
arf_submul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_ui(t, y); /* no need to free */
return arf_submul(z, x, t, prec, rnd);
}
static __inline__ int
arf_submul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
arf_init_set_si(t, y); /* no need to free */
return arf_submul(z, x, t, prec, rnd);
}
int arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd);
static __inline__ int
arf_submul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
if (!COEFF_IS_MPZ(*y))
return arf_submul_si(z, x, *y, prec, rnd);
else
return arf_submul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd);
}
int arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
static __inline__ int

164
arf/addmul.c Normal file
View file

@ -0,0 +1,164 @@
/*=============================================================================
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 "arf.h"
int
arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn, zn, tn, alloc;
mp_srcptr xptr, yptr, zptr;
mp_ptr tptr, tptr2;
fmpz_t texp;
long shift;
int tsgnbit, inexact;
ARF_MUL_TMP_DECL
if (arf_is_special(x) || arf_is_special(y) || arf_is_special(z))
{
if (arf_is_zero(z))
{
return arf_mul(z, x, y, prec, rnd);
}
else if (arf_is_finite(x) && arf_is_finite(y))
{
return arf_set_round(z, z, prec, rnd);
}
else
{
/* todo: speed up */
arf_t t;
arf_init(t);
arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
}
tsgnbit = ARF_SGNBIT(x) ^ ARF_SGNBIT(y);
ARF_GET_MPN_READONLY(xptr, xn, x);
ARF_GET_MPN_READONLY(yptr, yn, y);
ARF_GET_MPN_READONLY(zptr, zn, z);
fmpz_init(texp);
_fmpz_add2_fast(texp, ARF_EXPREF(x), ARF_EXPREF(y), 0);
shift = _fmpz_sub_small(ARF_EXPREF(z), texp);
alloc = tn = xn + yn;
ARF_MUL_TMP_ALLOC(tptr2, alloc)
tptr = tptr2;
ARF_MPN_MUL(tptr, xptr, xn, yptr, yn);
tn -= (tptr[0] == 0);
tptr += (tptr[0] == 0);
if (shift >= 0)
inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z),
tptr, tn, tsgnbit, shift, prec, rnd);
else
inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp,
zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd);
ARF_MUL_TMP_FREE(tptr2, alloc)
fmpz_clear(texp);
return inexact;
}
int
arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn, zn, tn, alloc;
mp_srcptr xptr, yptr, zptr;
mp_ptr tptr, tptr2;
fmpz_t texp, yexp;
long shift;
int tsgnbit, ysgnbit, inexact;
ARF_MUL_TMP_DECL
yn = FLINT_ABS(y->_mp_size);
if (arf_is_special(x) || yn == 0 || arf_is_special(z))
{
if (arf_is_zero(z))
{
return arf_mul_mpz(z, x, y, prec, rnd);
}
else if (arf_is_finite(x))
{
return arf_set_round(z, z, prec, rnd);
}
else
{
/* todo: speed up */
arf_t t;
arf_init(t);
arf_mul_mpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
}
ARF_GET_MPN_READONLY(xptr, xn, x);
yptr = y->_mp_d;
ysgnbit = (y->_mp_size < 0);
*yexp = yn * FLINT_BITS;
ARF_GET_MPN_READONLY(zptr, zn, z);
fmpz_init(texp);
tsgnbit = ARF_SGNBIT(x) ^ ysgnbit;
alloc = tn = xn + yn;
ARF_MUL_TMP_ALLOC(tptr2, alloc)
tptr = tptr2;
ARF_MPN_MUL(tptr, xptr, xn, yptr, yn);
shift = (tptr[tn - 1] == 0) * FLINT_BITS;
tn -= (tptr[tn - 1] == 0);
_fmpz_add2_fast(texp, ARF_EXPREF(x), yexp, -shift);
shift = _fmpz_sub_small(ARF_EXPREF(z), texp);
if (shift >= 0)
inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z),
tptr, tn, tsgnbit, shift, prec, rnd);
else
inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp,
zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd);
ARF_MUL_TMP_FREE(tptr2, alloc)
fmpz_clear(texp);
return inexact;
}

55
arf/mul_rnd_down.c
View file

@ -232,3 +232,58 @@ arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec)
}
}
int
arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn;
long fix, shift;
int sgnbit, inexact;
yn = FLINT_ABS(y->_mp_size);
xn = ARF_XSIZE(x);
xn >>= 1;
sgnbit = ARF_SGNBIT(x) ^ (y->_mp_size < 0);
/* Either operand is a special value. */
if (xn == 0 || yn == 0)
{
if (arf_is_finite(x))
{
arf_zero(z);
}
else
{
arf_t t;
arf_init_set_si(t, mpz_sgn(y));
arf_mul(z, x, t, prec, rnd);
arf_clear(t);
}
return 0;
}
else
{
mp_size_t zn, alloc;
mp_srcptr xptr, yptr;
mp_ptr tmp;
ARF_MUL_TMP_DECL
ARF_GET_MPN_READONLY(xptr, xn, x);
yptr = y->_mp_d;
alloc = zn = xn + yn;
ARF_MUL_TMP_ALLOC(tmp, alloc)
ARF_MPN_MUL(tmp, xptr, xn, yptr, yn);
shift = yn * FLINT_BITS - (tmp[zn - 1] == 0) * FLINT_BITS;
zn -= (tmp[zn - 1] == 0);
inexact = _arf_set_round_mpn(z, &fix, tmp, zn, sgnbit, prec, rnd);
_fmpz_add_fast(ARF_EXPREF(z), ARF_EXPREF(x), fix + shift);
ARF_MUL_TMP_FREE(tmp, alloc)
return inexact;
}
}

168
arf/submul.c Normal file
View file

@ -0,0 +1,168 @@
/*=============================================================================
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 "arf.h"
int
arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn, zn, tn, alloc;
mp_srcptr xptr, yptr, zptr;
mp_ptr tptr, tptr2;
fmpz_t texp;
long shift;
int tsgnbit, inexact;
ARF_MUL_TMP_DECL
if (arf_is_special(x) || arf_is_special(y) || arf_is_special(z))
{
if (arf_is_zero(z))
{
/* TODO: make more efficient */
arf_mul(z, x, y, ARF_PREC_EXACT, rnd);
return arf_neg_round(z, z, prec, rnd);
}
else if (arf_is_finite(x) && arf_is_finite(y))
{
return arf_set_round(z, z, prec, rnd);
}
else
{
/* todo: speed up */
arf_t t;
arf_init(t);
arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
}
tsgnbit = ARF_SGNBIT(x) ^ ARF_SGNBIT(y) ^ 1;
ARF_GET_MPN_READONLY(xptr, xn, x);
ARF_GET_MPN_READONLY(yptr, yn, y);
ARF_GET_MPN_READONLY(zptr, zn, z);
fmpz_init(texp);
_fmpz_add2_fast(texp, ARF_EXPREF(x), ARF_EXPREF(y), 0);
shift = _fmpz_sub_small(ARF_EXPREF(z), texp);
alloc = tn = xn + yn;
ARF_MUL_TMP_ALLOC(tptr2, alloc)
tptr = tptr2;
ARF_MPN_MUL(tptr, xptr, xn, yptr, yn);
tn -= (tptr[0] == 0);
tptr += (tptr[0] == 0);
if (shift >= 0)
inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z),
tptr, tn, tsgnbit, shift, prec, rnd);
else
inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp,
zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd);
ARF_MUL_TMP_FREE(tptr2, alloc)
fmpz_clear(texp);
return inexact;
}
int
arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn, zn, tn, alloc;
mp_srcptr xptr, yptr, zptr;
mp_ptr tptr, tptr2;
fmpz_t texp, yexp;
long shift;
int tsgnbit, ysgnbit, inexact;
ARF_MUL_TMP_DECL
yn = FLINT_ABS(y->_mp_size);
if (arf_is_special(x) || yn == 0 || arf_is_special(z))
{
if (arf_is_zero(z))
{
/* TODO: make more efficient */
arf_mul_mpz(z, x, y, ARF_PREC_EXACT, rnd);
return arf_neg_round(z, z, prec, rnd);
}
else if (arf_is_finite(x))
{
return arf_set_round(z, z, prec, rnd);
}
else
{
/* todo: speed up */
arf_t t;
arf_init(t);
arf_mul_mpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
}
ARF_GET_MPN_READONLY(xptr, xn, x);
yptr = y->_mp_d;
ysgnbit = (y->_mp_size > 0);
*yexp = yn * FLINT_BITS;
ARF_GET_MPN_READONLY(zptr, zn, z);
fmpz_init(texp);
tsgnbit = ARF_SGNBIT(x) ^ ysgnbit;
alloc = tn = xn + yn;
ARF_MUL_TMP_ALLOC(tptr2, alloc)
tptr = tptr2;
ARF_MPN_MUL(tptr, xptr, xn, yptr, yn);
shift = (tptr[tn - 1] == 0) * FLINT_BITS;
tn -= (tptr[tn - 1] == 0);
_fmpz_add2_fast(texp, ARF_EXPREF(x), yexp, -shift);
shift = _fmpz_sub_small(ARF_EXPREF(z), texp);
if (shift >= 0)
inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z),
tptr, tn, tsgnbit, shift, prec, rnd);
else
inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp,
zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd);
ARF_MUL_TMP_FREE(tptr2, alloc)
fmpz_clear(texp);
return inexact;
}

179
arf/test/t-addmul.c Normal file
View file

@ -0,0 +1,179 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_addmul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("addmul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, y, z, v;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(y);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
arf_randtest_special(y, state, 2000, 100);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(y)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 5))
{
case 0:
r1 = arf_addmul(z, x, y, prec, rnd);
r2 = arf_addmul_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); arf_print(y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 1:
r1 = arf_addmul(z, x, x, prec, rnd);
r2 = arf_addmul_naive(v, x, x, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL (aliasing 1)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 2:
r2 = arf_addmul_naive(v, v, v, prec, rnd);
r1 = arf_addmul(z, z, z, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing 2)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 3:
r2 = arf_addmul_naive(v, v, y, prec, rnd);
r1 = arf_addmul(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing 3)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = "); arf_print(y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_addmul_naive(v, x, v, prec, rnd);
r1 = arf_addmul(z, x, z, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL (aliasing 4)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(y);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

135
arf/test/t-addmul_fmpz.c Normal file
View file

@ -0,0 +1,135 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_addmul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_fmpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("addmul_fmpz....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
fmpz_t y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
fmpz_init(y);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
fmpz_randtest(y, state, 2000);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_addmul_fmpz(z, x, y, prec, rnd);
r2 = arf_addmul_fmpz_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); fmpz_print(y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_addmul_fmpz_naive(v, v, y, prec, rnd);
r1 = arf_addmul_fmpz(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = "); fmpz_print(y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
fmpz_clear(y);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

134
arf/test/t-addmul_si.c Normal file
View file

@ -0,0 +1,134 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
#include "long_extras.h"
int
arf_addmul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_si(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("addmul_si....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
long y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
y = z_randtest(state);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_addmul_si(z, x, y, prec, rnd);
r2 = arf_addmul_si_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = %ld", y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_addmul_si_naive(v, v, y, prec, rnd);
r1 = arf_addmul_si(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = %ld", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

133
arf/test/t-addmul_ui.c Normal file
View file

@ -0,0 +1,133 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_addmul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_ui(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_add(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("addmul_ui....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
ulong y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
y = n_randtest(state);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_addmul_ui(z, x, y, prec, rnd);
r2 = arf_addmul_ui_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = %lu", y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_addmul_ui_naive(v, v, y, prec, rnd);
r1 = arf_addmul_ui(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = %lu", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

124
arf/test/t-mul_fmpz.c Normal file
View file

@ -0,0 +1,124 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_mul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
arf_t t;
int r;
arf_init(t);
arf_set_fmpz(t, y);
r = arf_mul(z, x, t, prec, rnd);
arf_clear(t);
return r;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("mul_fmpz....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
arf_t x, z, v;
fmpz_t y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
fmpz_init(y);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 10);
fmpz_randtest(y, state, 2000);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10))
prec = ARF_PREC_EXACT;
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_mul_fmpz(z, x, y, prec, rnd);
r2 = arf_mul_fmpz_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); fmpz_print(y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_mul_fmpz_naive(v, x, y, prec, rnd);
r1 = arf_mul_fmpz(x, x, y, prec, rnd);
if (!arf_equal(x, v) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); fmpz_print(y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
fmpz_clear(y);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

123
arf/test/t-mul_si.c Normal file
View file

@ -0,0 +1,123 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
#include "long_extras.h"
int
arf_mul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
int r;
arf_init(t);
arf_set_si(t, y);
r = arf_mul(z, x, t, prec, rnd);
arf_clear(t);
return r;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("mul_si....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
arf_t x, z, v;
long y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 10);
y = z_randtest(state);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0)
prec = ARF_PREC_EXACT;
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_mul_si(z, x, y, prec, rnd);
r2 = arf_mul_si_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); printf("%ld", y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_mul_si_naive(v, x, y, prec, rnd);
r1 = arf_mul_si(x, x, y, prec, rnd);
if (!arf_equal(x, v) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); printf("%ld", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

123
arf/test/t-mul_ui.c Normal file
View file

@ -0,0 +1,123 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
#include "ulong_extras.h"
int
arf_mul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
int r;
arf_init(t);
arf_set_ui(t, y);
r = arf_mul(z, x, t, prec, rnd);
arf_clear(t);
return r;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("mul_ui....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
arf_t x, z, v;
ulong y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 10);
y = n_randtest(state);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0)
prec = ARF_PREC_EXACT;
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_mul_ui(z, x, y, prec, rnd);
r2 = arf_mul_ui_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); printf("%ld", y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_mul_ui_naive(v, x, y, prec, rnd);
r1 = arf_mul_ui(x, x, y, prec, rnd);
if (!arf_equal(x, v) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); printf("%ld", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

179
arf/test/t-submul.c Normal file
View file

@ -0,0 +1,179 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_submul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("submul....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, y, z, v;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(y);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
arf_randtest_special(y, state, 2000, 100);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(y)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 5))
{
case 0:
r1 = arf_submul(z, x, y, prec, rnd);
r2 = arf_submul_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); arf_print(y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 1:
r1 = arf_submul(z, x, x, prec, rnd);
r2 = arf_submul_naive(v, x, x, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL (aliasing 1)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 2:
r2 = arf_submul_naive(v, v, v, prec, rnd);
r1 = arf_submul(z, z, z, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing 2)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
case 3:
r2 = arf_submul_naive(v, v, y, prec, rnd);
r1 = arf_submul(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing 3)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = "); arf_print(y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_submul_naive(v, x, v, prec, rnd);
r1 = arf_submul(z, x, z, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL (aliasing 4)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(y);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

135
arf/test/t-submul_fmpz.c Normal file
View file

@ -0,0 +1,135 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_submul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_fmpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("submul_fmpz....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
fmpz_t y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
fmpz_init(y);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
fmpz_randtest(y, state, 2000);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_submul_fmpz(z, x, y, prec, rnd);
r2 = arf_submul_fmpz_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); fmpz_print(y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_submul_fmpz_naive(v, v, y, prec, rnd);
r1 = arf_submul_fmpz(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = "); fmpz_print(y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
fmpz_clear(y);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

134
arf/test/t-submul_si.c Normal file
View file

@ -0,0 +1,134 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
#include "long_extras.h"
int
arf_submul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_si(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("submul_si....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
long y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
y = z_randtest(state);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_submul_si(z, x, y, prec, rnd);
r2 = arf_submul_si_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = %ld", y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_submul_si_naive(v, v, y, prec, rnd);
r1 = arf_submul_si(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = %ld", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

133
arf/test/t-submul_ui.c Normal file
View file

@ -0,0 +1,133 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "arf.h"
int
arf_submul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd)
{
arf_t t;
int inexact;
arf_init(t);
arf_mul_ui(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN);
inexact = arf_sub(z, z, t, prec, rnd);
arf_clear(t);
return inexact;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("submul_ui....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arf_t x, z, v;
ulong y;
long prec, r1, r2;
arf_rnd_t rnd;
arf_init(x);
arf_init(z);
arf_init(v);
for (iter2 = 0; iter2 < 100; iter2++)
{
arf_randtest_special(x, state, 2000, 100);
y = n_randtest(state);
arf_randtest_special(z, state, 2000, 100);
arf_set(v, z);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(z)) < 10)
{
prec = ARF_PREC_EXACT;
}
switch (n_randint(state, 4))
{
case 0: rnd = ARF_RND_DOWN; break;
case 1: rnd = ARF_RND_UP; break;
case 2: rnd = ARF_RND_FLOOR; break;
default: rnd = ARF_RND_CEIL; break;
}
switch (n_randint(state, 2))
{
case 0:
r1 = arf_submul_ui(z, x, y, prec, rnd);
r2 = arf_submul_ui_naive(v, x, y, prec, rnd);
if (!arf_equal(z, v) || r1 != r2)
{
printf("FAIL!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("x = "); arf_print(x); printf("\n\n");
printf("y = %lu", y); printf("\n\n");
printf("z = "); arf_debug(z); printf("\n\n");
printf("v = "); arf_debug(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_submul_ui_naive(v, v, y, prec, rnd);
r1 = arf_submul_ui(z, z, y, prec, rnd);
if (!arf_equal(v, z) || r1 != r2)
{
printf("FAIL (aliasing)!\n");
printf("prec = %ld, rnd = %d\n\n", prec, rnd);
printf("y = %lu", y); printf("\n\n");
printf("v = "); arf_print(v); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
}
}
arf_clear(x);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}