Skip to content

Commit dd4199b

Browse files
Nikhil0250copybara-github
Nikhil0250
authored and
copybara-github
committed
Improve FastLog10()
PiperOrigin-RevId: 909989345
1 parent ef90e1b commit dd4199b

2 files changed

Lines changed: 154 additions & 103 deletions

File tree

‎hwy/contrib/math/fast_math-inl.h‎

Lines changed: 152 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -1614,7 +1614,8 @@ HWY_INLINE V FastLog2(D d, V x) {
16141614
* Fast approximation of log10(x).
16151615
*
16161616
* Valid Lane Types: float32, float64
1617-
* Max Relative Error: 0.008%
1617+
* Max Relative Error: 0.0012%
1618+
* Average Relative Error: 5.4e-6% for float32, 1.8e-7% for float64
16181619
* Valid Range: float32: (0, +FLT_MAX]
16191620
* float64: (0, +DBL_MAX]
16201621
*
@@ -1628,52 +1629,69 @@ HWY_INLINE V FastLog10(D d, V x) {
16281629
impl::FastLogRangeReduction<kHandleSubnormals>(d, x, y, exp);
16291630

16301631
constexpr size_t kLanes = HWY_MAX_LANES_D(D);
1631-
V b, c, d_coef;
1632+
V approx;
1633+
1634+
V a, b, c, d_val;
1635+
// Centering the approximation around y=1.0 by using z = y - 1.0 significantly
1636+
// improves accuracy for low-degree polynomials compared to approximating
1637+
// log(y) directly.
1638+
const V z = Sub(y, Set(d, static_cast<T>(1.0)));
16321639

16331640
if constexpr ((kLanes >= 4 && !HWY_HAVE_SCALABLE) ||
16341641
(HWY_HAVE_SCALABLE && sizeof(T) == 4 && detail::IsFull(d))) {
1642+
// --- Table Lookup ---
16351643
const auto scale = Set(d, static_cast<T>(11.3137085));
16361644
auto idx_i = ConvertInRangeTo(
16371645
RebindToSigned<D>(), MulAdd(y, scale, Set(d, static_cast<T>(-8.0))));
16381646

16391647
idx_i = Min(idx_i, Set(RebindToSigned<D>(), 7));
16401648

1649+
HWY_ALIGN static constexpr T arr_a[8] = {
1650+
static_cast<T>(0.3420769122219194),
1651+
static_cast<T>(0.24492300439548012),
1652+
static_cast<T>(0.18134323902060331),
1653+
static_cast<T>(0.137999443831595),
1654+
static_cast<T>(0.10743977319122217),
1655+
static_cast<T>(0.08527746618951795),
1656+
static_cast<T>(0.06881513239598655),
1657+
static_cast<T>(0.05634946779806663)};
1658+
16411659
HWY_ALIGN static constexpr T arr_b[8] = {
1642-
static_cast<T>(-1.00194730895928918),
1643-
static_cast<T>(-1.00042661239958708),
1644-
static_cast<T>(-1.0000255203465902),
1645-
static_cast<T>(-1),
1646-
static_cast<T>(-0.999929163668789478),
1647-
static_cast<T>(-0.999558969823431065),
1648-
static_cast<T>(-0.998743736501089163),
1649-
static_cast<T>(-0.997397894886509873)};
1650-
1651-
// Denominator coefficients are pre-multiplied by Ln(10) compared to FastLog
1652-
// to save a multiply instruction in the final step (preabsorption).
1660+
static_cast<T>(-0.1301481748440477),
1661+
static_cast<T>(-0.1906121071166758),
1662+
static_cast<T>(-0.21326579956586073),
1663+
static_cast<T>(-0.21718475171279292),
1664+
static_cast<T>(-0.21182605282145175),
1665+
static_cast<T>(-0.20205188075390862),
1666+
static_cast<T>(-0.19041912046596485),
1667+
static_cast<T>(-0.17830276936785788)};
1668+
16531669
HWY_ALIGN static constexpr T arr_c[8] = {
1654-
static_cast<T>(1.344377870361103122),
1655-
static_cast<T>(1.262218466064299438),
1656-
static_cast<T>(1.196453330732336395),
1657-
static_cast<T>(1.145230398439557540),
1658-
static_cast<T>(1.096932375640011115),
1659-
static_cast<T>(1.058095578246064594),
1660-
static_cast<T>(1.024424088002112043),
1661-
static_cast<T>(0.994880219820938994)};
1670+
static_cast<T>(0.44984736360963107),
1671+
static_cast<T>(0.4372109146505034),
1672+
static_cast<T>(0.4344592024931514),
1673+
static_cast<T>(0.4342992744270672),
1674+
static_cast<T>(0.4339565143878306),
1675+
static_cast<T>(0.4324981679587344),
1676+
static_cast<T>(0.4297421965958291),
1677+
static_cast<T>(0.4258045623390324)};
16621678

16631679
HWY_ALIGN static constexpr T arr_d[8] = {
1664-
static_cast<T>(1.008283402680355545),
1665-
static_cast<T>(1.058402359620750799),
1666-
static_cast<T>(1.109141922557689730),
1667-
static_cast<T>(1.157219973777604327),
1668-
static_cast<T>(1.210980553414537253),
1669-
static_cast<T>(1.261698563555383457),
1670-
static_cast<T>(1.312152899034920495),
1671-
static_cast<T>(1.362295845906852376)};
1680+
static_cast<T>(0.0010024790317717039),
1681+
static_cast<T>(0.00011601128161763332),
1682+
static_cast<T>(2.364622797584052e-06),
1683+
static_cast<T>(0.0),
1684+
static_cast<T>(7.886058205291105e-06),
1685+
static_cast<T>(8.148875978935532e-05),
1686+
static_cast<T>(0.00030040287702038374),
1687+
static_cast<T>(0.0007282733341565754)};
16721688

