const seed = 'The random seed'; //type=string
const hatching = 1; //min=0 max=2 step=1 (For screen, For plotter, None)
const centerPiece = 1; //min=0 max=2 step=1 (Ignore, Self (centered), Self (centered like circles in segments))
const shadesOfGrey = 6; //min=1 max=16 step=1 Fewer shades equals faster drawing
const radiusStep = 15; //min=5 max=40 step=1
const segmentsPerStep = 12; //min=8 max=20 step=1
const giveOrTake = 1; //min=0 max=4 step=1
const intraRadiusDelta = 2; //min=.5 max=5 step=.1
const intraMidOffsetP = .6; //min=0 max=4 step=.1
const border = 0; //min=0 max=50 step=1
const borderRadius = 0; //min=0 max=25 step=1
const penThickness = 1; //min=.1 max=10 step=.1

// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(hatching == 0? 1/shadesOfGrey: 1);

// Global code will be evaluated once.
init();
R.seed(seed);
const bales = hatching == 0? Bales(shadesOfGrey): [new Turtle(), new Turtle()];

const hatchings = (()=> {
    if(shadesOfGrey == 1) return [false];
    if(hatching == 0) return [[1,.15]];
    if(hatching == 2) return [false];
    return Array.from({length: shadesOfGrey}, (e,i) => [R.getAngle(), .15*(2*i+3)*penThickness]);
})();

const segments = [];
for(let r = centerPiece==0?0:.5*radiusStep, max = (100-border-(Math.SQRT2-1)*borderRadius)*Math.SQRT2 + radiusStep; r < max; r+=radiusStep) {
    const sc = R.getInt(segmentsPerStep - giveOrTake, segmentsPerStep + giveOrTake);
    const startOffset = R.getAngle();
    let sum = 0;
    const thisSegments = Array.from({length: sc}, (e, i, r = R.get(.3, 1)) => (sum+=r)? r:r).map(e => e * 2 * Math.PI / sum);
    for(let i = 0, sum = 0; i < thisSegments.length; sum+=thisSegments[i++]) {
        segments.push({
            startR: r,
            r: radiusStep,
            startA: sum + startOffset,
            a: thisSegments[i]
        });
    }
}

const borderElements = (() => {
    const topLeftBorderCorner = [];
    for(let i = 0, max = (.5*Math.PI*borderRadius+1)|0; i <= max; i++) {
        topLeftBorderCorner.push([
            borderRadius * -Math.cos(i*Math.PI*.5/max)-100+border+borderRadius,
            borderRadius * -Math.sin(i*Math.PI*.5/max)-100+border+borderRadius
        ]);
    }
    if(borderRadius == 0) topLeftBorderCorner.pop();
    const topBorder = [[[-150, -150], [-150, -100+border+borderRadius]].concat(topLeftBorderCorner)].map(c => c.concat(c.map(e => [-e[0], e[1]]).reverse())).pop().reverse();
    const bottomBorder = topBorder.map(v => [v[0], -v[1]]).reverse();
    const pol = new Polygons();
    drawP(pol, bales[0], topBorder);
    drawP(pol, bales[0], bottomBorder);
    const leftBorder = [[-150, -150], [border-100, -150], [border-100, 150], [-150, 150]];
    const rightBorder = leftBorder.map(v=>[-v[0], v[1]]).reverse();
    drawP(pol, bales[0], leftBorder);
    drawP(pol, bales[0], rightBorder);
    PT.drawTour(bales[0], topLeftBorderCorner.concat(
                          topLeftBorderCorner.map(v=>[-v[0],v[1]]).reverse(),
                          topLeftBorderCorner.map(v=>[-v[0],-v[1]]),
                          topLeftBorderCorner.map(v=>[v[0],-v[1]]).reverse(),
    ));
    return [topBorder, rightBorder, bottomBorder, leftBorder];
})();

// The walk function will be called until it returns false.
function walk(i) {
    if(i < segments.length) {
        drawSegment(segments[i]);
    } else {
        drawCenter(radiusStep/2)
    }
    return i < segments.length - (centerPiece > 0? 0: 1);
}

function drawCenter(radius) {
    const hatchPolygons = new Polygons();
    const outlinePolygons = new Polygons();
    const pts = PT.circle(radius).filter((e,i,a) => i < a.length*.75).concat([[0,-50], [-50,0], [0, 50], [50, 0]]);
    
    const center = centerPiece == 1? [0,0]: [[radius * ((1-intraMidOffsetP)*.5+intraMidOffsetP*R.get()), R.getAngle()]].map(r => [
        r[0] * Math.sin(r[1]),
        r[0] * Math.cos(r[1])
    ]).pop();
    
    const maxR = V.len(center) + radius + intraRadiusDelta;

    if(hatching < 2) drawP(hatchPolygons, bales[0], pts);
    drawP(outlinePolygons, bales[0], pts);
    if(hatching < 2) drawP(hatchPolygons, bales[0], pts.map(pt => V.scale(pt, -1)));
    drawP(outlinePolygons, bales[0], pts.map(pt => V.scale(pt, -1)));
    
    fillCircles(hatchPolygons, outlinePolygons, maxR, center);
}

