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, t, p): 23 24 p = p.astype(np.float32) 25 x = x.astype(np.float32) 26 27 m = ((x + 0.5/8) % (1/8)) - (0.5/8) 28 e = int((x - m) * 8) 29 30 y = ((((p[0]*m) + p[1])*m + p[2])*m + p[3])*m + p[4] 31 y = y * 2**(e // 8) * t[e % 8] 32 33 return y.astype(np.float32) 34 35def approx_exp2(): 36 37 x = np.arange(0, 1/8, step=1e-6) 38 p = np.polyfit(x, 2 ** x, 4) 39 t = [ 2**(i/8) for i in range(8) ] 40 41 x = np.arange(-10, 10, step=1e-3) 42 y = [ fast_exp2(x[i], t, p) for i in range(len(x)) ] 43 44 e = np.abs(y - 2**x) / (2 ** x) 45 46 print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, \n' 47 ' {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, '.format(*t)) 48 print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}'.format(*p)) 49 print('Max relative error: ', np.max(e)) 50 print('Max RMS error: ', np.sqrt(np.mean(e ** 2))) 51 52 if False: 53 fig, (ax1, ax2) = plt.subplots(2) 54 55 ax1.plot(x, 2**x, label='Reference') 56 ax1.plot(x, y, label='Approximation') 57 ax1.legend() 58 59 ax2.plot(x, e, label='Relative Error') 60 ax2.legend() 61 62 plt.show() 63 64 65def fast_log2(x, p): 66 67 p = p.astype(np.float32) 68 x = x.astype(np.float32) 69 70 (x, e) = np.frexp(x) 71 72 y = ((((p[0]*x) + p[1])*x + p[2])*x + p[3])*x + p[4] 73 74 return (e ) + y.astype(np.float32) 75 76def approx_log2(): 77 78 x = np.logspace(-1, 0, base=2, num=100) 79 p = np.polyfit(x, np.log2(x), 4) 80 81 x = np.logspace(-2, 5, num=10000) 82 y = [ fast_log2(x[i], p) for i in range(len(x)) ] 83 e = np.abs(y - np.log2(x)) 84 85 print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}' 86 .format(p[0], p[1], p[2], p[3], p[4])) 87 print('Max absolute error: ', np.max(e)) 88 print('Max RMS error: ', np.sqrt(np.mean(e ** 2))) 89 90 if False: 91 fig, (ax1, ax2) = plt.subplots(2) 92 93 ax1.plot(x, np.log2(x), label='Reference') 94 ax1.plot(x, y, label='Approximation') 95 ax1.legend() 96 97 ax2.plot(x, e, label = 'Absolute error') 98 ax2.legend() 99 100 plt.show() 101 102 103def table_db_q16(): 104 105 k = 10 * np.log10(2); 106 107 for i in range(32): 108 a = k * np.log2(np.ldexp(32 + i , -5)) - (i // 16) * (k/2); 109 b = k * np.log2(np.ldexp(32 + i+1, -5)) - (i // 16) * (k/2); 110 111 an = np.ldexp(a, 15) + 0.5 112 bn = np.ldexp(b - a, 15) + 0.5 113 print('{{ {:5d}, {:4d} }},' 114 .format(int(np.ldexp(a, 15) + 0.5), 115 int(np.ldexp(b - a, 15) + 0.5)), 116 end = ' ' if i % 4 < 3 else '\n') 117 118 119if __name__ == '__main__': 120 121 print('\n--- Approximation of 2^n ---') 122 approx_exp2() 123 124 print('\n--- Approximation of log2(n) ---') 125 approx_log2() 126 127 print('\n--- Table of fixed Q16 dB ---') 128 table_db_q16() 129 130 print('') 131