#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <errno.h>
#include <fenv.h>

const int DoubleExpLen=11;
const int DoubleMantissaLen=52;
const int DoubleBits = DoubleExpLen + DoubleMantissaLen + 1;

void print_double_in_binary (double a)
{
        uint64_t t;
        memcpy (&t, &a, sizeof(double));

        // sign
        printf ("%ld ", (t>>(DoubleBits-1))&1); // sign

        // exponent
        for (int i=(DoubleBits-2); i!=(DoubleBits-2)-DoubleExpLen; i--)
                printf ("%ld", (t>>i)&1);
        printf (" ");

        // mantissa
        for (int i=(DoubleBits-2)-DoubleExpLen; i>=0; i--)
                printf ("%ld", (t>>i)&1);
        printf (" 0x%016lx", t);
};

// https://en.cppreference.com/w/c/numeric/fenv/feround
void show_fe_current_rounding_direction(void)
{
        printf("current rounding direction:  ");
        switch (fegetround()) {
                case FE_TONEAREST:  printf ("FE_TONEAREST");  break;
                case FE_DOWNWARD:   printf ("FE_DOWNWARD");   break;
                case FE_UPWARD:     printf ("FE_UPWARD");     break;
                case FE_TOWARDZERO: printf ("FE_TOWARDZERO"); break;
                default:            printf ("unknown");
        };
        printf("\n");
}

int main(){
        double a;
        int i;

        printf ("0.1: "); print_double_in_binary (0.1); printf ("\n");
        printf ("0.2: "); print_double_in_binary (0.2); printf ("\n");
        printf ("0.3: "); print_double_in_binary (0.3); printf ("\n");
        a = 0.2;
        printf ("a = 0.2;\n");
        printf ("%f %e ", a, a); print_double_in_binary(a); printf ("\n");
        show_fe_current_rounding_direction();
#if 0        
        fesetround(FE_DOWNWARD);        // works OK
        fesetround(FE_UPWARD);          // not zero
        fesetround(FE_TONEAREST);       // not zero
        fesetround(FE_TOWARDZERO);      // works OK
#endif
        printf ("a += 0.1;\n");
        feclearexcept(FE_ALL_EXCEPT);
        a += 0.1;
        if(fetestexcept(FE_INEXACT))
                printf ("FE_INEXACT\n");
        printf ("%f %e ", a, a); print_double_in_binary(a); printf ("\n");

        printf ("a -= 0.3;\n");
        feclearexcept(FE_ALL_EXCEPT);
        a -= 0.3;
        if(fetestexcept(FE_INEXACT))
                printf ("FE_INEXACT\n");
        printf ("%f %e ", a, a); print_double_in_binary(a); printf ("\n");

        if (a==0.0)
                printf ("a==0\n");
        else
                printf ("a!=0\n");

        for (i = 0; a < 1.0; i++) 
        {
                printf ("i=%d ", i);
                printf ("%f %e ", a, a); print_double_in_binary(a); printf ("\n");
                a += a;
        };

        printf("stop. i=%d, a=%f %e\n", i, a, a);

        return 0;
}