function arcPoints(r, fromA, toA, res=0) {
    const pts = [];
    const angle = Math.abs(toA - fromA);
    for(let i = 0, steps = res==0?(1 + angle * r)|0: res; i <= steps; i++) {
        const f = fromA <= toA? fromA + angle * i/steps: fromA - angle * i/steps;
        pts.push([r * Math.sin(f), r * Math.cos(f)]);
    }
    return pts;
}

function drawSegment(segment) {
    const innerPts = arcPoints(segment.startR, segment.startA, segment.startA + segment.a);
    const outerPts = arcPoints(segment.startR + segment.r, segment.startA + segment.a, segment.startA);

    if(!innerPts.some(pt => -100 < pt[0] && pt[0] < 100 && -100 < pt[1] && pt[1] < 100)) return;
    
    const tourPolygons = new Polygons();
    const tour = tourPolygons.create();
    tour.addPoints(...innerPts.concat(outerPts))
    tour.addOutline();

    const cwDisc = [].concat(
        arcPoints(segment.startR, segment.startA, segment.startA + segment.a),
        arcPoints(160, segment.startA + segment.a, segment.startA + 2 * Math.PI, 8),
    );
    const ccwDisc = [].concat(
        arcPoints(segment.startR + segment.r, segment.startA - .1, segment.startA + segment.a + .1),
        arcPoints(160, segment.startA + segment.a + .1, segment.startA - .1, 3),
    );

    const polygons = new Polygons();
    const polygonsForCircleOutline = new Polygons();
    
    if(hatching < 2) drawP(polygons, bales[0], cwDisc);
    drawP(polygonsForCircleOutline, bales[0], cwDisc);
    if(hatching < 2) drawP(polygons, bales[0], ccwDisc);
    drawP(polygonsForCircleOutline, bales[0], ccwDisc);
    
    if((border > 0 || borderRadius > 0) && segment.startR + segment.r > 100 - border) {
        borderElements.forEach(be => {
            const t = tourPolygons.create();
            t.addPoints(...be);
            tourPolygons.draw(bales[0], t);
            if(hatching < 2) drawP(polygons, bales[0], be);
            drawP(polygonsForCircleOutline, bales[0], be);
        });
    }
    tourPolygons.draw(bales[0], tour);

    const center = [segment.startR + segment.r * ((1-intraMidOffsetP)*.5+intraMidOffsetP*R.get())].map(r => [
        r * Math.sin(segment.startA + segment.a/2),
        r * Math.cos(segment.startA + segment.a/2)
    ]).pop();
    
    const distance = (1+Math.max(V.len(V.sub(center, innerPts[0])), V.len(V.sub(center, outerPts[0]))))|0;
    
    fillCircles(polygons, polygonsForCircleOutline, distance + intraRadiusDelta, center);
}

function fillCircles(hatchPolygons, outlinePolygons, maxR, center) {
    for(let r = intraRadiusDelta/2; r <= maxR; r+=intraRadiusDelta) {
        const pts = PT.circle(r).map(pt => V.add(pt, center));
        
        if(hatching < 2) drawP(hatchPolygons, bales[R.getInt(1, bales.length)], pts, 
            [R.getInt(0, hatchings.length)].map(r => r == hatchings.length - 1 && hatching != 0? false: hatchings[R.getInt(0, hatchings.length)]).pop()
        );
        drawP(outlinePolygons, bales[0], pts, false, true);
    }
}

function drawP(polygons, turtle, pts, hatch = false, outline = false) {
    const p = polygons.create();
    p.addPoints(...pts);
    if(outline) p.addOutline();
    if(hatch !== false) p.addHatching(...hatch);
    polygons.draw(turtle, p);
}