1689+
a = Lookup8(d, arr_a, idx_i);
16731690
b = Lookup8(d, arr_b, idx_i);
16741691
c = Lookup8(d, arr_c, idx_i);
1675-
d_coef = Lookup8(d, arr_d, idx_i);
1692+
d_val = Lookup8(d, arr_d, idx_i);
16761693
} else {
1694+
// --- FALLBACK PATH: Blend Chain ---
16771695
const auto t0 = Set(d, static_cast<T>(0.7954951287634819));
16781696
const auto t1 = Set(d, static_cast<T>(0.8838834764038688));
16791697
const auto t2 = Set(d, static_cast<T>(0.9722718240442556));
@@ -1683,125 +1701,157 @@ HWY_INLINE V FastLog10(D d, V x) {
16831701
const auto t6 = Set(d, static_cast<T>(1.3258252146058032));
16841702

16851703
if constexpr (HWY_REGISTERS >= 32) {
1686-
auto b_low = Set(d, static_cast<T>(-1)); // idx 3
1687-
// Denominator coefficients pre-multiplied by Ln(10).
1688-
auto c_low = Set(d, static_cast<T>(1.145230398439557540));
1689-
auto d_low = Set(d, static_cast<T>(1.157219973777604327));
1704+
// Split into two parallel chains to reduce dependency latency.
1705+
// -- Chain 1: Indices 0 to 3 (Evaluated starting from t3 down to t0)
1706+
auto a_low = Set(d, static_cast<T>(0.137999443831595));
1707+
auto b_low = Set(d, static_cast<T>(-0.21718475171279292));
1708+
auto c_low = Set(d, static_cast<T>(0.4342992744270672));
1709+
auto d_low = Set(d, static_cast<T>(0.0));
16901710

16911711
auto mask = Lt(y, t2);
1712+
a_low =
1713+
IfThenElse(mask, Set(d, static_cast<T>(0.18134323902060331)), a_low);
16921714
b_low =
1693-
IfThenElse(mask, Set(d, static_cast<T>(-1.0000255203465902)), b_low);
1715+
IfThenElse(mask, Set(d, static_cast<T>(-0.21326579956586073)), b_low);
16941716
c_low =
1695-
IfThenElse(mask, Set(d, static_cast<T>(1.196453330732336395)), c_low);
1696-
d_low =
1697-
IfThenElse(mask, Set(d, static_cast<T>(1.109141922557689730)), d_low);
1717+
IfThenElse(mask, Set(d, static_cast<T>(0.4344592024931514)), c_low);
1718+
d_low = IfThenElse(mask, Set(d, static_cast<T>(2.364622797584052e-06)),
1719+
d_low);
16981720

16991721
mask = Lt(y, t1);
1722+
a_low =
1723+
IfThenElse(mask, Set(d, static_cast<T>(0.24492300439548012)), a_low);
17001724
b_low =
1701-
IfThenElse(mask, Set(d, static_cast<T>(-1.00042661239958708)), b_low);
1725+
IfThenElse(mask, Set(d, static_cast<T>(-0.1906121071166758)), b_low);
17021726
c_low =
1703-
IfThenElse(mask, Set(d, static_cast<T>(1.262218466064299438)), c_low);
1704-
d_low =
1705-
IfThenElse(mask, Set(d, static_cast<T>(1.058402359620750799)), d_low);
1727+
IfThenElse(mask, Set(d, static_cast<T>(0.4372109146505034)), c_low);
1728+
d_low = IfThenElse(mask, Set(d, static_cast<T>(0.00011601128161763332)),
1729+
d_low);
17061730

17071731
mask = Lt(y, t0);
1732+
a_low =
1733+
IfThenElse(mask, Set(d, static_cast<T>(0.3420769122219194)), a_low);
17081734
b_low =
1709-
IfThenElse(mask, Set(d, static_cast<T>(-1.00194730895928918)), b_low);
1735+
IfThenElse(mask, Set(d, static_cast<T>(-0.1301481748440477)), b_low);
17101736
c_low =
1711-
IfThenElse(mask, Set(d, static_cast<T>(1.344377870361103122)), c_low);
1712-
d_low =
1713-
IfThenElse(mask, Set(d, static_cast<T>(1.008283402680355545)), d_low);
1737+
IfThenElse(mask, Set(d, static_cast<T>(0.44984736360963107)), c_low);
1738+
d_low = IfThenElse(mask, Set(d, static_cast<T>(0.0010024790317717039)),
1739+
d_low);
17141740

1715-
auto b_high = Set(d, static_cast<T>(-0.997397894886509873)); // idx 7
1716-
// Denominator coefficients pre-multiplied by Ln(10).
1717-
auto c_high = Set(d, static_cast<T>(0.994880219820938994));
1718-
auto d_high = Set(d, static_cast<T>(1.362295845906852376));
1741+
// -- Chain 2: Indices 4 to 7 (Evaluated starting from t6 down to t4)
1742+
auto a_high = Set(d, static_cast<T>(0.05634946779806663));
1743+
auto b_high = Set(d, static_cast<T>(-0.17830276936785788));
1744+
auto c_high = Set(d, static_cast<T>(0.4258045623390324));
1745+
auto d_high = Set(d, static_cast<T>(0.0007282733341565754));
17191746

17201747
mask = Lt(y, t6);
1721-
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.998743736501089163)),
1748+
a_high =
1749+
IfThenElse(mask, Set(d, static_cast<T>(0.06881513239598655)), a_high);
1750+
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.19041912046596485)),
17221751
b_high);
1723-
c_high = IfThenElse(mask, Set(d, static_cast<T>(1.024424088002112043)),
1724-
c_high);
1725-
d_high = IfThenElse(mask, Set(d, static_cast<T>(1.312152899034920495)),
1752+
c_high =
1753+
IfThenElse(mask, Set(d, static_cast<T>(0.4297421965958291)), c_high);
1754+
d_high = IfThenElse(mask, Set(d, static_cast<T>(0.00030040287702038374)),
17261755
d_high);
17271756

17281757
mask = Lt(y, t5);
1729-
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.999558969823431065)),
1758+
a_high =
1759+
IfThenElse(mask, Set(d, static_cast<T>(0.08527746618951795)), a_high);
1760+
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.20205188075390862)),
17301761
b_high);
1731-
c_high = IfThenElse(mask, Set(d, static_cast<T>(1.058095578246064594)),
1732-
c_high);
1733-
d_high = IfThenElse(mask, Set(d, static_cast<T>(1.261698563555383457)),
1762+
c_high =
1763+
IfThenElse(mask, Set(d, static_cast<T>(0.4324981679587344)), c_high);
1764+
d_high = IfThenElse(mask, Set(d, static_cast<T>(8.148875978935532e-05)),
17341765
d_high);
17351766

17361767
mask = Lt(y, t4);
1737-
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.999929163668789478)),
1768+
a_high =
1769+
IfThenElse(mask, Set(d, static_cast<T>(0.10743977319122217)), a_high);
1770+
b_high = IfThenElse(mask, Set(d, static_cast<T>(-0.21182605282145175)),
17381771
b_high);
1739-
c_high = IfThenElse(mask, Set(d, static_cast<T>(1.096932375640011115)),
1740-
c_high);
1741-
d_high = IfThenElse(mask, Set(d, static_cast<T>(1.210980553414537253)),
1772+
c_high =
1773+
IfThenElse(mask, Set(d, static_cast<T>(0.4339565143878306)), c_high);
1774+
d_high = IfThenElse(mask, Set(d, static_cast<T>(7.886058205291105e-06)),
17421775
d_high);
17431776

1777+
// -- Merge the two chains
17441778
auto merge_mask = Lt(y, t3);
1779+
a = IfThenElse(merge_mask, a_low, a_high);
17451780
b = IfThenElse(merge_mask, b_low, b_high);
17461781
c = IfThenElse(merge_mask, c_low, c_high);
1747-
d_coef = IfThenElse(merge_mask, d_low, d_high);
1782+
d_val = IfThenElse(merge_mask, d_low, d_high);
17481783
} else {
1749-
b = Set(d, static_cast<T>(-0.997397894886509873));
1750-
// Denominator coefficients pre-multiplied by Ln(10).
1751-
c = Set(d, static_cast<T>(0.994880219820938994));
1752-
d_coef = Set(d, static_cast<T>(1.362295845906852376));
1784+
// Start with highest index (7)
1785+
a = Set(d, static_cast<T>(0.05634946779806663));
1786+
b = Set(d, static_cast<T>(-0.17830276936785788));
1787+
c = Set(d, static_cast<T>(0.4258045623390324));
1788+
d_val = Set(d, static_cast<T>(0.0007282733341565754));
17531789

1790+
// If y < t6 (idx 6)
17541791
auto mask = Lt(y, t6);
1755-
b = IfThenElse(mask, Set(d, static_cast<T>(-0.998743736501089163)), b);
1756-
c = IfThenElse(mask, Set(d, static_cast<T>(1.024424088002112043)), c);
1757-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.312152899034920495)),
1758-
d_coef);
1792+
a = IfThenElse(mask, Set(d, static_cast<T>(0.06881513239598655)), a);
1793+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.19041912046596485)), b);
1794+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4297421965958291)), c);
1795+
d_val = IfThenElse(mask, Set(d, static_cast<T>(0.00030040287702038374)),
1796+
d_val);
17591797

