jheriko

24th July 2003 14:04 UTC

General Cubic Solution
This is one for all the techies out there, a general cubic solution in avs code. I've tested it out by raytracing a folium of descartes in xy stretched out along the z axis. The result was (surprisingly) the correct shape with a reflection of itself (due to the cubic term this happens for the same reason as with a plane) the 'singular line' in the shape created a *really* ugly edge though so if you plan to make cubic DMs I would suggest avoiding shapes with singular points or lines. Here is the code from my DM that deals with solving the cubic (cx^3+dx^2+ex+f=0):

(init: i27=1/27;i3=1/3;twopi=acos(-1)*2;)

a=(3*c*e-sqr(d))/(3*sqr(c));
b=(2*sqr(d)*d-9*c*d*e+27*sqr(c)*f)/(27*sqr(c)*c);
tmp=sqr(b)+(4*a*sqr(a))*i27;
p3r=0.5*(-b+if(above(tmp,0),sqrt(tmp),0));
p3i=0.5*if(below(tmp,0),sqrt(tmp),0);
q3r=0.5*(b+if(above(tmp,0),sqrt(tmp),0));
q3i=0.5*if(below(tmp,0),sqrt(tmp),0);
pk=sqrt(sqr(p3r)+sqr(p3i));
qk=sqrt(sqr(q3r)+sqr(q3i));
p0=atan2(p3r,p3i);
q0=atan2(q3r,q3i);
pr1=pow(pk,i3)*cos(p0*twopi*i3);
pi1=pow(pk,i3)*sin(p0*twopi*i3);
pr2=pow(pk,i3)*cos(p0*2*twopi*i3);
pi2=pow(pk,i3)*sin(p0*2*twopi*i3);
pr3=pow(pk,i3)*cos(p0*3*twopi*i3);
pi3=pow(pk,i3)*sin(p0*3*twopi*i3);
qr1=pow(qk,i3)*cos(q0*twopi*i3);
qi1=pow(qk,i3)*sin(q0*twopi*i3);
qr2=pow(qk,i3)*cos(q0*2*twopi*i3);
qi2=pow(qk,i3)*sin(q0*2*twopi*i3);
qr3=pow(qk,i3)*cos(q0*3*twopi*i3);
qi3=pow(qk,i3)*sin(q0*3*twopi*i3);
t1=if(above(abs(pi1+qi1),0.01),10000,pr1+qr1);
t2=if(above(abs(pi2+qi2),0.01),10000,pr2+qr2);
t3=if(above(abs(pi3+qi3),0.01),10000,pr3+qr3);
t1=if(below(t1,0),10000,t1);
t2=if(below(t2,0),10000,t2);
t3=if(below(t3,0),10000,t3);
t=min(min(t1,t2),t3);

I haven't really optimised the code that much so that its 'easier' to see whats going on and so that I can go over a couple of useful things. Like the a and b variables - if you have a cubic x^3+px+q=0 then a and b in the code can be set a=p,b=q rather than setting c=1,d=0,e=p,f=q. The reason for this is that the first stage of solving a cubic is reducing it to this form.

Unlike when solving a quadratic we can't just ignore complex results from the square roots and cube roots in the intermediary stages since they may end up forming a purely real root at the end. As a result this code also contains a method for finding a cube root of a complex number.

prn=pow(pk,1/3)*cos(p0*n*twopi/3);
pin=pow(pk,1/3)*sin(p0*n*twopi/3);

If you replace 3 with another number and increment n from 1 to that number then you can find any root (4th root, 5th root etc..) of any complex number, not just that but you get a different root for each value of n which means that you can calculate all of the roots and not just one (real root) like sqrt(x) or pow(x,1/n) does.

Don't expect a quartic solution anytime soon, in algebraic notation its about 4 times bigger than a cubic.. in code that may translate to even more. No torii or helices just yet. :(


Hopefully someone other than me will find this code useful, or at least interesting.

Zana

25th July 2003 18:47 UTC

*shwoop*

that's way over my head....

UnConeD

25th July 2003 18:55 UTC

Wow ;). Very nice... unfortunately the more interesting shapes that I can think of are all described by a quartic, but it's good to have such code handy.

Anyway I went on an optimization spree:

