True. Reimplementing same thing for 64-bit precision (in Pascal) works okay.

I succeed with numeric Barstow method and avoid
monkey Ferrari formulas. Because evil Olympiad precision 1.e-09 I was forced to implement Barstow above BigGecimal in Java.
Who want good tests use MathCad. It's integration in case
rather accurate.
For those who inerested I write Barstow method from my prog:
shorten- also my method for cutting off long BigDecimals
public static BigDecimal Barstow(BigDecimal a[], int N, BigDecimal u0, BigDecimal v0, BigDecimal eps, int MaxIter, int MaxLen)[]

        {
            BigDecimal u, v;
            BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
            String str = null;
            int k, n,faza;
            u = u0;
            v = v0;
            BigDecimal[] b = new BigDecimal[N+1];
            BigDecimal[] c = new BigDecimal[N+1];
            BigDecimal rez[] = new BigDecimal[N+3];
            BigDecimal d, s, q,h;
            BigDecimal Z = new BigDecimal("0.0");
            BigDecimal bb1, bb2;
            h = new BigDecimal("0.0000000000000000000000000000000000000001");
            faza = 0;
            for (k = 0; k <= MaxIter; k++)
            {

                for (n = 0; n <= N; n++)
                {
                    if (n < 1) bb1 = Z; else bb1 = b[n - 1];
                    if (n < 2) bb2 = Z; else bb2 = b[n - 2];
                    b[n] = shorten(a[n].add(bb1.multiply(u)).add(bb2.multiply(v)), MaxLen);

                }
                for (n = 0; n <= N - 1; n++)
                {
                    if (n < 1) bb1 = Z; else bb1 = c[n - 1];
                    if (n < 2) bb2 = Z; else bb2 = c[n - 2];
                    c[n] = shorten((b[n].add(bb1.multiply(u))).add(bb2.multiply(v)), MaxLen);

                }
                d = shorten(((b[N].multiply(c[N - 3])).subtract(b[N - 1].multiply(c[N - 2]))), MaxLen);
                q = shorten((c[N - 2].multiply(c[N - 2])).subtract(c[N - 1].multiply(c[N - 3])), MaxLen);

                d = shorten(d.divide(q, BigDecimal.ROUND_FLOOR), MaxLen);
                s = shorten(((b[N - 1].multiply(c[N - 1])).subtract(b[N].multiply(c[N - 2]))), MaxLen);
                s = shorten(s.divide(q, BigDecimal.ROUND_FLOOR), MaxLen);
                if ((d.abs().compareTo(eps) < 0) && (s.abs().compareTo(eps) < 0))
                {
                    if (faza == 0) { faza = 1; MaxLen = MaxLen * 4; eps = eps.multiply(h); }
                    else break;
                }

                u = u.add(d);
                v = v.add(s);
                k = k;


            }
            rez[0] = u;
            rez[1] = v;
            for (n = 0; n <= N; n++) rez[n + 2] = b[n];
        return rez;

        }


Edited by author 27.09.2008 15:32