1798+
// If y < t5 (idx 5)
17601799
mask = Lt(y, t5);
1761-
b = IfThenElse(mask, Set(d, static_cast<T>(-0.999558969823431065)), b);
1762-
c = IfThenElse(mask, Set(d, static_cast<T>(1.058095578246064594)), c);
1763-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.261698563555383457)),
1764-
d_coef);
1800+
a = IfThenElse(mask, Set(d, static_cast<T>(0.08527746618951795)), a);
1801+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.20205188075390862)), b);
1802+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4324981679587344)), c);
1803+
d_val = IfThenElse(mask, Set(d, static_cast<T>(8.148875978935532e-05)),
1804+
d_val);
17651805

1806+
// If y < t4 (idx 4)
17661807
mask = Lt(y, t4);
1767-
b = IfThenElse(mask, Set(d, static_cast<T>(-0.999929163668789478)), b);
1768-
c = IfThenElse(mask, Set(d, static_cast<T>(1.096932375640011115)), c);
1769-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.210980553414537253)),
1770-
d_coef);
1808+
a = IfThenElse(mask, Set(d, static_cast<T>(0.10743977319122217)), a);
1809+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.21182605282145175)), b);
1810+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4339565143878306)), c);
1811+
d_val = IfThenElse(mask, Set(d, static_cast<T>(7.886058205291105e-06)),
1812+
d_val);
17711813

