summaryrefslogtreecommitdiff
path: root/lib/math/ancillary/log-overkill-constants.c
blob: 1af44806141f892ea4f73cdfb6dda94ef41b37b5 (plain)
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
/* cc -o log-overkill-constants log-overkill-constants.c -lmpfr */
#include <stdint.h>
#include <stdio.h>
#include <time.h>

#include <mpfr.h>

#define FLT64_TO_UINT64(f) (*((uint64_t *) ((char *) &f)))
#define UINT64_TO_FLT64(u) (*((double *) ((char *) &u)))

int
main(void)
{
        mpfr_t t1;
        mpfr_t t2;
        mpfr_t t3;
        double d = 0;
        uint64_t u = 0;
        size_t j = 0;
        long int k = 0;

        mpfr_init2(t1, 10000);
        mpfr_init2(t2, 10000);
        mpfr_init2(t3, 10000);

        /* C1 */
        for (k = -5; k >= -20; k -= 5) {
                printf(
                        "const C%ld : (uint64, uint64, uint64, uint64)[32] = [\n",
                        (k /
                         (
                                 -5)));

                for (j = 0; j < 32; ++j) {
                        mpfr_set_si(t1, 1, MPFR_RNDN);
                        mpfr_set_si(t2, k, MPFR_RNDN);
                        mpfr_exp2(t2, t2, MPFR_RNDN);
                        mpfr_mul_si(t2, t2, j, MPFR_RNDN);
                        mpfr_add(t1, t1, t2, MPFR_RNDN);

                        /* first, log(1 + ...) */
                        mpfr_log(t2, t1, MPFR_RNDN);
                        d = mpfr_get_d(t2, MPFR_RNDN);
                        u = FLT64_TO_UINT64(d);
                        printf("	(%#018lx, ", u);
                        mpfr_set_d(t3, d, MPFR_RNDN);
                        mpfr_sub(t2, t2, t3, MPFR_RNDN);
                        d = mpfr_get_d(t2, MPFR_RNDN);
                        u = FLT64_TO_UINT64(d);
                        printf("%#018lx, ", u);

                        /* now, 1/(1 + ...) */
                        mpfr_pow_si(t2, t1, -1, MPFR_RNDN);
                        d = mpfr_get_d(t2, MPFR_RNDN);
                        u = FLT64_TO_UINT64(d);
                        printf("    %#018lx, ", u);
                        mpfr_set_d(t3, d, MPFR_RNDN);
                        mpfr_sub(t2, t2, t3, MPFR_RNDN);
                        d = mpfr_get_d(t2, MPFR_RNDN);
                        u = FLT64_TO_UINT64(d);
                        printf("%#018lx), /* j = %zu */\n", u, j);
                }

                printf("]\n\n");
        }

        return 0;
}