src/libsphinxbase/fe/fixlog.c

   1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
   2 /* ====================================================================
   3  * Copyright (c) 2005 Carnegie Mellon University.  All rights
   4  * reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  *
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  *
  13  * 2. Redistributions in binary form must reproduce the above copyright
  14  *    notice, this list of conditions and the following disclaimer in
  15  *    the documentation and/or other materials provided with the
  16  *    distribution.
  17  *
  18  * This work was supported in part by funding from the Defense Advanced
  19  * Research Projects Agency and the National Science Foundation of the
  20  * United States of America, and the CMU Sphinx Speech Consortium.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
  23  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
  24  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  25  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
  26  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  27  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  28  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  29  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  30  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  31  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  32  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  33  *
  34  * ====================================================================
  35  *
  36  * File: fixlog.c
  37  *
  38  * Description: Fast approximate fixed-point logarithms
  39  *
  40  * Author: David Huggins-Daines <dhuggins@cs.cmu.edu>
  41  *
  42  */
  43
  44 #ifdef HAVE_CONFIG_H
  45 #include <config.h>
  46 #endif
  47
  48 #include "sphinxbase/prim_type.h"
  49 #include "sphinxbase/fixpoint.h"
  50
  51 #include "fe_internal.h"
  52
  53 /* Table of log2(x/64)*(1<<DEFAULT_RADIX) */
  54 /* perl -e 'for (0..63) {my $x = 1 + $_/64; print "\t(uint32)(", log($x)/log(2), "*(1<<DEFAULT_RADIX)),\n"}' */
  55 static uint32 logtable[] = {
  56     (uint32) (0 * (1 << DEFAULT_RADIX)),
  57     (uint32) (0.0223678130284545 * (1 << DEFAULT_RADIX)),
  58     (uint32) (0.0443941193584534 * (1 << DEFAULT_RADIX)),
  59     (uint32) (0.0660891904577724 * (1 << DEFAULT_RADIX)),
  60     (uint32) (0.0874628412503394 * (1 << DEFAULT_RADIX)),
  61     (uint32) (0.108524456778169 * (1 << DEFAULT_RADIX)),
  62     (uint32) (0.129283016944966 * (1 << DEFAULT_RADIX)),
  63     (uint32) (0.149747119504682 * (1 << DEFAULT_RADIX)),
  64     (uint32) (0.169925001442312 * (1 << DEFAULT_RADIX)),
  65     (uint32) (0.189824558880017 * (1 << DEFAULT_RADIX)),
  66     (uint32) (0.20945336562895 * (1 << DEFAULT_RADIX)),
  67     (uint32) (0.228818690495881 * (1 << DEFAULT_RADIX)),
  68     (uint32) (0.247927513443586 * (1 << DEFAULT_RADIX)),
  69     (uint32) (0.266786540694901 * (1 << DEFAULT_RADIX)),
  70     (uint32) (0.285402218862248 * (1 << DEFAULT_RADIX)),
  71     (uint32) (0.303780748177103 * (1 << DEFAULT_RADIX)),
  72     (uint32) (0.321928094887362 * (1 << DEFAULT_RADIX)),
  73     (uint32) (0.339850002884625 * (1 << DEFAULT_RADIX)),
  74     (uint32) (0.357552004618084 * (1 << DEFAULT_RADIX)),
  75     (uint32) (0.375039431346925 * (1 << DEFAULT_RADIX)),
  76     (uint32) (0.39231742277876 * (1 << DEFAULT_RADIX)),
  77     (uint32) (0.409390936137702 * (1 << DEFAULT_RADIX)),
  78     (uint32) (0.426264754702098 * (1 << DEFAULT_RADIX)),
  79     (uint32) (0.442943495848728 * (1 << DEFAULT_RADIX)),
  80     (uint32) (0.459431618637297 * (1 << DEFAULT_RADIX)),
  81     (uint32) (0.475733430966398 * (1 << DEFAULT_RADIX)),
  82     (uint32) (0.491853096329675 * (1 << DEFAULT_RADIX)),
  83     (uint32) (0.507794640198696 * (1 << DEFAULT_RADIX)),
  84     (uint32) (0.523561956057013 * (1 << DEFAULT_RADIX)),
  85     (uint32) (0.539158811108031 * (1 << DEFAULT_RADIX)),
  86     (uint32) (0.554588851677637 * (1 << DEFAULT_RADIX)),
  87     (uint32) (0.569855608330948 * (1 << DEFAULT_RADIX)),
  88     (uint32) (0.584962500721156 * (1 << DEFAULT_RADIX)),
  89     (uint32) (0.599912842187128 * (1 << DEFAULT_RADIX)),
  90     (uint32) (0.614709844115208 * (1 << DEFAULT_RADIX)),
  91     (uint32) (0.62935662007961 * (1 << DEFAULT_RADIX)),
  92     (uint32) (0.643856189774725 * (1 << DEFAULT_RADIX)),
  93     (uint32) (0.658211482751795 * (1 << DEFAULT_RADIX)),
  94     (uint32) (0.672425341971496 * (1 << DEFAULT_RADIX)),
  95     (uint32) (0.686500527183218 * (1 << DEFAULT_RADIX)),
  96     (uint32) (0.700439718141092 * (1 << DEFAULT_RADIX)),
  97     (uint32) (0.714245517666123 * (1 << DEFAULT_RADIX)),
  98     (uint32) (0.727920454563199 * (1 << DEFAULT_RADIX)),
  99     (uint32) (0.741466986401147 * (1 << DEFAULT_RADIX)),
 100     (uint32) (0.754887502163469 * (1 << DEFAULT_RADIX)),
 101     (uint32) (0.768184324776926 * (1 << DEFAULT_RADIX)),
 102     (uint32) (0.78135971352466 * (1 << DEFAULT_RADIX)),
 103     (uint32) (0.794415866350106 * (1 << DEFAULT_RADIX)),
 104     (uint32) (0.807354922057604 * (1 << DEFAULT_RADIX)),
 105     (uint32) (0.820178962415188 * (1 << DEFAULT_RADIX)),
 106     (uint32) (0.832890014164742 * (1 << DEFAULT_RADIX)),
 107     (uint32) (0.845490050944375 * (1 << DEFAULT_RADIX)),
 108     (uint32) (0.857980995127572 * (1 << DEFAULT_RADIX)),
 109     (uint32) (0.870364719583405 * (1 << DEFAULT_RADIX)),
 110     (uint32) (0.882643049361841 * (1 << DEFAULT_RADIX)),
 111     (uint32) (0.894817763307944 * (1 << DEFAULT_RADIX)),
 112     (uint32) (0.906890595608518 * (1 << DEFAULT_RADIX)),
 113     (uint32) (0.918863237274595 * (1 << DEFAULT_RADIX)),
 114     (uint32) (0.930737337562886 * (1 << DEFAULT_RADIX)),
 115     (uint32) (0.94251450533924 * (1 << DEFAULT_RADIX)),
 116     (uint32) (0.954196310386875 * (1 << DEFAULT_RADIX)),
 117     (uint32) (0.965784284662087 * (1 << DEFAULT_RADIX)),
 118     (uint32) (0.977279923499917 * (1 << DEFAULT_RADIX)),
 119     (uint32) (0.988684686772166 * (1 << DEFAULT_RADIX)),
 120 };
 121
 122 int32
 123 fixlog2(uint32 x)
 124 {
 125     uint32 y;
 126
 127     if (x == 0)
 128         return MIN_FIXLOG2;
 129
 130     /* Get the exponent. */
 131 #if defined(__GNUC__) && defined(__i386__)
 132   __asm__("bsrl %1, %0\n": "=r"(y): "g"(x):"cc");
 133     x <<= (31 - y);
 134 #elif defined(__ppc__)
 135   __asm__("cntlzw %0, %1\n": "=r"(y):"r"(x));
 136     x <<= y;
 137     y = 31 - y;
 138 #elif ((defined(__ARM_ARCH_5__) || defined(__ARM_ARCH_5T__) || \
 139         defined(__ARM_ARCH_5TE__)) && !defined(__thumb__))
 140   __asm__("clz %0, %1\n": "=r"(y):"r"(x));
 141     x <<= y;
 142     y = 31 - y;
 143 #else
 144     for (y = 31; y >= 0; --y) {
 145         if (x & 0x80000000)
 146             break;
 147         x <<= 1;
 148     }
 149 #endif
 150     y <<= DEFAULT_RADIX;
 151     /* Do a table lookup for the MSB of the mantissa. */
 152     x = (x >> 25) & 0x3f;
 153     return y + logtable[x];
 154 }
 155
 156 int
 157 fixlog(uint32 x)
 158 {
 159     int32 y;
 160
 161     if (x == 0)
 162         return MIN_FIXLOG;
 163
 164     y = fixlog2(x);
 165     return FIXMUL(y, FIXLN_2);
 166 }