/* I don't know much c programming, but i'm using MS visual c++ 6.0, and my problem is with basic c programming.

I have these variables:
*/
double k=.987645642, m=5.657897456;
double x=5.65465, y=54.98789456;
double d=2.654564, D=-54.65489478;
double f=-5.654572002404, F=-.56489798006405640;
double w=15, t, root, a, A;

/* It doesn't really matter what they represent, but I have an equation that uses them. I will split it up in four pieces, to represent it easier:

r=(2a+2Ak-2km+-root)/(2(k^2+1))
a=x+dt+ft^2/2
A=y+Dt+Ft^2/2
root=(-2a-2Ak+2km)^2-4(k^2+1)(a^2+A^2+m^2-w^2-2Am)

The equation is a function, r(t). I just want to find any value of t, where the root=0.

This means I can set the root to 0, and solve t. I did it with quickmath.com, giving 4 solutions. Now I don't care which, so I just use one for the example.*/

t=(2*D-2*d*k+sqrt((2*d*k-2*D)*(2*d*k-2*D)-8*(f*k-F)*(m-sqrt(k*k+1)*w+k*x-y)))/(2*(f*k-F));

a=x+d*t+f*t*t/2;
A=y+D*t+f*t*t/2;
root=(-2*a-2*A*k+2*k*m)*(-2*a-2*A*k+2*k*m)-4*(k*k+1)*(a*a+A*A+m*m-w*w-2*A*m);

/* now root should be 0, or close to 0. I get 0.149776, which is not too close to 0. Since double uses 17 digits precision, according to MS, it shouldn't be too much asking for 10 decimals being correct? */

/* How can I get 17 or 14 or as many digits precision as possible?! I don't want to use too much memory, and all these variables are necessary.

In the real programme I need this to be calculated every frame, along with other similar calculations.*/

If you don't want to use an MP library or even an ordinary numerical computing library, you'll have to take great care about how you code the calculations to avoid unnecessary loss of precision.

There are various rules about this: don't add two numbers of different magnitude, don't subtract two numbers of similar magnitude, remove unnecessary multiplications and run-time calculation of expressions involving only constants etc.

You may need to do a lot of work rearranging the algebra while taking into account these rules together with the domain of the constants as well as the variables and you should probably read Knuth's books or similar works for more in-depth analysis and tips.

One trick I know of to reduce multiplications in the calculation of polynomial expressions is Horner's method:

Instead of ax^n + bx^(n-1) + ... + z, calculate z + (...((ax + b)x +c)x + ...)x

which may help in the case of your quartic expression for root - or it may not! You could also check the derivatives at the four roots of "root" to see if "root" is flatter enough at some particular root to make it worthwhile switching.

Hmmm... Using your code exactly as it is and compiling with gcc gives:

1.194344337007608
-0.000000000001105

for t and root respectively. It seems the inaccuracy you're getting may not be coming from losses of precision due to the way your code is written and the nature of the calculation itself. I've no idea why your compiler would do so much worse - maybe you should read it's documentation.

Thanks for your replies. The problem seems to be in the code, but I had made a typo. My question is still how to solve the problem, since the root is too big.

I don't use ^ (powers) in my code, but I don't know how to do a square root without sqrt(). The variables also change, so it doesn't seem safe not to multiply and use subtraction with different values. And I'm not sure I can use maths multiplication, subtraction, and square roots.

So is it that I have to use an MP library to get somewhat accurate results?

----------------

I had typos in my code, resulting in the error. In my above code, somewhere close to the end, it says "A=y+D*t+f*t*t/2;" where it should be using F, ie "A=y+D*t+F*t*t/2;" is correct.

I have to use another example though, if you'd like to test it. I will paste the exact code here:

//start

double k=-1.125, m=287;
double x=140, y=106.922012;
double d=0, D=98.26169586;
double f=0, F=49.81;
double w=15, t, root, a, A;

t=(2*D-2*d*k+sqrt((2*d*k-2*D)*(2*d*k-2*D)-8*(f*k-F)*(m-sqrt(k*k+1)*w+k*x-y)))/(2*(f*k-F));

a=x+d*t+f*t*t/2;
A=y+D*t+F*t*t/2;
root=(-2*a-2*A*k+2*k*m)*(-2*a-2*A*k+2*k*m)-4*(k*k+1)*(a*a+A*A+m*m-w*w-2*A*m);

// END

From this I get:
t=-3.94546
root=0.0312500000000

This you can use for testing on quickmath.com. If you click solve under equations, then click advanced, and copy paste this:

