This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Is the maximum error bound in __mul really 1.001ULP?


Hi,

I've been trying to reproduce the computation that resulted in the
conclusion that the __mul implementation in mpa.c has a maximum error
bound of 1.001ULP.  I haven't succeeded in doing that and the best
estimate I have for maximum error bound is in the range of R^-1ULP,
i.e. 2^-24ULP.  I'm sure I must be doing something wrong, but I can't
figure out what it is:

The accurate multiplication result of two numbers (A, B) with
precision p and radix R is given by:

SUM{i = 2 to 2p} (T[i] * R ^ (E-i))

where each term T[i] is computed as:

up to pth term:  SUM{j = 1 to i - 1} (A[j] * B[i - j])
p+1th term to 2pth term: SUM{j = i - p to p} (A[j] * B[p - j])

The first term T[1] is reserved for carry, if any.

The __mul algorithm computes only the first p digits and two
additional digits as guard.  So the maximum absolute error in the
result is:

    SUM{i = 3 to p} ((p - i) * (R - 1)^2 * R^(E - p - i))

    = R^(E - p) * SUM{i = 3 to p} ((p - i + 1) * (R - 1)^2 * R^-i)

    = R^(E - p) * (R - 1)^2 * R^-3
      * SUM{i = 0 to p - 3} ((p - i - 2) * R^-i)

This is again bound by:

    2 * (p - 2) * R^(E - p) * (R - 1)^2 * R^-3

So the maximum error bound in ULP is the above divided by R^(E - p),
i.e.:

    2 * (p - 3) * (R - 1)^2 * R^-3

This is of the order of R^-1.  So the multiplication result should be
accurate up to ~2^-24ULP.

So what's wrong with this?

Siddhesh


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]