apps/codecs/libspeex/pseudofloat.h

   1 /* Copyright (C) 2005 Jean-Marc Valin */
   2 /**
   3    @file pseudofloat.h
   4    @brief Pseudo-floating point
   5  * This header file provides a lightweight floating point type for
   6  * use on fixed-point platforms when a large dynamic range is
   7  * required. The new type is not compatible with the 32-bit IEEE format,
   8  * it is not even remotely as accurate as 32-bit floats, and is not
   9  * even guaranteed to produce even remotely correct results for code
  10  * other than Speex. It makes all kinds of shortcuts that are acceptable
  11  * for Speex, but may not be acceptable for your application. You're
  12  * quite welcome to reuse this code and improve it, but don't assume
  13  * it works out of the box. Most likely, it doesn't.
  14  */
  15 /*
  16    Redistribution and use in source and binary forms, with or without
  17    modification, are permitted provided that the following conditions
  18    are met:
  19
  20    - Redistributions of source code must retain the above copyright
  21    notice, this list of conditions and the following disclaimer.
  22
  23    - Redistributions in binary form must reproduce the above copyright
  24    notice, this list of conditions and the following disclaimer in the
  25    documentation and/or other materials provided with the distribution.
  26
  27    - Neither the name of the Xiph.org Foundation nor the names of its
  28    contributors may be used to endorse or promote products derived from
  29    this software without specific prior written permission.
  30
  31    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  32    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  33    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  34    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
  35    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  36    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  37    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  38    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  39    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  40    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  41    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  42 */
  43
  44 #ifndef PSEUDOFLOAT_H
  45 #define PSEUDOFLOAT_H
  46
  47 #include "arch.h"
  48 #include "os_support.h"
  49 #include "math_approx.h"
  50 #include <math.h>
  51
  52 #ifdef FIXED_POINT
  53
  54 typedef struct {
  55    spx_int16_t m;
  56    spx_int16_t e;
  57 } spx_float_t;
  58
  59 static const spx_float_t FLOAT_ZERO = {0,0};
  60 static const spx_float_t FLOAT_ONE = {16384,-14};
  61 static const spx_float_t FLOAT_HALF = {16384,-15};
  62
  63 #ifdef MIN
  64 #undef MIN
  65 #endif
  66 #define MIN(a,b) ((a)<(b)?(a):(b))
  67 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x)
  68 {
  69    int e=0;
  70    int sign=0;
  71    if (x<0)
  72    {
  73       sign = 1;
  74       x = -x;
  75    }
  76    if (x==0)
  77    {
  78       spx_float_t r = {0,0};
  79       return r;
  80    }
  81    e = spx_ilog2(ABS32(x))-14;
  82    x = VSHR32(x, e);
  83    if (sign)
  84    {
  85       spx_float_t r;
  86       r.m = -x;
  87       r.e = e;
  88       return r;
  89    }
  90    else
  91    {
  92       spx_float_t r;
  93       r.m = x;
  94       r.e = e;
  95       return r;
  96    }
  97 }
  98
  99
 100 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b)
 101 {
 102    spx_float_t r;
 103    if (a.m==0)
 104       return b;
 105    else if (b.m==0)
 106       return a;
 107    if ((a).e > (b).e)
 108    {
 109       r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1));
 110       r.e = (a).e+1;
 111    }
 112    else
 113    {
 114       r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1));
 115       r.e = (b).e+1;
 116    }
 117    if (r.m>0)
 118    {
 119       if (r.m<16384)
 120       {
 121          r.m<<=1;
 122          r.e-=1;
 123       }
 124    } else {
 125       if (r.m>-16384)
 126       {
 127          r.m<<=1;
 128          r.e-=1;
 129       }
 130    }
 131    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 132    return r;
 133 }
 134
 135 static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b)
 136 {
 137    spx_float_t r;
 138    if (a.m==0)
 139       return b;
 140    else if (b.m==0)
 141       return a;
 142    if ((a).e > (b).e)
 143    {
 144       r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1));
 145       r.e = (a).e+1;
 146    }
 147    else
 148    {
 149       r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1);
 150       r.e = (b).e+1;
 151    }
 152    if (r.m>0)
 153    {
 154       if (r.m<16384)
 155       {
 156          r.m<<=1;
 157          r.e-=1;
 158       }
 159    } else {
 160       if (r.m>-16384)
 161       {
 162          r.m<<=1;
 163          r.e-=1;
 164       }
 165    }
 166    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 167    return r;
 168 }
 169
 170 static inline int FLOAT_LT(spx_float_t a, spx_float_t b)
 171 {
 172    if (a.m==0)
 173       return b.m>0;
 174    else if (b.m==0)
 175       return a.m<0;
 176    if ((a).e > (b).e)
 177       return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1));
 178    else
 179       return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1));
 180
 181 }
 182
 183 static inline int FLOAT_GT(spx_float_t a, spx_float_t b)
 184 {
 185    return FLOAT_LT(b,a);
 186 }
 187
 188 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b)
 189 {
 190    spx_float_t r;
 191    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
 192    r.e = (a).e+(b).e+15;
 193    if (r.m>0)
 194    {
 195       if (r.m<16384)
 196       {
 197          r.m<<=1;
 198          r.e-=1;
 199       }
 200    } else {
 201       if (r.m>-16384)
 202       {
 203          r.m<<=1;
 204          r.e-=1;
 205       }
 206    }
 207    /*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 208    return r;
 209 }
 210
 211 static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b)
 212 {
 213    spx_float_t r;
 214    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
 215    r.e = (a).e+(b).e+15;
 216    return r;
 217 }
 218
 219
 220 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b)
 221 {
 222    spx_float_t r;
 223    r.m = a.m;
 224    r.e = a.e+b;
 225    return r;
 226 }
 227
 228 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a)
 229 {
 230    if (a.e<0)
 231       return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e);
 232    else
 233       return a.m<<a.e;
 234 }
 235
 236 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a)
 237 {
 238    if (a.e<0)
 239       return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e;
 240    else
 241       return EXTEND32(a.m)<<a.e;
 242 }
 243
 244 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b)
 245 {
 246    return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15);
 247 }
 248
 249 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b)
 250 {
 251    int e1, e2;
 252    spx_float_t r;
 253    if (a==0 || b==0)
 254    {
 255       return FLOAT_ZERO;
 256    }
 257    e1 = spx_ilog2(ABS32(a));
 258    a = VSHR32(a, e1-14);
 259    e2 = spx_ilog2(ABS32(b));
 260    b = VSHR32(b, e2-14);
 261    r.m = MULT16_16_Q15(a,b);
 262    r.e = e1+e2-13;
 263    return r;
 264 }
 265
 266 /* Do NOT attempt to divide by a negative number */
 267 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b)
 268 {
 269    int e=0;
 270    spx_float_t r;
 271    if (a==0)
 272    {
 273       return FLOAT_ZERO;
 274    }
 275    e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15;
 276    a = VSHR32(a, e);
 277    if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15))
 278    {
 279       a >>= 1;
 280       e++;
 281    }
 282    r.m = DIV32_16(a,b.m);
 283    r.e = e-b.e;
 284    return r;
 285 }
 286
 287
 288 /* Do NOT attempt to divide by a negative number */
 289 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b)
 290 {
 291    int e0=0,e=0;
 292    spx_float_t r;
 293    if (a==0)
 294    {
 295       return FLOAT_ZERO;
 296    }
 297    if (b>32767)
 298    {
 299       e0 = spx_ilog2(b)-14;
 300       b = VSHR32(b, e0);
 301       e0 = -e0;
 302    }
 303    e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15;
 304    a = VSHR32(a, e);
 305    if (ABS32(a)>=SHL32(EXTEND32(b-1),15))
 306    {
 307       a >>= 1;
 308       e++;
 309    }
 310    e += e0;
 311    r.m = DIV32_16(a,b);
 312    r.e = e;
 313    return r;
 314 }
 315
 316 /* Do NOT attempt to divide by a negative number */
 317 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b)
 318 {
 319    int e=0;
 320    spx_int32_t num;
 321    spx_float_t r;
 322    if (b.m<=0)
 323    {
 324       speex_warning_int("Attempted to divide by", b.m);
 325       return FLOAT_ONE;
 326    }
 327    num = a.m;
 328    a.m = ABS16(a.m);
 329    while (a.m >= b.m)
 330    {
 331       e++;
 332       a.m >>= 1;
 333    }
 334    num = num << (15-e);
 335    r.m = DIV32_16(num,b.m);
 336    r.e = a.e-b.e-15+e;
 337    return r;
 338 }
 339
 340 static inline spx_float_t FLOAT_SQRT(spx_float_t a)
 341 {
 342    spx_float_t r;
 343    spx_int32_t m;
 344    m = SHL32(EXTEND32(a.m), 14);
 345    r.e = a.e - 14;
 346    if (r.e & 1)
 347    {
 348       r.e -= 1;
 349       m <<= 1;
 350    }
 351    r.e >>= 1;
 352    r.m = spx_sqrt(m);
 353    return r;
 354 }
 355
 356 #else
 357
 358 #define spx_float_t float
 359 #define FLOAT_ZERO 0.f
 360 #define FLOAT_ONE 1.f
 361 #define FLOAT_HALF 0.5f
 362 #define PSEUDOFLOAT(x) (x)
 363 #define FLOAT_MULT(a,b) ((a)*(b))
 364 #define FLOAT_AMULT(a,b) ((a)*(b))
 365 #define FLOAT_MUL32(a,b) ((a)*(b))
 366 #define FLOAT_DIV32(a,b) ((a)/(b))
 367 #define FLOAT_EXTRACT16(a) (a)
 368 #define FLOAT_EXTRACT32(a) (a)
 369 #define FLOAT_ADD(a,b) ((a)+(b))
 370 #define FLOAT_SUB(a,b) ((a)-(b))
 371 #define REALFLOAT(x) (x)
 372 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b))
 373 #define FLOAT_MUL32U(a,b) ((a)*(b))
 374 #define FLOAT_SHL(a,b) (a)
 375 #define FLOAT_LT(a,b) ((a)<(b))
 376 #define FLOAT_GT(a,b) ((a)>(b))
 377 #define FLOAT_DIVU(a,b) ((a)/(b))
 378 #define FLOAT_SQRT(a) (spx_sqrt(a))
 379
 380 #endif
 381
 382 #endif