[seqfan] Re: Reproducing A064690 (sign changes of a rational fraction iteration)

[Tim Peters]

[Sean A. Irvine <sairvin at gmail.com>]

> ...
> I am still somewhat concerned that the precision needed might be very high
> to be certain the numbers starting at 1081319 are correct, but perhaps the
> plots you have produced will be sufficiently convincing that there are no
> intervening terms?
>

The plot is suggestive, but not convincing on its own ;-)  What Kevin Ryde
posted later is convincing: using interval arithmetic to establish rigorous
bounds on the error.

I confirmed his results using a different interval package, the facilities
in the popular `mpmath` extension module for Python. That's floating-point
interval arithmetic, not the fixed-point Kevin used, but that shouldn't
make a material difference.

I'll give the code, to show how simple and fast this kind of thing _can_ be:

    from mpmath import iv
    iv.prec = 1000 # number of bits of precision
    one = iv.mpf(1)
    zero = iv.mpf(0)
    a = iv.mpf(0)
    b = iv.mpf(1)
    canned = []
    report = REPORT = 100_000
    i = 0

    def format(x):
        return f"({float(x.mid)} +- {float(x.delta)/2.0})"

    while True:
        if (a < zero) != (b < zero):
            print(i, format(a), format(b))
            canned.append(i)
        a, b = b, b - one/(b + one)
        i += 1
        if i >= report:
            print(f"{i=:_}", format(a))
            report += REPORT
        if zero in b:
            print("at", i, "interval includes 0", format(b))
            break

So it just continues until the most recent term contains 0 in its interval
of uncertainty. Then we can't know what sign the true value has. But,
before then, we can.

With the given thousand bits of precision, it gives up at iteration 3300370:

at 3300370 interval includes 0 (0.3806665240893587 +- 1.2334529875209603)