mrtab/rbjcs.c at main · brettyhale/mrtab · GitHub

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
/******************************************************************************/

/* RBJ.2.L2 : c(s) sequence: */

/* Copyright (c) 2020 Brett Hale.
 * distributed under BSD-2-Clause license terms. see: mrtab.c */

/******************************************************************************/

#include <inttypes.h>
#include <stdio.h>

#include <float.h>
#include <math.h>

#ifndef M_PI
#define M_PI (3.14159265358979323846264338327950288)
#endif

#if defined (QUADMATH)
#pragma GCC diagnostic ignored "-Wpedantic" /* (Q-suffix) */
#include <quadmath.h>
#endif


int main (void)
{
    unsigned int smax = (30), s; /* {0 .. smax} table: */
    double c[(30) + 1], r;

    /* bias summation terms such that: fp{c(s)} >= c(s) */

    r = M_PI * M_PI / 6.0;
    c[0] = nextafter(r, DBL_MAX);

    for (s = 1; s <= smax; s++)
    {
        double x = - 1.0 / (double) (s * s);
        r += nextafter(x, DBL_MAX);

        c[s] = (s + 1) * r; /* fp{r(s)} > r(s) */
    }

#if defined (QUADMATH)

    /* evaluate with quad precision to show that:
     * 0 <= (fp{c(s)} - c(s)) / c(s) < (s + 1) * (EPS) */

    __float128 cq[(30) + 1], rq, eq;

    rq = M_PIq * M_PIq / 6.0, eq = 0.0;
    cq[0] = rq;

    for (s = 1; s <= smax; s++)
    {
        __float128 y, t; /* Kahan summation: */

        y = - 1.0 / (__float128) (s * s) - eq;
        t = rq + y; eq = (t - rq) - y; rq = t;

        cq[s] = (s + 1) * rq;
    }

    for (s = 0; s <= smax; s++)
    {
        const char *const fmt = "s = %2u : rel = %s, "
            "rel / DBL_EPS = %s\n";

        __float128 rerr;
        char rbuf[128], ebuf[128];

        rerr = (c[s] - cq[s]) / cq[s];
        quadmath_snprintf(rbuf, sizeof(rbuf), "%.16Qe", rerr);

        rerr = rerr / DBL_EPSILON;
        quadmath_snprintf(ebuf, sizeof(ebuf), "%2.2Qf", rerr);

        fprintf(stdout, fmt, s, rbuf, ebuf);
    }

#endif

    fprintf(stdout, "\n    %.16e", c[0]);

    for (s = 1; s <= smax; s++)
        fprintf(stdout, (s % 3) ? ", %.16e" : ",\n    %.16e", c[s]);

    fprintf(stdout, "\n\n");

    return (0);
}

/******************************************************************************/