Skip to content

Instantly share code, notes, and snippets.

@fredrik-johansson
Created October 3, 2013 13:27
Show Gist options
  • Save fredrik-johansson/6809808 to your computer and use it in GitHub Desktop.
Save fredrik-johansson/6809808 to your computer and use it in GitHub Desktop.
example code for numerical integration
#include "fmpcb_calc.h"
#include "profiler.h"
int
sinx(fmpcb_ptr out, const fmpcb_t inp, void * params, long order, long prec)
{
int xlen = FLINT_MIN(2, order);
fmpcb_set(out, inp);
if (xlen > 1)
fmpcb_one(out + 1);
_fmpcb_poly_sin_series(out, out, xlen, order, prec);
return 0;
}
int
elliptic(fmpcb_ptr out, const fmpcb_t inp, void * params, long order, long prec)
{
fmpcb_ptr t;
t = _fmpcb_vec_init(order);
fmpcb_set(t, inp);
if (order > 1)
fmpcb_one(t + 1);
_fmpcb_poly_sin_series(t, t, FLINT_MIN(2, order), order, prec);
_fmpcb_poly_mullow(out, t, order, t, order, order, prec);
_fmpcb_vec_scalar_mul_2exp_si(t, out, order, -1);
fmpcb_sub_ui(t, t, 1, prec);
_fmpcb_vec_neg(t, t, order);
_fmpcb_poly_rsqrt_series(out, t, order, order, prec);
_fmpcb_vec_clear(t, order);
return 0;
}
int
bessel(fmpcb_ptr out, const fmpcb_t inp, void * params, long order, long prec)
{
fmpcb_ptr t;
fmpcb_t z;
ulong n;
t = _fmpcb_vec_init(order);
fmpcb_init(z);
fmpcb_set(t, inp);
if (order > 1)
fmpcb_one(t + 1);
n = 10;
fmprb_set_si(fmpcb_realref(z), 20);
fmprb_set_si(fmpcb_imagref(z), 10);
/* z sin(t) */
_fmpcb_poly_sin_series(out, t, FLINT_MIN(2, order), order, prec);
_fmpcb_vec_scalar_mul(out, out, order, z, prec);
/* t n */
_fmpcb_vec_scalar_mul_ui(t, t, FLINT_MIN(2, order), n, prec);
_fmpcb_poly_sub(out, t, FLINT_MIN(2, order), out, order, prec);
_fmpcb_poly_cos_series(out, out, order, order, prec);
_fmpcb_vec_clear(t, order);
fmpcb_clear(z);
return 0;
}
int main()
{
fmpcb_t r, s, a, b;
fmpr_t inr, outr;
long digits, prec;
fmpcb_init(r);
fmpcb_init(s);
fmpcb_init(a);
fmpcb_init(b);
fmpr_init(inr);
fmpr_init(outr);
/* fmprb_calc_verbose = 1; */
digits = 300;
prec = digits * 3.32193;
fmpr_set_d(inr, 0.125);
fmpr_set_d(outr, 1.0);
TIMEIT_ONCE_START
fmpcb_set_si(a, 0);
fmpcb_set_si(b, 100);
fmpcb_calc_integrate_taylor(r, sinx, NULL, a, b, inr, outr, prec, 1.1 * prec);
printf("\n");
fmpcb_printd(r, digits); printf("\n");
TIMEIT_ONCE_STOP
fmpr_set_d(inr, 0.2);
fmpr_set_d(outr, 0.5);
TIMEIT_ONCE_START
fmpcb_set_si(a, 0);
fmpcb_set_si(b, 6);
fmpcb_calc_integrate_taylor(r, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
fmpcb_set_si(a, 6);
fmprb_set_si(fmpcb_realref(b), 6);
fmprb_set_si(fmpcb_imagref(b), 6);
fmpcb_calc_integrate_taylor(s, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
fmpcb_add(r, r, s, prec);
printf("\n");
fmpcb_printd(r, digits); printf("\n");
TIMEIT_ONCE_STOP
fmpr_set_d(inr, 0.1);
fmpr_set_d(outr, 0.5);
TIMEIT_ONCE_START
fmpcb_set_si(a, 0);
fmpcb_const_pi(b, 3 * prec);
fmpcb_calc_integrate_taylor(r, bessel, NULL, a, b, inr, outr, prec, 3 * prec);
fmpcb_div(r, r, b, prec);
printf("\n");
fmpcb_printd(r, digits); printf("\n");
TIMEIT_ONCE_STOP
fmpcb_clear(r);
fmpcb_clear(s);
fmpcb_clear(a);
fmpcb_clear(b);
fmpr_clear(inr);
fmpr_clear(outr);
flint_cleanup();
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment