| .. _log_algorithm: |
| |
| ======================== |
| Log/Log10/Log2 Algorithm |
| ======================== |
| |
| .. default-role:: math |
| |
| In this short note, we will discuss in detail about the computation of |
| :math:`\log(x)` function, with double precision inputs, in particular, the range |
| reduction steps and error analysis. The algorithm is broken down into 2 main |
| phases as follow: |
| |
| 1. Fast phase: |
| |
| a. Range reduction |
| b. Polynomial approximation |
| c. Ziv's test |
| |
| 2. Accurate phase (if Ziv's test failed): |
| |
| a. Further range reduction |
| b. Polynomial approximation |
| |
| |
| Fast phase |
| ========== |
| |
| Range reduction |
| --------------- |
| |
| Let `x = 2^{e_x} (1 + m_x)` be a normalized double precision number, in which |
| `-1074 \leq e_x \leq 1022` and `0 \leq m_x < 1` such that |
| `2^{52} m_x \in \mathbb{Z}`. |
| |
| Then from the properties of logarithm: |
| |
| .. math:: |
| \log(x) &= \log\left( 2^{e_x} (1 + m_x) \right) \\ |
| &= \log\left( 2^{e_x} \right) + \log(1 + m_x) \\ |
| &= e_x \log(2) + \log(1 + m_x) |
| |
| the computation of `\log(x)` can be reduced to: |
| |
| 1. compute the product of `e_x` and `\log(2)`, |
| 2. compute `\log(1 + m_x)` for `0 \leq m_x < 1`, |
| 3. add step 1 and 2. |
| |
| To compute `\log(1 + m_x)` in step 2, we can reduce the range further by finding |
| `r > 0` such that: |
| |
| .. math:: |
| | r(1 + m_x) - 1 | < C \quad \quad \text{(R1)} |
| |
| for small `0 < C < 1`. Then if we let `u = r(1 + m_x) - 1`, `|u| < C`: |
| |
| .. math:: |
| \log(1 + m_x) &= \log \left( \frac{r (1 + m_x)}{r} \right) \\ |
| &= \log(r (1 + m_x) ) - \log(r) \\ |
| &= \log(1 + u) - \log(r) |
| |
| and step 2 can be computed with: |
| |
| a. extract `r` and `-\log(r)` from look-up tables, |
| b. compute the reduced argument `u = r(1 + m_x) - 1`, |
| c. compute `\log(1 + u)` by polynomial approximation or further range reduction, |
| d. add step a and step c results. |
| |
| |
| How to derive `r` |
| ----------------- |
| |
| For an efficient implementation, we would like to use the first `M` significant |
| bits of `m_x` to look up for `r`. In particular, we would like to find a value |
| of `r` that works for all `m_x` satisfying: |
| |
| .. math:: |
| k 2^{-M} \leq m_x < (k + 1) 2^{-M} \quad \text{for some} \quad |
| k = 0..2^{M} - 1. \quad\quad \text{(M1)} |
| |
| Let `r = 1 + s`, then `u` can be expressed in terms of `s` as: |
| |
| .. math:: |
| u &= r(1 + m_x) - 1 \\ |
| &= (1 + s)(1 + m_x) - 1 \\ |
| &= s m_x + s + m_x &\quad\quad \text{(U1)} \\ |
| &= s (1 + m_x) + m_x \\ |
| &= m_x (1 + s) + s. |
| |
| From the condition `\text{(R1)}`, `s` is bounded by: |
| |
| .. math:: |
| \frac{-C - m_x}{1 + m_x} < s < \frac{C - m_x}{1 + m_x} \quad\quad \text{(S1)}. |
| |
| Since our reduction constant `s` must work for all `m_x` in the interval |
| `I = \{ v: k 2^{-M} \leq v < (k + 1) 2^{-M} \}`, `s` is bounded by: |
| |
| .. math:: |
| \sup_{v \in I} \frac{-C - v}{1 + v} < s < \inf_{v \in I} \frac{C - v}{1 + v} |
| |
| For a fixed constant `|c| < 1`, let `f(v) = \frac{c - v}{1 + v}`, then its |
| derivative is: |
| |
| .. math:: |
| f'(v) = \frac{(-1)(1 + v) - (1)(c - v)}{(1 + v)^2} = \frac{-1 - c}{(1 + v)^2}. |
| |
| Since `|c| < 1`, `f'(v) < 0` for all `v \neq -1`, so: |
| |
| .. math:: |
| \sup_{v \in I} f(v) &= f \left( \inf\{ v: v \in I \} \right) |
| = f \left( k 2^{-M} \right) \\ |
| \inf_{v \in I} f(v) &= f \left( \sup\{ v: v \in I \} \right) |
| = f \left( (k + 1) 2^{-M} \right) |
| |
| Hence we have the following bound on `s`: |
| |
| .. math:: |
| \frac{-C - k 2^{-M}}{1 + k 2^{-M}} < s \leq |
| \frac{C - (k + 1) 2^{-M}}{1 + (k + 1) 2^{-M}}. \quad\quad \text{(S2)} |
| |
| In order for `s` to exist, we need that: |
| |
| .. math:: |
| \frac{C - (k + 1) 2^{-M}}{1 + (k + 1) 2^{-M}} > |
| \frac{-C - k 2^{-M}}{1 + k 2^{-M}} |
| |
| which is equivalent to: |
| |
| .. math:: |
| \quad\quad 2C - 2^{-M} + (2k + 1) 2^{-M} C > 0 |
| \iff C > \frac{2^{-M - 1}}{1 + (2k + 1) 2^{-M - 1}} \quad\quad \text{(C1)}. |
| |
| Consider the case `C = 2^{-N}`. Since `0 \leq k \leq 2^M - 1,` the right hand |
| side of `\text{(C1)}` is bounded by: |
| |
| .. math:: |
| 2^{-M - 1} > \frac{2^{-M - 1}}{1 + (2k + 1) 2^{-M - 1}} \geq |
| \frac{2^{-M - 1}}{1 + (2^{M + 1} - 1) 2^{-M - 1}} > 2^{-M - 2}. |
| |
| Hence, from `\text{(C1)}`, being an exact power of 2, `C = 2^{-N}` is bounded below |
| by: |
| |
| .. math:: |
| C = 2^{-N} \geq 2^{-M - 1}. |
| |
| To make the range reduction efficient, we will want to minimize `C` (maximize |
| `N`) while keeping the required precision of `s`(`r`) as low as possible. And |
| for that, we will consider the following two cases: `N = M + 1` and `N = M`. |
| |
| Case 1 - `N = M + 1` |
| ~~~~~~~~~~~~~~~~~~~~ |
| |
| When `N = M + 1`, `\text{(S2)}` becomes: |
| |
| .. math:: |
| \frac{-2^{-M - 1} - k 2^{-M}}{1 + k 2^{-M}} < s < |
| \frac{2^{-M - 1} - (k + 1) 2^{-M}}{1 + (k + 1) 2^{-M}}. |
| \quad\quad \text{(S2')} |
| |
| This is an interval of length: |
| |
| .. math:: |
| l &= \frac{2^{-M - 1} - (k + 1) 2^{-M}}{1 + (k + 1) 2^{-M}} - |
| \frac{-2^{-M - 1} - k 2^{-M}}{1 + k 2^{-M}} \\ |
| &= \frac{(2k + 1)2^{-2M - 1}}{(1 + k 2^{-M})(1 + (k + 1)2^{-M})} |
| \quad\quad \text{(L1)} |
| |
| As a function of `k`, the length `l` has its derivative with respect to `k`: |
| |
| .. math:: |
| \frac{dl}{dk} = |
| \frac{2^{2M + 1} - 2k(k + 1) - 1} |
| {2^{4M}(1 + k 2^{-M})^2 (1 + (k + 1) 2^{-M})^2} |
| |
| which is always positive for `0 \leq k \leq 2^M - 1`. So for all |
| `0 < k < 2^{-M}` (`k = 0` will be treated differently in edge cases), and for |
| `M > 2`, `l` is bounded below by: |
| |
| .. math:: |
| l > 2^{-2M}. |
| |
| It implies that we can always find `s` with `\operatorname{ulp}(s) = 2^{-2M}`. |
| And from `\text{(U1)}`, `u = s(1 + m_x) + m_x`, its `ulp` is: |
| |
| .. math:: |
| \operatorname{ulp}(u) &= \operatorname{ulp}(s) \cdot \operatorname{ulp}(m_x) \\ |
| &= 2^{-2M} \operatorname{ulp}(m_x). |
| |
| Since: |
| |
| .. math:: |
| |u| < C = 2^{-N} = 2^{-M - 1}, |
| |
| Its required precision is: |
| |
| .. math:: |
| \operatorname{prec}(u) &= \log_2(2^{-M-1} / \operatorname{ulp}(u)) \\ |
| &= \log_2(2^{M - 1} / \operatorname{ulp}(m_x)) \\ |
| &= M - 1 - \log_2(\operatorname{ulp}(m_x)). |
| |
| This means that in this case, we cannot restrict `u` to be exactly representable |
| in double precision for double precision input `x` with `M > 2`. Nonetheless, |
| for a reasonable value of `M`, we can have `u` exactly representable in double |
| precision for single precision input `x` (`\operatorname{ulp}(m_x) = 2^{-23}`) |
| such that `|u| < 2^{-M - 1}` using a look-up table of size `2^M`. |
| |
| A particular formula for `s` can be derived from `\text{(S2')}` by the midpoint |
| formula: |
| |
| .. math:: |
| s &= 2^{-2M} \operatorname{round}\left( 2^{2M} \cdot \operatorname{midpoint} |
| \left(-\frac{-2^{-M - 1} - k2^{-M}}{1 + k 2^{-M}}, |
| \frac{2^{-M-1} - (k + 1)2^{-M}}{1 + (k + 1) 2^{-M}}\right) \right) \\ |
| &= 2^{-2M} \operatorname{round}\left( 2^{2M} \cdot \frac{1}{2} \left( |
| \frac{-2^{-M - 1} - k2^{-M}}{1 + k 2^{-M}} + |
| \frac{2^{-M - 1} + (k + 1)2^{-M}}{1 + (k + 1) 2^{-M}} |
| \right) \right) \\ |
| &= 2^{-2M} \operatorname{round}\left( \frac{ |
| - \left(k + \frac{1}{2} \right) \left(2^M - k - \frac{1}{2} \right) } |
| {(1 + k 2^{-N})(1 + (k + 1) 2^{-N})} \right) \\ |
| &= - 2^{-2M} \operatorname{round}\left( \frac{ |
| \left(k + \frac{1}{2} \right) \left(2^M - k - \frac{1}{2} \right) } |
| {(1 + k 2^{-N})(1 + (k + 1) 2^{-N})} \right) \quad\quad \text{(S3)} |
| |
| The corresponding range and formula for `r = 1 + s` are: |
| |
| .. math:: |
| \frac{1 - 2^{-M - 1}}{1 + k 2^{-M}} < r \leq |
| \frac{1 + 2^{-M - 1}}{1 + (k + 1) 2^{-M}} |
| |
| .. math:: |
| r &= 2^{-2M} \operatorname{round}\left( 2^{2M} \cdot |
| \operatorname{midpoint}\left( \frac{1 - 2^{-M - 1}}{1 + k 2^{-M}}, |
| \frac{1 + 2^{-M - 1}}{1 + (k + 1) 2^{-M}}\right) \right) \\ |
| &= 2^{-2M} \operatorname{round}\left( 2^{2M} \cdot \frac{1}{2} \left( |
| \frac{1 + 2^{-M-1}}{1 + (k + 1) 2^{-M}} + \frac{1 - 2^{-M-1}}{1 + k 2^{-M}} |
| \right) \right) \\ |
| &= 2^{-2M} \operatorname{round}\left( 2^{2M} \cdot \frac{ |
| 1 + \left(k + \frac{1}{2} \right) 2^{-M} - 2^{-2M-2} }{(1 + k 2^{-M}) |
| (1 + (k + 1) 2^{-M})} \right) |
| |
| Case 1 - `N = M` |
| ~~~~~~~~~~~~~~~~ |
| |
| When `N = M`, `\text{(S2)}` becomes: |
| |
| .. math:: |
| \frac{-(k + 1)2^{-M}}{1 + k 2^{-M}} < s < \frac{-k 2^{-M}}{1 + (k + 1) 2^{-M}} |
| \quad\quad \text{(S2")} |
| |
| This is an interval of length: |
| |
| .. math:: |
| l &= \frac{- k 2^{-M}}{1 + (k + 1) 2^{-M}} - |
| \frac{- (k + 1) 2^{-M}}{1 + k 2^{-M}} \\ |
| &= \frac{2^{-M} (1 + (2k + 1) 2^{-M})}{(1 + k 2^{-M})(1 + (k + 1)2^{-M})} |
| \quad\quad \text{(L1')} |
| |
| As a function of `k`, its derivative with respect to `k`: |
| |
| .. math:: |
| \frac{dl}{dk} = |
| -\frac{2^{-2M}(k(k + 1)2^{-M + 1} + 2^{-M} + 2k + 1)} |
| {(1 + k 2^{-M})^2 (1 + (k + 1) 2^{-M})^2} |
| |
| which is always negative for `0 \leq k \leq 2^M - 1`. So for `M > 1`, `l` is |
| bounded below by: |
| |
| .. math:: |
| l > \frac{2^{-M - 1} (3 - 2^{-M})}{2 - 2^{-M}} > 2^{-M - 1}. |
| |
| It implies that we can always find `s` with `\operatorname{ulp}(s) = 2^{-M-1}`. |
| And from `\text{(U1)}`, `u = s(1 + m_x) + m_x`, its `ulp` is: |
| |
| .. math:: |
| \operatorname{ulp}(u) &= \operatorname{ulp}(s) \cdot \operatorname{ulp}(m_x) \\ |
| &= 2^{-M - 1} \operatorname{ulp}(m_x). |
| |
| Since: |
| |
| .. math:: |
| |u| < C = 2^{-N} = 2^{-M}, |
| |
| Its required precision is: |
| |
| .. math:: |
| \operatorname{prec}(u) &= \log_2(2^{-M} / \operatorname{ulp}(u)) \\ |
| &= \log_2(2 / \operatorname{ulp}(m_x)) \\ |
| &= 1 - \log_2(\operatorname{ulp}(m_x)). |
| |
| Hence, for double precision `x`, `\operatorname{ulp}(m_x) = 2^{-52}`, and the |
| precision needed for `u` is `\operatorname{prec}(u) = 53`, i.e., `u` can be |
| exactly representable in double precision. And in this case, `s` can be |
| derived from `\text{(S2")}` by the midpoint formula: |
| |
| .. math:: |
| s &= 2^{-M - 1} \operatorname{round}\left( 2^{M + 1} \cdot |
| \operatorname{midpoint} \left(-\frac{-(k + 1)2^{-M}}{1 + k 2^{-M}}, |
| \frac{-k2^{-M}}{1 + (k + 1) 2^{-M}}\right) \right) \\ |
| &= 2^{-M - 1} \operatorname{round}\left( 2^{M + 1} \cdot \frac{1}{2} \left( |
| \frac{-(k + 1)2^{-M}}{1 + k 2^{-M}} + \frac{-k2^{-M}}{1 + (k + 1) 2^{-M}} |
| \right) \right) \\ |
| &= -2^{-M - 1} \operatorname{round}\left( \frac{ |
| (2k + 1) + (2k^2 + 2k + 1) 2^{-M} } |
| {(1 + k 2^{-N})(1 + (k + 1) 2^{-N})} \right) \quad\quad \text{(S3')} |
| |
| The corresponding range and formula for `r = 1 + s` are: |
| |
| .. math:: |
| \frac{1 - 2^{-M}}{1 + k 2^{-M}} < r \leq \frac{1 + 2^{-M}}{1 + (k + 1) 2^{-M}} |
| |
| .. math:: |
| r &= 2^{-M-1} \operatorname{round}\left( 2^{M + 1} \cdot |
| \operatorname{midpoint}\left( \frac{1 - 2^{-M}}{1 + k 2^{-M}}, |
| \frac{1 + 2^{-M}}{1 + (k + 1) 2^{-M}}\right) \right) \\ |
| &= 2^{-M-1} \operatorname{round}\left( 2^{M + 1} \cdot \frac{1}{2} \left( |
| \frac{1 + 2^{-M}}{1 + (k + 1) 2^{-M}} + \frac{1 - 2^{-M}}{1 + k 2^{-M}} |
| \right) \right) \\ |
| &= 2^{-M - 1} \operatorname{round}\left( 2^{M + 1} \cdot \frac{ |
| 1 + \left(k + \frac{1}{2} \right) 2^{-M} - 2^{-2M-1} }{(1 + k 2^{-M}) |
| (1 + (k + 1) 2^{-M})} \right) |
| |
| Edge cases |
| ---------- |
| |
| 1. When `k = 0`, notice that: |
| |
| .. math:: |
| 0 = k 2^{-N} \leq m_x < (k + 1) 2^{-N} = 2^{-N} = C, |
| |
| so we can simply choose `r = 1` so that `\log(r) = 0` is exact, then `u = m_x`. |
| This will help reduce the accumulated errors when `m_x` is close to 0 while |
| maintaining the range reduction output's requirements. |
| |
| 2. When `k = 2^{N} - 1`, `\text{(S2)}` becomes: |
| |
| .. math:: |
| -\frac{1}{2} - \frac{C - 2^{-M-1}}{2 - 2^{-M}} <> s \leq |
| -\frac{1}{2} + \frac{C}{2}. |
| |
| so when `C > 2^{-M - 1}` is a power of 2, we can always choose: |
| |
| .. math:: |
| s = -\frac{1}{2}, \quad \text{i.e.} \quad r = \frac{1}{2}. |
| |
| This reduction works well to avoid catastrophic cancellation happening when |
| `e_x = -1`. |
| |
| This also works when `C = 2^{-M - 1}` if we relax the condition on `u` to |
| `|u| \leq C = 2^{-M-1}`. |
| |
| Intermediate precision, and Ziv's test |
| -------------------------------------- |
| |
| In the fast phase, we want extra precision while performant, so we use |
| double-double precision for most intermediate computation steps, and employ Ziv |
| test to see if the result is accurate or not. In our case, the Ziv's test can |
| be described as follow: |
| |
| 1. Let `re = re.hi + re.lo` be the double-double output of the fast phase |
| computation. |
| 2. Let `err` be an estimated upper bound of the errors of `re`. |
| 3. If `\circ(re.hi + (re.lo - err)) == \circ(re.hi + (r.lo + err))` then the |
| result is correctly rounded to double precision for the current rounding mode |
| `\circ`. Otherwise, the accurate phase with extra precision is needed. |
| |
| For an easy and cheap estimation of the error bound `err`, since the range |
| reduction step described above is accurate, the errors of the result: |
| |
| .. math:: |
| \log(x) &= e_x \log(2) - \log(r) + \log(1 + u) \\ |
| &\approx e_x \log(2) - \log(r) + u P(u) |
| |
| come from 2 parts: |
| |
| 1. the look-up part: `e_x \log(2) - \log(r)` |
| 2. the polynomial approximation part: `u P(u)` |
| |
| The errors of the first part can be computed with a single `\operatorname{fma}` |
| operation: |
| |
| .. math:: |
| err_1 = \operatorname{fma}(e_x, err(\log(2)), err(\log(r))), |
| |
| and then combining with the errors of the second part for another |
| `\operatorname{fma}` operation: |
| |
| .. math:: |
| err = \operatorname{fma}(u, err(P), err_1) |
| |
| |
| Accurate phase |
| ============== |
| |
| Extending range reduction |
| ------------------------- |
| |
| Since the output `u = r(1 + m_x) - 1` of the fast phase's range reduction |
| is computed exactly, we can apply further range reduction steps by |
| using the following formula: |
| |
| .. math:: |
| u_{i + 1} = r_i(1 + u_i) - 1 = u_i \cdot r_i + (r_i - 1), |
| |
| where `|u_i| < 2^{-N_i}` and `u_0 = u` is representable in double precision. |
| |
| Let `s_i = r_i - 1`, then we can rewrite it as: |
| |
| .. math:: |
| u_{i + 1} &= (1 + s_i)(1 + u_i) - 1 \\ |
| &= s_i u_i + u_i + s_i \\ |
| &= u_i (1 + s_i) + s_i |
| &= s_i (1 + u_i) + u_i. |
| |
| Then the bound on `u_{i + 1}` is translated to `s_i` as: |
| |
| .. math:: |
| \frac{-2^{-N_{i + 1}} - u_i}{1 + u_i} < s_i < \frac{2^{-N_{i + 1}} - u_i}{1 + u_i}. |
| |
| Let say we divide the interval `[0, 2^-{N_i})` into `2^{M_i}` subintervals |
| evenly and use the index `k` such that: |
| |
| .. math:: |
| k 2^{-N_i - M_i} \leq u_i < (k + 1) 2^{-N_i - M_i}, |
| |
| to look-up for the reduction constant `s_{i, k}`. In other word, `k` is given |
| by the formula: |
| |
| .. math:: |
| k = \left\lfloor 2^{N_i + M_i} u_i \right\rfloor |
| |
| Notice that our reduction constant `s_{i, k}` must work for all `u_i` in the |
| interval `I = \{ v: k 2^{-N_i - M_i} \leq v < (k + 1) 2^{-N_i - M_i} \}`, |
| so it is bounded by: |
| |
| .. math:: |
| \sup_{v \in I} \frac{-2^{-N_{i + 1}} - v}{1 + v} < s_{i, k} < \inf_{v \in I} \frac{2^{-N_{i + 1}} - v}{1 + v} |
| |
| For a fixed constant `|C| < 1`, let `f(v) = \frac{C - v}{1 + v}`, then its derivative |
| is: |
| |
| .. math:: |
| f'(v) = \frac{(-1)(1 + v) - (1)(C - v)}{(1 + v)^2} = \frac{-1 - C}{(1 + v)^2}. |
| |
| Since `|C| < 1`, `f'(v) < 0` for all `v \neq -1`, so: |
| |
| .. math:: |
| \sup_{v \in I} f(v) &= f \left( \inf\{ v: v \in I \} \right) |
| = f \left( k 2^{-N_i - M_i} \right) \\ |
| \inf_{v \in I} f(v) &= f \left( \sup\{ v: v \in I \} \right) |
| = f \left( (k + 1) 2^{-N_i - M_i} \right) |
| |
| Hence we have the following bound on `s_{i, k}`: |
| |
| .. math:: |
| \frac{-2^{-N_{i + 1}} - k 2^{-N_i - M_i}}{1 + k 2^{-N_i - M_i}} < s_{i, k} |
| \leq \frac{2^{-N_{i + 1}} - (k + 1) 2^{-N_i - M_i}}{1 + (k + 1) 2^{-N_i - M_i}} |
| |
| This interval is of length: |
| |
| .. math:: |
| l &= \frac{2^{-N_{i + 1}} - (k + 1) 2^{-N_i - M_i}}{1 + (k + 1) 2^{-N_i - M_i}} - |
| \frac{-2^{-N_{i + 1}} - k 2^{-N_i - M_i}}{1 + k 2^{-N_i - M_i}} \\ |
| &= \frac{2^{-N_{i + 1} + 1} - 2^{-N_i - M_i} + (2k + 1) 2^{-N_{i + 1} - N_i - M_i}} |
| {(1 + k 2^{-N_i - M_i})(1 + (k + 1) 2^{-N_i -M_i})} |
| |
| So in order to be able to find `s_{i, k}`, we need that: |
| |
| .. math:: |
| 2^{-N_{i + 1} + 1} - 2^{-N_i - M_i} + (2k + 1) 2^{-N_{i + 1} - N_i - M_i} > 0 |
| |
| This give us the following bound on `N_{i + 1}`: |
| |
| .. math:: |
| N_{i + 1} \leq N_i + M_i + 1. |
| |
| To make the range reduction effective, we will want to maximize `N_{i + 1}`, so |
| let consider the two cases: `N_{i + 1} = N_i + M_i + 1` and |
| `N_{i + 1} = N_i + M_i`. |
| |
| |
| |
| The optimal choice to balance between maximizing `N_{i + 1}` and minimizing the |
| precision needed for `s_{i, k}` is: |
| |
| .. math:: |
| N_{i + 1} = N_i + M_i, |
| |
| and in this case, the optimal `\operatorname{ulp}(s_{i, k})` is: |
| |
| .. math:: |
| \operatorname{ulp}(s_{i, k}) = 2^{-N_i - M_i} |
| |
| and the corresponding `\operatorname{ulp}(u_{i + 1})` is: |
| |
| .. math:: |
| \operatorname{ulp}(u_{i + 1}) &= \operatorname{ulp}(u_i) \operatorname{ulp}(s_{i, k}) \\ |
| &= \operatorname{ulp}(u_i) \cdot 2^{-N_i - M_i} \\ |
| &= \operatorname{ulp}(u_0) \cdot 2^{-N_0 - M_0} \cdot 2^{-N_0 - M_0 - M_1} \cdots 2^{-N_0 - M_0 - M_1 - \cdots - M_i} \\ |
| &= 2^{-N_0 - 53} \cdot 2^{-N_0 - M_0} \cdot 2^{-N_0 - M_0 - M_1} \cdots 2^{-N_0 - M_0 - M_1 - \cdots - M_i} |
| |
| Since `|u_{i + 1}| < 2^{-N_{i + 1}} = 2^{-N_0 - M_1 - ... -M_i}`, the precision |
| of `u_{i + 1}` is: |
| |
| .. math:: |
| \operatorname{prec}(u_{i + 1}) &= (N_0 + 53) + (N_0 + M_0) + \cdots + |
| (N_0 + M_0 + \cdots + M_i) - (N_0 + M_0 + \cdots + M_i) \\ |
| &= (i + 1) N_0 + i M_0 + (i - 1) M_1 + \cdots + M_{i - 1} + 53 |
| |
| If we choose to have the same `M_0 = M_1 = \cdots = M_i = M`, this can be |
| simplified to: |
| |
| .. math:: |
| \operatorname{prec}(u_{i + 1}) = (i + 1) N_0 + \frac{i(i + 1)}{2} \cdot M + 53. |
| |
| We summarize the precision analysis for extending the range reduction in the |
| table below: |
| |
| +-------+-----+-----------+------------+--------------+-----------------+-------------------+ |
| | `N_0` | `M` | No. steps | Table size | Output bound | ulp(`s_{i, k}`) | prec(`u_{i + 1}`) | |
| +-------+-----+-----------+------------+--------------+-----------------+-------------------+ |
| | 7 | 4 | 1 | 32 | `2^{-11}` | `2^{-12}` | 60 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 2 | 64 | `2^{-15}` | `2^{-16}` | 71 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 3 | 96 | `2^{-19}` | `2^{-20}` | 86 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 4 | 128 | `2^{-23}` | `2^{-24}` | 105 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 5 | 160 | `2^{-27}` | `2^{-28}` | 128 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 6 | 192 | `2^{-31}` | `2^{-32}` | 155 | |
| | +-----+-----------+------------+--------------+-----------------+-------------------+ |
| | | 5 | 3 | 192 | `2^{-22}` | `2^{-23}` | 89 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 4 | 256 | `2^{-27}` | `2^{-28}` | 111 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 5 | 320 | `2^{-32}` | `2^{-33}` | 138 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 6 | 384 | `2^{-37}` | `2^{-38}` | 170 | |
| | +-----+-----------+------------+--------------+-----------------+-------------------+ |
| | | 6 | 3 | 384 | `2^{-25}` | `2^{-26}` | 92 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 4 | 512 | `2^{-31}` | `2^{-32}` | 117 | |
| | +-----+-----------+------------+--------------+-----------------+-------------------+ |
| | | 7 | 1 | 256 | `2^{-24}` | `2^{-15}` | 60 | |
| | | +-----------+------------+--------------+-----------------+-------------------+ |
| | | | 2 | 512 | `2^{-21}` | `2^{-22}` | 74 | |
| +-------+-----+-----------+------------+--------------+-----------------+-------------------+ |
| |
| where: |
| |
| - Number of steps = `i + 1` |
| - Table size = `(i + 1) 2^{M + 1}` |
| - Output bound = `2^{-N_{i + 1}} = 2^{-N_0 - (i + 1) M}` |
| - `\operatorname{ulp}(s_{i, k}) = 2^{-N_{i + 1} - 1}` |
| - `\operatorname{prec}(u_{i + 1}) = (i + 1) N_0 + \frac{i(i + 1)}{2} \cdot M + 53` |
