pseudofloat.h

00001 /* Copyright (C) 2005 Jean-Marc Valin */
00002 /**
00003    @file pseudofloat.h
00004    @brief Pseudo-floating point
00005  * This header file provides a lightweight floating point type for
00006  * use on fixed-point platforms when a large dynamic range is 
00007  * required. The new type is not compatible with the 32-bit IEEE format,
00008  * it is not even remotely as accurate as 32-bit floats, and is not
00009  * even guaranteed to produce even remotely correct results for code
00010  * other than Speex. It makes all kinds of shortcuts that are acceptable
00011  * for Speex, but may not be acceptable for your application. You're
00012  * quite welcome to reuse this code and improve it, but don't assume
00013  * it works out of the box. Most likely, it doesn't.
00014  */
00015 /*
00016    Redistribution and use in source and binary forms, with or without
00017    modification, are permitted provided that the following conditions
00018    are met:
00019    
00020    - Redistributions of source code must retain the above copyright
00021    notice, this list of conditions and the following disclaimer.
00022    
00023    - Redistributions in binary form must reproduce the above copyright
00024    notice, this list of conditions and the following disclaimer in the
00025    documentation and/or other materials provided with the distribution.
00026    
00027    - Neither the name of the Xiph.org Foundation nor the names of its
00028    contributors may be used to endorse or promote products derived from
00029    this software without specific prior written permission.
00030    
00031    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00032    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00033    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00034    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
00035    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00036    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00037    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00038    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00039    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00040    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00041    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00042 */
00043 
00044 #ifndef PSEUDOFLOAT_H
00045 #define PSEUDOFLOAT_H
00046 
00047 #include "misc.h"
00048 #include "math_approx.h"
00049 #include <math.h>
00050 
00051 #ifdef FIXED_POINT
00052 
00053 typedef struct {
00054    spx_int16_t m;
00055    spx_int16_t e;
00056 } spx_float_t;
00057 
00058 static const spx_float_t FLOAT_ZERO = {0,0};
00059 static const spx_float_t FLOAT_ONE = {16384,-14};
00060 static const spx_float_t FLOAT_HALF = {16384,-15};
00061 
00062 #define MIN(a,b) ((a)<(b)?(a):(b))
00063 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x)
00064 {
00065    int e=0;
00066    int sign=0;
00067    if (x<0)
00068    {
00069       sign = 1;
00070       x = -x;
00071    }
00072    if (x==0)
00073    {
00074       spx_float_t r = {0,0};
00075       return r;
00076    }
00077    e = spx_ilog2(ABS32(x))-14;
00078    x = VSHR32(x, e);
00079    if (sign)
00080    {
00081       spx_float_t r;
00082       r.m = -x;
00083       r.e = e;
00084       return r;
00085    }
00086    else      
00087    {
00088       spx_float_t r;
00089       r.m = x;
00090       r.e = e;
00091       return r;
00092    }
00093 }
00094 
00095 
00096 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b)
00097 {
00098    spx_float_t r;
00099    if (a.m==0)
00100       return b;
00101    else if (b.m==0)
00102       return a;
00103    if ((a).e > (b).e) 
00104    {
00105       r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1));
00106       r.e = (a).e+1;
00107    }
00108    else 
00109    {
00110       r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1));
00111       r.e = (b).e+1;
00112    }
00113    if (r.m>0)
00114    {
00115       if (r.m<16384)
00116       {
00117          r.m<<=1;
00118          r.e-=1;
00119       }
00120    } else {
00121       if (r.m>-16384)
00122       {
00123          r.m<<=1;
00124          r.e-=1;
00125       }
00126    }
00127    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
00128    return r;
00129 }
00130 
00131 static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b)
00132 {
00133    spx_float_t r;
00134    if (a.m==0)
00135       return b;
00136    else if (b.m==0)
00137       return a;
00138    if ((a).e > (b).e)
00139    {
00140       r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1));
00141       r.e = (a).e+1;
00142    }
00143    else 
00144    {
00145       r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1);
00146       r.e = (b).e+1;
00147    }
00148    if (r.m>0)
00149    {
00150       if (r.m<16384)
00151       {
00152          r.m<<=1;
00153          r.e-=1;
00154       }
00155    } else {
00156       if (r.m>-16384)
00157       {
00158          r.m<<=1;
00159          r.e-=1;
00160       }
00161    }
00162    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
00163    return r;
00164 }
00165 
00166 static inline int FLOAT_LT(spx_float_t a, spx_float_t b)
00167 {
00168    if (a.m==0)
00169       return b.m>0;
00170    else if (b.m==0)
00171       return a.m<0;   
00172    if ((a).e > (b).e)
00173       return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1));
00174    else 
00175       return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1));
00176 
00177 }
00178 
00179 static inline int FLOAT_GT(spx_float_t a, spx_float_t b)
00180 {
00181    return FLOAT_LT(b,a);
00182 }
00183 
00184 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b)
00185 {
00186    spx_float_t r;
00187    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
00188    r.e = (a).e+(b).e+15;
00189    if (r.m>0)
00190    {
00191       if (r.m<16384)
00192       {
00193          r.m<<=1;
00194          r.e-=1;
00195       }
00196    } else {
00197       if (r.m>-16384)
00198       {
00199          r.m<<=1;
00200          r.e-=1;
00201       }
00202    }
00203    /*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
00204    return r;   
00205 }
00206 
00207 static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b)
00208 {
00209    spx_float_t r;
00210    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
00211    r.e = (a).e+(b).e+15;
00212    return r;   
00213 }
00214 
00215 
00216 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b)
00217 {
00218    spx_float_t r;
00219    r.m = a.m;
00220    r.e = a.e+b;
00221    return r;
00222 }
00223 
00224 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a)
00225 {
00226    if (a.e<0)
00227       return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e);
00228    else
00229       return a.m<<a.e;
00230 }
00231 
00232 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a)
00233 {
00234    if (a.e<0)
00235       return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e;
00236    else
00237       return EXTEND32(a.m)<<a.e;
00238 }
00239 
00240 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b)
00241 {
00242    return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15);
00243 }
00244 
00245 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b)
00246 {
00247    int e1, e2;
00248    spx_float_t r;
00249    if (a==0 || b==0)
00250    {
00251       return FLOAT_ZERO;
00252    }
00253    e1 = spx_ilog2(ABS32(a));
00254    a = VSHR32(a, e1-14);
00255    e2 = spx_ilog2(ABS32(b));
00256    b = VSHR32(b, e2-14);
00257    r.m = MULT16_16_Q15(a,b);
00258    r.e = e1+e2-13;
00259    return r;
00260 }
00261 
00262 /* Do NOT attempt to divide by a negative number */
00263 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b)
00264 {
00265    int e=0;
00266    spx_float_t r;
00267    if (a==0)
00268    {
00269       return FLOAT_ZERO;
00270    }
00271    e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15;
00272    a = VSHR32(a, e);
00273    if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15))
00274    {
00275       a >>= 1;
00276       e++;
00277    }
00278    r.m = DIV32_16(a,b.m);
00279    r.e = e-b.e;
00280    return r;
00281 }
00282 
00283 
00284 /* Do NOT attempt to divide by a negative number */
00285 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b)
00286 {
00287    int e0=0,e=0;
00288    spx_float_t r;
00289    if (a==0)
00290    {
00291       return FLOAT_ZERO;
00292    }
00293    if (b>32767)
00294    {
00295       e0 = spx_ilog2(b)-14;
00296       b = VSHR32(b, e0);
00297       e0 = -e0;
00298    }
00299    e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15;
00300    a = VSHR32(a, e);
00301    if (ABS32(a)>=SHL32(EXTEND32(b-1),15))
00302    {
00303       a >>= 1;
00304       e++;
00305    }
00306    e += e0;
00307    r.m = DIV32_16(a,b);
00308    r.e = e;
00309    return r;
00310 }
00311 
00312 /* Do NOT attempt to divide by a negative number */
00313 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b)
00314 {
00315    int e=0;
00316    spx_int32_t num;
00317    spx_float_t r;
00318    if (b.m<=0)
00319    {
00320       speex_warning_int("Attempted to divide by", b.m);
00321       return FLOAT_ONE;
00322    }
00323    num = a.m;
00324    a.m = ABS16(a.m);
00325    while (a.m >= b.m)
00326    {
00327       e++;
00328       a.m >>= 1;
00329    }
00330    num = num << (15-e);
00331    r.m = DIV32_16(num,b.m);
00332    r.e = a.e-b.e-15+e;
00333    return r;
00334 }
00335 
00336 static inline spx_float_t FLOAT_SQRT(spx_float_t a)
00337 {
00338    spx_float_t r;
00339    spx_int32_t m;
00340    m = SHL32(EXTEND32(a.m), 14);
00341    r.e = a.e - 14;
00342    if (r.e & 1)
00343    {
00344       r.e -= 1;
00345       m <<= 1;
00346    }
00347    r.e >>= 1;
00348    r.m = spx_sqrt(m);
00349    return r;
00350 }
00351 
00352 #else
00353 
00354 #define spx_float_t float
00355 #define FLOAT_ZERO 0.f
00356 #define FLOAT_ONE 1.f
00357 #define FLOAT_HALF 0.5f
00358 #define PSEUDOFLOAT(x) (x)
00359 #define FLOAT_MULT(a,b) ((a)*(b))
00360 #define FLOAT_AMULT(a,b) ((a)*(b))
00361 #define FLOAT_MUL32(a,b) ((a)*(b))
00362 #define FLOAT_DIV32(a,b) ((a)/(b))
00363 #define FLOAT_EXTRACT16(a) (a)
00364 #define FLOAT_EXTRACT32(a) (a)
00365 #define FLOAT_ADD(a,b) ((a)+(b))
00366 #define FLOAT_SUB(a,b) ((a)-(b))
00367 #define REALFLOAT(x) (x)
00368 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b))
00369 #define FLOAT_MUL32U(a,b) ((a)*(b))
00370 #define FLOAT_SHL(a,b) (a)
00371 #define FLOAT_LT(a,b) ((a)<(b))
00372 #define FLOAT_GT(a,b) ((a)>(b))
00373 #define FLOAT_DIVU(a,b) ((a)/(b))
00374 #define FLOAT_SQRT(a) (spx_sqrt(a))
00375 
00376 #endif
00377 
00378 #endif