1814+
// If y < t3 (idx 3)
17721815
mask = Lt(y, t3);
1773-
b = IfThenElse(mask, Set(d, static_cast<T>(-1)), b);
1774-
c = IfThenElse(mask, Set(d, static_cast<T>(1.145230398439557540)), c);
1775-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.157219973777604327)),
1776-
d_coef);
1816+
a = IfThenElse(mask, Set(d, static_cast<T>(0.137999443831595)), a);
1817+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.21718475171279292)), b);
1818+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4342992744270672)), c);
1819+
d_val = IfThenElse(mask, Set(d, static_cast<T>(0.0)), d_val);
17771820

1821+
// If y < t2 (idx 2)
17781822
mask = Lt(y, t2);
1779-
b = IfThenElse(mask, Set(d, static_cast<T>(-1.0000255203465902)), b);
1780-
c = IfThenElse(mask, Set(d, static_cast<T>(1.196453330732336395)), c);
1781-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.109141922557689730)),
1782-
d_coef);
1823+
a = IfThenElse(mask, Set(d, static_cast<T>(0.18134323902060331)), a);
1824+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.21326579956586073)), b);
1825+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4344592024931514)), c);
1826+
d_val = IfThenElse(mask, Set(d, static_cast<T>(2.364622797584052e-06)),
1827+
d_val);
17831828

