src/lib/dl/ext/speex/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 #define MIN(a,b) ((a)<(b)?(a):(b))
  64 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x)
  65 {
  66    int e=0;
  67    int sign=0;
  68    if (x<0)
  69    {
  70       sign = 1;
  71       x = -x;
  72    }
  73    if (x==0)
  74    {
  75       spx_float_t r = {0,0};
  76       return r;
  77    }
  78    e = spx_ilog2(ABS32(x))-14;
  79    x = VSHR32(x, e);
  80    if (sign)
  81    {
  82       spx_float_t r;
  83       r.m = -x;
  84       r.e = e;
  85       return r;
  86    }
  87    else
  88    {
  89       spx_float_t r;
  90       r.m = x;
  91       r.e = e;
  92       return r;
  93    }
  94 }
  95
  96
  97 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b)
  98 {
  99    spx_float_t r;
 100    if (a.m==0)
 101       return b;
 102    else if (b.m==0)
 103       return a;
 104    if ((a).e > (b).e)
 105    {
 106       r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1));
 107       r.e = (a).e+1;
 108    }
 109    else
 110    {
 111       r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1));
 112       r.e = (b).e+1;
 113    }
 114    if (r.m>0)
 115    {
 116       if (r.m<16384)
 117       {
 118          r.m<<=1;
 119          r.e-=1;
 120       }
 121    } else {
 122       if (r.m>-16384)
 123       {
 124          r.m<<=1;
 125          r.e-=1;
 126       }
 127    }
 128    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 129    return r;
 130 }
 131
 132 static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b)
 133 {
 134    spx_float_t r;
 135    if (a.m==0)
 136       return b;
 137    else if (b.m==0)
 138       return a;
 139    if ((a).e > (b).e)
 140    {
 141       r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1));
 142       r.e = (a).e+1;
 143    }
 144    else
 145    {
 146       r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1);
 147       r.e = (b).e+1;
 148    }
 149    if (r.m>0)
 150    {
 151       if (r.m<16384)
 152       {
 153          r.m<<=1;
 154          r.e-=1;
 155       }
 156    } else {
 157       if (r.m>-16384)
 158       {
 159          r.m<<=1;
 160          r.e-=1;
 161       }
 162    }
 163    /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 164    return r;
 165 }
 166
 167 static inline int FLOAT_LT(spx_float_t a, spx_float_t b)
 168 {
 169    if (a.m==0)
 170       return b.m>0;
 171    else if (b.m==0)
 172       return a.m<0;
 173    if ((a).e > (b).e)
 174       return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1));
 175    else
 176       return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1));
 177
 178 }
 179
 180 static inline int FLOAT_GT(spx_float_t a, spx_float_t b)
 181 {
 182    return FLOAT_LT(b,a);
 183 }
 184
 185 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b)
 186 {
 187    spx_float_t r;
 188    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
 189    r.e = (a).e+(b).e+15;
 190    if (r.m>0)
 191    {
 192       if (r.m<16384)
 193       {
 194          r.m<<=1;
 195          r.e-=1;
 196       }
 197    } else {
 198       if (r.m>-16384)
 199       {
 200          r.m<<=1;
 201          r.e-=1;
 202       }
 203    }
 204    /*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
 205    return r;
 206 }
 207
 208 static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b)
 209 {
 210    spx_float_t r;
 211    r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
 212    r.e = (a).e+(b).e+15;
 213    return r;
 214 }
 215
 216
 217 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b)
 218 {
 219    spx_float_t r;
 220    r.m = a.m;
 221    r.e = a.e+b;
 222    return r;
 223 }
 224
 225 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a)
 226 {
 227    if (a.e<0)
 228       return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e);
 229    else
 230       return a.m<<a.e;
 231 }
 232
 233 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a)
 234 {
 235    if (a.e<0)
 236       return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e;
 237    else
 238       return EXTEND32(a.m)<<a.e;
 239 }
 240
 241 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b)
 242 {
 243    return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15);
 244 }
 245
 246 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b)
 247 {
 248    int e1, e2;
 249    spx_float_t r;
 250    if (a==0 || b==0)
 251    {
 252       return FLOAT_ZERO;
 253    }
 254    e1 = spx_ilog2(ABS32(a));
 255    a = VSHR32(a, e1-14);
 256    e2 = spx_ilog2(ABS32(b));
 257    b = VSHR32(b, e2-14);
 258    r.m = MULT16_16_Q15(a,b);
 259    r.e = e1+e2-13;
 260    return r;
 261 }
 262
 263 /* Do NOT attempt to divide by a negative number */
 264 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b)
 265 {
 266    int e=0;
 267    spx_float_t r;
 268    if (a==0)
 269    {
 270       return FLOAT_ZERO;
 271    }
 272    e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15;
 273    a = VSHR32(a, e);
 274    if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15))
 275    {
 276       a >>= 1;
 277       e++;
 278    }
 279    r.m = DIV32_16(a,b.m);
 280    r.e = e-b.e;
 281    return r;
 282 }
 283
 284
 285 /* Do NOT attempt to divide by a negative number */
 286 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b)
 287 {
 288    int e0=0,e=0;
 289    spx_float_t r;
 290    if (a==0)
 291    {
 292       return FLOAT_ZERO;
 293    }
 294    if (b>32767)
 295    {
 296       e0 = spx_ilog2(b)-14;
 297       b = VSHR32(b, e0);
 298       e0 = -e0;
 299    }
 300    e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15;
 301    a = VSHR32(a, e);
 302    if (ABS32(a)>=SHL32(EXTEND32(b-1),15))
 303    {
 304       a >>= 1;
 305       e++;
 306    }
 307    e += e0;
 308    r.m = DIV32_16(a,b);
 309    r.e = e;
 310    return r;
 311 }
 312
 313 /* Do NOT attempt to divide by a negative number */
 314 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b)
 315 {
 316    int e=0;
 317    spx_int32_t num;
 318    spx_float_t r;
 319    if (b.m<=0)
 320    {
 321       speex_warning_int("Attempted to divide by", b.m);
 322       return FLOAT_ONE;
 323    }
 324    num = a.m;
 325    a.m = ABS16(a.m);
 326    while (a.m >= b.m)
 327    {
 328       e++;
 329       a.m >>= 1;
 330    }
 331    num = num << (15-e);
 332    r.m = DIV32_16(num,b.m);
 333    r.e = a.e-b.e-15+e;
 334    return r;
 335 }
 336
 337 static inline spx_float_t FLOAT_SQRT(spx_float_t a)
 338 {
 339    spx_float_t r;
 340    spx_int32_t m;
 341    m = SHL32(EXTEND32(a.m), 14);
 342    r.e = a.e - 14;
 343    if (r.e & 1)
 344    {
 345       r.e -= 1;
 346       m <<= 1;
 347    }
 348    r.e >>= 1;
 349    r.m = spx_sqrt(m);
 350    return r;
 351 }
 352
 353 #else
 354
 355 #define spx_float_t float
 356 #define FLOAT_ZERO 0.f
 357 #define FLOAT_ONE 1.f
 358 #define FLOAT_HALF 0.5f
 359 #define PSEUDOFLOAT(x) (x)
 360 #define FLOAT_MULT(a,b) ((a)*(b))
 361 #define FLOAT_AMULT(a,b) ((a)*(b))
 362 #define FLOAT_MUL32(a,b) ((a)*(b))
 363 #define FLOAT_DIV32(a,b) ((a)/(b))
 364 #define FLOAT_EXTRACT16(a) (a)
 365 #define FLOAT_EXTRACT32(a) (a)
 366 #define FLOAT_ADD(a,b) ((a)+(b))
 367 #define FLOAT_SUB(a,b) ((a)-(b))
 368 #define REALFLOAT(x) (x)
 369 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b))
 370 #define FLOAT_MUL32U(a,b) ((a)*(b))
 371 #define FLOAT_SHL(a,b) (a)
 372 #define FLOAT_LT(a,b) ((a)<(b))
 373 #define FLOAT_GT(a,b) ((a)>(b))
 374 #define FLOAT_DIVU(a,b) ((a)/(b))
 375 #define FLOAT_SQRT(a) (spx_sqrt(a))
 376
 377 #endif
 378
 379 #endif