const char hz_gl4_curve_to_sdf_fragment_shader[] = 
{
    "#version 450 core\n"
    "#extension GL_ARB_separate_shader_objects : enable\n"
    "\n"
    "layout (location = 0) in vec2 f_pixel_coord;\n"
    "layout (location = 1) in flat vec2 f_bezier_v1;\n"
    "layout (location = 2) in flat vec2 f_bezier_v2;\n"
    "layout (location = 3) in flat vec2 f_bezier_c1;\n"
    "layout (location = 4) in flat vec2 f_bezier_c2;\n"
    "layout (location = 5) in flat int f_bezier_type;\n"
    "\n"
    "layout (std140, binding = 0) uniform UboData {\n"
    "    mat4 view_matrix;\n"
    "    float max_sdf_distance;\n"
    "};\n"
    "\n"
    "float dot2( in vec2 v ) { return dot(v,v); }\n"
    "#define saturate(x) clamp(x,0.,1.)\n"
    "\n"
    "float cross2(in vec2 a, in vec2 b) {\n"
    "    return a.x*b.y - a.y*b.x;\n"
    "}\n"
    "\n"
    "vec2 B1(in vec2 P0, in vec2 P1, in float t) {\n"
    "    return P0 + t * (P1 - P0);\n"
    "}\n"
    "\n"
    "float ud_line(vec2 P0, vec2 P1, vec2 P)\n"
    "{\n"
    "    // P0 + t * (P1 - P0)\n"
    "    vec2 dv = P1 - P0;\n"
    "    float t = saturate(dot(P-P0,dv)/dot2(dv));\n"
    "    return distance(P,(P0+t*dv));\n"
    "}\n"
    "\n"
    "float sd_line(in vec2 P0, in vec2 P1, in vec2 P) {\n"
    "    // P0 + t * (P1 - P0)\n"
    "    vec2 dv = P1 - P0;\n"
    "    float t = saturate(dot(P-P0,dv)/dot2(dv));\n"
    "    vec2 z = (P0+t*dv) - P;\n"
    "    return sign(cross2(dv,z)) * sqrt(dot2(z));\n"
    "}\n"
    "\n"
    "vec2 B2(in vec2 P0, in vec2 p1, in vec2 p2, in float t) {\n"
    "    return t*t*p2 + 2.0*t*p1 + P0;\n"
    "}\n"
    "\n"
    "vec2 B3(in vec2 P0, in vec2 p1, in vec2 p2, in vec2 p3, in float t) {\n"
    "    return t*t*t*p3 + 3.0*t*t*p2 + 3.0*t*p1 + P0;\n"
    "}\n"
    "\n"
    "// Taken from scholarius at https://www.shadertoy.com/view/wts3RX\n"
    "#define NEWTON_ITER 2\n"
    "#define HALLEY_ITER 0\n"
    "\n"
    "float cbrt( float x ) {\n"
    "    if (sign(x) == 0.0) return 0.0;\n"
    "    \n"
    "    float y = sign(x) * uintBitsToFloat( floatBitsToUint( abs(x) ) / 3u + 0x2a514067u );\n"
    "        \n"
    "    \n"
    "    for( int i = 0; i < NEWTON_ITER; ++i )\n"
    "        y = ( 2. * y + x / ( y * y ) ) * .333333333;\n"
    "\n"
    "    for( int i = 0; i < HALLEY_ITER; ++i )\n"
    "    {\n"
    "        float y3 = y * y * y;\n"
    "        y *= ( y3 + 2. * x ) / ( 2. * y3 + x );\n"
    "    }\n"
    "    \n"
    "    return y;\n"
    "}\n"
    "\n"
    "\n"
    "float dist_cubic_inflexion_to_extrema(float p, float q) {\n"
    "    return pow(q/2.0,2.0) + pow(p/3.0,3.0);\n"
    "}\n"
    "\n"
    "#define FLT_EPSILON 1.e-7\n"
    "#define PI 3.14159\n"
    "\n"
    "// Complex cube root\n"
    "vec3 complex_cbrt(vec2 c) {\n"
    "    float m = pow(dot(c,c),1.0/6.0); // magnitude\n"
    "    float t = atan(c.y/c.x);// theta original angle\n"
    "    float t0 = (t+2.0*PI*0.)/3.0;\n"
    "    float t1 = (t+2.0*PI*1.0)/3.0;\n"
    "    float t2 = (t+2.0*PI*2.0)/3.0;\n"
    "    vec2 z0 = vec2(cos(t0),sin(t0))*m;\n"
    "    vec2 z1 = vec2(cos(t1),sin(t1))*m;\n"
    "    vec2 z2 = vec2(cos(t2),sin(t2))*m;\n"
    "    return vec3(z0[0], z1[0], z2[0]);\n"
    "}\n"
    "\n"
    "vec2 complex_conjugate(vec2 z) {\n"
    "    return vec2(z[0], -z[1]);\n"
    "}\n"
    "\n"
    "// float ud_bezier2(vec2 P, vec2 P0, vec2 P1, vec2 P2)\n"
    "// {\n"
    "//     float res = 0.0;\n"
    "    \n"
    "//     vec2 p0 = P - P0;\n"
    "//     vec2 p1 = P1 - P0;\n"
    "//     vec2 p2 = P2 - 2.0*P1 + P0;\n"
    "\n"
    "//     // declare coefficients for equation of form t^3 + at^2 + bt + c = 0\n"
    "//     float a = dot2(p2);\n"
    "//     float b = 3.0*dot(p1,p2);\n"
    "//     float c = dot(2.0*p1,p1)-dot(p2,p0);\n"
    "//     float d = -dot(p1,p0);\n"
    "\n"
    "//     // declare y, p and q for depressed cubic form y^3 + py + q = 0 where:\n"
    "//     // t = y - b/(3a)\n"
    "//     // p = c/a - (b^2)/(3a^2) \n"
    "//     // q = (2b^3)/(27a^3) - (bc)/(3a^2) + d/a\n"
    "//     float p = c/a - (b*b)/(3.0*a*a);\n"
    "//     float q = (2.0*b*b*b)/(27.0*a*a*a) - (b*c)/(3.0*a*a) + d/a;\n"
    "\n"
    "//     // compute the difference between maxima and inflection point, \n"
    "//     // this tells us how many roots the cubic has.\n"
    "//     float z = pow(q/2.0,2.0) + pow(p/3.0,3.0);\n"
    "//     float h = b/(3.0*a); // x-delta of inflextion point\n"
    "\n"
    "//     if (z >= 0.0) {\n"
    "//         // 1 root\n"
    "//         z = sqrt(z);\n"
    "//         vec2 x = vec2(z,-z)-q/2.0;\n"
    "//         vec2 uv = sign(x)*pow(abs(x),vec2(1.0/3.0));\n"
    "//         float t = clamp(uv.x+uv.y-h,0.0,1.0);\n"
    "//         vec2 v = B2(P0,p1,p2,t) - P;\n"
    "//         res = dot2(v);\n"
    "//     } else {\n"
    "//         // 3 roots\n"
    "//         // Formula based on the paper on Joukowsky maps and their relation\n"
    "//         // to the cubic equation.\n"
    "//         // https://maa.tandfonline.com/doi/full/10.1080/00029890.2019.1528814#.YRBwTi295B1\n"
    "//         // This solves the three-root case using an equilateral triangle's vertices\n"
    "//         // projected onto the x-axis.\n"
    "//         // h is the point of inflexion on the x-axis.\n"
    "//         // h = b/(3a)\n"
    "//         // tx = h+2rcos(theta)\n"
    "//         // ty = h+2rcos(theta+(2pi)/3)\n"
    "//         // only the first two of these roots are required.\n"
    "//         float r = sqrt(p/-3.0), k = acos(q/(-2.0*r*r*r))/3.0;\n"
    "//         float n = cos(k);\n"
    "//         float m = cos(k+2.09439510239);\n"
    "//         vec2 t = clamp(vec2(n,m)*2.0*r-h,0.0,1.0);\n"
    "//         // t = clamp(t+vec2(h,-h), 0.0, 1.0);\n"
    "//         vec2 vx = B2(P0,p1,p2,t.x) - P, vy = B2(P0,p1,p2,t.y) - P;\n"
    "//         float xd = dot2(vx), yd = dot2(vy);\n"
    "//         if (xd<yd) res=xd; else res=yd;\n"
    "//     }\n"
    "    \n"
    "//     return sqrt(res);\n"
    "// }\n"
    "\n"
    "float udBezier( in vec2 pos, in vec2 A, in vec2 B, in vec2 C )\n"
    "{    \n"
    "    vec2 a = B - A;\n"
    "    vec2 b = A - 2.0*B + C;\n"
    "    vec2 c = a * 2.0;\n"
    "    vec2 d = A - pos;\n"
    "\n"
    "    float kk = 1.0/dot(b,b);\n"
    "    float kx = kk * dot(a,b);\n"
    "    float ky = kk * (2.0*dot(a,a)+dot(d,b))/3.0;\n"
    "    float kz = kk * dot(d,a);      \n"
    "\n"
    "    float res = 0.0;\n"
    "\n"
    "    float p  = ky - kx*kx;\n"
    "    float q  = kx*(2.0*kx*kx - 3.0*ky) + kz;\n"
    "    float p3 = p*p*p;\n"
    "    float q2 = q*q;\n"
    "    float h  = q2 + 4.0*p3;\n"
    "\n"
    "    if( h>=0.0 ) \n"
    "    {   // 1 root\n"
    "        h = sqrt(h);\n"
    "        vec2 x = (vec2(h,-h)-q)/2.0;\n"
    "\n"
    "        #if 0\n"
    "        // When p≈0 and p<0, h-q has catastrophic cancelation. So, we do\n"
    "        // h=√(q²+4p³)=q·√(1+4p³/q²)=q·√(1+w) instead. Now we approximate\n"
    "        // √ by a linear Taylor expansion into h≈q(1+½w) so that the q's\n"
    "        // cancel each other in h-q. Expanding and simplifying further we\n"
    "        // get x=vec2(p³/q,-p³/q-q). And using a second degree Taylor\n"
    "        // expansion instead: x=vec2(k,-k-q) with k=(1-p³/q²)·p³/q\n"
    "        if( abs(p)<0.001 )\n"
    "        {\n"
    "            float k = p3/q;              // linear approx\n"
    "          //float k = (1.0-p3/q2)*p3/q;  // quadratic approx \n"
    "            x = vec2(k,-k-q);  \n"
    "        }\n"
    "        #endif\n"
    "\n"
    "        vec2 uv = sign(x)*pow(abs(x), vec2(1.0/3.0));\n"
    "        float t = clamp( uv.x+uv.y-kx, 0.0, 1.0 );\n"
    "        vec2  q = d+(c+b*t)*t;\n"
    "        res = dot2(q);\n"
    "    }\n"
    "    else \n"
    "    {   // 3 roots\n"
    "        float z = sqrt(-p);\n"
    "        float v = acos(q/(p*z*2.0))/3.0;\n"
    "        float m = cos(v);\n"
    "        float n = sin(v)*1.732050808;\n"
    "        vec3  t = clamp( vec3(m+m,-n-m,n-m)*z-kx, 0.0, 1.0 );\n"
    "        vec2  qx=d+(c+b*t.x)*t.x; float dx=dot2(qx);\n"
    "        vec2  qy=d+(c+b*t.y)*t.y; float dy=dot2(qy);\n"
    "        if( dx<dy ) res=dx; else res=dy;\n"
    "    }\n"
    "    \n"
    "    return sqrt( res );\n"
    "}\n"
    "\n"
    "layout (location = 0) out vec4 outColor;\n"
    "\n"
    "void main() {\n"
    "    float t = max_sdf_distance;\n"
    "\n"
    "    switch (f_bezier_type) {\n"
    "        case 2: { //line\n"
    "            t = ud_line(f_bezier_v1, f_bezier_v2, f_pixel_coord);\n"
    "            break;\n"
    "        }\n"
    "\n"
    "        case 3: { //quadratic bezier\n"
    "            t = udBezier(f_pixel_coord,f_bezier_v1, f_bezier_c1, f_bezier_v2);\n"
    "            break;\n"
    "        }\n"
    "\n"
    "        case 4: { // cubic bezier\n"
    "            break;\n"
    "        }\n"
    "    }\n"
    "\n"
    "    if (t > max_sdf_distance) // discard fragments outside of curve max distance\n"
    "        discard;\n"
    "\n"
    "    t = 0.5 - (t / max_sdf_distance / 2.0);\n"
    "    outColor = vec4(t,t,t,t);\n"
    "}\n"
};
