flint-sys 0.9.0

Bindings to the FLINT C library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

/*
    Copyright (C) 2010,2011 Fredrik Johansson

    This file is part of FLINT.

    FLINT is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
*/

#include "test_helpers.h"
#include "gmpcompat.h"
#include "mpn_extras.h"
#include "fmpq.h"
#include "arith.h"

void numerical_test(fmpq_t res, slong n, double ans)
{
    const double tol = 1e-13;
    double err;

    mpq_t tmp;
    mpq_init(tmp);

    arith_harmonic_number(res, n);
    fmpq_get_mpq(tmp, res);
    err = mpq_get_d(tmp) - ans;
    err = FLINT_ABS(err);

    if (err > tol)
    {
        flint_printf("FAIL: %wd %.16f %.16f\n", n, mpq_get_d(tmp), ans);
        fflush(stdout);
        flint_abort();
    }

    mpq_clear(tmp);
}

void
mpq_harmonic_balanced(mpq_t res, slong a, slong b)
{
    slong k;
    mpq_t t;

    mpq_init(t);

    if (b - a < 50)
    {
        flint_mpq_set_ui(res, 0, UWORD(1));
        for (k = a; k <= b; k++)
        {
            flint_mpq_set_ui(t, UWORD(1), k);
            mpq_add(res, res, t);
        }
    }
    else
    {
        mpq_harmonic_balanced(res, a, (a+b)/2);
        mpq_harmonic_balanced(t, (a+b)/2+1, b);
        mpq_add(res, res, t);
    }

    mpq_clear(t);
}

TEST_FUNCTION_START(arith_harmonic_number, state)
{
    slong i;
    mpq_t x, y;
    fmpq_t t;


    fmpq_init(t);
    mpq_init(x);
    mpq_init(y);

    for (i = -2; i < 1000; i++)
    {
        mpq_harmonic_balanced(x, 1, i);
        arith_harmonic_number(t, i);
        fmpq_get_mpq(y, t);

        if (!mpq_equal(x, y))
        {
            flint_printf("FAIL: %wd\n", i);
            fflush(stdout);
            flint_abort();
        }
    }

    numerical_test(t, 1000, 7.4854708605503449127);
    numerical_test(t, 1001, 7.4864698615493459117);
    numerical_test(t, 1002, 7.4874678655413618797);
    numerical_test(t, 1003, 7.4884648745144426375);

    numerical_test(t, 10000, 9.7876060360443822642);
    numerical_test(t, 10001, 9.7877060260453821642);
    numerical_test(t, 10002, 9.7878060060493813643);
    numerical_test(t, 10003, 9.7879059760583786652);
    numerical_test(t, 10004, 9.7880059360743722677);

    numerical_test(t, 20000, 10.480728217229327573);
    numerical_test(t, 30000, 10.886184992119899362);
    numerical_test(t, 40000, 11.173862897945522882);
    numerical_test(t, 50000, 11.397003949278482638);
    numerical_test(t, 60000, 11.579323839415955783);
    numerical_test(t, 70000, 11.733473328773164956);
    numerical_test(t, 80000, 11.867003828544530692);
    numerical_test(t, 90000, 11.984786169759202469);

    numerical_test(t, 100000, 12.090146129863427947);
    numerical_test(t, 100001, 12.090156129763428947);
    numerical_test(t, 100002, 12.090166129563432947);
    numerical_test(t, 100003, 12.090176129263441947);
    numerical_test(t, 100004, 12.090186128863457946);

    numerical_test(t, 300000, 13.188755085205611713);
    numerical_test(t, 500000, 13.699580042305528322);
    numerical_test(t, 700000, 14.036051993212618803);
    numerical_test(t, 900000, 14.287366262763433338);

    mpq_clear(x);
    mpq_clear(y);
    fmpq_clear(t);

    TEST_FUNCTION_END(state);
}