Keskustelu: Ohjelmointikysymykset: C: Ison tarkkuuden lukujen ohjelman käännös

#include "fmprb.h"
#include "fmpcb.h"
#include "fmprb_poly.h"
#include "zeta.h"
#include "profiler.h"

void log_zeta(fmprb_struct * z, long n, long prec)
{
    fmpcb_struct * w;
    fmprb_struct * t;
    fmpcb_t s, a;
    long i;

    t = _fmprb_vec_init(n);
    w = _fmpcb_vec_init(n);

    fmpcb_init(s);
    fmpcb_init(a);

    fmpcb_zero(s);
    fmpcb_one(a);

    zeta_series(w, s, a, 0, n, prec);

    for (i = 0; i < n; i++)
        fmprb_neg(t + i, fmpcb_realref(w + i));

    _fmprb_poly_log_series(z, t, n, prec);

    fmpcb_clear(s);
    fmpcb_clear(a);

    _fmpcb_vec_clear(w, n);
    _fmprb_vec_clear(t, n);
}

void log_gamma(fmprb_struct * z, long n, long prec)
{
    long i;

    if (n > 0) fmprb_zero(z);
    if (n > 1) fmprb_const_euler(z + 1, prec);
    if (n > 2) zeta_ui_vec(z + 2, 2, n - 2, prec);

    for (i = 2; i < n; i++)
        fmprb_div_ui(z + i, z + i, i, prec);

    for (i = 1; i < n; i += 2)
        fmprb_neg(z + i, z + i);

    for (i = 0; i < n; i++)
        fmprb_mul_2exp_si(z + i, z + i, -i);
}

void log_1s(fmprb_struct * z, long n, long prec)
{
    long i;

    if (n > 0) fmprb_zero(z);

    if (n > 1)
    {
        fmprb_const_pi(z + 1, prec);
        fmprb_log(z + 1, z + 1, prec);
        fmprb_mul_2exp_si(z + 1, z + 1, -1);
        fmprb_add_ui(z + 1, z + 1, 1, prec);
    }

    for (i = 2; i < n; i++)
    {
        fmprb_one(z + i);
        fmprb_div_ui(z + i, z + i, i, prec);
    }

    _fmprb_vec_neg(z, z, n);
}

int main()
{
    fmprb_struct *z, *t, *a;
    timeit_t t0;

    long i, n, prec;

    n = 100 + 1;
    prec = 6 * n;  /* must choose wisely */

    z = _fmprb_vec_init(n);
    t = _fmprb_vec_init(n);
    a = _fmprb_vec_init(n);

    timeit_start(t0);
    log_zeta(z, n, prec);
    timeit_stop(t0);
    printf("zeta: %ld ms\n", t0->cpu);

    timeit_start(t0);
    log_gamma(t, n, prec);
    timeit_stop(t0);
    printf("gamma: %ld ms\n", t0->cpu);

    _fmprb_vec_add(z, z, t, n, prec);
    log_1s(t, n, prec);
    _fmprb_vec_add(z, z, t, n, prec);

    for (i = 0; i < n; i++)
        fmprb_set_si(t + i, (i == 0) ? 0 : -1);

    timeit_start(t0);
    _fmprb_poly_compose_series(a, z, n, t, n, n, prec);
    timeit_stop(t0);
    printf("composition: %ld ms\n", t0->cpu);

    for (i = 0; i < n; i++)
    {
        printf("%ld: ", i); fmprb_printd(a + i, 15);
        printf("\n");
    }

    _fmprb_vec_clear(z, n);
    _fmprb_vec_clear(t, n);
    _fmprb_vec_clear(a, n);
}