https://cohost.org/tomforsyth/post/982199-polynomial-interpola [noscript] log insign up tomforsyth Tom Forsyth @tomforsyth GPU designer and gfx coder * he/him * eelpi.gotdns.org/ Gfx coder and chip designer. Worked at Oculus, Valve, RAD, Muckyfoot, 3Dlabs, now back at Intel. Blade2, Larrabee, TF2, VR and many other atrocities. log in tomforsyth Tom Forsyth@tomforsyth2/8/2023, 5:13 AM --------------------------------------------------------------------- Polynomial interpolation [Polynomila] Context: this is part of a series of reposts of some of my more popular blog posts. It was originally posted in June 2013 here http:/ /tomforsyth1000.github.io/blog.wiki.html# %5B%5BPolynomial%20interpolation%5D%5D. This version has been edited slightly to update for the intervening years. I was reminded of this one recently when I saw this video by Stand Up Maths [https://www.youtube.com/watch?v=F_43oTnTXiw] / Matt Parker [ https://mastodon.gamedev.place/@[email protected]/109819371402062543 ]. He uses some heavy tools, but I remembered that you don't need a maths degree to figure out an approximating polynomial for any function, and indeed it's fast enough you can do it in a game at runtime - it doesn't need to be precalculated - so it's enormously useful for all sorts of game-related things like graphics and gameplay. Anyway, on to the post. In work recently [2013, PC VR at Oculus] I've been dealing with fitting polynomial curves to set of points. Normally for this job you'd reach for Mathematica or some other heavy-math package and use the right incantations. It would do Magic Math Stuff and spit out the right answer. The problem if you didn't really understand the question, and so you don't really understand the solution, and certainly it's hard to do things like implement at runtime in your own code, or explain why it doesn't work in some cases. So I do like to at least try to do it myself with my high-school math skills before reaching for the Math Wand. --------------------------------------------------------------------- The question here is - I have a function, and a bunch of data points, and they're not that badly-controlled, but I do need to approximate them with a nice simple polynomial I can put in something like a pixel shader. To keep it concrete, say it's a cubic polynomial in this example. I'm quite happy to stick down a few control points myself and do the noise-smoothing approximation process by eye (humans are very good at this!), but how does that help me find the right numbers to put in the shader code? In a shader, a cubic polynomial looks like this: y = A + B*x + C*x^2 + D*x^3; Actually, if you want the ultimate efficiency, you should write it this way to express its fundamental "fused multiply-add" nature. To us they're the same, but to IEEE754 rules they're not EXACTLY the same, so it sometime pays to write it out very specifically. y = A + x*( B + x*( C + x*( D ) ) ); OK... but that doesn't help me. I have four points (p1,q1) (p2,q2) (p3,q3) (p4,q4) that I want this curve to pass through. How do I find A, B, C and D? Intuition says there should be a solution - we have four constraints (the four values of q) and four constants (A, B, C, D). There should be a solution! I have seen many people throw their hands up at this stage of problems like it and say "I'm not a mathematician - this is too hard for me". Except it's not - it's not actually difficult maths. And in fact there's plenty of mathematicians that will make a mess of this too, precisely because they know too much maths - so they'll start to use a higher-order sledgehammer to crack this fairly simple nut and soon you'll be #include-ing all sorts of libraries with legally intimidating licenses. So now I've let the cat out of the bag and told you it doesn't require anything resembling graduate maths, let's smash our foreheads into it like the dirty stubborn hackers we are. The truly stubborn way to approach this is to write out four equations and use high-school math to solve them simultaneously by adding and subtracting stuff, so: q1 = A + p1*( B + p1*( C + p1*( D ) ) ); q2 = A + p2*( B + p2*( C + p2*( D ) ) ); q3 = A + p3*( B + p3*( C + p3*( D ) ) ); q4 = A + p4*( B + p4*( C + p4*( D ) ) ); So we can... um... multiply everything by... er... oh that's all going to be mess I think. Maybe stubborn isn't all you need. This looks like the wrong path. Let's see if we can solve some simpler cases - see if that helps us. Well... what if p1 was zero? Oh, then life gets much simpler, because from the first equation: q1 = A + p1*( B + p1*( C + p1*( D ) ) ); = A; Well that was easy! Except p1 probably isn't zero. But what if we reformatted the polynomial so it was all relative to zero, i.e. rewrite it: y = A + (x-p1)*( B + (x-p1)*( C + (x-p1)*( D ) ) ); So by using (x-p1) instead of x we can make all those extra terms vanish at x==p1. But this isn't the same A as before because we changed things. We could probably multiply it all out and find the old A from the new A, once we've found the new B, C and D. And we haven't found them yet. This is still looking ugly. Though - that's a neat trick multiplying things you don't like by (x-p1). What if we did it for all the other P values? What if we wrote the polynomial like this: y = f1 * (x-p2) * (x-p3) * (x-p4) + f2 * (x-p1) * (x-p3) * (x-p4) + f3 * (x-p1) * (x-p2) * (x-p4) + f4 * (x-p1) * (x-p2) * (x-p3) ; Clearly this is still a cubic polynomial, and I can theoretically convert it back to the one I want to put in the shader as long as I find out f1-f4. The nice thing is that when you set x=p1 and y=q1, the second third and fourth lines are zero and you just get a fairly simple thing to solve: q1 = f1 * (p1-p2) * (p1-p3) * (p1-p4); Therefore: f1 = q1 / ( (p1-p2) * (p1-p3) * (p1-p4) ); Best of all, it's symmetrical - I don't have to really figure out the other four lines, they're just obvious extensions of the first: f2 = q2 / ( (p2-p1) * (p2-p3) * (p2-p4) ); f3 = q3 / ( (p3-p1) * (p3-p2) * (p3-p4) ); f4 = q4 / ( (p4-p1) * (p4-p2) * (p4-p3) ); Symmetry often hints that we're on the right tracks. OK, so now I have f1-f4, how do I get A,B,C,D? Well, now we do have to multiply out and gather terms together for constants, x, x^2, x^3. But it's really not that hard, and you can still use lots of symmetry to help yourself. Taking each line of the new polynomial separately: f1 * (x-p2) * (x-p3) * (x-p4) = -f1*p2*p3*p4 + x * ( f1*(p2*p3+p3*p4+p4*p2) ) + x^2 * (-f1*(p2+p3+p4)) + x^3 * (f1) f2 * (x-p1) * (x-p3) * (x-p4) = -f2*p1*p3*p4 + x * ( f2*(p1*p3+p3*p4+p4*p1) ) + x^2 * (-f2*(p1+p3+p4)) + x^3 * (f2) f3 * (x-p1) * (x-p2) * (x-p4) = -f3*p1*p2*p4 + x * ( f3*(p1*p2+p2*p4+p4*p1) ) + x^2 * (-f3*(p1+p2+p4)) + x^3 * (f3) f4 * (x-p1) * (x-p2) * (x-p3) = -f4*p1*p2*p3 + x * ( f4*(p1*p2+p2*p3+p3*p1) ) + x^2 * (-f4*(p1+p2+p3)) + x^3 * (f4) ...and then gather vertically by powers of x to get: A = -( f1*p2*p3*p4 + f2*p1*p3*p4 + f3*p1*p2*p4 + f4*p1*p2*p3 ); B = f1*(p2*p3+p3*p4+p4*p2) + f2*(p1*p3+p3*p4+p4*p1) + f3*(p1*p2+p2*p4+p4*p1) + f4*(p1*p2+p2*p3+p3*p1); C = -( f1*(p2+p3+p4) + f2*(p1+p3+p4) + f3*(p1+p2+p4) + f4*(p1+p2+p3) ); D = f1 + f2 + f3 +f4; Hey - that's actually pretty elegant! And of course if you want to get perf-crazy you can gather some terms together and so on, but sometimes it really is fine to just let the compiler do that for you. There's always the danger that you fail to Finish Your Derviations, but in this case I don't see one (there is an obviously elegant way to express all this symbolically, but it's a bit dodgy because when you actually go to evaluate it you end up dividing by zero unless you pre-cancel terms, and if you allow those sort of shenanigans you can accidentally prove that 1==2) Also, what are the failure points? Well, the obvious ones are the divides when calculating f1-f4. They fail when any of the P values are the same as any of the other P values, e.g. p1==p2, and you divide by zero and get infinities all over the place. That's a pretty reasonable failure - you can't have a simple polynomial having two different values at the same input value. There's also clear numerical instability if two P values are close to each other - the polynomial is going to go pretty bonkers - it might actually go through the points you want, but its behaviour between them is going to go up and down all over the place. Other than that, it seems pretty well-behaved. Now of course the mathematically-inclined will point out that these are simply the products of some matrix multiplies or start (ab)using Kronecker deltas. That's all fine, but that's not the point of my post. Of course you folks can solve it, and many harder problems as well - that wasn't in doubt. The point is that for simple things like this, hackers like me don't have to go asking the maths gods (or their programs) for help - we can solve this on our own by applying some basic problem-solving skills. The key is knowing just how far your maths knowledge will spread, and the trick is to develop the intuition to know, as I did, that this was well-constrained enough that there probably was a simple solution, and to at least explore towards a solution without prematurely invoking post-high-school maths. This intuition can often be wrong - there are some really simple-seeming problems that are indeed evil bastards that need heavy number-crunching. But in this case the solution is so fast and so elegant we could perform it millions of times a frame, and that opens up all sorts of interesting applications. Whereas if you handed it off to Maple, you'd of course get the right answer, but it would be far slower, and when it failed it would not be particularly obvious why it failed or what to do about it. #maths#math#polynomial#interpolation#gamedev --------------------------------------------------------------------- 0 comments You must log in to comment. * (c) 2023 anti software software club llc * thanks for using cohost Legal * Terms of Use * Privacy Notice * Community Guidelines About * @staff * SupportCredits * cohost status * cohost on twitter * ASSC on twitter * Careers