Different floating point behavior

output

type: float
iv1: -0x1p+0, -0x1.ccccccp-1
Interval!float(-1, -0.9, 2, 0, -nan, 0.118267, 2.00563, LinearFun!float(nan, nan, nan), LinearFun!float(nan, nan, nan))
iv.lx -1.0000000000 (-0x1p+0)
hat.a 0.1182672605 (0x1.e46c36p-4)
hat.slope 1.1826723814 (0x1.2ec39ep+0)
hat.y -0.8999999762 (-0x1.ccccccp-1)
hat(p1) -0.0000000075 (-0x1p-27) //-0
hat(p2) 0.1182672605 (0x1.e46c36p-4)
---
type: double
iv1: -0x1p+0, -0x1.ccccccccccccdp-1
Interval!double(-1, -0.9, 2, 0, -nan, 0.118267, 2.00562, LinearFun!double(nan, nan, nan), LinearFun!double(nan, nan, nan))
iv.lx -1.0000000000 (-0x1p+0)
hat.a 0.1182672100 (0x1.e46c28723a042p-4)
hat.slope 1.1826721000 (0x1.2ec399476442ap+0)
hat.y -0.9000000000 (-0x1.ccccccccccccdp-1)
hat(p1) 0.0000000000 (0x1p-56) // +0
hat(p2) 0.1182672100 (0x1.e46c28723a042p-4)
---
type: real
iv1: -0x8p-3, -0xe.666666666666666p-4
Interval!real(-1, -0.9, 2, 0, -nan, 0.118267, 2.00562, LinearFun!real(nan, nan, nan), LinearFun!real(nan, nan, nan))
iv.lx -1.0000000000 (-0x8p-3)
hat.a 0.1182672100 (0xf.23614391d021968p-7)
hat.slope 1.1826721000 (0x9.761cca3b2214fdfp-3)
hat.y -0.9000000000 (-0xe.666666666666666p-4)
hat(p1) -0.0000000000 (-0x8p-70) //-0
hat(p2) 0.1182672100 (0xf.23614391d021968p-7)

test.d

void transform(S)(S c, S f0, S f1, ref S t0, ref S t1)
{
    import std.math: pow, exp, copysign;
    // for c=0 no transformations are applied
    t0 = (c == 0) ? f0 : copysign(S(1), c) * exp(c * f0);
    t1 = (c == 0) ? f1 : c * t0 * f1;
}

Interval!S transformToInterval(S)(in S l, in S r, in S c,
            in S lf0, in S lf1,
            in S rf0, in S rf1)
{
    S lt0 = void, lt1 = void, rt0 = void, rt1 = void;
    transform(c, lf0, lf1, lt0, lt1);
    transform(c, rf0, rf1, rt0, rt1);
    return Interval!S(l, r, c, lt0, lt1, rt0, rt1);
}

LinearFun!S secant(S)(in S xl, in S xr, in S yl, in S yr)
{
    auto slope = (yr - yl) / (xr - xl);
    // y (aka x0) is defined to be the maximal point of the boundary
    if (yl >= yr)
        return LinearFun!S(slope, xl, yl);
    else
        return LinearFun!S(slope, xr, yr);
}

void determineSqueezeAndHat(S)(ref Interval!S iv)
{
    auto sec = secant(iv.lx, iv.rx, iv.ltx, iv.rtx);
    auto t_l = LinearFun!S(iv.lx, iv.ltx, iv.lt1x);
    auto t_r = LinearFun!S(iv.rx, iv.rtx, iv.rt1x);

    iv.squeeze = iv.ltx < iv.rtx ? t_l : t_r;
    iv.hat = sec;
}

struct LinearFun(S)
{
    /// direction and steepness
    S slope; // aka beta

    S _y;

    S _a;

    S a() @property const
    {
        return _a;
    }

    this(in S slope, in S _y, in S _a)
    {
        this.slope = slope;
        this._y = _y;
        this._a = _a;
    }

    /// call the linear function with x
    S opCall(in S x) const
    {
        return _a + slope * (x - _y);
    }
}

struct Interval(S)
{

    /// positions of the interval
    immutable S lx;
    S rx;

    /// T_c family of the interval
    immutable S c;

    /// transformed values
    immutable S ltx, lt1x;
    S rtx, rt1x;

    LinearFun!S hat;
    LinearFun!S squeeze;

    /// only supported constructor
    this (S lx, S rx, S c,
          S ltx, S lt1x,
          S rtx, S rt1x)
    {
        this.lx = lx;
        this.rx = rx;
        this.c = c;

        this.ltx  = ltx;
        this.lt1x = lt1x;

        this.rtx  = rtx;
        this.rt1x = rt1x;
    }
}

void main(string[] args)
{
    import std.stdio;
    import std.math : log;
    import std.meta : AliasSeq;
    import std.range : lockstep;
    import std.range : dropOne, save;
    foreach (S; AliasSeq!(float, double, real))
    {
        auto f0 = (S x) => log(1 - x^^4);
        auto f1 = (S x) => -4 * x^^3 / (1 - x^^4);
        auto it = (S l, S r, S c) => transformToInterval!S(l, r, c, f0(l), f1(l),
                                                                f0(r), f1(r));

        writefln("type: %s", S.stringof);
        S p1 = -1.0;
        S p2 = -0.9;

        writefln("iv1: %a, %a", p1, p2);

        S c = 2;
        auto iv = it(p1, p2, c);
        writeln(iv);

        determineSqueezeAndHat!S(iv);

        writefln("iv.lx %.10f (%1$a)", iv.lx);
        writefln("hat.a %.10f (%1$a)", iv.hat.a);
        writefln("hat.slope %.10f (%1$a)", iv.hat.slope);
        writefln("hat.y %.10f (%1$a)", iv.hat._y);
        writefln("hat(p1) %.10f (%1$a)", iv.hat(iv.lx)); // bum!!
        writefln("hat(p2) %.10f (%1$a)", iv.hat(iv.rx));
        writeln("---");
    }
}