{"id":"fb4d519e95cdbd83","slug":"cubic-spline","trashed":false,"description":"","likes":34,"publish_level":"public","forks":0,"fork_of":null,"has_importers":true,"update_time":"2021-02-04T15:48:13.581Z","first_public_version":null,"paused_version":null,"publish_time":"2019-10-21T00:50:11.167Z","publish_version":1977,"latest_version":1977,"thumbnail":"7af80c742971e17dfeceabdc5fab7d891b59cb0543db18d8ad9b8f1e9d60521b","default_thumbnail":"7af80c742971e17dfeceabdc5fab7d891b59cb0543db18d8ad9b8f1e9d60521b","roles":[],"sharing":null,"owner":{"id":"9b8a8e177a49059a","avatar_url":"https://avatars.observableusercontent.com/avatar/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","login":"jrus","name":"Jacob Rus","bio":"geodetic computer","home_url":"","type":"team","tier":"starter_2024"},"creator":{"id":"ddd7235babc79ce0","avatar_url":"https://avatars.observableusercontent.com/avatar/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","login":"jrus","name":"Jacob Rus","bio":"geodetic computer","home_url":"","tier":"pro"},"authors":[{"id":"ddd7235babc79ce0","avatar_url":"https://avatars.observableusercontent.com/avatar/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","name":"Jacob Rus","login":"jrus","bio":"geodetic computer","home_url":"","tier":"pro","approved":true,"description":""}],"collections":[],"files":[],"comments":[{"id":"3c2b3beacc76f323","content":"Looks like this static import isn't used anywhere and can be removed.","node_id":1530,"create_time":"2021-12-10T21:28:42.826Z","update_time":null,"resolved":false,"user":{"id":"07362516b5994994","avatar_url":"https://avatars.observableusercontent.com/avatar/0d7defa821f38094c03bad23b9b360a5364e6e97e21fc238c39ddc48db7994ad","login":"mootari","name":"Fabian Iwand","bio":"Web dev and tinkerer.","home_url":"https://mootari.de/","tier":"pro"}},{"id":"7ec5e004c9ba1fba","content":"Typo: \"and returns a return a function\"","node_id":1829,"create_time":"2021-12-12T12:20:28.864Z","update_time":null,"resolved":false,"user":{"id":"07362516b5994994","avatar_url":"https://avatars.observableusercontent.com/avatar/0d7defa821f38094c03bad23b9b360a5364e6e97e21fc238c39ddc48db7994ad","login":"mootari","name":"Fabian Iwand","bio":"Web dev and tinkerer.","home_url":"https://mootari.de/","tier":"pro"}}],"commenting_lock":null,"suggestion_from":null,"suggestions_to":[],"version":1977,"title":"Cubic Spline","license":"isc","copyright":"Copyright © 2019–2021 Jacob Rus","nodes":[{"id":0,"value":"Title = md`# Cubic Spline`","pinned":false,"mode":"js","data":null,"name":null},{"id":1951,"value":"md`\nNothing here is fundamentally new, but I tried to pick notation and write the algebraic explanation in a way that I find a bit clearer and less cluttered than other explanations I have seen in books or online. YMMV.\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":938,"value":"Thumbnail = {\n const svgwidth = width; // Math.min(width, 640);\n const svgheight = svgwidth * 5/8\n const points = [[.25, .4], [.4, .5], [.5, .55], [.65, .4], [.95, .7]];\n \n const svg = d3.select(DOM.svg(svgwidth, svgheight))\n .attr('style', `background-color: #fbfcfd;`);\n \n const colors = [\n \"#380432\", \"#061a47\", \"#581341\", \"#032f6c\", \"#7c2250\", \"#024983\",\n \"#a43452\", \"#02678c\", \"#c84c57\", \"#088590\", \"#e16c5d\", \"#07a48a\",\n \"#f3906e\", \"#63bd8d\", \"#fdb774\", \"#a7d46c\", \"#fdde88\", \"#dfe976\",\n \"black\"]\n \n const spread = [\n +0.45, -0.45, +0.40, -0.40, +0.35, -0.35,\n +0.30, -0.30, +0.25, -0.25, +0.20, -0.20,\n +0.15, -0.15, +0.10, -0.10, +0.05, -0.05, 0]\n \n const paths = new Array(spread.length);\n \n const curve = svg.append(\"g\")\n .attr(\"fill\", \"none\")\n .selectAll(\"path\")\n .data(spread)\n .enter().append(\"path\")\n .attr(\"stroke\", (p, i) => colors[i])\n .attr(\"stroke-width\", (p, i) => ((i == colors.length - 1) ? 2.5 : 1) * svgwidth*1e-3 + .5 );\n\n const colorpoints = svg.append(\"g\")\n .selectAll(\"circle\")\n .data(spread)\n .enter().append(\"circle\")\n .attr(\"r\", 2e-3*svgwidth + .5)\n .attr(\"fill\", (k,i) => colors[i]);\n \n const dragged = function dragged(d, i) {\n const xraw = d3.event.x / svgwidth + d3.event.subject.x;\n const yraw = d3.event.y / svgheight + d3.event.subject.y;\n \n const prevx = ((i > 0) ? points[i-1][0] : 0) + .01;\n const nextx = ((i < (points.length - 1)) ? points[i+1][0] : 1) - .01;\n \n d[0] = Math.max(Math.min(xraw, nextx), prevx);\n d[1] = Math.max(Math.min(yraw, .99), .01);\n update();\n }\n \n const point = svg.append(\"g\")\n .attr(\"cursor\", \"move\")\n .attr(\"pointer-events\", \"all\")\n .attr(\"stroke\", \"transparent\")\n .attr(\"stroke-width\", 20)\n .selectAll(\"circle\")\n .data(points)\n .enter().append(\"circle\")\n .attr(\"r\", svgwidth*4e-3 + 1)\n .attr(\"fill\", 'black')\n .call(d3.drag()\n .subject(([x, y]) => ({x, y}))\n .on(\"drag\", dragged));\n\n const update = function update() {\n for (let i = 0; i < spread.length; i++) {\n let t = points.map(x => x[0]);\n let y = points.map(x => x[1]);\n y[1] += spread[i];\n \n let spline = CyclicCubicSpline(t, y)\n\n const bezpts = draw_hermite_spline(spline.knots, spline.values, spline.dvalues)\n .map((x, i) => !(i % 2) ? x*svgwidth : x*svgheight); // scale\n paths[i] = bezpts_to_svgpath(bezpts);\n\n// let m = cyclic_cubic_spline(t, y)\n// const n = t.length;\n// t = [t[n-1]-1, ...t, t[0] + 1];\n// y = [y[n-1], ...y, y[0]];\n// m = [m[n-1], ...m, m[0]];\n\n// const bezpts = draw_hermite_spline(t, y, m)\n// .map((x, i) => !(i % 2) ? x*svgwidth : x*svgheight); // scale\n// paths[i] = bezpts_to_svgpath(bezpts);\n }\n \n curve\n .attr(\"d\", (p, i) => paths[i]);\n point\n .attr(\"cx\", ([x,y]) => x * svgwidth)\n .attr(\"cy\", ([x,y]) => y * svgheight);\n colorpoints\n .attr(\"cx\", (k) => points[1][0] * svgwidth)\n .attr(\"cy\", (k) => (points[1][1] + k) * svgheight);\n\n }\n\n update(); \n \n return html`\n
\n ${svg.node()}\n
\n A cyclic cubic spline. Changing the value at one point affects the spline globally.\n
\n
\n
`\n}","pinned":false,"mode":"js","data":null,"name":null},{"id":1813,"value":"ThumbnailNatural = {\n const svgwidth = width; // Math.min(width, 640);\n const svgheight = svgwidth * 5/8\n const points = [[.25, .55], [.4, .55], [.5, .5], [.66, .35], [.95, .7]];\n \n const svg = d3.select(DOM.svg(svgwidth, svgheight))\n .attr('style', `background-color: #fbfcfd;`);\n \n const colors = [\n \"#380432\", \"#061a47\", \"#581341\", \"#032f6c\", \"#7c2250\", \"#024983\",\n \"#a43452\", \"#02678c\", \"#c84c57\", \"#088590\", \"#e16c5d\", \"#07a48a\",\n \"#f3906e\", \"#63bd8d\", \"#fdb774\", \"#a7d46c\", \"#fdde88\", \"#dfe976\",\n \"black\"]\n \n const spread = [\n +0.45, -0.45, +0.40, -0.40, +0.35, -0.35,\n +0.30, -0.30, +0.25, -0.25, +0.20, -0.20,\n +0.15, -0.15, +0.10, -0.10, +0.05, -0.05, 0]\n \n const paths = new Array(spread.length);\n \n const curve = svg.append(\"g\")\n .attr(\"fill\", \"none\")\n .selectAll(\"path\")\n .data(spread)\n .enter().append(\"path\")\n .attr(\"stroke\", (p, i) => colors[i])\n .attr(\"stroke-width\", (p, i) => ((i == colors.length - 1) ? 2.5 : 1) * svgwidth*1e-3 + .5 );\n\n const colorpoints = svg.append(\"g\")\n .selectAll(\"circle\")\n .data(spread)\n .enter().append(\"circle\")\n .attr(\"r\", 2e-3*svgwidth + .5)\n .attr(\"fill\", (k,i) => colors[i]);\n \n const dragged = function dragged(d, i) {\n const xraw = d3.event.x / svgwidth + d3.event.subject.x;\n const yraw = d3.event.y / svgheight + d3.event.subject.y;\n \n const prevx = ((i > 0) ? points[i-1][0] : 0) + .01;\n const nextx = ((i < (points.length - 1)) ? points[i+1][0] : 1) - .01;\n \n d[0] = Math.max(Math.min(xraw, nextx), prevx);\n d[1] = Math.max(Math.min(yraw, .99), .01);\n update();\n }\n \n const point = svg.append(\"g\")\n .attr(\"cursor\", \"move\")\n .attr(\"pointer-events\", \"all\")\n .attr(\"stroke\", \"transparent\")\n .attr(\"stroke-width\", 20)\n .selectAll(\"circle\")\n .data(points)\n .enter().append(\"circle\")\n .attr(\"r\", svgwidth*4e-3 + 1)\n .attr(\"fill\", 'black')\n .call(d3.drag()\n .subject(([x, y]) => ({x, y}))\n .on(\"drag\", dragged));\n\n const update = function update() {\n for (let i = 0; i < spread.length; i++) {\n let t = points.map(x => x[0]);\n let y = points.map(x => x[1]);\n y[2] += spread[i];\n let spline = NaturalCubicSpline(t, y)\n\n const bezpts = draw_hermite_spline(spline.knots, spline.values, spline.dvalues)\n .map((x, i) => !(i % 2) ? x*svgwidth : x*svgheight); // scale\n paths[i] = bezpts_to_svgpath(bezpts);\n }\n \n curve\n .attr(\"d\", (p, i) => paths[i]);\n point\n .attr(\"cx\", ([x,y]) => x * svgwidth)\n .attr(\"cy\", ([x,y]) => y * svgheight);\n colorpoints\n .attr(\"cx\", (k) => points[2][0] * svgwidth)\n .attr(\"cy\", (k) => (points[2][1] + k) * svgheight);\n\n }\n\n update(); \n \n return html`\n
\n ${svg.node()}\n
\n A natural cubic spline with linear extrapolation. Try moving the black points with the mouse.\n
\n
\n
`\n}","pinned":false,"mode":"js","data":null,"name":null},{"id":1122,"value":"Setup = md`\n## Problem setup\n\nImagine we start with a list of pairs ${$`(t_i, y_i)`}, with ${$`i`} ranging from ${$`0`} to ${$`n`} and ${$`t_i`} in increasing order so that ${$`t_i < t_{i+1}`}, and then we want to fit a smooth interpolating function ${$`f(t)`}\nthrough all of them, covering the interval ${$`t \\in [t_0, t_n]`}.\nThe ${$`t_i`} are called knots, and the ${$`y_i`} are called values.\n\nWe can build our interpolating function out of several pieces which [smoothly][] fit together (${$`C^2`} continuity), by first breaking our domain up into ${$`n`} subintervals ${$`[t_i, t_{i+1}]`}.\nOn each subinterval we can construct one cubic polynomial segment ${$`f_i(t)`}, under the interpolation constraint that ${$`f_i(t_{i}) = y_i`} and ${$`f_i(t_{i+1}) = y_{i+1}`} and also under the smoothness constraint that the first and second derivatives match across each knot: ${$`f_{i-1}'(t_i) = f_i'(t_i)`} and ${$`f_{i-1}''(t_i) = f_i''(t_i)`}.\n\n [smoothly]: https://en.wikipedia.org/wiki/Smoothness\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":162,"value":"Affine = md`\nFirst some preliminaries...\n\n## Affine complement\n\nI often find it convenient to write polynomials on the interval ${$`[0, 1]`} using what I call the [affine][] complement, ${$`\\bt = 1-t`}.\nNote that ${$`1 = \\bt + t`}, which makes it easy to elevate the degree of any term.\nIf we normalize our polynomials so that every term has the same total degree in terms of ${$`\\bt`} and ${$`t`}, the resulting polynomial basis is a rescaled version of the [Bernstein][] basis.\nThis form makes it very easy to evaluate polynomials at ${$`0`} or at ${$`1`}, because for the function ${$`f(t) = \\bt`}, we have ${$`f(0) = 1`}, ${$`f(1) = 0`}.\nThis means that any polynomial term ${$`\\bt^a`} evaluates to ${$`1`} at the left end of the interval, any polynomial term ${$`t^b`} evaluates to ${$`1`} at the right end, and any mixed term ${$`\\bt^at^b`} evaluates to ${$`0`} at both ends of the interval.\n\nThe derivative is ${$`\\frac{\\mathrm{d}}{\\mathrm{d}t}\\bt = –1`}, so we can easily differentiate any product of powers of ${$`t`} and ${$`\\bt`},\n\n${$$`\\tfrac{\\mathrm{d}}{\\mathrm{d}t}\\bt^{\\,m}t^n = n\\bt^{\\,m}t^{n-1} - m\\bt^{\\,m-1}t^n`}\n\n [affine]: https://en.wikipedia.org/wiki/Affine_combination\n [Bernstein]: https://en.wikipedia.org/wiki/Bernstein_polynomial\n\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":188,"value":"Hermite = md`\n## Hermite basis\n\nGiven any cubic polynomial ${$`p(t)`} defined on the interval, we can express it in terms of its values and first derivatives at the endpoints, ${$`\\{p_0, p_1, p'_0, p'_1\\}`}.\nThis is called the [Hermite][] basis.\n\n [Hermite]: https://en.wikipedia.org/wiki/Cubic_Hermite_spline\n\n${$$`\\begin{aligned}\np(t) &= p_0(\\bt^{\\,3} + 3\\bt^{\\,2}t) + p'_0(\\bt^{\\,2}t)\n ${(width < 500) ? `\\\\\\\\ &\\\\quad` : \"\"}\n + p_1(3\\bt t^2 + t^3) + p'_1(-\\bt t^2) \\\\[.2em]\n &= p_0\\bt^{\\,3} + (3p_0 + p'_0)\\bt^{\\,2}t\n ${(width < 500) ? `\\\\\\\\ &\\\\quad` : \"\"}\n + (3p_1 - p'_1)\\bt t^2 + p_1t^3\n\\end{aligned}`}\n\nIf we take derivatives of this polynomial we get:\n\n${responsive_katex(invalidation, $$`\n\\begin{aligned}\np'(t) &= p'_0\\bt^{\\,2} - (2p'_0 + 2p'_1 - 6(p_1 - p_0))\\bt t + p'_1t^2 \\\\[.2em]\n\np''(t) &= -(4p'_0 + 2p'_1 - 6(p_1 - p_0))\\bt\n ${(width < 500) ? '\\\\\\\\ &\\\\quad' : ''}\n + (2p'_0 + 4p'_1 - 6(p_1 - p_0))t \\\\[.2em]\n\np'''(t) &= 6(p'_0 + p'_1) -12(p_1 - p_0) \\\\[.2em]\n\\end{aligned}`)}\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":157,"value":"Slope = md`\n## Fitting together a spline\n\nRemember, we have a list of knots ${$`t_i`} and values ${$`y_i`} through which we want to fit a smooth interpolating function ${$`f(t)`}, built out of several cubic polynomial pieces ${$`f_i(t)`}, under constraints ${$`f_{i-1}'(t_i) = f_i'(t_i)`} and ${$`f_{i-1}''(t_i) = f_i''(t_i)`}.\n\nWe can write ${$`s_i = 1/(t_{i+1} - t_i)`} and ${$`r_i = (y_{i+1} - y_i)/(t_{i+1} - t_i)`} to simplify our equations below.\n\nSince we start already knowing the values at the knots, we only need to choose a first derivative at each knot, which we will call ${$`m_i = f_{i-1}'(t_i) = f_{i}'(t_i)`}, in order to fully specify the cubic polynomial ${$`f_{i}(t)`} on each subinterval.\n\nOur constraint on choosing these ${$`m_i`} is that we require the second derivatives to be equal across each knot ${$`f_{i-1}''(t_i) = f_{i}''(t_i)`}.\nWe can look at the second derivative of a polynomial in Hermite basis described above, remembering that when we re-scale the domain to a width of ${$`1/s_i`} from a width of ${$`1`}, we must scale first and second derivatives proportionally.\n\n${$$`\\begin{aligned}\n& (2s_{i-1}m_{i-1} + 4s_{i-1}m_{i}) - 6s_{i-1}r_{i-1}\n${(width < 550) ? `\\\\\\\\ & \\\\qquad` : \"\"}\n= - (4s_im_i + 2s_im_{i+1}) + 6s_ir_i\n\\end{aligned}`}\n\nMoving the ${$`m_i`} which we want to solve for to the left, and dividing everything by ${$`2`}, we get\n\n${$$`\\begin{aligned}\n& s_{i-1}m_{i-1} + (2s_{i-1} + 2s_i)m_i + s_im_{i+1}\n${(width < 550) ? `\\\\\\\\ & \\\\qquad` : \"\"}\n= 3s_{i-1}r_{i-1} + 3s_ir_i\n\\end{aligned}`}\n\nThis system of equations turns into a [tridiagonal matrix][tridiag], which can be solved in linear time using the Thomas algorithm. This algorithm is provably numerically stable in this case because the matrix is [diagonally dominant][].\n\n [tridiag]: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm\n [diagonally dominant]: https://en.wikipedia.org/wiki/Diagonally_dominant_matrix\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":419,"value":"Boundary = md`\n## Boundary conditions\n\nOne slight oversimplification above is that we have only treated intermediate knots, without discussing how to treat the endpoints.\nThere are a few common approaches.\nWe can treat the domain as [periodic][], in which case there are ${$`n+1`} subintervals, and the matrix is cyclic tridiagonal with matching numbers of variables and constraints:\n\n [periodic]: https://en.wikipedia.org/wiki/Periodic_function\n\n${responsive_katex(invalidation, $$`\n\\def\\dd{\\raisebox{-.4em}{$\\ddots$}}\n\n\\begin{aligned}\n&\\begin{pmatrix}\n 2s_{n} + 2s_0 & s_0 & & s_n \\\\\n s_0 & 2s_0 + 2s_1 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & s_{n-1} \\\\\n s_n & & s_{n-1} & 2s_{n-1} + 2s_n\n\\end{pmatrix}\n\\begin{pmatrix}\n m_0 \\\\\n m_1 \\\\\n \\vdots \\\\\n m_n\n\\end{pmatrix}\n\n${(width < 600) ? `\\\\\\\\ &\\\\qquad` : \"\"}\n= \n\n3\\begin{pmatrix}\n s_{n}r_{n} + s_0r_0 \\\\[.3em]\n s_{0}r_{0} + s_1r_1 \\\\\n \\vdots \\\\\n s_{n-1}r_{n-1} + s_nr_n\n\\end{pmatrix}\n\n\\end{aligned}`)}\n\n(In the matrices on this page, empty spaces indicate ${$`0`} entries.)\n\nOtherwise, we have ${$`2`} more variables than constraints, and we must choose an additional ${$`2`} constraints to apply. If we explicitly set the derivatives at the endpoints ${$`m_0`} and ${$`m_n`}, for instance setting them to ${$`0`}, then this is called a “clamped” spline.\nIf we instead force the third derivatives to match across the first interior knots ${$`f_0'''(t_1) = f_1'''(t_1)`} and ${$`f_{n-2}'''(t_{n-1}) = f_{n-1}'''(t_{n-1})`} then this is called the “not-a-knot” condition, because it forces the cubic polynomials to be the same, ${$`f_0(t) = f_1(t)`} and ${$`f_{n-2}(t) = f_{n-1}(t)`}.\nFinally if we set the second derivatives to zero at the endpoints ${$`f_0''(t_0) = 0`} and ${$`f_{n-1}''(t_{n}) = 0`} then this is called a “natural” spline:\n\n${responsive_katex(invalidation, $$`\n\\def\\dd{\\raisebox{-.4em}{$\\ddots$}}\n\n\\begin{aligned}\n&\\begin{pmatrix}\n 2s_0 & s_0 \\\\\n s_0 & 2s_0 + 2s_1 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & s_{n-1} \\\\\n & & s_{n-1} & 2s_{n-1}\n\\end{pmatrix}\n\\begin{pmatrix}\n m_0 \\\\\n m_1 \\\\\n \\vdots \\\\\n m_n\n\\end{pmatrix}\n${(width < 470) ? `\\\\\\\\ &\\\\qquad` : \"\"}\n= \n\n3\\begin{pmatrix}\n s_0r_0 \\\\[.3em]\n s_{0}r_{0} + s_1r_1 \\\\\n \\vdots \\\\\n s_{n-1}r_{n-1}\n\\end{pmatrix}\n\n\\end{aligned}`)}\n\nOne notable thing about the natural spline is that because the second derivative is zero at the far endpoints, it can be smoothly extrapolated by a linear function at each end.\n\nIf the knots are equally spaced, all of the ${$`s_i`} are the same and can be factored out and the equations become somewhat simpler.\nFor example the matrix equation for finding the slopes of the natural spline can be written:\n\n${responsive_katex(invalidation, $$`\n\\def\\dd{\\raisebox{-.4em}{$\\ddots$}}\n\n\\begin{aligned}\n&\\begin{pmatrix}\n 2 & 1 \\\\\n 1 & 4 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & 1 \\\\\n & & 1 & 4 & 1 \\\\\n & & & 1 & 2\n\\end{pmatrix}\n\\begin{pmatrix}\n m_0 \\\\\n m_1 \\\\\n \\vdots \\\\\n m_{n-1} \\\\\n m_n\n\\end{pmatrix}\n${(width < 430) ? `\\\\\\\\ &\\\\qquad` : \"\"}\n= \n\n\\frac{3}{t_1 - t_0}\\begin{pmatrix}\n y_1 - y_0 \\\\[.3em]\n y_2 - y_0 \\\\\n \\vdots \\\\\n y_n - y_{n-2} \\\\\n y_n - y_{n-1}\n\\end{pmatrix}\n\n\\end{aligned}`)}\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":464,"value":"Tridiagonal = md`\n## Solving a tridiagonal system\n\nGiven a matrix equation ${$`Ax = b`}, for symmetric tridiagonal matrix ${$`A`}, we can solve for ${$`x`} by performing two passes of Gaussian elimination, conceptually factoring ${$`A`} into a product of a lower-triangular and upper-triangular matrix.\nThis is called the Thomas algorithm.\nTo reduce the number of division operations required, we can use a variant algorithm (cf. Quarteroni, Sacco and Saleri (2007) [*Numerical Mathematics*][qss], pp. 93–96) based on a factorization ${$`A = LDL^T`}:\n\n [qss]: https://www.springer.com/us/book/9783540346586\n\n${responsive_katex(invalidation, $$`\n\\def\\dd{\\raisebox{-.4em}{$\\ddots$}}\n\n\\begin{aligned}\n\n& \\begin{pmatrix}\n 2s_0 & s_0 \\\\\n s_0 & 2s_0 + 2s_1 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & s_{n-1} \\\\\n & & s_{n-1} & 2s_{n-1}\n\\end{pmatrix}\n\\\\[3.5em]\n= &\n\n\\begin{pmatrix}\n d_0^{-1} \\\\\n s_0 & d_1^{-1} \\\\[.8em]\n & \\ddots & \\ddots \\\\\n & & s_{n-1} & d_n^{-1}\n\\end{pmatrix}\n\\begin{pmatrix}\n d_0 \\\\\n & d_1 \\\\[.8em]\n & & \\ddots \\\\\n & & & d_n\n\\end{pmatrix}\n\\begin{pmatrix}\n d_0^{-1} & s_0 \\\\\n & d_1^{-1} & \\dd \\\\[.8em]\n & & \\ddots & s_{n-1} \\\\\n & & & d_n^{-1}\n\\end{pmatrix}\n\n\\end{aligned}\n`)}\n\nThe ${$`d_i`} can be found recursively starting from the top, with ${$`1`} division, ${$`3`} multiplications, and ${$`2`} additions per row.\n\nThen to solve ${$`Ax = b`}, we can first solve for ${$`y`} in ${$`Ly = b`}, and then solve for ${$`x`} in ${$`Ux = y`}, where ${$`U = DL^T`}.\nEach of these is also a simple recurrence, first forward through the rows, then backward.\nEach step requires only ${$`2`} multiplications and ${$`1`} addition.\n\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":547,"value":"Cyclic = md`\n#### Cyclic tridiagonal system\n\nIn the cyclic case ${$`Cx = b`}, we can first break the matrix down into the sum of a tridiagonal matrix and a [rank][] ${$`1`} correction, ${$`C = A + uv^T`}, and then find ${$`x`} using a method based on the [Sherman–Morrison formula][sm].\n\n [rank]: https://en.wikipedia.org/wiki/Rank_(linear_algebra)\n [sm]: https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula\n\n${responsive_katex(invalidation, $$`\n\\def\\dd{\\raisebox{-.4em}{$\\ddots$}}\n\n\\begin{aligned}\n\n&\n\\begin{pmatrix}\n 2s_{n} + 2s_0 & s_0 & & s_n \\\\\n s_0 & 2s_0 + 2s_1 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & s_{n-1} \\\\\n s_n & & s_{n-1} & 2s_{n-1} + 2s_n\n\\end{pmatrix}\n\\\\ = &\n\n\\begin{pmatrix}\n s_{n} + 2s_0 & s_0 \\\\\n s_0 & 2s_0 + 2s_1 & \\dd \\\\[.8em]\n & \\ddots & \\ddots & s_{n-1} \\\\\n & & s_{n-1} & 2s_{n-1} + s_n\n\\end{pmatrix}\n+\n\\begin{pmatrix}\n s_n \\\\\n \\\\[.8em]\n \\\\\n s_n\n\\end{pmatrix}\n\\begin{pmatrix}\n 1 & & & 1\n\\end{pmatrix}\n\n\\end{aligned}`)}\n\nWe can find ${$`x`} by solving ${$`Aw = b`} and ${$`Az = u`}, then computing ${$`x = w - \\dfrac{v^Tw}{1 + v^Tz}z`}.\n\nThis works because:\n\n${responsive_katex(invalidation, $$`\n\\begin{aligned}\nCx &= \\left(A + u v^T \\right)\\left(w - \\dfrac{v^T w}{1 + v^T z}z\\right)\\\\[.8em]\n &= Aw + (uv^T)w\n ${(width < 380) ? `\\\\\\\\ &\\\\quad` : \"\"}\n - \\left( Az + (uv^T)z \\right)\\dfrac{ v^T w}{1 + v^T z}\\\\[.8em]\n &= b + uv^Tw - \\dfrac{u\\left( 1 + v^Tz \\right) v^T w}{1 + v^T z}\\\\\n &= b\n\\end{aligned}`)}\n\nWe can use the same factorization of ${$`A`} for solving both of these matrix equations ${$`Aw = b`} and ${$`Az = u`}, making this cyclic algorithm nearly as fast as the non-cyclic version.\n\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":1328,"value":"Parametric = md`\n## Splines as parametric curves\n 760 ? -220 : Math.min(-70, 540-width)}px;\">\nThe discussion above is all about functions ${$`f : \\mathbb{R} \\to \\mathbb{R}`}, but the word [*spline*][spline] originally refers to a thin strip of wood or metal used to design planar curves, especially for shipbuilding.\n\nIn the computer graphics and computer-aided design literature, there are many ways of constructing piecewise parametric curves, often called “splines”, through a sequence of control points, including by [modeling an idealization of a metal strip](https://levien.com/phd/phd.html). Most of these methods have nothing to do with the cubic spline described earlier in this paper, but it is also possible to use our cubic spline to make a piecewise parametric curve through an ordered list of points ${$`\\mathbf{p}_i`} in the affine plane or in space.\n\nThis is done by choosing a corresponding parameter value ${$`t_i`} for each point; representing the point using some coordinate system, e.g. ${$`(p_{x, i}, p_{y, i})`}; and then finding an interpolating cubic spline (in the sense of this paper) for each coordinate in terms of the parameter, ${$`\\left(f_{x}(t), f_{y}(t)\\right)`}. The simplest choice of parameter values for the control points is just counting numbers ${$`[t_i] = [0, 1, 2, 3, \\dots]`}, but sometimes parameter values are chosen so that their intervals are proportional to the distance between points in the Euclidean plane, ${$`t_{i+1} - t_i = \\left\\|\\mathbf{p}_{i+1} - \\mathbf{p}_i\\right\\|`}, called a chord-length parametrization.\n\n [spline]: https://en.wikipedia.org/wiki/Flat_spline\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":1130,"value":"Code = md`## Code`","pinned":false,"mode":"js","data":null,"name":null},{"id":863,"value":"// Takes in two arrays of knots and values, and returns a\n// Float64Array of the slope at each knot.\nnatural_cubic_spline = function natural_cubic_spline(knots, values) {\n const\n t = knots, y = values, n = t.length - 1,\n s = new Float64Array(n+1), // horizontal scales\n d = new Float64Array(n+1), // diagonal matrix\n m = new Float64Array(n+1); // result vector: slope at each knot\n for (let i = 0; i < n; i++) {\n s[i] = 1 / (t[i+1] - t[i]); // Note s[n] == 0\n // Setup pass; at this point 'm' represents the right-hand side\n m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }\n d[0] = 0.5 / s[0];\n m[0] = d[0] * m[0];\n for (let i = 1; i <= n; i++) {\n d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);\n // Forward pass; at this point 'm' represents Um.\n m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }\n for (let i = n-1; i >= 0; i--) {\n m[i] -= d[i]*s[i]*m[i+1]; }\n return m;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":815,"value":"// Takes in an array of values at knots t = [0, 1, 2, ..., n]\n// and returns a Float64Array of the slope at each knot.\nnatural_cubic_spline_equi = function natural_cubic_spline_equi(values) {\n const\n y = values, n = y.length - 1,\n d = new Float64Array(n+1), // diagonal matrix\n m = new Float64Array(n+1); // result vector: slope at each knot\n for (let i = 1; i < n; i++)\n m[i] = 3*(y[i+1] - y[i-1]);\n d[0] = 0.5, m[0] = 1.5*(y[1] - y[0]), m[n] = 3*(y[n] - y[n-1]);\n for (let i = 1; i < n; i++)\n d[i] = 1 / (4 - d[i-1]), m[i] = (m[i] - m[i-1]) * d[i];\n d[n] = 1 / (2 - d[n-1]), m[n] = (m[n] - m[n-1]) * d[n];\n for (let i = n-1; i >= 0; i--)\n m[i] -= d[i]*m[i+1];\n return m;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":725,"value":"// Takes in two arrays of knots and values, and returns a\n// Float64Array of the slope at each knot.\n// Assumes that the knots are in the periodic interval [0, 1].\ncyclic_cubic_spline = function cyclic_cubic_spline(knots, values) {\n const\n t = knots, y = values, n = t.length - 1,\n s = new Float64Array(n+1), // horizontal scales\n d = new Float64Array(n+1), // diagonal matrix\n z = new Float64Array(n+1), // correction vector\n m = new Float64Array(n+1); // result vector: slope at each knot\n for (let i = 0; i < n; i++) {\n s[i] = 1 / (t[i+1] - t[i]);\n m[i] += (m[i+1] += 3*s[i]*s[i]*(y[i+1] - y[i])); }\n s[n] = 1 / (t[0] + 1 - t[n]);\n const bn = 3*s[n]*s[n]*(y[0] - y[n]);\n m[0] += bn; m[n] += bn;\n d[0] = 1 / (2*s[0] + s[n]);\n m[0] = d[0] * m[0];\n z[0] = d[0] * s[n];\n for (let i = 1; i < n; i++) {\n d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);\n m[i] = d[i]*(m[i] - s[i-1]*m[i-1]);\n z[i] = d[i]*(0 - s[i-1]*z[i-1]); }\n d[n] = 1 / (s[n] + 2*s[n-1] - s[n-1]*s[n-1]*d[n-1]);\n m[n] = d[n]*(m[n] - s[n-1]*m[n-1]);\n z[n] = d[n]*(s[n] - s[n-1]*z[n-1]);\n for (let i = n-1; i >= 0; i--) {\n m[i] -= d[i]*s[i]*m[i+1];\n z[i] -= d[i]*s[i]*z[i+1]; }\n const zscale = (m[0] + m[n]) / (1 + z[0] + z[n]);\n for (let i = 0; i <= n; i++) {\n m[i] -= zscale * z[i]; }\n return m;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":777,"value":"// Takes a list of knots, values, and derivatives, and constructs\n// a cubic polynomial in Bernstein–Bézier basis on each subinterval,\n// output as a Float64Array.\ndraw_hermite_spline = function draw_hermite_spline (knots, values, dvalues) {\n const t = knots, y = values, m = dvalues, n = t.length - 1;\n const bezpts = new Float64Array(8*n);\n for (let i = 0; i < n; i++) {\n const h = (t[i+1] - t[i]) * 0.3333333333333333;\n bezpts.set([\n t[i], y[i], t[i] + h, y[i] + m[i]*h,\n t[i+1] - h, y[i+1] - m[i+1]*h, t[i+1], y[i+1]\n ], 8*i); }\n return bezpts;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":1518,"value":"herm_to_cheb3 = function herm_to_cheb3(y0, y1, m0, m1) {\n return [\n 0.5*(y0 + y1 + 0.125*(m0 - m1)),\n 0.03125*(18*(y1 - y0) - m0 - m1 ),\n 0.0625*(m1 - m0),\n 0.03125*(2*(y0 - y1) + m0 + m1)];\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":4,"value":"// Given a list of knots and a target point, finds the subinterval\n// in which the target point falls. Points outside the interval are \n// considered to be in the first or last interval.\nbinary_search = function binary_search(knots, target) {\n let n = knots.length - 1, low = 0, half;\n while ((half = (n/2)|0) > 0) {\n low += half * (knots[low + half] <= target);\n n -= half; }\n return low;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":1720,"value":"md`\nWe can make an entirely self-contained version with no dependencies:\n`","pinned":false,"mode":"js","data":null,"name":null},{"id":1575,"value":"// Given an array of knots and an array of values,\n// return a function which evaluates a natural cubic spline\n// extrapolated by linear functions on either end.\nNaturalCubicSpline = function NaturalCubicSpline(knots, values) {\n const\n n = knots.length - 1,\n t = new Float64Array(n+3), // knots\n y = new Float64Array(n+3), // values\n s = new Float64Array(n+2), // horizontal scales\n m = new Float64Array(n+3), // slope at each knot\n d = new Float64Array(n+2), // diagonal matrix\n coeffs = new Float64Array(4*(n+2)); // Chebyshev coefficients\n t.set(knots, 1), y.set(values, 1);\n \n // Natural cubic spline algorithm\n for (let i = 1; i < n+1; i++) {\n s[i] = 1 / (t[i+1] - t[i]);\n m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }\n d[1] = 0.5 / s[1];\n m[1] = d[1] * m[1];\n for (let i = 2; i <= n+1; i++) {\n d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);\n m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }\n for (let i = n; i >= 1; i--) {\n m[i] -= d[i]*s[i]*m[i+1]; }\n\n // Linear extrapolation\n t[0] = t[1] - 1; t[n+2] = t[n+1] + 1;\n y[0] = y[1] - m[1]; y[n+2] = y[n+1] + m[n+1];\n s[0] = s[n+1] = 1;\n m[0] = m[1]; m[n+2] = m[n+1];\n \n // Set up Chebyshev coefficients\n for (let i = 0; i < n+2; i++) {\n const h = t[i+1] - t[i];\n const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;\n coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));\n coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);\n coeffs[j+2] = 0.0625*(m1 - m0);\n coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }\n\n return Object.assign(\n function eval_spline(u) {\n // Binary search\n let nn = t.length - 1, low = 0, half;\n while ((half = (nn/2)|0) > 0) {\n low += half * (t[low + half] <= u);\n nn -= half; }\n const i = low, j = 4*i;\n u = (2*u - t[i] - t[i+1]) * s[i]; // scale to [–1, 1]\n // Clenshaw's method.\n let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;\n b1 = coeffs[j+1] + 2 * u * b2 - b1;\n return coeffs[j] + u * b1 - b2; },\n { knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":1829,"value":"// Takes in two arrays of knots and values, and returns a\n// return a function which evaluates a cyclic cubic spline.\n// Assumes that the knots are in the periodic interval [0, 1],\n// or optionally some other interval.\nCyclicCubicSpline = function CyclicCubicSpline(knots, values, domain=[0,1]) {\n const\n n = knots.length - 1,\n scale = 1 / (domain[1] - domain[0]),\n t = new Float64Array(n+3), // knots\n y = new Float64Array(n+3), // values\n s = new Float64Array(n+2), // horizontal scales\n m = new Float64Array(n+3), // slope at each knot\n d = new Float64Array(n+2), // diagonal matrix\n z = new Float64Array(n+2), // correction vector\n coeffs = new Float64Array(4*(n+2)); // Chebyshev coefficients\n y.set(values, 1);\n \n for (let i = n; i >= 0; i--) {\n t[i+1] = (knots[i] - domain[0]) * scale; }\n t[0] = t[n+1] - 1; t[n+2] = t[1] + 1;\n y[0] = y[n+1]; y[n+2] = y[1];\n \n for (let i = 1; i < n+1; i++) {\n s[i] = 1 / (t[i+1] - t[i]);\n m[i] += (m[i+1] += 3*s[i]*s[i]*(y[i+1] - y[i])); }\n s[n+1] = 1 / (t[1] + 1 - t[n+1]);\n const bn = 3*s[n+1]*s[n+1]*(y[1] - y[n+1]);\n m[1] += bn; m[n+1] += bn;\n d[1] = 1 / (2*s[1] + s[n+1]);\n m[1] = d[1] * m[1];\n z[1] = d[1] * s[n+1];\n for (let i = 2; i < n+1; i++) {\n d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);\n m[i] = d[i]*(m[i] - s[i-1]*m[i-1]);\n z[i] = d[i]*(0 - s[i-1]*z[i-1]); }\n d[n+1] = 1 / (s[n+1] + 2*s[n] - s[n]*s[n]*d[n]);\n m[n+1] = d[n+1]*(m[n+1] - s[n]*m[n]);\n z[n+1] = d[n+1]*(s[n+1] - s[n]*z[n]);\n for (let i = n; i >= 1; i--) {\n m[i] -= d[i]*s[i]*m[i+1];\n z[i] -= d[i]*s[i]*z[i+1]; }\n const zscale = (m[1] + m[n+1]) / (1 + z[1] + z[n+1]);\n for (let i = 1; i <= n+1; i++) {\n m[i] -= zscale * z[i]; }\n \n s[0] = s[n+1] = 1 / (t[1] - t[0]);\n m[0] = m[n+1]; m[n+2] = m[1];\n \n // Set up Chebyshev coefficients\n for (let i = 0; i < n+2; i++) {\n const h = t[i+1] - t[i];\n const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;\n coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));\n coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);\n coeffs[j+2] = 0.0625*(m1 - m0);\n coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }\n\n return Object.assign(\n function eval_spline(u) {\n u = ((u - domain[0]) * scale) % 1 // scale interval to [0, 1]\n // Binary search\n let nn = t.length - 1, low = 0, half;\n while ((half = (nn/2)|0) > 0) {\n low += half * (t[low + half] <= u);\n nn -= half; }\n const i = low, j = 4*i;\n u = (2*u - t[i] - t[i+1]) * s[i]; // scale subinterval to [–1, 1]\n // Clenshaw's method.\n let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;\n b1 = coeffs[j+1] + 2 * u * b2 - b1;\n return coeffs[j] + u * b1 - b2; },\n { knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":1923,"value":"// Takes in two arrays of knots and values, and returns a\n// Float64Array of cubic Bézier segments\nparametric_NCS = function parametric_NCS(points, chordal=true) {\n const\n n = points.length - 1,\n y0 = new Float64Array(n+1),\n y1 = new Float64Array(n+1),\n h = new Float64Array(n), // chord lengths\n s = new Float64Array(n+1), // scale = 1/h\n d = new Float64Array(n+1), // diagonal matrix\n m0 = new Float64Array(n+1), // result vectors: parametric\n m1 = new Float64Array(n+1), // slope at each knot \n bezpts = new Float64Array(8*n);\n for (let i = 0; i <= n; i++) {\n [y0[i], y1[i]] = points[i]; }\n \n if (chordal) {\n for (let i = 0; i < n; i++) {\n const dy0 = y0[i+1] - y0[i];\n const dy1 = y1[i+1] - y1[i];\n h[i] = Math.sqrt(dy0*dy0 + dy1*dy1)\n s[i] = 1 / h[i]; } } // Note s[n] == 0\n else {\n h.fill(1, 0, n);\n s.fill(1, 0, n); }\n \n for (let i = 0; i < n; i++) {\n // Setup pass; at this point 'm' represents the right-hand side\n m0[i] += (m0[i+1] = 3*s[i]*s[i]*(y0[i+1] - y0[i])); \n m1[i] += (m1[i+1] = 3*s[i]*s[i]*(y1[i+1] - y1[i])); }\n d[0] = 0.5 / s[0];\n m0[0] = d[0] * m0[0];\n m1[0] = d[0] * m1[0];\n for (let i = 1; i <= n; i++) {\n d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);\n // Forward pass; at this point 'm' represents Um.\n m0[i] = d[i]*(m0[i] - s[i-1]*m0[i-1]); \n m1[i] = d[i]*(m1[i] - s[i-1]*m1[i-1]); }\n for (let i = n-1; i >= 0; i--) {\n m0[i] -= d[i]*s[i]*m0[i+1];\n m1[i] -= d[i]*s[i]*m1[i+1]; }\n \n for (let i = 0; i < n; i++) {\n const hi = h[i] * 0.3333333333333333;\n bezpts.set([\n y0[i], y1[i], y0[i] + hi*m0[i], y1[i] + hi*m1[i],\n y0[i+1] - hi*m0[i+1], y1[i+1] - hi*m1[i+1], y0[i+1], y1[i+1],\n ], 8*i); }\n return bezpts;\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":939,"value":"d3 = require(\"d3@5\")","pinned":true,"mode":"js","data":null,"name":null},{"id":687,"value":"import {bezeval, bezpts_to_svgpath} from '@jrus/bezplot'","pinned":true,"mode":"js","data":null,"name":null},{"id":1530,"value":"import {evaluate as chebeval} from '@jrus/cheb'","pinned":true,"mode":"js","data":null,"name":null},{"id":57,"value":"import {responsive_katex, math_css, $ as $tex} with {$css as $css} from \"@jrus/misc\"","pinned":true,"mode":"js","data":null,"name":null},{"id":54,"value":"$css = html`${math_css}`","pinned":false,"mode":"js","data":null,"name":null},{"id":1501,"value":"$ = $tex.options({\n macros: { \"\\\\bt\": \"\\\\bar{t}\" }\n});","pinned":true,"mode":"js","data":null,"name":null},{"id":1509,"value":"$$ = $tex.options({\n displayMode: true,\n macros: { \"\\\\bt\": \"\\\\bar{t}\" }\n});","pinned":true,"mode":"js","data":null,"name":null},{"id":884,"value":"// version from D3\nfunction controlPoints(x) {\n var i, n = x.length - 1, m, a = new Array(n), b = new Array(n), r = new Array(n);\n a[0] = 0, b[0] = 2, r[0] = x[0] + 2 * x[1];\n for (i = 1; i < n - 1; i++) a[i] = 1, b[i] = 4, r[i] = 4 * x[i] + 2 * x[i + 1];\n a[n - 1] = 2, b[n - 1] = 7, r[n - 1] = 8 * x[n - 1] + x[n];\n for (i = 1; i < n; ++i) m = a[i] / b[i - 1], b[i] -= m, r[i] -= m * r[i - 1];\n a[n - 1] = r[n - 1] / b[n - 1];\n for (i = n - 2; i >= 0; --i) a[i] = (r[i] - a[i + 1]) / b[i];\n b[n - 1] = (x[n] + a[n - 1]) / 2;\n for (i = 0; i < n - 1; ++i) b[i] = 2 * x[i + 1] - a[i + 1];\n return [a, b];\n}","pinned":true,"mode":"js","data":null,"name":null},{"id":888,"value":"// Alternative to D3 version\ncontrol_points = function control_points(y) {\n var i, n = y.length - 1, d = new Array(n+1), m = new Array(n+1);\n for (i = 1; i < n; i++) m[i] = (y[i+1] - y[i-1]);\n m[0] = (d[0] = .5)*(y[1] - y[0]), m[n] = (y[n] - y[n-1]);\n for (i = 1; i < n; i++) m[i] = (m[i] - m[i-1]) * (d[i] = 1 / (4 - d[i-1]));\n m[n] = (m[n] - m[n-1]) * (d[n] = 1 / (2 - d[n-1]));\n for (i = n-1; i >= 0; i--) m[i] -= d[i]*m[i+1];\n for (i = 0; i < n; i++) d[i] = y[i+1] - m[i+1], m[i] = y[i] + m[i];\n d.pop(), m.pop();\n return [m, d];\n}","pinned":true,"mode":"js","data":null,"name":null}],"resolutions":[],"schedule":null,"last_view_time":null}