11//
22// code for doing exponentials and logarithms
33// Copyright 2020 Total Spectrum Software Inc.
4- // MIT licensed
4+ // SPDX-License-Identifier: MIT
55//
66
77#define CONST_E 2.71828182846f
@@ -18,83 +18,6 @@ typedef union ForU {
1818static inline unsigned __asuint (float f ) { ForU u ; u .f = f ; return u .u ; }
1919static inline float __asfloat (unsigned d ) { ForU u ; u .u = d ; return u .f ; }
2020
21- #ifdef __P2__
22-
23- // calculate log2(mant) as a 5.27 fixed point number
24- unsigned __builtin_qlog (unsigned mant )
25- {
26- unsigned r ;
27- __asm {
28- qlog mant
29- getqx r
30- };
31- return r ;
32- }
33-
34- // calculate 2^mant, where mant is a 5.27 fixed point number
35- unsigned __builtin_qexp (unsigned mant )
36- {
37- unsigned r ;
38- __asm {
39- qexp mant
40- getqx r
41- };
42- return r ;
43- }
44- #else
45- // calculate log base 2 of an unsigned num as a 5.27 fixed point number
46- // uses the lookup table in ROM
47- // algorithm comes from Propeller manual
48- unsigned __builtin_qlog (unsigned num )
49- {
50- unsigned exp ;
51- unsigned r ;
52- unsigned r2 ;
53- unsigned frac ;
54- exp = __builtin_clz (num );
55- num = num << exp ; // left justify
56- // we will look up based on 11 bits, so the bottom 16 bits are the fraction
57- frac = num & 0xffff ;
58- num = (num >> (30 - 11 )) & 0xffe ; // 30 because we will multiply by 2 for word access
59- r = * ((unsigned short * )(0xC000 + num ));
60- r2 = * ((unsigned short * )(0xC002 + num ));
61- r = r + (((r2 - r )* frac ) >> 16 );
62- exp = exp ^ 0x1f ;
63- r = r | (exp << 16 );
64-
65- r = r << (27 - 16 );
66- return r ;
67- }
68- // calculate 2^num, where "num" is a 5.27 fixed point num
69- // uses the lookup table in ROM
70- // algorithm comes from Propeller manual
71- unsigned __builtin_qexp (unsigned orig_num )
72- {
73- unsigned exp ;
74- unsigned r , r2 ;
75- unsigned frac ;
76- unsigned num = orig_num ;
77-
78- exp = (num >> 27 );
79- // the original value is 5.27
80- // we will look up based on 11 bits, so the bottom 16 bits are the fraction
81- frac = num & 0xffff ;
82- num = (num >> (26 - 11 )) & 0x0ffe ;
83- //__builtin_printf("...num=0x%08x lookup=%08x exp=%d\n", orig_num, num, exp);
84- r = * ((unsigned short * )(0xD000 + num ));
85- r2 = * ((unsigned short * )(0xD002 + num ));
86- //__builtin_printf("...r1=0x%08x r2=%08x frac=%04x\n", r, r2, frac);
87- r = r + ((frac * (r2 - r )) >> 16 );
88- r = r << 15 ; // shift into 30..15
89- r |= 0x80000000 ;
90- //printf("...r=0x%08x\n", r);
91- exp ^= 0x1f ;
92- r = r >> exp ;
93- return r ;
94- }
95-
96- #endif
97-
9821//
9922// floating point: calculate log base 2 of a float number
10023//
@@ -165,7 +88,7 @@ float __builtin_log2f(float x)
16588 }
16689 }
16790#endif
168- r = __builtin_qlog (mant );
91+ r = _qlog (mant );
16992 // at this point r is a 5.27 fixed point number giving log2(mant)
17093 // note that mant was 1.23 by construction, so the upper 5 bits
17194 // is "24", i.e. we have 0xbnnnnnn
@@ -213,7 +136,7 @@ float __builtin_exp2f(float x)
213136 //printf("... u=%08x\n", u);
214137 u = u <<11 ; // 5.27
215138 u |= (16 <<27 );
216- r = __builtin_qexp (u ); // as a 16.16 number
139+ r = _qexp (u ); // as a 16.16 number
217140 //printf("... r=%08x\n", r);
218141
219142 if (n >= 0 ) {
@@ -259,7 +182,7 @@ float __builtin_exp2f(float x)
259182 }
260183 u |= (24 <<27 );
261184 //printf("..u=%x (input) ", u);
262- r = __builtin_qexp (u );
185+ r = _qexp (u );
263186 //printf("..u=%x r=%x\n", u, r);
264187 // round:
265188 r = (r + 1 )>>1 ;
0 commit comments