k=-1.125
m=287
x=140
y=106.922012
d=0
D=98.26169586
f=0
F=49.81
w=15
t=(2*D-2*d*k+sqrt((2*d*k-2*D)*(2*d*k-2*D)-8*(f*k-F)*(m-sqrt(k*k+1)*w+k*x-y)))/(2*(f*k-F))
a=x+d*t+f*t*t/2
A=y+D*t+F*t*t/2
V=(-2*a-2*A*k+2*k*m)*(-2*a-2*A*k+2*k*m)-4*(k*k+1)*(a*a+A*A+m*m-w*w-2*A*m)

Click solve button, and get:
t=-3.94564,
V=0 (this is the root, I had to change it here because I can't use variables with multiple letters)

Yes, I didn't type in your code exactly as you wrote it - I corrected for your typo before I ran the program, so the result I got was sound.

With your new set of constants I got:

-4.018061558096647
0.000000000002217

So the root I got is less accurate than yours but the evaluation of "root" turns out much better! It seems the problem in your program is occurring in the expression for the variable "root" (this time at least). So have you tried rearranging that calculation in your program to see if it varies the accuracy? Perhaps taking it step by step - you could then use printf to see exactly where the divergence is happening - or perhaps expanding it in t, collecting terms and writing it in Horner's method form.

The trouble is of course that whatever you do, it's efficacy will likely vary with the numerical parameters, which means you would have to do all the hard work I mentioned if you wanted to be sure of consistent accuracy over the whole domain. Using an MP library would be an easy, brute force solution (probably) to this problem but (depending on your program requirements) maybe you could just test for when t is in the vicinity of the roots rather than when "root" is in the vicinity of zero or use some other tricks.

This sort of problem is really the domain of the math library or game programmer or CS student/graduate and no doubt there is a sensible solution that doesn't involve having to use an MP library. Since I'm not primarily a programmer, all I know is a few basic rules and beyond that I always rely on the power of libraries and tools such as gmp, gsl, maxima, pari/gp etc.

You could spend hours analysing your particular calculations, the domain and the behaviour of the compiler and arrive at a satisfactory solution or you could use an MP library but I am pretty sure there must be some general solution to this sort of problem that I don't happen to know about.

Note that quickmath.com does not appear to be using high precision arithmetic to calculate it's numerical results and if it gave V=0, it most likely evaluated it algebraically. The root found numerically by your C program was closer than the answer you got from quickmath. (When I checked the accuracy of the C code using Maxima and pari/gp, they both agreed on the result).

I'm making a physics engine and I do need both the root and t to be as accurate as possible. Because it is a simulation game and I need it to behave similarly on each computer. Now if it is quite innacurate it will have different error ratios on different framerates.

Then it is quite strange that you got another value for t, because both quickmath and my program got exactly the same. I don't want to rely on the order of calculations, because I will need different values to be correct at different times, and I might need to calculate back and forth.

I'm not sure however that it really is the inaccuracy of the numbers that is the problem. I might just try a library or do some more experimenting.

My question then is, could the compiler be as inaccurate as losing more than a dousin digits from a not extremely complex calculation?

Thank you for your time.

Exactly the same? - your program's -3.94546 was closer than quickmath's -3.94564 unless that was a typo. Anyway...

No - not that inaccurate with that calculation and those numbers. In fact the calculation of "root" (V) is quite insensitive to variation about the accurate root shown in the screenshot.

Your first example was wildly inaccurate too (assuming what you got for t was fairly close that time), so you're right - your compiled program is doing something very odd somewhere in this bit:

a=x+d*t+f*t*t/2;
A=y+D*t+F*t*t/2;
root=(-2*a-2*A*k+2*k*m)*(-2*a-2*A*k+2*k*m)-4*(k*k+1)*(a*a+A*A+m*m-w*w-2*A*m);

You need to find out where exactly it is going so wrong.

Corrections:

The Mathematica at quickmath, Maxima and pari/gp all seem to be using an approximation algorithm to find the roots when you enter them numerically, rather than plugging the numbers into the algebraic forms of the roots they've already calculated. I couldn't understand at first why Maxima was giving very different answers: one when I typed in the algebraic expression for one of the roots it had found and then substituted numbers into it, and another different answer when I substituted the numbers straight into the set of (algebraically found) roots.

It seems in the latter case that instead of doing the obvious and easy substitution, it goes all the way back to the original equation and then solves it numerically by some approximation method.

The roots it finds by approximation are in fact much less accurate: -4.01806... is really much better than -3.94564... and when I said the expression was insensitive about the root, I must have been asleep or looking at a different equation because it is actually highly sensitive - as is obvious given the mutual proximity of the roots and the sizes of the coefficients.

Now although there is some rule breaking going on in the C program (especially in the second example), gcc and the glibc math library still managed to get the results at a useful level of accuracy - even with that second example - and if your compiler or math library can't do as well, you should change them. (I don't think it's likely to be the machine - mine is an AMD (an athlon-xp) too).