Re: Accuracy in float and double (in C)

From: Peter Spellucci (spellucci_at_fb04373.mathematik.tu-darmstadt.de)
Date: 02/25/05 

Date: Fri, 25 Feb 2005 17:12:27 +0000 (UTC)

In article <ekjgkytxje1j@legacy>,
 carvegio@gmail.com (Steve) writes:
>Hello, I am trying to solve this problem in c and I am having a great
>deal of problem trying to figure out what is going on.
>
>I was hoping someone could please explain it to me.
>
>It involves a three term recurrence x[i]=A*x[i-1]+B*x[i-2]
>with x[0]=x[1]=...x[n]=1/3 (0.33333r).
>
>A is (double)(100001+1), B is (double)(-100001)
>
>When I work this out using float I get perfectly accurate answers
>but when I change type to double, my answers are horribly wrong
>and get worse and worse as we progress, and similarily with long
>double.
>
>I don't understand why, I know floating point can not represent
>1/3 accurately. But I would of thought that float would have the
>same problem storing it as float, and cause bad results.
>But this is not the case.
>
>Also if we change the recurrence to
>x[i]=(A-1)*x[i-1]+B*x[i-2]+x[i-1] we get perfect answers, and I
>dont understand how this is so different to what we had before,
> just taking one part of A outside.
>
>Hope someone can help
>
>Thanks
>

this is quite easy;
you have a linear recurrence which is unstable.
the general solution is

   x[n]=A+B*(100001)^n.
by selecting x[0]=x[1]=c you get theoretically x[i]=c for all i.
but in evaluating this numerically there will be roundoff errors
and these will be different in different arithmetic.
in float the roundoffs cancel out and also if you rewrite the recurrence
as x[i]=100001*x[i-1]-100001*x[i-2]+x[i-1];
but in double and in original form you get , because rouding is different for th
e
two multiplications, a blowup of 100001^i*eps (in this case
eps .= 10^{-19} because the data will be held in registers) in step i
and this explains what you see.
hth
peter