#define PI 3.1415926535897932384626433

#define sq(x) dot(x, x)

vec2 sic(float t)
{
    return vec2(cos(t), sin(t));
}
vec2 perp(vec2 v)
{
    return (v * vec2(1, -1)).yx;
}

float kp = 33., 
kq = 44.,
r1 = 0.25, 
r2 = 0.125, 
r3 = 0.065;// + 0.3*(sin(iTime)+1.)/2.;//Change kp and kq!


vec3 torusKnot(float t)
{
    vec2 sicXY = sic(kp * t);
    vec2 sicRZ = r2 * sic(kq * t);
    
    return vec3((sicRZ.x + r1)*sicXY, sicRZ.y);
}
vec3 torusKnotDerivative(float t )
{
    vec2 sicXY = sic(kp * t);
    vec2 sicRZ = r2 * sic(kq * t);
    
    vec2 dSicXY = kp * vec2(-1, 1) * sicXY.yx;
    vec2 dSicRZ = kq * vec2(-1, 1) * sicRZ.yx;
    
    return vec3(dSicRZ.x*sicXY + dSicXY*(sicRZ.x + r1), dSicRZ.y);
}

float torusKnotSqDistance(float t, vec3 p )
{
    return sq(torusKnot(t) - p);
}
float torusKnotSqDistanceDerivative(float t, vec3 p )
{
    return 2.*dot(torusKnot(t) - p, torusKnotDerivative(t));
}


vec3 torusKnotSqDistanceMinimumInside(vec3 p /*, TorusKnotParameters tkp*/ )
{
    //all hard coded numbers depend on the 'r's
    
    int sections = int(2.*max(kq,kp));
    float sectionLength = 2.*PI/float(sections);
    
    float lerningRate = 1.5/(max(kq,kp));
    const int maxIterations = 50;//50
    
    float minDist;
    float bestT;
    
    for(int j = 0; j < sections; j++)
    {
        float t = sectionLength * float(j);
        
        for(int i = 0; i < maxIterations; i++)
        {
            float dt = torusKnotSqDistanceDerivative(t, p);
            dt *= lerningRate;
            
            if(abs(dt) < 0.001)
            {
                float sqDist = torusKnotSqDistance(t, p);
        
                if(sqDist <= sq(r3))
                {
                    return vec3(t, sqDist, 1.);
                }
                
                break;
            }
            
            t -= dt;
            
            
            if(t != clamp(t, sectionLength * (float(j)-1.), sectionLength * (float(j)+1.)))
            {
                break;
            }
        }
    }
    
    return vec3(0.);
}

int gcd(ivec2 v)
{
    while(v.x != v.y)
    {
        if(v.x > v.y)
            v.x -= v.y;
        else
            v.y -= v.x;
    }
    
    return v.x;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 table = vec2(5, 3);
    
    vec2 uv = table*fragCoord/iResolution.xy;
    ivec2 kpkq = ivec2(uv) + 1;
    kpkq = kpkq.yx;
    uv = mod(uv, 1.);
    uv -= 0.5;
    
    //uv.x = table.x*fragCoord.x/iResolution.x
    uv.x *= iResolution.x*table.y/(table.x*iResolution.y);
	
    const float clickSize = 0.;
    if(iMouse.z < 0. && length(iMouse.xy + iMouse.zw) <= clickSize)//ivec2(iMouse.xy*table/iResolution.xy) == ivec2(-iMouse.zw*table/iResolution.xy)
    {
		uv = (fragCoord - iResolution.xy/2.)/iResolution.y;
        kpkq = ivec2(table*iMouse.xy/iResolution.xy) + 1;
        kpkq = kpkq.yx;
    }

    vec2 angles = /*false &&*/ iMouse.z > 0. ? PI*(2.*iMouse.xy/iResolution.xy - 1.) : 1.5*iTime*vec2(0.7, 1);//PI*vec2(0.25, 0.25)*(sic(iTime) + vec2(0, 1))
    vec2 sic0 = sic(angles[0]);
    vec2 sic1 = sic(angles[1]);
    
    vec3 f = vec3(sic1.x * sic0, sic1.y);
    vec3 r = vec3(perp(sic0), 0);
    vec3 u = cross(f, r);
    
    ////////////////////////////////////The torus knot parameters//////////////////////////////////////////
    //float kp = 33., kq = 44., r1 = 0.25, r2 = 0.125, r3 = 0.065;// + 0.3*(sin(iTime)+1.)/2.;//Change kp and kq!
    // p and q are flipped. The parameters: 'lerningRate', 'maxIterations' and 'sections' also may need to be changed.
    
    //ivec2 kpkq = ivec2(10.*iMouse.xy/iResolution.xy) + 1;
    
    kpkq /= gcd(kpkq);
    
    kp = float(kpkq.x);
    kq = float(kpkq.y);
    
    
    //TorusKnotParameters tkp = TorusKnotParameters(kp, kq, r1, r2, r3);
    
    vec3 p = uv.x*r + uv.y*f;
    
    vec3 res = torusKnotSqDistanceMinimumInside(p);
    
    vec3 col = vec3(0);
    
    
    if(res[2] > 0.5)
    {
        uv.x *= table.x*iResolution.y/(iResolution.x*table.y);
        
        p = (uv.x*r + uv.y*f);//square of size 1x1x1
        
        col = 2.*p + 0.5;
        
        float len = length(col);
        col = normalize(col);
        
        //col = vec3(vec2(kpkq)/(table+1.), 0);
        
        col *= smoothstep(r3, r3 - 1./iResolution.y, sqrt(res[1]));
    }

    // Output to screen
    fragColor = vec4(col,1.0);
}