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

add sub

Browse source
This commit is contained in:
Fredrik Johansson 2014-05-02 17:59:34 +02:00
parent 81055070cf
commit 3a16b8587d
3 changed files with 312 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

42
arf.h
View file

@ -777,6 +777,46 @@ _fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
r3 += r2 < t1; \
} while (0)
#define ARF_MPN_MUL(_z, _x, _xn, _y, _yn) \
if ((_xn) == (_yn)) \
{ \
if ((_xn) == 1) \
{ \
umul_ppmm((_z)[1], (_z)[0], (_x)[0], (_y)[0]); \
} \
else if ((_xn) == 2) \
{ \
mp_limb_t __x1, __x0, __y1, __y0; \
__x0 = (_x)[0]; \
__x1 = (_x)[1]; \
__y0 = (_y)[0]; \
__y1 = (_y)[1]; \
nn_mul_2x2((_z)[3], (_z)[2], (_z)[1], (_z)[0], __x1, __x0, __y1, __y0); \
} \
else if ((_x) == (_y)) \
{ \
mpn_sqr((_z), (_x), (_xn)); \
} \
else \
{ \
mpn_mul_n((_z), (_x), (_y), (_xn)); \
} \
} \
else if ((_xn) > (_yn)) \
{ \
if ((_yn) == 1) \
(_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_x), (_xn), (_y)[0]); \
else \
mpn_mul((_z), (_x), (_xn), (_y), (_yn)); \
} \
else \
{ \
if ((_xn) == 1) \
(_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_y), (_yn), (_x)[0]); \
else \
mpn_mul((_z), (_y), (_yn), (_x), (_xn)); \
}
#define ARF_MUL_STACK_ALLOC 40
#define ARF_MUL_TLS_ALLOC 1000
@ -872,6 +912,8 @@ int _arf_add_mpn(arf_t z, mp_srcptr xp, mp_size_t xn, int xsgnbit,
int arf_add(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
int arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
#ifdef __cplusplus
}
#endif

88
arf/sub.c Normal file
View file

@ -0,0 +1,88 @@
/*=============================================================================
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_sub_special(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
{
if (arf_is_zero(x))
{
if (arf_is_zero(y))
{
arf_zero(z);
return 0;
}
else
return arf_neg_round(z, y, prec, rnd);
}
else if (arf_is_zero(y))
{
return arf_set_round(z, x, prec, rnd);
}
else if (arf_is_nan(x) || arf_is_nan(y)
|| (arf_is_pos_inf(x) && arf_is_pos_inf(y))
|| (arf_is_neg_inf(x) && arf_is_neg_inf(y)))
{
arf_nan(z);
return 0;
}
else if (arf_is_special(x))
{
arf_set(z, x);
return 0;
}
else
{
arf_neg(z, y);
return 0;
}
}
int
arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd)
{
mp_size_t xn, yn;
mp_srcptr xptr, yptr;
long shift;
if (arf_is_special(x) || arf_is_special(y))
{
return arf_sub_special(z, x, y, prec, rnd);
}
shift = _fmpz_sub_small(ARF_EXPREF(x), ARF_EXPREF(y));
ARF_GET_MPN_READONLY(xptr, xn, x);
ARF_GET_MPN_READONLY(yptr, yn, y);
if (shift >= 0)
return _arf_add_mpn(z, xptr, xn, ARF_SGNBIT(x), ARF_EXPREF(x),
yptr, yn, ARF_SGNBIT(y) ^ 1, shift, prec, rnd);
else
return _arf_add_mpn(z, yptr, yn, ARF_SGNBIT(y) ^ 1, ARF_EXPREF(y),
xptr, xn, ARF_SGNBIT(x), -shift, prec, rnd);
}

182
arf/test/t-sub.c Normal file
View file

@ -0,0 +1,182 @@
/*=============================================================================
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_sub_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
{
fmpr_t a, b;
long r;
fmpr_init(a);
fmpr_init(b);
arf_get_fmpr(a, x);
arf_get_fmpr(b, y);
r = fmpr_sub(a, a, b, prec, rnd);
arf_set_fmpr(z, a);
fmpr_clear(a);
fmpr_clear(b);
return (r == FMPR_RESULT_EXACT) ? 0 : 1;
}
int main()
{
long iter, iter2;
flint_rand_t state;
printf("sub....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; 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);
prec = 2 + n_randint(state, 2000);
if (n_randint(state, 10) == 0 &&
fmpz_bits(ARF_EXPREF(x)) < 10 &&
fmpz_bits(ARF_EXPREF(y)) < 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_sub(z, x, y, prec, rnd);
r2 = arf_sub_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_sub(z, x, x, prec, rnd);
r2 = arf_sub_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_sub_naive(v, x, x, prec, rnd);
r1 = arf_sub(x, x, x, prec, rnd);
if (!arf_equal(v, x) || r1 != r2)
{
printf("FAIL (aliasing 2)!\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 3:
r2 = arf_sub_naive(v, x, y, prec, rnd);
r1 = arf_sub(x, x, y, prec, rnd);
if (!arf_equal(x, v) || r1 != r2)
{
printf("FAIL (aliasing 3)!\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("v = "); arf_print(v); printf("\n\n");
printf("r1 = %ld, r2 = %ld\n", r1, r2);
abort();
}
break;
default:
r2 = arf_sub_naive(v, y, x, prec, rnd);
r1 = arf_sub(x, y, x, prec, rnd);
if (!arf_equal(x, 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("y = "); arf_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(y);
arf_clear(z);
arf_clear(v);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}