x=sqr(c);y=sqr(d);
a=(c*e-.333333*y)/x;
b=((.0740740*y-.333333*c*e)*d+x*f)/(x*c);
tmp=sqr(b)+(4*a*sqr(a))*.0370370;
x=above(tmp,0);y=sqrt(tmp);
p3r=if(x,y,0)-b;
p3i=if(x,0,y);
q3r=if(x,y,0)+b;
q3i=if(x,0,y);
pk=pow((sqr(p3r)+sqr(p3i))*.25,.166667);
qk=pow((sqr(q3r)+sqr(q3i))*.25,.166667);
p0=atan2(p3r,p3i)*2.09440;
q0=atan2(q3r,q3i)*2.09440;
pr1=pk*cos(p0);pi1=pk*sin(p0);
pr2=pk*cos(p0*2);pi2=pk*sin(p0*2);
pr3=pk*cos(p0*3);pi3=pk*sin(p0*3);
qr1=qk*cos(q0);qi1=qk*sin(q0);
qr2=qk*cos(q0*2);qi2=qk*sin(q0*2);
qr3=qk*cos(q0*3);qi3=qk*sin(q0*3);
t1=if(above(abs(pi1+qi1),0.0001),10000,pr1+qr1);
t2=if(above(abs(pi2+qi2),0.0001),10000,pr2+qr2);
t3=if(above(abs(pi3+qi3),0.0001),10000,pr3+qr3);
t1=if(below(t1,0),10000,t1);
t2=if(below(t2,0),10000,t2);
t3=if(below(t3,0),10000,t3);
t=min(min(t1,t2),t3);

Non-obvious stuff:
- Got rid of two sqrt's (merged them into pow())
- Replaced below(tmp,0) by bnot(above(tmp,0)). This is a minor accuracy loss, but it shouldn't matter that much really.

I replaced constant subexpressions by a single constant (w/ 6 significant digits). Experience with the 'variable LED display' thingy shows that this is about as much accuracy as you can expect from AVS.

I didn't optimize p0*2, p0*3 (and same for q0) because sometimes simple statements like this are faster than the extra work required for storing it in a variable.

Oh and of course for the people who hadn't realized, the cubic solver stops before the t1 part, and the solutions are (pr1+qr1, pi1+qi1), (pr2+qr2, pi2+qi2), (pr3+qr3, pi3+qi3). The first t1,t2,t3-statements are discarding the imaginaries, while the second series is discarding the negative reals (for raytracing).


Oh and Jheriko, do you have a mathematical description of this method? Certain parts puzzle me... like why you multiply p0,q0 by twopi, if they're already in radians (due to atan2).

anubis2003

25th July 2003 19:15 UTC

Can someone who understands this tell me what it would be useful in doing? Is it just to make complex-shaped 3D DMs???

UnConeD

25th July 2003 19:46 UTC

anubis2003, it solves a cubic equation, i.e. an equation of the form "cx^3+dx^2+ex+f=0" where x is unknown and c,d,e,f are given constants.

For a linear equation ax+b=0, the result is x=-b/a. For a quadratic equation ax^2+bx+c, if you set D=b^2-4*a*c the solutions are x=(-b ± sqrt(D))/2a

It is possible to write out a general solution to a quartic (4th degree) equation, but it would be very difficult and complicated.

I believe it has been proven (well something along these lines, Jheriko can elaborate) that there is impossible to find a general solution for any equation with degree 5 or higher, but (numerical) methods exist to solve them.

Cubic equations will be rare for the stuff used in AVS, but if you ever run into one, you can now solve it using this code ;).

anubis2003

25th July 2003 19:58 UTC

Yeah, I knew what the Cubic solution was, but what I meant by my question is what would be the use in AVS? Just to make some more complex DMs?

dirkdeftly

25th July 2003 21:50 UTC

Originally posted by UnConeD
- Replaced below(tmp,0) by bnot(above(tmp,0)). This is a minor accuracy loss, but it shouldn't matter that much really.

wtf? that's one extra function call in order to lose accuracy...what the hell is the point?

jheriko

25th July 2003 22:07 UTC

I'll go through the entire solution step by step here:

the first two lines a=blah,b=blah are reducing the cubic cy^3+dy^2+ey+f to x^3+ax+b by substituting y=x-d/3c.. something which i just realised i forgot to undo at the end of the solution :P so my code is wrong! but this is a small error and is easily fixed by adding d/3c to the solution. if you wish to use the cubic code this will be important! it also explains why i couldn't move the camera in my test... which i had originally put down to the presence of an extruded singularity.

the idea behind solving this is that we subsitute x=p+q into the reduced cubic noting that

(p+q)^3=p^3+q^3+3pq(p+q)

so we can see that b=-p^3-q^3 and a = -3pq

