1#!/usr/bin/env python3
2#
3# Copyright 2022 Google LLC
4#
5# Licensed under the Apache License, Version 2.0 (the "License");
6# you may not use this file except in compliance with the License.
7# You may obtain a copy of the License at
8#
9#     http://www.apache.org/licenses/LICENSE-2.0
10#
11# Unless required by applicable law or agreed to in writing, software
12# distributed under the License is distributed on an "AS IS" BASIS,
13# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14# See the License for the specific language governing permissions and
15# limitations under the License.
16#
17
18import numpy as np
19import matplotlib.pyplot as plt
20
21
22def fast_exp2(x, p):
23
24    p = p.astype(np.float32)
25    x = x.astype(np.float32)
26
27    y = (((((p[0]*x) + p[1])*x + p[2])*x + p[3])*x + p[4])*x + 1
28
29    return np.power(y.astype(np.float32), 16)
30
31def approx_exp2():
32
33    x = np.arange(-8, 8, step=1e-3)
34
35    p = np.polyfit(x, ((2 ** (x/16)) - 1) / x, 4)
36    y = [ fast_exp2(x[i], p) for i in range(len(x)) ]
37    e = np.abs(y - 2**x) / (2 ** x)
38
39    print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}'
40                .format(p[0], p[1], p[2], p[3], p[4]))
41    print('Max relative error: ', np.max(e))
42    print('Max RMS error: ', np.sqrt(np.mean(e ** 2)))
43
44    if False:
45        fig, (ax1, ax2) = plt.subplots(2)
46
47        ax1.plot(x, 2**x, label='Reference')
48        ax1.plot(x, y, label='Approximation')
49        ax1.legend()
50
51        ax2.plot(x, e, label='Relative Error')
52        ax2.legend()
53
54        plt.show()
55
56
57def fast_log2(x, p):
58
59    p = p.astype(np.float32)
60    x = x.astype(np.float32)
61
62    (x, e) = np.frexp(x)
63
64    y = ((((p[0]*x) + p[1])*x + p[2])*x + p[3])*x + p[4]
65
66    return (e ) + y.astype(np.float32)
67
68def approx_log2():
69
70    x = np.logspace(-1, 0, base=2, num=100)
71    p = np.polyfit(x, np.log2(x), 4)
72
73    x = np.logspace(-2, 5, num=10000)
74    y = [ fast_log2(x[i], p) for i in range(len(x)) ]
75    e = np.abs(y - np.log2(x))
76
77    print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}'
78                .format(p[0], p[1], p[2], p[3], p[4]))
79    print('Max absolute error: ', np.max(e))
80    print('Max RMS error: ', np.sqrt(np.mean(e ** 2)))
81
82    if False:
83        fig, (ax1, ax2) = plt.subplots(2)
84
85        ax1.plot(x, np.log2(x),  label='Reference')
86        ax1.plot(x, y, label='Approximation')
87        ax1.legend()
88
89        ax2.plot(x, e, label = 'Absolute error')
90        ax2.legend()
91
92        plt.show()
93
94
95def table_db_q16():
96
97    k = 10 * np.log10(2);
98
99    for i in range(32):
100        a = k * np.log2(np.ldexp(32 + i  , -5)) - (i // 16) * (k/2);
101        b = k * np.log2(np.ldexp(32 + i+1, -5)) - (i // 16) * (k/2);
102
103        an = np.ldexp(a, 15) + 0.5
104        bn = np.ldexp(b - a, 15) + 0.5
105        print('{{ {:5d}, {:4d} }},'
106            .format(int(np.ldexp(a, 15) + 0.5),
107                    int(np.ldexp(b - a, 15) + 0.5)),
108            end = ' ' if i % 4 < 3 else '\n')
109
110
111if __name__ == '__main__':
112
113    print('\n--- Approximation of 2^n ---')
114    approx_exp2()
115
116    print('\n--- Approximation of log2(n) ---')
117    approx_log2()
118
119    print('\n--- Table of fixed Q16 dB ---')
120    table_db_q16()
121
122    print('')
123