diff --git a/floating_point/Makefile b/floating_point/Makefile index 6716405..97f0f61 100644 --- a/floating_point/Makefile +++ b/floating_point/Makefile @@ -1,6 +1,6 @@ CC := gcc CFLAGS := -Wall -Wextra -lm -O2 -std=c11 -LDLIBS := -lm +LDLIBS := -lmpfr -lgmp -lm TARGET := bf16_test SRCS := bf16.c bf16_test.c diff --git a/floating_point/README.md b/floating_point/README.md new file mode 100644 index 0000000..d67887f --- /dev/null +++ b/floating_point/README.md @@ -0,0 +1,29 @@ +# BF16 +## Test Method +### Dependencies +Libraries: +* MPFR (GNU Multiple Precision Floating-Point Reliable Library) +* GMP (GNU Multiple Precision Arithmetic Library) + +```bash +sudo apt install build-essential libmpfr-dev libgmp-dev +``` +### Running Tests +``` +make test +``` + +Produces three log files: + +``` +.range1.log -> [0.140625, pi/2] From 0x3E10 to 0x3FC9 +.range2.log -> [pi/2, +largest] From 0x3FC9 to 0x7F7F +.range3.log -> [-smallest, -largest] From 0x8080 to 0xFF7F +``` + +where `+largest/-largest` is the largest positive/negative normal BF16, that is `+largest=0x7F7F` and `-largest=0xFF7F`. + +Each line in log files contains: +``` +,,, +``` \ No newline at end of file diff --git a/floating_point/bf16_test.c b/floating_point/bf16_test.c index 6dca346..bb11512 100644 --- a/floating_point/bf16_test.c +++ b/floating_point/bf16_test.c @@ -1,67 +1,159 @@ +/* bf16_test.c + * + * MPFR-based validation harness for bf16_sin. + * Uses mpfr_sin (200-bit) as oracle, then rounds MPFR -> float -> BF16 + * and compares to your bf16_sin result. + * + * Requirements: + * - libmpfr, libgmp available at link time + * + * Header of logs: + * absolute_error, input_bits (hex), ref_bits, approx_bits + */ + #include #include -#include "bf16.h" #include +#include +#include -int main() +#include "bf16.h" + +int main(void) { - FILE *f1 = fopen(".range1.log", "w"); - FILE *f2 = fopen(".range2.log", "w"); - FILE *f3 = fopen(".range3.log", "w"); - - // x between [pi/2, 0] - bf16_t i = (bf16_t){.bits = 0b0011111111001001}; // Initial - int j = 0, record_k, error_count = 0; - float difference, max_difference = 0; - for (j; j < 450; j++){ - float glibc_sin = bf16_to_fp32(fp32_to_bf16(sinf(bf16_to_fp32(i)))); - float payne = bf16_to_fp32(bf16_sin(i, &record_k)); - difference = fabs(payne-glibc_sin); - if (difference > 0) { - if (difference > max_difference) max_difference = difference; - error_count++; - } - fprintf(f1, "%f\n",fabs(payne-glibc_sin)); - i.bits--; - } - printf("Test range [0, pi/2]:\nMax difference = %f, Total %d numbers difference\n\n", max_difference, error_count); - - // x larger than pi/2, need range reduction - i = (bf16_t){.bits = 0b0011111111001010}; // Initial - j = 0; - for (j; j < 16310; j++){ - float glibc_sin = bf16_to_fp32(fp32_to_bf16(sinf(bf16_to_fp32(i)))); - float payne = bf16_to_fp32(bf16_sin(i, &record_k)); - difference = fabs(payne-glibc_sin); - if (difference > 0) { - if (difference > max_difference) max_difference = difference; - error_count++; - } - fprintf(f2, "%f\n",fabs(payne-glibc_sin)); - i.bits++; - } - printf("Test range [pi/2, inf):\nMax difference = %f, Total %d numbers difference\n\n", max_difference, error_count); - - // x is negative, need range reduction - i = (bf16_t){.bits = 0x8080}; // Initial value, smallest negative normal BF16 - j = 0; - for (j; j < 32512; j++){ - float glibc_sin = bf16_to_fp32(fp32_to_bf16(sinf(bf16_to_fp32(i)))); - float payne = bf16_to_fp32(bf16_sin(i, &record_k)); - difference = fabs(payne-glibc_sin); - if (difference > 0) { - if (difference > max_difference) max_difference = difference; - error_count++; - } - fprintf(f3, "%f\n",fabs(payne-glibc_sin)); - i.bits++; - } - printf("Test range [-smallest, -inf):\nMax difference = %f," - "Total %d numbers difference\n\n", max_difference, error_count); - - fclose(f1); - fclose(f2); - fclose(f3); - return 0; -} + mpfr_t x_mp, y_mp; + mpfr_init2(x_mp, 200); /* 200-bit precision */ + mpfr_init2(y_mp, 200); + + FILE *f1 = fopen(".range1.log", "w"); + FILE *f2 = fopen(".range2.log", "w"); + FILE *f3 = fopen(".range3.log", "w"); + + if (!f1 || !f2 || !f3) { + fprintf(stderr, "Failed to open log files\n"); + return 1; + } + + bf16_t i; + int record_k; + + double global_max_abs = 0.0; + unsigned long global_diff_count = 0; + + /* ---------------------------------------------------------- + * Range 1: [0.140625, pi/2] + * From 0x3E10 (0.140625) to 0x3FC9 (1.5703125, rounded pi/2) + * ----------------------------------------------------------*/ + double max_abs = 0.0; + unsigned long diff_count = 0; + + for (i.bits = 0x3E10; i.bits < 0x3FC9; i.bits++) { + float x = bf16_to_fp32(i); + + mpfr_set_flt(x_mp, x, MPFR_RNDN); + mpfr_sin(y_mp, x_mp, MPFR_RNDN); + float mpfr_sinf = mpfr_get_flt(y_mp, MPFR_RNDN); + + /* Round MPFR result to BF16 properly by converting MPFR->float then float->bf16 */ + bf16_t ref_bf = fp32_to_bf16(mpfr_sinf); + float ref_rounded_f = bf16_to_fp32(ref_bf); + + bf16_t approx_bf = bf16_sin(i, &record_k); + float approx_f = bf16_to_fp32(approx_bf); + + float abs_err = fabsf(approx_f - ref_rounded_f); + + if (abs_err > max_abs) max_abs = abs_err; + if (abs_err > 0.0f) diff_count++; + + fprintf(f1, "%f, 0x%04x, 0x%04x, 0x%04x\n", + abs_err, + (unsigned) i.bits, (unsigned) ref_bf.bits, (unsigned) approx_bf.bits); + } + + printf("Test range [0.140625, pi/2]:\nMax difference = %.9g, Total %lu numbers difference.\n\n", + max_abs, diff_count); + if (max_abs > global_max_abs) global_max_abs = max_abs; + global_diff_count += diff_count; + + /* ------------------------------------------------------------- + * Range 2: [pi/2, +largest] + * From 0x3FC9 (rounded pi/2) to 0x7F7F(Largest positive normal) + * Do range reduction from [pi/2, +largeset] to [0, pi/2] + * -------------------------------------------------------------*/ + max_abs = 0.0; + diff_count = 0; + + for (i.bits = 0x3FC9; i.bits < 0x7F7F; i.bits++) { + float x = bf16_to_fp32(i); + mpfr_set_flt(x_mp, x, MPFR_RNDN); + mpfr_sin(y_mp, x_mp, MPFR_RNDN); + float mpfr_sinf = mpfr_get_flt(y_mp, MPFR_RNDN); + + bf16_t ref_bf = fp32_to_bf16(mpfr_sinf); + float ref_rounded_f = bf16_to_fp32(ref_bf); + + bf16_t approx_bf = bf16_sin(i, &record_k); + float approx_f = bf16_to_fp32(approx_bf); + + float abs_err = fabsf(approx_f - ref_rounded_f); + + if (abs_err > max_abs) max_abs = abs_err; + if (abs_err > 0.0f) diff_count++; + + fprintf(f2, "%f, 0x%04x, 0x%04x, 0x%04x\n", + abs_err, + (unsigned) i.bits, (unsigned) ref_bf.bits, (unsigned) approx_bf.bits); + } + + printf("Test range [pi/2, +largest]:\nMax difference = %.9g, Total %lu numbers difference.\n\n", + max_abs, diff_count); + if (max_abs > global_max_abs) global_max_abs = max_abs; + global_diff_count += diff_count; + + /* -------------------------------------------------------------------------- + * Range 3: [-smallest, -largest] + * From 0x8080 (smallest negative normal) to 0xFF7F (largest negative normal) + * Do range reduction from [-smallest, -largest] to [0, pi/2] + * --------------------------------------------------------------------------*/ + max_abs = 0.0; + diff_count = 0; + + for (i.bits = 0x8080; i.bits < 0xFF7F; i.bits++) { + float x = bf16_to_fp32(i); + + mpfr_set_flt(x_mp, x, MPFR_RNDN); + mpfr_sin(y_mp, x_mp, MPFR_RNDN); + float mpfr_sinf = mpfr_get_flt(y_mp, MPFR_RNDN); + + bf16_t ref_bf = fp32_to_bf16(mpfr_sinf); + float ref_rounded_f = bf16_to_fp32(ref_bf); + + bf16_t approx_bf = bf16_sin(i, &record_k); + float approx_f = bf16_to_fp32(approx_bf); + + float abs_err = fabsf(approx_f - ref_rounded_f); + + if (abs_err > max_abs) max_abs = abs_err; + if (abs_err > 0.0f) diff_count++; + + fprintf(f3, "%f, 0x%04x, 0x%04x, 0x%04x\n", + abs_err, + (unsigned) i.bits, (unsigned) ref_bf.bits, (unsigned) approx_bf.bits); + } + + printf("Test range [-smallest, -largest]:\nMax difference = %.9g, Total %lu numbers difference.\n", + max_abs, diff_count); + if (max_abs > global_max_abs) global_max_abs = max_abs; + global_diff_count += diff_count; + + fclose(f1); + fclose(f2); + fclose(f3); + + mpfr_clear(x_mp); + mpfr_clear(y_mp); + + return 0; +}