function init() {
    ///////////////////////////////////////////////////////
    // Vector functions - Created by Jurgen Westerhof 2024
    // https://turtletoy.net/turtle/d068ad6040
    ///////////////////////////////////////////////////////
    class Vector {
        static add  (a,b) { return a.map((v,i)=>v+b[i]); }
        static sub  (a,b) { return a.map((v,i)=>v-b[i]); }
        static mul  (a,b) { return a.map((v,i)=>v*b[i]); }
        static div  (a,b) { return a.map((v,i)=>v/b[i]); }
        static scale(a,s) { return a.map(v=>v*s); }
    
        static det(m)                { return m.length == 1? m[0][0]: m.length == 2 ? m[0][0]*m[1][1]-m[0][1]*m[1][0]: m[0].reduce((r,e,i) => r+(-1)**(i+2)*e*this.det(m.slice(1).map(c => c.filter((_,j) => i != j))),0); }
        static angle(a)              { return Math.PI - Math.atan2(a[1], -a[0]); } //compatible with turtletoy heading
        static rot2d(angle)          { return [[Math.cos(angle), -Math.sin(angle)], [Math.sin(angle), Math.cos(angle)]]; }
        static rot3d(yaw,pitch,roll) { return [[Math.cos(yaw)*Math.cos(pitch), Math.cos(yaw)*Math.sin(pitch)*Math.sin(roll)-Math.sin(yaw)*Math.cos(roll), Math.cos(yaw)*Math.sin(pitch)*Math.cos(roll)+Math.sin(yaw)*Math.sin(roll)],[Math.sin(yaw)*Math.cos(pitch), Math.sin(yaw)*Math.sin(pitch)*Math.sin(roll)+Math.cos(yaw)*Math.cos(roll), Math.sin(yaw)*Math.sin(pitch)*Math.cos(roll)-Math.cos(yaw)*Math.sin(roll)],[-Math.sin(pitch), Math.cos(pitch)*Math.sin(roll), Math.cos(pitch)*Math.cos(roll)]]; }
        static trans(matrix,a)       { return a.map((v,i) => a.reduce((acc, cur, ci) => acc + cur * matrix[ci][i], 0)); }
        //Mirror vector a in a ray through [0,0] with direction mirror
        static mirror2d(a,mirror)    { return [Math.atan2(...mirror)].map(angle => this.trans(this.rot2d(angle), this.mul([-1,1], this.trans(this.rot2d(-angle), a)))).pop(); }

        static equals(a,b)   { return !a.some((e, i) => e != b[i]); }
        static approx(a,b,p) { return this.len(this.sub(a,b)) < (p === undefined? .001: p); }
        static norm  (a)     { return this.scale(a,1/this.len(a)); }
        static len   (a)     { return Math.hypot(...a); }
        static lenSq (a)     { return a.reduce((a,c)=>a+c**2,0); }
        static lerp  (a,b,t) { return a.map((v, i) => v*(1-t) + b[i]*t); }
        static dist  (a,b)   { return Math.hypot(...this.sub(a,b)); }
        
        static dot  (a,b)   { return a.reduce((a,c,i) => a+c*b[i], 0); }
        static cross(...ab) { return ab[0].map((e, i) => ab.map(v => v.filter((ee, ii) => ii != i))).map((m,i) => (i%2==0?-1:1)*this.det(m)); }
        
        static clamp(a,min,max) { return a.map((e,i) => Math.min(Math.max(e, min[i]), max[i])) };
        static rotateClamp(a,min,max) { return a.map((e,i) => {const d = max[i]-min[i];if(d == 0) return min[i];while(e < min[i]) { e+=d; }while(e > max[i]) { e-=d; }return e;});
        }
    }
    this.V = Vector;
    
    class Intersection2D {
        //a-start, a-direction, b-start, b-direction
        //returns false on no intersection or [[intersection:x,y], scalar a-direction, scalar b-direction
        static info(as, ad, bs, bd) { const d = V.sub(bs, as), det = -V.det([bd, ad]); if(det === 0) return false; const res = [V.det([d, bd]) / det, V.det([d, ad]) / det]; return [V.add(as, V.scale(ad, res[0])), ...res]; }
        static ray(a, b, c, d) { return this.info(a, b, c, d); }
        static segment(a,b,c,d, inclusiveStart = true, inclusiveEnd = true) { const i = this.info(a, V.sub(b, a), c, V.sub(d, c)); return i === false? false: ( (inclusiveStart? 0<=i[1] && 0<=i[2]: 0<i[1] && 0<i[2]) && (inclusiveEnd?   i[1]<=1 && i[2]<=1: i[1]<1 && i[2]<1) )?i[0]:false;}
        static tour(tour, segmentStart, segmentDirection) { return tour.map((e, i, a) => [i, this.info(e, V.sub(a[(i+1)%a.length], e), segmentStart, segmentDirection)]).filter(e => e[1] !== false && 0 <= e[1][1] && e[1][1] <= 1).filter(e => 0 <= e[1][2]).map(e => ({position: e[1][0],tourIndex: e[0],tourSegmentPortion: e[1][1],segmentPortion: e[1][2],}));}
        static inside(tour, pt) { return tour.map((e,i,a) => this.segment(e, a[(i+1)%a.length], pt, [Number.MAX_SAFE_INTEGER, 0], true, false)).filter(e => e !== false).length % 2 == 1; }
        static circles(centerA, radiusA, centerB, radiusB) {const result = {intersect_count: 0,intersect_occurs: true,one_is_in_other: false,are_equal: false,point_1: [null, null],point_2: [null, null],};const dx = centerB[0] - centerA[0];const dy = centerB[1] - centerA[1];const dist = Math.hypot(dy, dx);if (dist > radiusA + radiusB) {result.intersect_occurs = false;}if (dist < Math.abs(radiusA - radiusB) && !N.approx(dist, Math.abs(radiusA - radiusB))) {result.intersect_occurs = false;result.one_is_in_other = true;}if (V.approx(centerA, centerB) && radiusA === radiusB) {result.are_equal = true;}if (result.intersect_occurs) {const centroid = (radiusA**2 - radiusB**2 + dist * dist) / (2.0 * dist);const x2 = centerA[0] + (dx * centroid) / dist;const y2 = centerA[1] + (dy * centroid) / dist;const prec = 10000;const h = (Math.round(radiusA**2 * prec)/prec - Math.round(centroid**2 * prec)/prec)**.5;const rx = -dy * (h / dist);const ry = dx * (h / dist);result.point_1 = [x2 + rx, y2 + ry];result.point_2 = [x2 - rx, y2 - ry];if (result.are_equal) {result.intersect_count = Infinity;} else if (V.equals(result.point_1, result.point_2)) {result.intersect_count = 1;} else {result.intersect_count = 2;}}return result;}
    }
    this.Intersection = Intersection2D;
    
    class PathTools {
        static bezier(p1, cp1, cp2, p2, steps = null) {steps = (steps === null? Math.max(3, (V.len(V.sub(cp1, p1)) + V.len(V.sub(cp2, cp1)) + V.len(V.sub(p2, cp2))) | 0): steps) - 1; return Array.from({length: steps + 1}).map((v, i, a, f = i/steps) => [[V.lerp(p1, cp1, f),V.lerp(cp1, cp2, f),V.lerp(cp2, p2, f)]].map(v => V.lerp(V.lerp(v[0], v[1], f), V.lerp(v[1], v[2], f), f))[0]);}
        // https://stackoverflow.com/questions/18655135/divide-bezier-curve-into-two-equal-halves#18681336
        static splitBezier(p1, cp1, cp2, p2, t=.5) {const e = V.lerp(p1, cp1, t);const f = V.lerp(cp1, cp2, t);const g = V.lerp(cp2, p2, t);const h = V.lerp(e, f, t);const j = V.lerp(f, g, t);const k = V.lerp(h, j, t);return [[p1, e, h, k], [k, j, g, p2]];}
        static circular(radius,verticeCount,rotation=0) {return Array.from({length: verticeCount}).map((e,i,a,f=i*2*Math.PI/verticeCount+rotation) => [radius*Math.cos(f),radius*Math.sin(f)])}
        static circle(r){return this.circular(r,Math.max(12, r*2*Math.PI|0));}
        static arc(radius, extend = 2 * Math.PI, clockWiseStart = 0, steps = null, includeLast = false) { return [steps == null? (radius*extend+1)|0: steps].map(steps => Array.from({length: steps}).map((v, i, a) => [radius * Math.cos(clockWiseStart + extend*i/(a.length-(includeLast?1:0))), radius * Math.sin(clockWiseStart + extend*i/(a.length-(includeLast?1:0)))])).pop(); }
        static draw(turtle, path) {path.forEach((pt, i) => turtle[i==0?'jump':'goto'](pt));}
        static drawTour(turtle, path) {this.draw(turtle, path.concat([path[0]]));}
        static drawPoint(turtle, pt, r = .1) {this.drawTour(turtle, this.circle(r).map(e => V.add(e, pt)));}
        static drawArrow(turtle, s, d, width = 6, length = 3) {turtle.jump(s);const arrowHeadBase = V.add(s,d);turtle.goto(arrowHeadBase);turtle.goto(V.add(arrowHeadBase, V.trans(V.rot2d(-V.angle(d)), [-length, width/2])));turtle.jump(V.add(arrowHeadBase, V.trans(V.rot2d(-V.angle(d)), [-length, -width/2])));turtle.goto(arrowHeadBase);}
        static circlesTangents(c1_center, c1_radius, c2_center, c2_radius, internal = false) {let middle_circle = [V.scale(V.sub(c1_center, c2_center), .5)].map(hwp => [V.add(c2_center, hwp), V.len(hwp)]).pop();if(!internal && c1_radius == c2_radius) {let target = V.sub(c2_center, c1_center);let scaledTarget = V.scale(target, c1_radius/V.len(target));let partResult = [V.add(c1_center, V.trans(V.rot2d(Math.PI/2), scaledTarget)),V.add(c1_center, V.trans(V.rot2d(Math.PI/-2), scaledTarget))];return [partResult,partResult.map(pt => V.add(pt, target))]}let swap = !internal && c2_radius > c1_radius;if(swap) {let t = [[...c1_center], c1_radius];c1_center = c2_center;c1_radius = c2_radius;c2_center = t[0];c2_radius = t[1];}let internal_waypoints = Intersection.circles(c1_center, c1_radius + (internal?c2_radius:-c2_radius), ...middle_circle);if(!internal_waypoints.intersect_occurs) return [];const circlePointAtDirection2 = (circle_center, radius, direction) => V.add(circle_center, V.scale(direction, radius/V.len(direction)));const result = [[circlePointAtDirection2(c1_center, c1_radius, V.sub(internal_waypoints.point_1, c1_center)),circlePointAtDirection2(c1_center, c1_radius, V.sub(internal_waypoints.point_2, c1_center))],[circlePointAtDirection2(c2_center, c2_radius, internal?V.sub(c1_center, internal_waypoints.point_1):V.sub(internal_waypoints.point_1, c1_center)),circlePointAtDirection2(c2_center, c2_radius, internal?V.sub(c1_center, internal_waypoints.point_2):V.sub(internal_waypoints.point_2, c1_center))]];return swap? [[result[1][1],result[1][0]],[result[0][1],result[0][0]]]: result;}
        static vectors(path) {return Array.from({length: path.length - 1}).map((e, i) => V.sub(path[i+1], path[i]));}
        static path(vectors) {return vectors.reduce((a,c) => a.length==0?[c]:[...a, V.add(c, a[a.length-1])], []);}
        static redistributeLinear(path, length = .5) {const p = path.map(pt => [...pt]);const result = [[...p[0]]];let pointer = 1;doneAll: while(pointer < p.length) {let l = length;while(pointer < p.length) {const distance = V.len(V.sub(p[pointer], p[pointer - 1]));if(distance < l) {l-=distance;pointer++;continue;}if(distance == l) {if(pointer < p.length - 1) result.push([...p[pointer]]);pointer++;break doneAll;}if(l < distance) {const newPoint = V.lerp(p[pointer-1], p[pointer], l/distance);if(pointer < p.length - 1) result.push([...newPoint]);p[pointer - 1] = newPoint;break;}}}result.push(p.pop());return result;}
        static length(path) { return this.lengths(path).reduce((c, a) => a + c, 0); }
        static lengths(path) { return path.map((e, i, a) => V.len(V.sub(e, a[(i+1)%a.length]))).filter((e, i, a) => i < a.length - 1); }
        static intersectInfoRay(path, origin, direction) {const vectors = this.vectors(path);const ri = vectors.map((e, i) => [i, Intersection.info(origin, direction, path[i], e)]).filter(e => 0 <= e[1][2] && e[1][2] <= 1 && 0 < e[1][1]).sort(e => e[1][1]);if(ri.length == 0) return false;const hit = ri[0];const lengths = this.lengths(path);const length = lengths.reduce((a, c) => a + c, 0);let l = 0;for(let i = 0; i < hit[0]; i++) {l += lengths[i];}return [hit[1][0], (l + (lengths[hit[0]] * hit[1][2])) / length, hit[1][1]];}
        static lerp(path, part) {if(part < 0 || 1 < part) throw new Error('Range of part is 0 to 1, got ' + path);const lengths = this.lengths(path);const length = lengths.reduce((a, c) => a + c, 0);let l = length * part;for(let i = 0; i < lengths.length; i++) {if(lengths[i] < l) {l-=lengths[i];continue;}return V.lerp(path[i], path[i+1], l / V.len(V.sub(path[i+1], path[i])));}return [...path[path.length - 1]];}
    }
    
    this.PT = PathTools;
    
    class Complex {
        static add(a,b)     { return V.add(a,b); }
        static sub(a,b)     { return V.sub(a,b); }
        static scale(a,s)   { return V.scale(a,s); }
        static mult(a,b)    { return [a[0]*b[0]-a[1]*b[1],a[0]*b[1]+a[1]*b[0]]; }
        static sqrt(a)      { return [[Math.hypot(...a)**.5, Math.atan2(...a.reverse()) / 2]].map(ra => [ra[0]*Math.cos(ra[1]), ra[0]*Math.sin(ra[1])]).pop(); }
    }
    this.C = Complex;
    
    class Numbers {
        static approx(a,b,p)            { return Math.abs(a-b) < (p === undefined? .001: p); }
        static clamp(a, min, max)       { return Math.min(Math.max(a, min), max); }
        static rotateClamp(a, min, max) { if(min == max) return min;while (a < min) { a+=(max-min); }while (a > max) { a-=(max-min); }return a;}
        static lerp(a, b, t)            { return V.lerp([a], [b], t); };
    }
    this.N = Numbers;

    class Matrix {
        static bayer(order) { return [...Array(1<<order)].map((_,y,a) => { const g = (k=order,x)=>k--&&4*g(k,x)|2*(x>>k)+3*(y>>k&1)&3; return a.map(g); }); }
        static rotate(m) { return m[0].map((e, i) => m.map(r => r[i]).reverse()); }
        static rotateCCW(m) { return m[0].map((e, i) => m.map(r => r[r.length-1-i])); }
        static add(a,b) { return a.map((e, c) => e.map((e, r) => a[c][r] + b[c][r])); }
        static sub(a,b) { return a.map((e, c) => e.map((e, r) => a[c][r] - b[c][r])); }
        static multiply(a,b) { return Array.from({length: b.length}, (e,resCol) => Array.from({length: a[0].length}, (e,resRow) => b[resCol].reduce((acc, c, bRow) => acc + a[bRow][resRow] * b[resCol][bRow], 0)));}
        static scale(a,s) { return a.map((e, c) => e.map((e, r) => a[c][r] * s)); }
        static random(c,r,fillFn = Math.random) { return Array.from({length: c}, (e,i) => Array.from({length: r}, e => fillFn(c, r))); }
        static identity(d) { return Array.from({length: d}, (e,c) => Array.from({length: d}, (e, r) => c==r?1:0 )); }
        static log(m, name, logFn = console.log) { if(name != undefined) logFn(name); if(m === undefined || (typeof m == 'object' && (m[0] === undefined || m[0][0] === undefined))) { return logFn(`Failed to log matrix:`, m); } logFn(m[0].map((e,r) => m.map((e,c) => m[c][r]).join(', ')).join('\n')); }
        static invert(m) { let _A = m.map(col => col.map(cell => cell));/*clone matrix*/let temp;const N = _A.length;const E = Array.from({length: N}, (e,i) => Array.from({length: _A[0].length}, (e,j) => i==j?1:0));for (let k = 0; k < N; k++) {temp = _A[k][k];for (let j = 0; j < N; j++) {_A[k][j] /= temp;E[k][j] /= temp;}for (let i = k + 1; i < N; i++) {temp = _A[i][k];for (let j = 0; j < N; j++) {_A[i][j] -= _A[k][j] * temp;E[i][j] -= E[k][j] * temp;}}}for (let k = N - 1; k > 0; k--) {for (let i = k - 1; i >= 0; i--) {temp = _A[i][k];for (let j = 0; j < N; j++) {_A[i][j] -= _A[k][j] * temp;E[i][j] -= E[k][j] * temp;}}}return E; }
        static determinant(m) { return m.length == 1 ?m[0][0] :m.length == 2 ? m[0][0]*m[1][1]-m[0][1]*m[1][0] :m[0].reduce((r,e,i) => r+(-1)**(i+2)*e*this.determinant(m.slice(1).map(c => c.filter((_,j) => i != j))), 0)}
        static flip(m) { return Array.from({length: m[0].length}, (_, r) => Array.from({length: m.length}, (e, c) => m[c][r])); }
        static sum(m) { return m.reduce((a, c) => a + c.reduce((aa, cc) => aa + cc, 0), 0); }
    }
    this.M = Matrix;
    
    class Algorithms {
        static nthTriangular(n) { return ((n * n) + n) / 2; }
    }
    this.A = Algorithms;
    
    class Random {
        static #apply(seed) {
            // Seedable random number generator by David Bau: http://davidbau.com/archives/2010/01/30/random_seeds_coded_hints_and_quintillions.html
            !function(a,b,c,d,e,f,g,h,i){function j(a){var b,c=a.length,e=this,f=0,g=e.i=e.j=0,h=e.S=[];for(c||(a=[c++]);d>f;)h[f]=f++;for(f=0;d>f;f++)h[f]=h[g=s&g+a[f%c]+(b=h[f])],h[g]=b;(e.g=function(a){for(var b,c=0,f=e.i,g=e.j,h=e.S;a--;)b=h[f=s&f+1],c=c*d+h[s&(h[f]=h[g=s&g+b])+(h[g]=b)];return e.i=f,e.j=g,c})(d)}function k(a,b){var c,d=[],e=typeof a;if(b&&"object"==e)for(c in a)try{d.push(k(a[c],b-1))}catch(f){}return d.length?d:"string"==e?a:a+"\0"}function l(a,b){for(var c,d=a+"",e=0;e<d.length;)b[s&e]=s&(c^=19*b[s&e])+d.charCodeAt(e++);return n(b)}function m(c){try{return o?n(o.randomBytes(d)):(a.crypto.getRandomValues(c=new Uint8Array(d)),n(c))}catch(e){return[+new Date,a,(c=a.navigator)&&c.plugins,a.screen,n(b)]}}function n(a){return String.fromCharCode.apply(0,a)}var o,p=c.pow(d,e),q=c.pow(2,f),r=2*q,s=d-1,t=c["seed"+i]=function(a,f,g){var h=[];f=1==f?{entropy:!0}:f||{};var o=l(k(f.entropy?[a,n(b)]:null==a?m():a,3),h),s=new j(h);return l(n(s.S),b),(f.pass||g||function(a,b,d){return d?(c[i]=a,b):a})(function(){for(var a=s.g(e),b=p,c=0;q>a;)a=(a+c)*d,b*=d,c=s.g(1);for(;a>=r;)a/=2,b/=2,c>>>=1;return(a+c)/b},o,"global"in f?f.global:this==c)};if(l(c[i](),b),g&&g.exports){g.exports=t;try{o=require("crypto")}catch(u){}}else h&&h.amd&&h(function(){return t})}(this,[],Math,256,6,52,"object"==typeof module&&module,"function"==typeof define&&define,"random");
            Math.seedrandom(seed);
        }
        static seedRandom() { this.#apply(new Date().getMilliseconds()); }
        static seedDaily() { this.#apply(new Date().toDateString()); }
        static seed(seed) { this.#apply(seed); }
        static getInt(min, max) { if(max == undefined) {max = min + 1; min = 0; } const [mi, ma] = [Math.min(min, max), Math.max(min, max)]; return (mi + Math.random() * (ma - mi)) | 0;}
        static get(min, max) {if(min == undefined) {return Math.random();}if(max == undefined) {max = min;min = 0;}const [mi, ma] = [Math.min(min, max), Math.max(min, max)];return mi + Math.random() * (ma - mi);}
        static getAngle(l = 1) { return l * this.get(0, 2*Math.PI); }
        // Standard Normal variate using Box-Muller transform.
        static getGaussian(mean=.5, stdev=.1) {const u = 1 - this.get(); /* Converting [0,1) to (0,1] */const v = this.get();const z = ( -2.0 * Math.log( u ) )**.5 * Math.cos( 2.0 * Math.PI * v );/* Transform to the desired mean and standard deviation: */return z * stdev + mean;}
        static skew(value, skew = 0) { /*skew values (from 0 to 1) by a skew from -1 to 1, respectively right and left skewed (resp more values to left or to right), 0 is not skewed*/return (skew < 0)? value - (this.skew(1-value, -skew) - (1-value)): Math.pow(value, 1-Math.abs(skew));}
        static getNormalDistributed(skew = 0) { /*skew values (from 0 to 1) by a skew from -1 to 1, respectively right and left skewed (resp more values to left or to right), 0 is not skewed*/let v = -1;while(v < 0 || 1 <= v) { v = this.getGaussian(.5, .1) };return this.skew(v, skew);}
    }
    this.R = Random;
}

////////////////////////////////////////////////////////////////
// Bale utility code - Created by Jurgen Westerhof 2022
// https://turtletoy.net/turtle/beb59d67ae
// Abusing the opacity, usage:
//      Canvas.setpenopacity(1/paletteSize);
//      const bales = Bales(paletteSize); // Bales(count, includeFullTransparent = true, turtleClass = null)
// Then use bales[x] wherever you would use a turtle object to 'draw'
// in 'color' x (i.e Polygon hatching with a bale object and .15 interspacing)
//      bales[x].jump(0,0);
//      bales[x].goto(40,0);
////////////////////////////////////////////////////////////////
function Bale(n, turtleClass = null) {class Bale {constructor(n, turtleClass = null) { this.turtles = Array.apply(null,{length: n}).map(i => turtleClass == null? new Turtle(): new turtleClass()); }back(e)         { this.turtles.forEach(t => t.back(e)); return this; }backward(e)     { this.turtles.forEach(t => t.backward(e)); return this; }bk(e)           { this.turtles.forEach(t => t.bk(e)); return this; }fd(e)           { this.turtles.forEach(t => t.fd(e)); return this; }forward(e)      { this.turtles.forEach(t => t.forward(e)); return this; }left(e)         { this.turtles.forEach(t => t.left(e)); return this; }lt(e)           { this.turtles.forEach(t => t.lt(e)); return this; }right(e)        { this.turtles.forEach(t => t.right(e)); return this; }rt(e)           { this.turtles.forEach(t => t.rt(e)); return this; }seth(e)         { this.turtles.forEach(t => t.seth(e)); return this; }setheading(e)   { this.turtles.forEach(t => t.setheading(e)); return this; }setx(e)         { this.turtles.forEach(t => t.setx(e)); return this; }sety(e)         { this.turtles.forEach(t => t.sety(e)); return this; }setpos(x, y)        { this.turtles.forEach(t => t.setpos(x, y)); return this; }setposition(x, y)   { this.turtles.forEach(t => t.setposition(x, y)); return this; }toradians(e)    { this.turtles.forEach(t => t.toradians(e)); return this; }degrees(e)      { this.turtles.forEach(t => t.degrees(e)); return this; }goto(x, y)      { this.turtles.forEach(t => t.goto(x, y)); return this; }jmp(x, y)       { this.turtles.forEach(t => t.jmp(x, y)); return this; }jump(x, y)      { this.turtles.forEach(t => t.jump(x, y)); return this; }circle(radius, extent, steps) { this.turtles.map(t => t.circle(radius, extent, steps)); return this; }clone()         { let b = new Bale(this.turtles.length); this.turtles.forEach((t, k) => b.turtles[k] = t.clone()); return b; }h()             { return this.turtles.length == 0? null: this.turtles[0].h(); }heading()       { return this.turtles.length == 0? null: this.turtles[0].heading(); }home()          { this.turtles.forEach(t => t.home()); return this; }isdown()        { return this.turtles.length == 0? null: this.turtles[0].isdown(); }pos()           { return this.turtles.length == 0? null: this.turtles[0].pos(); }position()      { return this.turtles.length == 0? null: this.turtles[0].position(); }pd()            { this.turtles.forEach(t => t.pd()); return this; }pendown()       { this.turtles.forEach(t => t.pendown()); return this; }penup()         { this.turtles.forEach(t => t.penup()); return this; }pu()            { this.turtles.forEach(t => t.pu()); return this; }down()          { this.turtles.forEach(t => t.down()); return this; }up()            { this.turtles.forEach(t => t.up()); return this; }radians()       { this.turtles.forEach(t => t.radians()); return this; }x()             { return this.turtles.length == 0? null: this.turtles[0].x(); }xcor()          { return this.turtles.length == 0? null: this.turtles[0].xcor(); }y()             { return this.turtles.length == 0? null: this.turtles[0].y(); }ycor()          { return this.turtles.length == 0? null: this.turtles[0].ycor(); }set(key, value) { this.turtles.forEach(i => i[key] = value); return this; }get(key) { return this.turtles.length == 0? null: this.turtles[0][key]; }}return new Bale(n, turtleClass);}
function Bales(count, includeFullTransparent = true, turtleClass = null) { if(count == 1) return [new Bale(1, turtleClass)]; const getExponent = (base, target) => Math.log(target) / Math.log(base); const baleSize = count - (includeFullTransparent?1:0); const n = Array.apply(null,{length: baleSize}).map((v,k) => Math.round(getExponent(1 - 1/count, 1 - (count - k == count?.99:(baleSize - k)/baleSize)))); if(includeFullTransparent) n.push(0); return n.map(i => new Bale(i, turtleClass));}

////////////////////////////////////////////////////////////////
// Polygon Clipping utility code - Created by Reinder Nijhoff 2019
// (Polygon binning by Lionel Lemarie 2021) https://turtletoy.net/turtle/95f33bd383
// (Delegated Hatching by Jurgen Westerhof 2024) https://turtletoy.net/turtle/d068ad6040
// (Deferred Polygon Drawing by Jurgen Westerhof 2024) https://turtletoy.net/turtle/6f3d2bc0b5
// https://turtletoy.net/turtle/a5befa1f8d
//
// const polygons = new Polygons();
// const p = polygons.create();
// polygons.draw(turtle, p);
// polygons.list();
// polygons.startDeferSession();
// polygons.stopDeferring();
// polygons.finalizeDeferSession(turtle);
//
// p.addPoints(...[[x,y],]);
// p.addSegments(...[[x,y],]);
// p.addOutline();
// p.addHatching(angle, distance); OR p.addHatching(HatchObject); where HatchObject has a method 'hatch(PolygonClass, thisPolygonInstance)'
// p.inside([x,y]);
// p.boolean(polygon, diff = true);
// p.segment_intersect([x,y], [x,y], [x,y], [x,y]);
////////////////////////////////////////////////////////////////
function Polygons(){const t=[],s=25,e=Array.from({length:s**2},t=>[]),n=class{constructor(){this.cp=[],this.dp=[],this.aabb=[]}addPoints(...t){let s=1e5,e=-1e5,n=1e5,h=-1e5;(this.cp=[...this.cp,...t]).forEach(t=>{s=Math.min(s,t[0]),e=Math.max(e,t[0]),n=Math.min(n,t[1]),h=Math.max(h,t[1])}),this.aabb=[s,n,e,h]}addSegments(...t){t.forEach(t=>this.dp.push(t))}addOutline(){for(let t=0,s=this.cp.length;t<s;t++)this.dp.push(this.cp[t],this.cp[(t+1)%s])}draw(t){for(let s=0,e=this.dp.length;s<e;s+=2)t.jump(this.dp[s]),t.goto(this.dp[s+1])}addHatching(t, s) {if(typeof t == 'object') return t.hatch(n, this);const e=new n;e.cp.push([-1e5,-1e5],[1e5,-1e5],[1e5,1e5],[-1e5,1e5]);const h=Math.sin(t)*s,o=Math.cos(t)*s,a=200*Math.sin(t),i=200*Math.cos(t);for(let t=.5;t<150/s;t++) {e.dp.push([h*t+i,o*t-a],[h*t-i,o*t+a]);e.dp.push([-h*t+i,-o*t-a],[-h*t-i,-o*t+a]);}e.boolean(this,!1);this.dp=[...this.dp,...e.dp]}inside(t){let s=0;for(let e=0,n=this.cp.length;e<n;e++)this.segment_intersect(t,[.1,-1e3],this.cp[e],this.cp[(e+1)%n])&&s++;return 1&s}boolean(t,s=!0){const e=[];for(let n=0,h=this.dp.length;n<h;n+=2){const h=this.dp[n],o=this.dp[n+1],a=[];for(let s=0,e=t.cp.length;s<e;s++){const n=this.segment_intersect(h,o,t.cp[s],t.cp[(s+1)%e]);!1!==n&&a.push(n)}if(0===a.length)s===!t.inside(h)&&e.push(h,o);else{a.push(h,o);const n=o[0]-h[0],i=o[1]-h[1];a.sort((t,s)=>(t[0]-h[0])*n+(t[1]-h[1])*i-(s[0]-h[0])*n-(s[1]-h[1])*i);for(let n=0;n<a.length-1;n++)(a[n][0]-a[n+1][0])**2+(a[n][1]-a[n+1][1])**2>=.001&&s===!t.inside([(a[n][0]+a[n+1][0])/2,(a[n][1]+a[n+1][1])/2])&&e.push(a[n],a[n+1])}}return(this.dp=e).length>0}segment_intersect(t,s,e,n){const h=(n[1]-e[1])*(s[0]-t[0])-(n[0]-e[0])*(s[1]-t[1]);if(0===h)return!1;const o=((n[0]-e[0])*(t[1]-e[1])-(n[1]-e[1])*(t[0]-e[0]))/h,a=((s[0]-t[0])*(t[1]-e[1])-(s[1]-t[1])*(t[0]-e[0]))/h;return o>=0&&o<=1&&a>=0&&a<=1&&[t[0]+o*(s[0]-t[0]),t[1]+o*(s[1]-t[1])]}};const y=function(n,j=[]){const h={},o=200/s;for(var a=0;a<s;a++){const c=a*o-100,r=[0,c,200,c+o];if(!(n[3]<r[1]||n[1]>r[3]))for(var i=0;i<s;i++){const c=i*o-100;r[0]=c,r[2]=c+o,n[0]>r[2]||n[2]<r[0]||e[i+a*s].forEach(s=>{const e=t[s];n[3]<e.aabb[1]||n[1]>e.aabb[3]||n[0]>e.aabb[2]||n[2]<e.aabb[0]||j.includes(s)||(h[s]=1)})}}return Array.from(Object.keys(h),s=>t[s])};return{list:()=>t,create:()=>new n,draw:(n,h,o=!0)=>{rpl=y(h.aabb, this.dei === undefined? []: Array.from({length: t.length - this.dei}).map((e, i) => this.dsi + i));for(let t=0;t<rpl.length&&h.boolean(rpl[t]);t++);const td=n.isdown();if(this.dsi!==undefined&&this.dei===undefined)n.pu();h.draw(n),o&&function(n){t.push(n);const h=t.length-1,o=200/s;e.forEach((t,e)=>{const a=e%s*o-100,i=(e/s|0)*o-100,c=[a,i,a+o,i+o];c[3]<n.aabb[1]||c[1]>n.aabb[3]||c[0]>n.aabb[2]||c[2]<n.aabb[0]||t.push(h)})}(h);if(td)n.pd();},startDeferSession:()=>{if(this.dei!==undefined)throw new Error('Finalize deferring before starting new session');this.dsi=t.length;},stopDeferring:()=>{if(this.dsi === undefined)throw new Error('Start deferring before stopping');this.dei=t.length;},finalizeDeferSession:(n)=>{if(this.dei===undefined)throw new Error('Stop deferring before finalizing');for(let i=this.dsi;i<this.dei;i++) {rpl = y(t[i].aabb,Array.from({length:this.dei-this.dsi+1}).map((e,j)=>i+j));for(let j=0;j<rpl.length&&t[i].boolean(rpl[j]);j++);t[i].draw(n);}this.dsi=undefined;this.dei=undefined;}}}