we can find the values of p and q from this using several methods my favorite one is to square p^3+q^3 and p^3-q^3 giving us p^6+2p^3 q^3+q^6 and p^6-2p^3 q^3+q^6

subtracting the expansions

b^2-(p^3-q^3)^2 = 4p^3q^3

finding -4p^3q^3 in terms of a gives us:

(p^3-q^3)=sqrt(b^2+(4a^3)/27)

adding subtracting our (p^3+q^3) and (p^3-q^3) gives us 2p^3 or -2q^3 and since we now know what p^3+q^3 and p^3-q^3 are in terms of a and b we can then find what 2^p^3 and -2q^3 are... the rest of the solution is fairly trivial maths: +/- and cube roots and results in an expression for (p+q) too big and ugly for me to be bothered to type up here.

now, you're still wondering where sin cos and 2pi come into it all, well the problem is that when solving a cubic you may have to use complex numbers in the intermediary stages even if you have real roots so since avs can't handle complex numbers you have to calculate the square roots and cube roots of them yourself. i did this using de moivre's theorem, a very simple consequence to euler's theorem: e^ix=cos(x)-isin(x) which is k^x=|k|^x*(cos(k*arg(k)+isin(k*arg(k)), the fraction of two pi is used to get all of three of the cube roots rather than just one by cycling k around the origin by 1/3 of a circle each time.


it is impossible to solve polynomials of degree higher than 4 so qunitics are out of the question, however i can't really elaborate much on why because the maths behind it is *VERY* advanced, in fact the guy who first proved this invented a whole branch of mathematics in doing so and some of his work is *still* not understood by the worlds leading mathematicians more than 100 years after his death.

edit: i mean general polynomials higher than degree 4 cant be solved in terms of radicals (+,-,*,/ and ^) since some higher degree polynomials can be solved analytically. there is a solution to the general quintic (degree 5) that I know of but it involves what is called a 'hypergeometric' function and would be of no use to avs since like sin and cosine are the sum of infinite series hypergeometrics are products of such series.

UnConeD

25th July 2003 23:03 UTC

Atero: if you had bothered to read the code, rather than shouting bloody murder because you spotted what you thought was an error and you are always right, you would've seen that I did not put in bnot() at all. I replaced below(tmp,0) by above(tmp,0) (which is calculated once) and switched the places of the true/false clauses in an if statement. I wrote bnot() to make it clear why there was a loss in accuracy.
Plus, it doesn't really matter if stuff is a function call... to AVS, a*b is no more or less a function call than atan2(a,b).

Jheriko: I know about the whole handyness of writing complex numbers in polar form. The n-roots of a number are the corner points of an n-gon.

Your typography is a bit messed up (see this Penny Arcade comic), but the correct expression is I believe:

k^x=|k|^x*(cos(arg(k)*x)+isin(arg(k)*x))

So if atan2 already returns the angle in radians (on a scale of 2π ) why multiply by 2π again?

jheriko

26th July 2003 04:13 UTC

Yeah... I neglected the *x's in de Moivre.

Those '*' should really be '+' and the '/3' should apply to both the pi multiples and the argument... damn another mistake... I really should have tested this code more thoroughly than I did.. :/ I apologize for unleashing this buggy confusing mess upon you. I'm going to have to check this code over very thoroughly now. Maybe I can't move the camera because my code is just a pile of cack. :(

Strange though that it draws the correct shape when the camera is on the xy origin line. After fixing this mistake I still can't move the camera, the shape doesn't become such a huge mess as it did before though.

I've tested the code against my calculator's in built cubic solver, the way in which it calculates cube roots is correct but the solutions come out wrong. I can't find flaw in the method and its one which I've converted sucessfully to qbasic, c++ and the languages of two calculators so I must have made some error in converting it to code. :(

I will fix this! The only good thing so far is that the code isn't completely wrong since it does draw my cubic correctly (well it looks that way) from 0,0,z. I'll post back here when its fixed. I need to sleep now.

Typography better?

UnConeD

26th July 2003 04:33 UTC

If you construct your own cubic with known solutions (x-a)(x-b)(x-c) and solve it again, you can easily verify using an equal() or below(abs(x-a),epsilon).
Or you could use the AVS variable displayer...

I believe the cube roots of p,q should be calculated like this:
- Calculate D = sqrt(sqr(Re)+sqr(Im)), R = atan2(Im,Re)
- Take the cube-root of D = d
- Divide R by 3 = r

p = d*(cos(r)+i*sin(r))
p' = d*(cos(r*2)+i*sin(r*2))
p'' = d*(cos(r*3)+i*sin(r*3))

No π anywhere...

jheriko

26th July 2003 16:38 UTC

You need pi to calculate the second and third cube roots. Consider your method above, do the three roots describe a regular polygon centred on the origin? You need to use 2pi to rotate around a fraction of a circle in radians... there is no other way that I can think of.

I am still having real trouble identifying whats wrong with the solver code, I've attached my test preset which tests the calculated roots of

x^3-2x+1=0

They should be -1.618, .618 and 1 the code gets the 1 but not the other two.

The cube root code is where the three roots 'split' but this code uses the correct method according to equivalent code on my calculator. :/

jheriko

26th July 2003 17:01 UTC

Here is an 'AVS proof' that my cube root method works fine.

UnConeD

26th July 2003 19:15 UTC

Yeah indeed, mine didn't form a regular n-gon. My fault for doing 5AM math stuff ;).

Zevensoft

27th July 2003 11:53 UTC

Why is it when digging in my square garden I end up with cubic roots :igor:

jheriko

27th July 2003 20:28 UTC

Because digging adds a third dimension (up/down) to your square garden so that you can now define sections of cubes. :P

jheriko

28th July 2003 15:37 UTC

I have derived the quartic solution so once I can get this cubic working I can proceed to a quartic solver... I've started a fresh cubic solver from scratch using a different method. For those of you interested in the solution to the general quartic the equation that I have derived when all the substitutions are removed covers a little more than a single side of A4 paper, and that with writing square roots and divisions on single lines to save on the space (e.g. x= sqrt(x/2)) i was hoping to fit it one side for my ease of use. Porting it to code will be a mammoth task.

I have just right now thought of somethings that could optimise the solution of any polynomial including quadratics which is that the product of the roots for a polynomial ax^n+bx^n-1...nx+k = k/a, also the sum of the roots is equal to -n/a, this means that you can calculate the last root by subtracting the others from -n/a.

jheriko

28th July 2003 20:57 UTC

damn.. too late for edit :(

n/a should read b/a.

jheriko

28th July 2003 21:56 UTC

I have coded in to avs another method for solving the cubic, it doesn't work either which is infuriating. I've attached two presets, one using the new method and one using the old, they both 'look' correct on the xy origin, moving the camera along z causes no problems but moving it in x or y creates a bizarre effect. Now that I've coded it twice and check my method through by calculator and mathematica I'm convinced that there must be something else wrong with the preset other than the cubic solver part, but I can't see any flaws myself. I hate to say it but I'm going to give up on this.

anubis2003

1st August 2003 02:14 UTC

Now if someone figures out the quartic solver and gets it working, then we would be able to have DM torus's, right? And with those, couldn't we make DM tunnels that have nice smooth turns? Just intersect two cylinders with a torus, and you got a curved tunnel.

UnConeD

1st August 2003 02:57 UTC

I'm busy with a different approach... numerical solving.

I've tried Newton-Raphson, but it sucks for the limited amount of iterations that you can do in AVS. I've also tried Halley Rational/Irrational iteration, which converge really great if you have a good initial guess, but are bad for certain cases where f'(x) and/or f''(x) are near zero. Sometimes they also converge to the wrong root. I thought of using Muller-iteration too, but it would flip on quartics with only two roots (unless you're interested in complex roots).

Going to try Bairstow's method now.

anubis2003

1st August 2003 03:18 UTC

Never heard of any of that.:p Hope it works though - I would love to see some sweet DMs.

jheriko

1st August 2003 13:15 UTC

My theory on the complete failure of my cubic solver is that like when raytracing planes 'negative images' (possibly two of them?) are created, due to the sign preserving properties of x^3, and due to the fact that cubics are not planar these image(s) intersect the actual surface which could result in the sort of messy nastyness that is seen when the camera moves from the origin. If that is the case then a quartic (like a quadratic) should not create such images and so a straight analytical solution should work. If anyone can be bothered to convert this behemoth into code here is a mathematica generated solution for x^4+ax^2+bx+c=0, the reduced quartic obtained from dy^4+ey^3+fy^2+gy+h=0 when the substitution y=x-(e/4d) is made. Each x-> (blah) is a seperate root.

anubis2003

1st August 2003 13:18 UTC

forget something?

jheriko

1st August 2003 13:36 UTC

Damn it! It didn't post my attachment or tell me that it was too big. :(
Not my fault! I sat there and watched it take its time uploading my attachment only to have it not get posted. Grrr....


I'll just have to post it in two parts.

Interestingly 2 color bitmap produced a smaller filesize than jpg by about 20kb and png by abut 120kb.

jheriko

1st August 2003 13:40 UTC

...

UnConeD

1st August 2003 13:54 UTC

What about a 2-bit png though? :)

I don't really think the 'negative image' has anything to do with the cubicness. Remember that this is only caused because we solve for a ray that extends in both directions, not just in the looking direction.
So a cubic should still only result in 3 solutions, some of which might be negative/complex and should be discarded.

For a quadratic it worked exactly the same as with a plane: solve and discard the negatives/complexes, so there's no reason it should be different from a cubic/quartic.

UnConeD

1st August 2003 15:39 UTC

I can't seem to get Bairstow's method to do anything remotely useful in AVS. It's guaranteed to give a correct result eventually, but when that happens is very dubious, and seems only to occur after about 20-30 iterations. In C++ code that's nothing, but in AVS you can only get away with 10 max, and less in a high-detail DM.

If anyone's interested, here's my code. It starts with a quartic a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0, and factors it into (x^2 - r*x - s)*(k4*x^2 + k3*x + k2) = 0. You can then solve these two quadratics the regular way.

Remember that this is a numerical approximation, not an exact result! The iteration loop is marked as such, and can be repeated to increase the accuracy. You'll need at least 5-6 to get a somewhat acceptable result, but you'll probably need (a lot) more.

I'm now experimenting with a hybrid Bairstow/Halley approach, which uses a low-iteration bairstow approach to get a good initial guess and then uses Halley to refine it. It seems to work better, but it's still not perfect.

The main problem is that most numerical methods deal with getting a really accurate result. In AVS, you only need a somewhat accurate result, but calculated really fast. No such methods exist, at least none that guarantee a somewhat accurate result everywhere.

£ Setup Quartic a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0;
a4=1;a3=2;a2=3;a1=4;a0=5;
£ Initial guess for r,s (tweak until it's right for your problem);
r=1;s=1;

£ Bairstow Method Iteration (repeat this block as much as you need)
k3=a3+r*a4;k2=a2+r*k3+s*a4;k1=a1+r*k2+s*k3;k0=a0+r*k1+s*k2;
c3=k3+r*a4;c2=k2+r*c3+s*a4;c1=k1+r*c2+s*c3;
x=1/(sqr(c2)-c1*c3);
r=r+(k0*c3-k1*c2)*x;s=s+(k1*c1-k0*c2)*x;
£ End Bairstow Iteration (stop repeat here);

£ Solve the quadratics (reals only, imaginaries are set to 5000);
disc=sqr(r)*.5+2*s;x=above(disc,0);r=r*.5;disc=sqrt(disc);
k0=if(x,r-disc,5000);k1=if(x,r+disc,5000);
disc=sqr(k3)-4*k2*a4;x=above(disc,0);y=.5/a4;disc=sqrt(disc);
k2=if(x,(-k3-disc)*y,5000);k3=if(x,(-k3+disc)*y,5000);

£ The solutions are now in k0,k1,k2,k3;

And here's the code for Halley iteration of a quartic root. Given an approximation k, it calculates a better approximation:

£ Halley Iteration;
fk0=(((a4*k+a3)*k+a2)*k+a1)*k+a0;
fk1=((4*a4*k+3*a3)*k+2*a2)*k+a1;
fk2=(12*a4*k+6*a3)*k+2*a2;
k=k-((fk0*fk1)/(sqr(fk1)-fk0*fk2*.5));

fk0 is the value at k, fk1/fk2 are the first/second derivatives at k. You can use this technique for any function really, but if the second derivative is too complicated, you're better off with Newton-Raphson because it only uses the first derivative (though Newton-Raphson only converges with order 2, while Halley has order 3).

Zevensoft

2nd August 2003 03:04 UTC

Tell me if this idea is right:

jheriko

2nd August 2003 13:20 UTC

Not sure what the diagram is supposed to represent but you do usually want the smallest positve root, sometimes the largest positve root is useful too though (like for making a transparent cylinder).

snoozin_hermit

27th August 2003 00:17 UTC

stuff is great for those with 'nuff education to know that 2+2=4 but I only went to grammer school and all i know is 1+1=2

S-uper_T-oast

28th August 2003 03:30 UTC

Stupid ass n00b.
1st: DON'T REVIVE OLD THREADS!
2nd: I am assuming your first language is english so use it correctly, this isn't some loser AIM chatroom were you are in a rush to type everything out.

snoozin_hermit

2nd September 2003 04:51 UTC

nwr