1829+
// If y < t1 (idx 1)
17841830
mask = Lt(y, t1);
1785-
b = IfThenElse(mask, Set(d, static_cast<T>(-1.00042661239958708)), b);
1786-
c = IfThenElse(mask, Set(d, static_cast<T>(1.262218466064299438)), c);
1787-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.058402359620750799)),
1788-
d_coef);
1831+
a = IfThenElse(mask, Set(d, static_cast<T>(0.24492300439548012)), a);
1832+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.1906121071166758)), b);
1833+
c = IfThenElse(mask, Set(d, static_cast<T>(0.4372109146505034)), c);
1834+
d_val = IfThenElse(mask, Set(d, static_cast<T>(0.00011601128161763332)),
1835+
d_val);
17891836

1837+
// If y < t0 (idx 0)
17901838
mask = Lt(y, t0);
1791-
b = IfThenElse(mask, Set(d, static_cast<T>(-1.00194730895928918)), b);
1792-
c = IfThenElse(mask, Set(d, static_cast<T>(1.344377870361103122)), c);
1793-
d_coef = IfThenElse(mask, Set(d, static_cast<T>(1.008283402680355545)),
1794-
d_coef);
1839+
a = IfThenElse(mask, Set(d, static_cast<T>(0.3420769122219194)), a);
1840+
b = IfThenElse(mask, Set(d, static_cast<T>(-0.1301481748440477)), b);
1841+
c = IfThenElse(mask, Set(d, static_cast<T>(0.44984736360963107)), c);
1842+
d_val = IfThenElse(mask, Set(d, static_cast<T>(0.0010024790317717039)),
1843+
d_val);
17951844
}
17961845
}
1797-
1798-
auto num = Add(y, b);
1799-
auto den = MulAdd(c, y, d_coef);
1800-
auto approx = Div(num, den);
1846+
// Math: approx = (a*z + b)*z^2 + (c*z + d_val)
1847+
const auto z2 = Mul(z, z);
1848+
const auto pab = MulAdd(a, z, b);
1849+
const auto pcd = MulAdd(c, z, d_val);
1850+
approx = MulAdd(pab, z2, pcd);
18011851

18021852
const auto kLog10_2 = Set(d, static_cast<T>(0.3010299956639812)); // log10(2)
18031853
// Computes exp * log10(2) + approx. Since approx was scaled by 1/Ln(10)
1804-
// via the pre-scaled denominator, this yields the correct log10 result
1854+
// via the pre-scaled coefficients, this yields the correct log10 result
18051855
// using a single MulAdd instruction.
18061856
return MulAdd(exp, kLog10_2, approx);
18071857
}
@@ -1810,7 +1860,8 @@ HWY_INLINE V FastLog10(D d, V x) {
18101860
* Fast approximation of log(1 + x).
18111861
*
18121862
* Valid Lane Types: float32, float64
1813-
* Max Relative Error: 0.0081%
1863+
* Max Relative Error: 0.0012%
1864+
* Average Relative Error: 0.00013% for float32, 0.000039% for float64
18141865
* Valid Range: float32: [-1 + epsilon, +FLT_MAX]
18151866
* float64: [-1 + epsilon, +DBL_MAX]
18161867
*

‎hwy/contrib/math/math_test.cc‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -427,7 +427,7 @@ struct TestFastLog2 {
427427
struct TestFastLog10 {
428428
template <class T, class D>
429429
HWY_NOINLINE void operator()(T, D d) {
430-
const double max_relative_error = 0.000082;
430+
const double max_relative_error = 1.15E-5;
431431
const uint64_t samples = 1000000;
432432
if (sizeof(T) == 4) {
433433
TestMathRelative<T, D>("FastLog10", std::log10, CallFastLog10, d,
@@ -452,7 +452,7 @@ struct TestFastLog10 {
452452
struct TestFastLog1p {
453453
template <class T, class D>
454454
HWY_NOINLINE void operator()(T, D d) {
455-
const double max_relative_error = 0.00009;
455+
const double max_relative_error = 1.15E-5;
456456
const uint64_t samples = 1000000;
457457
if (sizeof(T) == 4) {
458458
TestMathRelative<T, D>("FastLog1p", std::log1p, CallFastLog1p, d,

0 commit comments

Comments
 (0)