Lill's method

// Mathologer's video: https://www.youtube.com/watch?v=IUC-8P0zXe8

​

​

//////// Examples ////////

​

// 1) x - 3

// let coefs = [1, -3]; // highest order coefficient comes first

​

// 2) (x-1)(x+2)(x-3) = x^3 - 2x^2 - 5x + 6

let coefs = [1, -2, -5, 6];

​

// 3) x(x-1) = x^2 - x

// let coefs = [1, -1, 0];

​

​

let unit = 30; // 1.0 = unit*px

function setup() {

  createCanvas(windowWidth, windowHeight);

  background(0);

  stroke(255);

  noLoop();

}

​

​

let state = 0;    // 0 - instructions, 1 - drawing

let sol_angle = 0;

​

function draw() {

    switch(state)

    {

      case 0:

        textSize(16);

        fill(255);

        strokeWeight(0);

        text("Lill's method for finding roots of polynomials.\n\n\

Polynomial ax^n + bx^(n-1) +... is represented by line segments at right angles with lengths a, b ...\n\n\

Use your mouse to adjust the incidence angle of the laser.\n\n\

If the laser hits green circle, negative tangent of angle solves the polynomial.\n\n\

Use mouse wheel to zoom in/out.\n\n\

Change polynomial coefficients at the top of source code.", 100, 100);

        fill(0, 200, 200);

        textSize(32);

​

        text("Press <ENTER>", 100, 400);

        break;

      case 1:

        translate(width/2, height/2);

        scale(1, -1);

​

        ////////////////// Draw the polynomial segments ///////////////////

        let x = 0;

        let y = 0;

        let xnew = 0;

        let ynew = 0;

        let phi = 0;

​

        let points = [[0, 0]];  // store corner points

        for (let i=0; i<coefs.length; i++)

        {

          ax = cos(phi * PI / 180);

          ay = sin(phi * PI / 180);

​

          strokeWeight(3);

          stroke(255, 0, 102);

​

          if (i%2 == 0) // horizontal segment

          {

            xnew = x + ax*coefs[i]*unit;

            line(x, y, xnew, ynew);

          }

          if (i%2 != 0) // vertical segment

          {

            ynew = y + ay*coefs[i]*unit;

            line(x, y, xnew, ynew);

          }

          stroke(255);

​

          x = xnew;

          y = ynew;

          phi += 90;

          points.push([x, y]);

        }

​

        strokeWeight(0);

        fill(0, 255, 0);

        circle(x, y, 10);

        strokeWeight(3);

        ///////////////////////////////////////////////////////////

​

        let angle = sol_angle * PI / 180;

        let direction = createVector(1, tan(angle));

​

        let xc = 0; // laser coordinates

        let yc = 0;

​

        let x0, y0; // coordinates where the laser hits the next segment

​

        for (let i=1; i<points.length-1; i++)

        {

          // extrapolate laser path and calculate where it hits the next segment

          if (i%2 != 0)

          {

            x0 = points[i][0];

            y0 = tan(angle) * (x0 - xc) + yc;

          }

​

          else

          {

            y0 = points[i][1];

            x0 = (y0 - yc)/tan(angle) + xc;

          }

​

          // bookkeeping to figure out whether to _refract_ left or right when hitting imaginary line

          let left = true;

          let a_deg = angle * PI / 180;

​

          if (i%2 != 0) // vertical segment

          {

            let dif = x0 - xc;

​

            if (dif > 0)

            {

              if (0 < a_deg && a_deg < 90){left = false;}

            }

            if (dif < 0)

            {

            if (180 < a_deg && a_deg < 270){left = false;}

            }

​

            // laser lands on the segment, _reflect_ instead of _refract_

            if ((points[i][0] < x0 && x0 < points[i+1][0]) || (points[i+1][0] < x0 && x0 < points[i][0]))

              left = 1-left;

          }

​

          else // horizontal segment

          {

            let dif = y0 - yc;

            if (dif > 0)

            {

              if (90 < a_deg && a_deg < 180){left = false;}

            }

            if (dif < 0)

            {

              if (270 < a_deg && a_deg < 360){left = false;}

            }

​

            if ((points[i][1] < y0 && y0 < points[i+1][1]) || (points[i+1][1] < y0 && y0 < points[i][1]))

              left = 1-left;

          }

​

          if (left == true)

            direction.rotate(PI/2);

          else

            direction.rotate(-PI/2);

​

          stroke(0, 255, 255);

          line(xc, yc, x0, y0);

          stroke(255);

​

          xc = x0;

          yc = y0;

          angle = direction.heading();

        }

​

        let px_last = points[points.length-1][0];

        let py_last = points[points.length-1][1];

​

        let error = sqrt((px_last-xc)*(px_last-xc) + (py_last-yc)*(py_last-yc)) / unit;

        error = round(error*100)/100;

​

        let solution = -round(tan(PI/180*sol_angle) * 100)/100;

​

        scale(1, -1);

        textSize(height/30);

        fill(255);

        strokeWeight(0);

​

        text("x = -tan(angle) = " + solution, -width/2+width/20, -height/2+height/20);

        text("error = " + error, -width/2+width/3, -height/2+height/20);

        text("Polynomial coefficients: [" + coefs + "]", 0, -height/2+height/20);

        break;

    }

}

​

​

function keyPressed() {

  if (keyCode == ENTER)

  {

    state = 1 - state;

    background(0);

    redraw();

  }

}

​

​

function mouseWheel(delta) {

  unit += 5*event.delta/53;

  if (unit <= 0)

    unit = 0;

  background(0);

  redraw();

}

​

​

function mouseMoved() {

  let mx = mouseX - width/2;

  let my = height/2 - mouseY;

  sol_angle = atan(my/mx)*180/PI;

​

  background(0);

  redraw();

}