1/* gpt2.rs's constants_are_their_values and scalar_functions_are_close (references in double, computed
2 * here without libm), then a golden digest from the Rust over 300 000 inputs of each function. */
3#include "mesh/ref/gpt2_scalar.h"
4
5#include "mesh/ftz.h"
6
7#include <stdio.h>
8
9static int failures;
10
11static double dabs(double x) { return x < 0 ? -x : x; }
12
13/* exp by range reduction to |r| <= ln2/2 and a long Taylor series. */
14static double dexp(double x) {
15    const double ln2 = 0.6931471805599453;
16    int k = (int)(x / ln2 + (x < 0 ? -0.5 : 0.5));
17    double r = x - k * ln2, term = 1, sum = 1;
18    for (int n = 1; n < 30; n++) term *= r / n, sum += term;
19    while (k > 0) sum *= 2, k--;
20    while (k < 0) sum /= 2, k++;
21    return sum;
22}
23
24static double dtanh(double x) { return 1 - 2 / (dexp(2 * x) + 1); }
25
26static double dsqrt(double x) {
27    double y = x > 1 ? x : 1;
28    for (int i = 0; i < 200; i++) y = 0.5 * (y + x / y);
29    return y;
30}
31
32static uint64_t fh(uint64_t h, uint64_t w) { return (h ^ w) * 0x100000001b3ull; }
33
34static uint32_t mix(uint32_t x) {
35    x ^= x >> 16;
36    x *= 0x7feb352du;
37    x ^= x >> 15;
38    x *= 0x846ca68bu;
39    return x ^ (x >> 16);
40}
41
42static uint32_t gen(uint32_t m, uint32_t i) {
43    uint32_t h = mix((i * 0x9e3779b9u) ^ (m * 0x85ebca6bu));
44    uint32_t h2 = mix(h ^ 0x27d4eb2fu);
45    if (m == 0) return h;
46    if (m == 1) return (h & 0x80000000u) | ((h2 % 41) << 23) | (h & 0x007fffffu);
47    return (h & 0x80000000u) | ((50 + h2 % 31) << 23) | (h & 0x007fffffu);
48}
49
50int main(void) {
51    /* constants_are_their_values */
52    if (mesh_ftz_from_bits(MESH_GPT2_LN_EPS_BITS) != 1e-5f || mesh_ftz_from_bits(MESH_GPT2_GELU_C_BITS) != 0.7978846f ||
53        mesh_ftz_from_bits(MESH_GPT2_GELU_A_BITS) != 0.044715f || mesh_ftz_from_bits(MESH_GPT2_LOG2E_BITS) != 1.442695f ||
54        mesh_ftz_from_bits(MESH_GPT2_LN2_HI_BITS) != 0.69314575f ||
55        mesh_ftz_from_bits(MESH_GPT2_LN2_LO_BITS) != 1.4286068e-6f || mesh_ftz_from_bits(MESH_GPT2_ATTN_SCALE_BITS) != 0.125f)
56        failures++;
57    const float want[11] = {1.0f / 3628800.0f, 1.0f / 362880.0f, 1.0f / 40320.0f, 1.0f / 5040.0f, 1.0f / 720.0f,
58                            1.0f / 120.0f,     1.0f / 24.0f,     1.0f / 6.0f,     0.5f,           1.0f, 1.0f};
59    for (int i = 0; i < 11; i++)
60        if (mesh_ftz_from_bits(mesh_gpt2_exp_poly_coeffs[i]) != want[i]) failures++;
61
62    /* scalar_functions_are_close */
63    for (int i = -400; i <= 400; i++) {
64        float x = (float)i / 8.0f;
65        if (x > -87.0f && x < 60.0f) {
66            double e = mesh_gpt2_exp_poly(x);
67            if (dabs(e / dexp(x) - 1.0) >= 4e-7) fprintf(stderr, "exp %g: %g\n", x, e), failures++;
68        }
69        double t = mesh_gpt2_tanh_poly(x / 8.0f);
70        if (dabs(t - dtanh((double)x / 8.0)) >= 3e-7) fprintf(stderr, "tanh %g\n", x), failures++;
71    }
72    const float as[7] = {1e-5f, 0.01f, 0.5f, 1.0f, 2.0f, 1234.5f, 3e7f};
73    for (int i = 0; i < 7; i++) {
74        double s = mesh_gpt2_sqrt_newton(as[i]);
75        if (dabs(s / dsqrt(as[i]) - 1.0) >= 2e-7) fprintf(stderr, "sqrt %g\n", as[i]), failures++;
76    }
77    if (mesh_gpt2_exp_poly(0.0f) != 1.0f || mesh_gpt2_exp_poly(-100.0f) != 0.0f) failures++;
78
79    /* golden: gold gpt2 "scalars" */
80    uint64_t h = 0xcbf29ce484222325ull;
81    for (uint32_t i = 0; i < 300000; i++) {
82        float x = mesh_ftz_from_bits(gen(i % 3, i));
83        float y = ((float)i - 150000.0f) / 1000.0f;
84        float v[7] = {mesh_gpt2_exp_poly(x),  mesh_gpt2_exp_poly(y), mesh_gpt2_tanh_poly(x), mesh_gpt2_tanh_poly(y),
85                      mesh_gpt2_gelu_tanh(x), mesh_gpt2_gelu_tanh(y),
86                      mesh_gpt2_sqrt_newton(mesh_ftz_from_bits(gen(i % 3, i) & 0x7fffffffu))};
87        for (int j = 0; j < 7; j++) h = fh(h, mesh_ftz_bits(v[j]));
88    }
89    if (h != 0xcd832a1d2b3bc7b0ull) fprintf(stderr, "scalars %#018llx\n", (unsigned long long)h), failures++;
90    return failures != 0;
91}