1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
|
/*
* Copyright (c) 2010 Broadcom Corporation
*
* Permission to use, copy, modify, and/or distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
* SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
* OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
* CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
#include "phy_qmath.h"
/*
* Description: This function make 16 bit unsigned multiplication.
* To fit the output into 16 bits the 32 bit multiplication result is right
* shifted by 16 bits.
*/
u16 qm_mulu16(u16 op1, u16 op2)
{
return (u16) (((u32) op1 * (u32) op2) >> 16);
}
/*
* Description: This function make 16 bit multiplication and return the result
* in 16 bits. To fit the multiplication result into 16 bits the multiplication
* result is right shifted by 15 bits. Right shifting 15 bits instead of 16 bits
* is done to remove the extra sign bit formed due to the multiplication.
* When both the 16bit inputs are 0x8000 then the output is saturated to
* 0x7fffffff.
*/
s16 qm_muls16(s16 op1, s16 op2)
{
s32 result;
if (op1 == (s16) 0x8000 && op2 == (s16) 0x8000)
result = 0x7fffffff;
else
result = ((s32) (op1) * (s32) (op2));
return (s16) (result >> 15);
}
/*
* Description: This function add two 32 bit numbers and return the 32bit
* result. If the result overflow 32 bits, the output will be saturated to
* 32bits.
*/
s32 qm_add32(s32 op1, s32 op2)
{
s32 result;
result = op1 + op2;
if (op1 < 0 && op2 < 0 && result > 0)
result = 0x80000000;
else if (op1 > 0 && op2 > 0 && result < 0)
result = 0x7fffffff;
return result;
}
/*
* Description: This function add two 16 bit numbers and return the 16bit
* result. If the result overflow 16 bits, the output will be saturated to
* 16bits.
*/
s16 qm_add16(s16 op1, s16 op2)
{
s16 result;
s32 temp = (s32) op1 + (s32) op2;
if (temp > (s32) 0x7fff)
result = (s16) 0x7fff;
else if (temp < (s32) 0xffff8000)
result = (s16) 0xffff8000;
else
result = (s16) temp;
return result;
}
/*
* Description: This function make 16 bit subtraction and return the 16bit
* result. If the result overflow 16 bits, the output will be saturated to
* 16bits.
*/
s16 qm_sub16(s16 op1, s16 op2)
{
s16 result;
s32 temp = (s32) op1 - (s32) op2;
if (temp > (s32) 0x7fff)
result = (s16) 0x7fff;
else if (temp < (s32) 0xffff8000)
result = (s16) 0xffff8000;
else
result = (s16) temp;
return result;
}
/*
* Description: This function make a 32 bit saturated left shift when the
* specified shift is +ve. This function will make a 32 bit right shift when
* the specified shift is -ve. This function return the result after shifting
* operation.
*/
s32 qm_shl32(s32 op, int shift)
{
int i;
s32 result;
result = op;
if (shift > 31)
shift = 31;
else if (shift < -31)
shift = -31;
if (shift >= 0) {
for (i = 0; i < shift; i++)
result = qm_add32(result, result);
} else {
result = result >> (-shift);
}
return result;
}
/*
* Description: This function make a 16 bit saturated left shift when the
* specified shift is +ve. This function will make a 16 bit right shift when
* the specified shift is -ve. This function return the result after shifting
* operation.
*/
s16 qm_shl16(s16 op, int shift)
{
int i;
s16 result;
result = op;
if (shift > 15)
shift = 15;
else if (shift < -15)
shift = -15;
if (shift > 0) {
for (i = 0; i < shift; i++)
result = qm_add16(result, result);
} else {
result = result >> (-shift);
}
return result;
}
/*
* Description: This function make a 16 bit right shift when shift is +ve.
* This function make a 16 bit saturated left shift when shift is -ve. This
* function return the result of the shift operation.
*/
s16 qm_shr16(s16 op, int shift)
{
return qm_shl16(op, -shift);
}
/*
* Description: This function return the number of redundant sign bits in a
* 32 bit number. Example: qm_norm32(0x00000080) = 23
*/
s16 qm_norm32(s32 op)
{
u16 u16extraSignBits;
if (op == 0) {
return 31;
} else {
u16extraSignBits = 0;
while ((op >> 31) == (op >> 30)) {
u16extraSignBits++;
op = op << 1;
}
}
return u16extraSignBits;
}
/* This table is log2(1+(i/32)) where i=[0:1:32], in q.15 format */
static const s16 log_table[] = {
0,
1455,
2866,
4236,
5568,
6863,
8124,
9352,
10549,
11716,
12855,
13968,
15055,
16117,
17156,
18173,
19168,
20143,
21098,
22034,
22952,
23852,
24736,
25604,
26455,
27292,
28114,
28922,
29717,
30498,
31267,
32024,
32767
};
#define LOG_TABLE_SIZE 32 /* log_table size */
#define LOG2_LOG_TABLE_SIZE 5 /* log2(log_table size) */
#define Q_LOG_TABLE 15 /* qformat of log_table */
#define LOG10_2 19728 /* log10(2) in q.16 */
/*
* Description:
* This routine takes the input number N and its q format qN and compute
* the log10(N). This routine first normalizes the input no N. Then N is in
* mag*(2^x) format. mag is any number in the range 2^30-(2^31 - 1).
* Then log2(mag * 2^x) = log2(mag) + x is computed. From that
* log10(mag * 2^x) = log2(mag * 2^x) * log10(2) is computed.
* This routine looks the log2 value in the table considering
* LOG2_LOG_TABLE_SIZE+1 MSBs. As the MSB is always 1, only next
* LOG2_OF_LOG_TABLE_SIZE MSBs are used for table lookup. Next 16 MSBs are used
* for interpolation.
* Inputs:
* N - number to which log10 has to be found.
* qN - q format of N
* log10N - address where log10(N) will be written.
* qLog10N - address where log10N qformat will be written.
* Note/Problem:
* For accurate results input should be in normalized or near normalized form.
*/
void qm_log10(s32 N, s16 qN, s16 *log10N, s16 *qLog10N)
{
s16 s16norm, s16tableIndex, s16errorApproximation;
u16 u16offset;
s32 s32log;
/* normalize the N. */
s16norm = qm_norm32(N);
N = N << s16norm;
/* The qformat of N after normalization.
* -30 is added to treat the no as between 1.0 to 2.0
* i.e. after adding the -30 to the qformat the decimal point will be
* just rigtht of the MSB. (i.e. after sign bit and 1st MSB). i.e.
* at the right side of 30th bit.
*/
qN = qN + s16norm - 30;
/* take the table index as the LOG2_OF_LOG_TABLE_SIZE bits right of the
* MSB */
s16tableIndex = (s16) (N >> (32 - (2 + LOG2_LOG_TABLE_SIZE)));
/* remove the MSB. the MSB is always 1 after normalization. */
s16tableIndex =
s16tableIndex & (s16) ((1 << LOG2_LOG_TABLE_SIZE) - 1);
/* remove the (1+LOG2_OF_LOG_TABLE_SIZE) MSBs in the N. */
N = N & ((1 << (32 - (2 + LOG2_LOG_TABLE_SIZE))) - 1);
/* take the offset as the 16 MSBS after table index.
*/
u16offset = (u16) (N >> (32 - (2 + LOG2_LOG_TABLE_SIZE + 16)));
/* look the log value in the table. */
s32log = log_table[s16tableIndex]; /* q.15 format */
/* interpolate using the offset. q.15 format. */
s16errorApproximation = (s16) qm_mulu16(u16offset,
(u16) (log_table[s16tableIndex + 1] -
log_table[s16tableIndex]));
/* q.15 format */
s32log = qm_add16((s16) s32log, s16errorApproximation);
/* adjust for the qformat of the N as
* log2(mag * 2^x) = log2(mag) + x
*/
s32log = qm_add32(s32log, ((s32) -qN) << 15); /* q.15 format */
/* normalize the result. */
s16norm = qm_norm32(s32log);
/* bring all the important bits into lower 16 bits */
/* q.15+s16norm-16 format */
s32log = qm_shl32(s32log, s16norm - 16);
/* compute the log10(N) by multiplying log2(N) with log10(2).
* as log10(mag * 2^x) = log2(mag * 2^x) * log10(2)
* log10N in q.15+s16norm-16+1 (LOG10_2 is in q.16)
*/
*log10N = qm_muls16((s16) s32log, (s16) LOG10_2);
/* write the q format of the result. */
*qLog10N = 15 + s16norm - 16 + 1;
return;
}
|
