#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;
}

struct TorusKnotParameters
{
    float kp;
    float kq;
    float r1;
    float r2;
    float r3;
};

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

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

struct Ray
{
  vec3 ro;
  vec3 rd;
};

vec3 torusKnotSqDistanceMinimumInside(vec3 p, TorusKnotParameters tkp)
{
    int sections = 1*int(max(tkp.kq,tkp.kp));
    float sectionLength = 2.*PI/float(sections);
    
    float lerningRate = 1./(max(tkp.kq,tkp.kp)-0.75);
    const int maxIterations = 900;//50
    
    float minDist;
    float bestT;
    
    for(int j = 0; j < sections; j++)
    {
        float t = sectionLength * float(j);
        bool failed = false;
        
        for(int i = 0; i < maxIterations; i++)
        {
            float dt = torusKnotSqDistanceDerivative(t, p, tkp);
            
            if(abs(lerningRate*dt) < 0.003)break;
            
            t -= lerningRate*dt;
            
            
            if(t != clamp(t, sectionLength * (float(j)-1.), sectionLength * (float(j)+1.)))
            {
                failed = true;
                break;
            }
        }
        if(failed)continue;
        
        float sqDist = torusKnotSqDistance(t, p, tkp);
        
        if(sqDist <= sq(tkp.r3))
        {
            return vec3(t, sqDist, 1.);
        }
    }
    
    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 uv = 5.*fragCoord/iResolution.xy;
    ivec2 kpkq = ivec2(uv) + 1;
    uv = mod(uv, 1.);
    uv -= 0.5;
    uv.x *= iResolution.x/iResolution.y;
	*/
    
    
    vec2 uv = (fragCoord - iResolution.xy/2.)/iResolution.y;

    vec2 angles = /*false &&*/ iMouse.z > 0.5 ? PI*(2.*iMouse.xy/iResolution.xy - 1.) : 1.5*iTime*vec2(0.3, 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 = 3., kq = 4., r1 = 0.25, r2 = 0.125, r3 = 0.05;// + 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, tkp);
    
    vec3 col = vec3(0);
    
    
    if(res[2] > 0.5)
    {
        uv = fragCoord/iResolution.xy - 0.5;//[-0.5, 0.5] 
        p = (uv.x*r + uv.y*f);//square of size 1x1x1
        
        col = 2.*p + 0.5;
        
        float len = length(col);
        col = normalize(col);
        
        col *= smoothstep(tkp.r3, tkp.r3 - 1./iResolution.y, sqrt(res[1]));
    }

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