This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Is the maximum error bound in __mul really 1.001ULP?
- From: Siddhesh Poyarekar <siddhesh at redhat dot com>
- To: libc-alpha at sourceware dot org
- Date: Fri, 19 Apr 2013 16:47:59 +0530
- Subject: 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