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