apps/plugins/pdbox/PDa/intern/rsqrt~.c

   1 #ifdef ROCKBOX
   2 #include "plugin.h"
   3 #include "../../pdbox.h"
   4 #endif
   5
   6 #include "../src/m_pd.h"
   7 #include <../src/m_fixed.h>
   8
   9 /* sigrsqrt - reciprocal square root good to 8 mantissa bits  */
  10
  11 #define DUMTAB1SIZE 256
  12 #define DUMTAB2SIZE 1024
  13
  14 static float rsqrt_exptab[DUMTAB1SIZE], rsqrt_mantissatab[DUMTAB2SIZE];
  15
  16 static void init_rsqrt(void)
  17 {
  18     int i;
  19     for (i = 0; i < DUMTAB1SIZE; i++)
  20     {
  21 #ifdef ROCKBOX
  22         union f2i f2i;
  23         f2i.i = (i ? (i == DUMTAB1SIZE-1 ? DUMTAB1SIZE-2 : i) : 1)<< 23;
  24         rsqrt_exptab[i] = 1./sqrt(f2i.f);
  25 #else /* ROCKBOX */
  26         float f;
  27         long l = (i ? (i == DUMTAB1SIZE-1 ? DUMTAB1SIZE-2 : i) : 1)<< 23;
  28         *(long *)(&f) = l;
  29         rsqrt_exptab[i] = 1./sqrt(f);
  30 #endif /* ROCKBOX */
  31     }
  32     for (i = 0; i < DUMTAB2SIZE; i++)
  33     {
  34         float f = 1 + (1./DUMTAB2SIZE) * i;
  35         rsqrt_mantissatab[i] = 1./sqrt(f);
  36     }
  37 }
  38
  39     /* these are used in externs like "bonk" */
  40
  41 float q8_rsqrt(float f)
  42 {
  43 #ifdef ROCKBOX
  44     union f2i f2i;
  45     if(f < 0.0)
  46         return 0.0;
  47     else
  48     {
  49         f2i.f = f;
  50         return (rsqrt_exptab[(f2i.i >> 23) & 0xff] *
  51             rsqrt_mantissatab[(f2i.i >> 13) & 0x3ff]);
  52     }
  53 #else /* ROCKBOX */
  54     long l = *(long *)(&f);
  55     if (f < 0) return (0);
  56     else return (rsqrt_exptab[(l >> 23) & 0xff] *
  57             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
  58 #endif /* ROCKBOX */
  59 }
  60
  61 float q8_sqrt(float f)
  62 {
  63 #ifdef ROCKBOX
  64     union f2i f2i;
  65     if(f < 0.0)
  66         return 0.0;
  67     else
  68     {
  69         f2i.f = f;
  70         return (f * rsqrt_exptab[(f2i.i >> 23) & 0xff] *
  71             rsqrt_mantissatab[(f2i.i >> 13) & 0x3ff]);
  72     }
  73 #else /* ROCKBOX */
  74     long l = *(long *)(&f);
  75     if (f < 0) return (0);
  76     else return (f * rsqrt_exptab[(l >> 23) & 0xff] *
  77             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
  78 #endif /* ROCKBOX */
  79 }
  80
  81     /* the old names are OK unless we're in IRIX N32 */
  82
  83 #ifndef N32
  84 float qsqrt(float f) {return (q8_sqrt(f)); }
  85 float qrsqrt(float f) {return (q8_rsqrt(f)); }
  86 #endif
  87
  88
  89
  90 typedef struct sigrsqrt
  91 {
  92     t_object x_obj;
  93     float x_f;
  94 } t_sigrsqrt;
  95
  96 static t_class *sigrsqrt_class;
  97
  98 static void *sigrsqrt_new(void)
  99 {
 100     t_sigrsqrt *x = (t_sigrsqrt *)pd_new(sigrsqrt_class);
 101     outlet_new(&x->x_obj, gensym("signal"));
 102     x->x_f = 0;
 103     return (x);
 104 }
 105
 106 static t_int *sigrsqrt_perform(t_int *w)
 107 {
 108     float *in = *(t_float **)(w+1), *out = *(t_float **)(w+2);
 109     t_int n = *(t_int *)(w+3);
 110     while (n--)
 111     {
 112         float f = *in;
 113         long l = *(long *)(in++);
 114         if (f < 0) *out++ = 0;
 115         else
 116         {
 117             float g = rsqrt_exptab[(l >> 23) & 0xff] *
 118                 rsqrt_mantissatab[(l >> 13) & 0x3ff];
 119             *out++ = 1.5 * g - 0.5 * g * g * g * f;
 120         }
 121     }
 122     return (w + 4);
 123 }
 124
 125 static void sigrsqrt_dsp(t_sigrsqrt *x, t_signal **sp)
 126 {
 127 #ifdef ROCKBOX
 128     (void) x;
 129 #endif
 130     dsp_add(sigrsqrt_perform, 3, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
 131 }
 132
 133 void rsqrt_tilde_setup(void)
 134 {
 135     init_rsqrt();
 136     sigrsqrt_class = class_new(gensym("rsqrt~"), (t_newmethod)sigrsqrt_new, 0,
 137         sizeof(t_sigrsqrt), 0, 0);
 138             /* an old name for it: */
 139     class_addcreator(sigrsqrt_new, gensym("q8_rsqrt~"), 0);
 140     CLASS_MAINSIGNALIN(sigrsqrt_class, t_sigrsqrt, x_f);
 141     class_addmethod(sigrsqrt_class, (t_method)sigrsqrt_dsp, gensym("dsp"), 0);
 142 }
 143