Hi, I'm writing a shader which is using emulated double precision (stored as vec2 with high part as vec2.x and low part as vec2.y) as described here https://www.thasler.com/blog/blog/glsl-part2-emu which also fortunately had the functions for computing addition and multiplication with emulated doubles but not for subtraction and division which I also need. In the article the DSFUN90 library is mentioned as the source for the emulated double maths so I followed that link but I couldn't find where the author got their functions from. Afterwards I tried writing my own subtraction and division functions with no success. I imagine that subtraction is relatively simple but I have no idea how division would work.

Here's the addition and multiplication functions from the article (slightly modified):

vec2 dAdd (vec2 doubleA, vec2 doubleB){

    vec2 doubleResult;
    float t1, t2, e;

    t1 = doubleA.x + doubleB.x;
    e = t1 - doubleA.x;
    t2 = ((doubleB.x - e) + (doubleA.x - (t1 - e))) + doubleA.y + doubleB.y;

    doubleResult.x = t1 + t2;
    doubleResult.y = t2 - (doubleResult.x - t1);
    return doubleResult;

}

vec2 dMul (vec2 doubleA, vec2 doubleB){

    vec2 doubleResult;
    float c11, c21, c2, e, t1, t2;
    float a1, a2, b1, b2, cona, conb, split = 8193.0; //what is the significance of the split value?

    cona = doubleA.x * split;
    conb = doubleB.x * split;
    a1 = cona - (cona - doubleA.x);
    b1 = conb - (conb - doubleB.x);
    a2 = doubleA.x - a1;
    b2 = doubleB.x - b1;

    c11 = doubleA.x * doubleB.x;
    c21 = a2 * b2 + (a2 * b1 + (a1 * b2 + (a1 * b1 - c11)));

    c2 = doubleA.x * doubleB.y + doubleA.y * doubleB.x;

    t1 = c11 + c2;
    e = t1 - c11;
    t2 = doubleA.y * doubleB.y + ((c2 - e) + (c11 - (t1 - e))) + c21;

    doubleResult.x = t1 + t2;
    doubleResult.y = t2 - (doubleResult.x - t1);

    return doubleResult;
    
}

Can anyone explain what would be necessary for the subtraction and division functions?

return dAdd(doubleA, vec2(-doubleB.x, doubleB.y);

Oh, yeah that is really simple, but why wouldn't the low part of doubleB also be negative?

Apologies for bumping this thread but I have another question, how could I go about implementing an emulated double precision exp function?

Start by implementing normal exp function without exp.