mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
code in progress for sin/cos
This commit is contained in:
Fredrik Johansson
parent
855389cdde
commit
696cf0b7e3
8 changed files with 3253 additions and 1 deletions
28
arb.h
28
arb.h
|
|
@ -942,6 +942,34 @@ void arb_exp_arf_bb(arb_t z, const arf_t x, long prec, int minus_one);
|
|||
int _arb_get_mpn_fixed_mod_log2(mp_ptr w, fmpz_t q, mp_limb_t * error,
|
||||
const arf_t x, mp_size_t wn);
|
||||
|
||||
/* sin/cos implementation */
|
||||
|
||||
/* only goes up to (pi/4) * 256 */
|
||||
#define ARB_SIN_COS_TAB1_NUM 203
|
||||
#define ARB_SIN_COS_TAB1_BITS 8
|
||||
#define ARB_SIN_COS_TAB1_PREC 512
|
||||
#define ARB_SIN_COS_TAB1_LIMBS (ARB_SIN_COS_TAB1_PREC / FLINT_BITS)
|
||||
|
||||
/* only goes up to log(2) * 32 */
|
||||
#define ARB_SIN_COS_TAB21_NUM 26
|
||||
#define ARB_SIN_COS_TAB21_BITS 5
|
||||
#define ARB_SIN_COS_TAB22_NUM (1 << ARB_SIN_COS_TAB22_BITS)
|
||||
#define ARB_SIN_COS_TAB22_BITS 5
|
||||
#define ARB_SIN_COS_TAB2_PREC 2560
|
||||
#define ARB_SIN_COS_TAB2_LIMBS (ARB_SIN_COS_TAB2_PREC / FLINT_BITS)
|
||||
|
||||
#define ARB_PI4_TAB_LIMBS (4608 / FLINT_BITS)
|
||||
extern const mp_limb_t arb_pi4_tab[ARB_PI4_TAB_LIMBS];
|
||||
|
||||
void _arb_sin_cos_taylor_naive(mp_ptr ysin, mp_ptr ycos, mp_limb_t * error,
|
||||
mp_srcptr x, mp_size_t xn, ulong N);
|
||||
|
||||
void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos,
|
||||
mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N);
|
||||
|
||||
int _arb_get_mpn_fixed_mod_pi4(mp_ptr w, fmpz_t q, int * octant,
|
||||
mp_limb_t * error, const arf_t x, mp_size_t wn);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
114
arb/get_mpn_fixed_mod_pi4.c
Normal file
114
arb/get_mpn_fixed_mod_pi4.c
Normal file
|
|
@ -0,0 +1,114 @@
|
|||
/*=============================================================================
|
||||
|
||||
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"
|
||||
|
||||
#define TMP_ALLOC_LIMBS(__n) TMP_ALLOC((__n) * sizeof(mp_limb_t))
|
||||
|
||||
int _arf_get_integer_mpn(mp_ptr y, mp_srcptr x, mp_size_t xn, long exp);
|
||||
|
||||
int _arf_set_mpn_fixed(arf_t z, mp_srcptr xp, mp_size_t xn, mp_size_t fixn, int negative, long prec);
|
||||
|
||||
int
|
||||
_arb_get_mpn_fixed_mod_pi4(mp_ptr w, fmpz_t q, int * octant,
|
||||
mp_limb_t * error, const arf_t x, mp_size_t wn)
|
||||
{
|
||||
mp_srcptr xp;
|
||||
mp_size_t xn;
|
||||
long exp;
|
||||
|
||||
ARF_GET_MPN_READONLY(xp, xn, x);
|
||||
exp = ARF_EXP(x);
|
||||
|
||||
if (exp <= -1)
|
||||
{
|
||||
flint_mpn_zero(w, wn);
|
||||
*error = _arf_get_integer_mpn(w, xp, xn, exp + wn * FLINT_BITS);
|
||||
*octant = 0;
|
||||
if (q != NULL)
|
||||
fmpz_zero(q);
|
||||
return 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
mp_ptr qp, rp, np;
|
||||
mp_srcptr dp;
|
||||
mp_size_t qn, rn, nn, dn, tn, alloc;
|
||||
TMP_INIT;
|
||||
|
||||
tn = ((exp + 2) + FLINT_BITS - 1) / FLINT_BITS;
|
||||
|
||||
dn = wn + tn; /* denominator */
|
||||
nn = wn + 2 * tn; /* numerator */
|
||||
qn = nn - dn + 1; /* quotient */
|
||||
rn = dn; /* remainder */
|
||||
|
||||
if (dn > ARB_PI4_TAB_LIMBS)
|
||||
return 0;
|
||||
|
||||
TMP_START;
|
||||
|
||||
alloc = qn + rn + nn;
|
||||
qp = TMP_ALLOC_LIMBS(alloc);
|
||||
rp = qp + qn;
|
||||
np = rp + rn;
|
||||
|
||||
dp = arb_pi4_tab + ARB_PI4_TAB_LIMBS - dn;
|
||||
|
||||
flint_mpn_zero(np, nn);
|
||||
_arf_get_integer_mpn(np, xp, xn, exp + dn * FLINT_BITS);
|
||||
|
||||
mpn_tdiv_qr(qp, rp, 0, np, nn, dp, dn);
|
||||
|
||||
*octant = qp[0] % 8;
|
||||
|
||||
if (*octant % 2 == 0)
|
||||
{
|
||||
flint_mpn_copyi(w, rp + tn, wn);
|
||||
*error = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
mpn_sub_n(w, dp + tn, rp + tn, wn);
|
||||
*error = 3;
|
||||
}
|
||||
|
||||
if (q != NULL)
|
||||
{
|
||||
/* read the exponent */
|
||||
while (qn > 1 && qp[qn-1] == 0)
|
||||
qn--;
|
||||
|
||||
if (qn == 1)
|
||||
fmpz_set_ui(q, qp[0]);
|
||||
else
|
||||
fmpz_set_mpn_large(q, qp, qn, 0);
|
||||
}
|
||||
|
||||
TMP_END;
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
|
||||
2559
arb/sin_cos_tab.c
Normal file
2559
arb/sin_cos_tab.c
Normal file
File diff suppressed because it is too large
Load diff
97
arb/sin_cos_taylor_naive.c
Normal file
97
arb/sin_cos_taylor_naive.c
Normal file
|
|
@ -0,0 +1,97 @@
|
|||
/*=============================================================================
|
||||
|
||||
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"
|
||||
|
||||
void
|
||||
_arb_sin_cos_taylor_naive(mp_ptr ysin, mp_ptr ycos, mp_limb_t * error,
|
||||
mp_srcptr x, mp_size_t xn, ulong N)
|
||||
{
|
||||
ulong k;
|
||||
mp_ptr s, s2, t, u, v;
|
||||
mp_size_t nn = xn + 1;
|
||||
|
||||
if (N == 0)
|
||||
{
|
||||
flint_mpn_zero(ysin, xn);
|
||||
flint_mpn_zero(ycos, xn);
|
||||
error[0] = 0;
|
||||
return;
|
||||
}
|
||||
|
||||
s = flint_malloc(sizeof(mp_limb_t) * (nn + 1));
|
||||
s2 = flint_malloc(sizeof(mp_limb_t) * (nn + 1));
|
||||
t = flint_malloc(sizeof(mp_limb_t) * nn);
|
||||
v = flint_malloc(sizeof(mp_limb_t) * nn);
|
||||
u = flint_malloc(sizeof(mp_limb_t) * 2 * nn);
|
||||
|
||||
/* s = 1 */
|
||||
flint_mpn_zero(s, nn);
|
||||
s[nn] = 1;
|
||||
/* s2 = 0 */
|
||||
flint_mpn_zero(s2, nn + 1);
|
||||
|
||||
/* t = v = x */
|
||||
flint_mpn_zero(t, nn);
|
||||
flint_mpn_copyi(t + 1, x, xn);
|
||||
flint_mpn_copyi(v, t, nn);
|
||||
|
||||
for (k = 1; k < 2 * N; k++)
|
||||
{
|
||||
if (k % 4 == 0)
|
||||
s[nn] += mpn_add_n(s, s, t, nn);
|
||||
else if (k % 4 == 1)
|
||||
s2[nn] += mpn_add_n(s2, s2, t, nn);
|
||||
else if (k % 4 == 2)
|
||||
s[nn] -= mpn_sub_n(s, s, t, nn);
|
||||
else
|
||||
s2[nn] -= mpn_sub_n(s2, s2, t, nn);
|
||||
|
||||
/* t = t * x / (k + 1) */
|
||||
mpn_mul_n(u, t, v, nn);
|
||||
flint_mpn_copyi(t, u + nn, nn);
|
||||
mpn_divrem_1(t, 0, t, nn, k + 1);
|
||||
}
|
||||
|
||||
if (s[nn] != 0)
|
||||
{
|
||||
flint_mpn_store(ycos, xn, LIMB_ONES);
|
||||
flint_mpn_copyi(ysin, s2 + 1, xn);
|
||||
}
|
||||
else
|
||||
{
|
||||
flint_mpn_copyi(ycos, s + 1, xn);
|
||||
flint_mpn_copyi(ysin, s2 + 1, xn);
|
||||
}
|
||||
|
||||
error[0] = 2;
|
||||
|
||||
flint_free(s);
|
||||
flint_free(s2);
|
||||
flint_free(t);
|
||||
flint_free(u);
|
||||
flint_free(v);
|
||||
}
|
||||
|
||||
178
arb/sin_cos_taylor_rs.c
Normal file
178
arb/sin_cos_taylor_rs.c
Normal file
|
|
@ -0,0 +1,178 @@
|
|||
/*=============================================================================
|
||||
|
||||
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"
|
||||
|
||||
#define TMP_ALLOC_LIMBS(size) TMP_ALLOC((size) * sizeof(mp_limb_t))
|
||||
|
||||
#define FACTORIAL_TAB_SIZE 288
|
||||
|
||||
extern const mp_limb_t factorial_tab_numer[FACTORIAL_TAB_SIZE];
|
||||
extern const mp_limb_t factorial_tab_denom[FACTORIAL_TAB_SIZE];
|
||||
|
||||
void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos,
|
||||
mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N)
|
||||
{
|
||||
mp_ptr s, t, xpow;
|
||||
mp_limb_t new_denom, old_denom, c;
|
||||
long power, k, m;
|
||||
int cosorsin;
|
||||
|
||||
TMP_INIT;
|
||||
TMP_START;
|
||||
|
||||
if (2 * N >= FACTORIAL_TAB_SIZE - 1)
|
||||
{
|
||||
printf("_arb_sin_cos_taylor_rs: N too large!\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
if (N <= 1)
|
||||
{
|
||||
if (N == 0)
|
||||
{
|
||||
flint_mpn_zero(ysin, xn);
|
||||
flint_mpn_zero(ycos, xn);
|
||||
error[0] = 0;
|
||||
}
|
||||
else if (N == 1)
|
||||
{
|
||||
flint_mpn_store(ycos, xn, LIMB_ONES);
|
||||
flint_mpn_copyi(ysin, x, xn);
|
||||
error[0] = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Choose m ~= sqrt(num_terms) (m must be even, >= 2) */
|
||||
m = 2;
|
||||
while (m * m < N)
|
||||
m += 2;
|
||||
|
||||
/* todo: merge allocations */
|
||||
xpow = TMP_ALLOC_LIMBS((m + 1) * xn);
|
||||
s = TMP_ALLOC_LIMBS(xn + 2);
|
||||
t = TMP_ALLOC_LIMBS(2 * xn + 2); /* todo: 1 limb too much? */
|
||||
|
||||
/* higher index ---> */
|
||||
/* | ---xn--- | */
|
||||
/* xpow = | <temp> | x^m | x^(m-1) | ... | x^2 | x | */
|
||||
|
||||
#define XPOW_WRITE(__k) (xpow + (m - (__k)) * xn)
|
||||
#define XPOW_READ(__k) (xpow + (m - (__k) + 1) * xn)
|
||||
|
||||
mpn_sqr(XPOW_WRITE(1), x, xn);
|
||||
mpn_sqr(XPOW_WRITE(2), XPOW_READ(1), xn);
|
||||
|
||||
for (k = 4; k <= m; k += 2)
|
||||
{
|
||||
mpn_mul_n(XPOW_WRITE(k - 1), XPOW_READ(k / 2), XPOW_READ(k / 2 - 1), xn);
|
||||
mpn_sqr(XPOW_WRITE(k), XPOW_READ(k / 2), xn);
|
||||
}
|
||||
|
||||
for (cosorsin = 0; cosorsin < 2; cosorsin++)
|
||||
{
|
||||
flint_mpn_zero(s, xn + 1);
|
||||
|
||||
/* todo: skip one nonscalar multiplication (use x^m)
|
||||
when starting on x^0 */
|
||||
power = (N - 1) % m;
|
||||
|
||||
for (k = N - 1; k >= 0; k--)
|
||||
{
|
||||
c = factorial_tab_numer[2 * k + cosorsin];
|
||||
new_denom = factorial_tab_denom[2 * k + cosorsin];
|
||||
old_denom = factorial_tab_denom[2 * k + cosorsin + 2];
|
||||
|
||||
/* change denominators */
|
||||
if (new_denom != old_denom && k < N - 1)
|
||||
{
|
||||
if (k % 2 == 0)
|
||||
/* mpn_neg_n(s, s, xn + 1); */
|
||||
s[xn] += old_denom;
|
||||
|
||||
mpn_divrem_1(s, 0, s, xn + 1, old_denom);
|
||||
|
||||
if (k % 2 == 0)
|
||||
s[xn] -= 1;
|
||||
/* mpn_neg_n(s, s, xn + 1); */
|
||||
}
|
||||
|
||||
if (power == 0)
|
||||
{
|
||||
/* add c * x^0 -- only top limb is affected */
|
||||
if (k % 2 == 0)
|
||||
s[xn] += c;
|
||||
else
|
||||
s[xn] -= c;
|
||||
|
||||
/* Outer polynomial evaluation: multiply by x^m */
|
||||
if (k != 0)
|
||||
{
|
||||
mpn_mul(t, s, xn + 1, XPOW_READ(m), xn);
|
||||
flint_mpn_copyi(s, t + xn, xn + 1);
|
||||
}
|
||||
|
||||
power = m - 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (k % 2 == 0)
|
||||
s[xn] += mpn_addmul_1(s, XPOW_READ(power), xn, c);
|
||||
else
|
||||
s[xn] -= mpn_submul_1(s, XPOW_READ(power), xn, c);
|
||||
|
||||
power--;
|
||||
}
|
||||
}
|
||||
|
||||
/* finally divide by denominator */
|
||||
if (cosorsin == 0)
|
||||
{
|
||||
mpn_divrem_1(t, 0, s, xn + 1, factorial_tab_denom[0]);
|
||||
|
||||
/* perturb down to a number < 1 if necessary. note that this
|
||||
does not invalidate the error bound: 1 - ulp is either
|
||||
1 ulp too small or must be closer to the exact value */
|
||||
if (t[xn] == 0)
|
||||
flint_mpn_copyi(ycos, t, xn);
|
||||
else
|
||||
flint_mpn_store(ycos, xn, LIMB_ONES);
|
||||
}
|
||||
else
|
||||
{
|
||||
mpn_divrem_1(s, 0, s, xn + 1, factorial_tab_denom[0]);
|
||||
mpn_mul(t, s, xn + 1, x, xn);
|
||||
flint_mpn_copyi(ysin, t + xn, xn);
|
||||
}
|
||||
}
|
||||
|
||||
/* error bound (ulp) */
|
||||
error[0] = 2;
|
||||
}
|
||||
|
||||
TMP_END;
|
||||
}
|
||||
|
||||
|
|
@ -32,7 +32,7 @@ int main()
|
|||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("arb_get_mpn_fixed_mod_log2....");
|
||||
printf("get_mpn_fixed_mod_log2....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
|
|
|||
135
arb/test/t-get_mpn_fixed_mod_pi4.c
Normal file
135
arb/test/t-get_mpn_fixed_mod_pi4.c
Normal 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) 2014 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "arb.h"
|
||||
|
||||
int _arf_set_mpn_fixed(arf_t z, mp_srcptr xp, mp_size_t xn, mp_size_t fixn, int negative, long prec);
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("get_mpn_fixed_mod_pi4....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
/* _flint_rand_init_gmp(state); */
|
||||
|
||||
for (iter = 0; iter < 100000; iter++)
|
||||
{
|
||||
arf_t x;
|
||||
int octant;
|
||||
fmpz_t q;
|
||||
mp_ptr w;
|
||||
arb_t wb, t, u;
|
||||
mp_size_t wn;
|
||||
long prec, prec2;
|
||||
int success;
|
||||
mp_limb_t error;
|
||||
|
||||
prec = 2 + n_randint(state, 10000);
|
||||
wn = 1 + n_randint(state, 200);
|
||||
prec2 = FLINT_MAX(prec, wn * FLINT_BITS) + 100;
|
||||
|
||||
arf_init(x);
|
||||
arb_init(wb);
|
||||
arb_init(t);
|
||||
arb_init(u);
|
||||
fmpz_init(q);
|
||||
w = flint_malloc(sizeof(mp_limb_t) * wn);
|
||||
|
||||
arf_randtest(x, state, prec, 14);
|
||||
|
||||
/* this should generate numbers close to multiples of pi/4 */
|
||||
if (n_randint(state, 4) == 0)
|
||||
{
|
||||
arb_const_pi(t, prec);
|
||||
arb_mul_2exp_si(t, t, -2);
|
||||
fmpz_randtest(q, state, 200);
|
||||
arb_mul_fmpz(t, t, q, prec);
|
||||
arf_add(x, x, arb_midref(t), prec, ARF_RND_DOWN);
|
||||
}
|
||||
|
||||
arf_abs(x, x);
|
||||
|
||||
success = _arb_get_mpn_fixed_mod_pi4(w, q, &octant, &error, x, wn);
|
||||
|
||||
if (success)
|
||||
{
|
||||
/* could round differently */
|
||||
if (fmpz_fdiv_ui(q, 8) != octant)
|
||||
{
|
||||
printf("bad octant\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
_arf_set_mpn_fixed(arb_midref(wb), w, wn, wn, 0, FLINT_BITS * wn);
|
||||
mag_set_ui_2exp_si(arb_radref(wb), error, -FLINT_BITS * wn);
|
||||
|
||||
arb_const_pi(u, prec2);
|
||||
arb_mul_2exp_si(u, u, -2);
|
||||
arb_set(t, wb);
|
||||
if (octant % 2 == 1)
|
||||
arb_sub(t, u, t, prec2);
|
||||
arb_addmul_fmpz(t, u, q, prec2);
|
||||
|
||||
if (!arb_contains_arf(t, x))
|
||||
{
|
||||
printf("FAIL (containment)\n");
|
||||
printf("x = "); arf_printd(x, 50); printf("\n\n");
|
||||
printf("q = "); fmpz_print(q); printf("\n\n");
|
||||
printf("w = "); arb_printd(wb, 50); printf("\n\n");
|
||||
printf("t = "); arb_printd(t, 50); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_const_pi(t, prec2);
|
||||
arb_mul_2exp_si(t, t, -2);
|
||||
|
||||
if (arf_sgn(arb_midref(wb)) < 0 ||
|
||||
arf_cmp(arb_midref(wb), arb_midref(t)) >= 0)
|
||||
{
|
||||
printf("FAIL (expected 0 <= w < pi/4)\n");
|
||||
printf("x = "); arf_printd(x, 50); printf("\n\n");
|
||||
printf("w = "); arb_printd(wb, 50); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
flint_free(w);
|
||||
fmpz_clear(q);
|
||||
arf_clear(x);
|
||||
arb_clear(wb);
|
||||
arb_clear(t);
|
||||
arb_clear(u);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
141
arb/test/t-sin_cos_taylor_rf.c
Normal file
141
arb/test/t-sin_cos_taylor_rf.c
Normal file
|
|
@ -0,0 +1,141 @@
|
|||
/*=============================================================================
|
||||
|
||||
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"
|
||||
#include "mpn_extras.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("sin_cos_taylor_rf....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
_flint_rand_init_gmp(state);
|
||||
|
||||
for (iter = 0; iter < 100000; iter++)
|
||||
{
|
||||
mp_ptr x, y1s, y1c, y2s, y2c, t;
|
||||
mp_limb_t err1, err2;
|
||||
ulong N;
|
||||
mp_size_t xn;
|
||||
int cmp, result;
|
||||
|
||||
N = n_randint(state, 144 - 1);
|
||||
xn = 1 + n_randint(state, 20);
|
||||
|
||||
x = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
y1s = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
y1c = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
y2s = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
y2c = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
t = flint_malloc(sizeof(mp_limb_t) * xn);
|
||||
|
||||
flint_mpn_rrandom(x, state->gmp_state, xn);
|
||||
flint_mpn_rrandom(y1s, state->gmp_state, xn);
|
||||
flint_mpn_rrandom(y1c, state->gmp_state, xn);
|
||||
flint_mpn_rrandom(y2s, state->gmp_state, xn);
|
||||
flint_mpn_rrandom(y2c, state->gmp_state, xn);
|
||||
x[xn - 1] &= (LIMB_ONES >> 4);
|
||||
|
||||
_arb_sin_cos_taylor_naive(y1s, y1c, &err1, x, xn, N);
|
||||
_arb_sin_cos_taylor_rs(y2s, y2c, &err2, x, xn, N);
|
||||
|
||||
cmp = mpn_cmp(y1s, y2s, xn);
|
||||
|
||||
if (cmp == 0)
|
||||
{
|
||||
result = 1;
|
||||
}
|
||||
else if (cmp > 0)
|
||||
{
|
||||
mpn_sub_n(t, y1s, y2s, xn);
|
||||
result = flint_mpn_zero_p(t + 1, xn - 1) && (t[0] <= err2);
|
||||
}
|
||||
else
|
||||
{
|
||||
mpn_sub_n(t, y2s, y1s, xn);
|
||||
result = flint_mpn_zero_p(t + 1, xn - 1) && (t[0] <= err2);
|
||||
}
|
||||
|
||||
if (!result)
|
||||
{
|
||||
printf("FAIL\n");
|
||||
printf("N = %ld xn = %ld\n", N, xn);
|
||||
printf("x =");
|
||||
flint_mpn_debug(x, xn);
|
||||
printf("y1s =");
|
||||
flint_mpn_debug(y1s, xn);
|
||||
printf("y2s =");
|
||||
flint_mpn_debug(y2s, xn);
|
||||
abort();
|
||||
}
|
||||
|
||||
cmp = mpn_cmp(y1c, y2c, xn);
|
||||
|
||||
if (cmp == 0)
|
||||
{
|
||||
result = 1;
|
||||
}
|
||||
else if (cmp > 0)
|
||||
{
|
||||
mpn_sub_n(t, y1c, y2c, xn);
|
||||
result = flint_mpn_zero_p(t + 1, xn - 1) && (t[0] <= err2);
|
||||
}
|
||||
else
|
||||
{
|
||||
mpn_sub_n(t, y2c, y1c, xn);
|
||||
result = flint_mpn_zero_p(t + 1, xn - 1) && (t[0] <= err2);
|
||||
}
|
||||
|
||||
if (!result)
|
||||
{
|
||||
printf("FAIL\n");
|
||||
printf("N = %ld xn = %ld\n", N, xn);
|
||||
printf("x =");
|
||||
flint_mpn_debug(x, xn);
|
||||
printf("y1c =");
|
||||
flint_mpn_debug(y1c, xn);
|
||||
printf("y2c =");
|
||||
flint_mpn_debug(y2c, xn);
|
||||
abort();
|
||||
}
|
||||
|
||||
flint_free(x);
|
||||
flint_free(y1s);
|
||||
flint_free(y1c);
|
||||
flint_free(y2s);
|
||||
flint_free(y2c);
|
||||
flint_free(t);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
Loading…
Add table
Reference in a new issue