// --- Managing discontinuity when computing hardware derivatives ( -> Jacobian, fwidth )
// See also: computing Jacobian to cancel distortions https://www.shadertoy.com/view/WlByRW

void mainImage(out vec4 O,  vec2 u ){

    vec2 R = iResolution.xy,  I, z, 
         U = (u - .5*R) / R.y,                           // normalized coordinates
         D = vec2(6,-1);

    U*=2.; U-=.5; U/=dot(U,U); U+=.5;                    // distortion
    float a = atan(U.y, U.x)/6.283;   
    U = log(length(U)) +  a * D - .3*iTime;              // spiraling
                       // a cause a discontinuity. unseen in spiral because param fits.
    vec2 dFx = dFdx(U), dFy = dFdy(U); // but derivatives (fwidth,dFdx,dFdy) will get a +1 or -1 jump
    if (u.x/R.x<.5)                                      // --- left: the naive way
        z = fwidth(U);
    else {                                               // --- right: eliminating the jump trough discontinuity
        float dax = dFdx(a), day = dFdy(a),
               dx = abs(dax) > .5 ? sign(dax) : 0.,      // detect the jump
               dy = abs(day) > .5 ? sign(day) : 0.;      // ( jump is +- 1 since atan/2PI, 0.5 threshold is large )
        dFx -= D*dx, dFy -= D*dy,                        // eliminates in derivatives
        z = abs(dFx)+abs(dFy);                           // recomputes fwidth manually
/* shorter:
    dax = dFdx(a), day = dFdy(a);
    dFX = dFdx(U) - ( abs(dax) > .5 ? D*sign(dax) : R-R ),
    dFY = dFdy(U) - ( abs(day) > .5 ? D*sign(day) : R-R );  
    -> Jacobian, det, width, length...
*/     
    }
                                                         // --- display
    U.x+=.5; I = round(U); U = fract(U);                 // show distorted checker  
    O = vec4(.5*U,0,0); 
    
    if (iMouse.z <= 0.)                                  // draw antialiased circles
        O += smoothstep(3.*length(z),0., abs(length(U*2.-1.) - .5)-.01 ); 
    else {                                               // draw derivatives
#if 1                                                    //   fwidth   
        O += vec4(z,0,0);                                  
     // O += length(z);              
#else                                                    //   Jacobian   
        mat2 M = inverse(transpose(mat2(dFx, dFy)))/R.y; // Jacobian to go back to screen
        O  = vec4( .5+ M);
     // O += max( 0., -30.*determinant(M)*R.y );
#endif
    }
    if (abs(u.x-R.x/2.)<2.) O.b ++;                      // separator
}