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

   1 #include <m_pd.h>
   2 #include <m_fixed.h>
   3
   4 /* sigrsqrt - reciprocal square root good to 8 mantissa bits  */
   5
   6 #define DUMTAB1SIZE 256
   7 #define DUMTAB2SIZE 1024
   8
   9 static float rsqrt_exptab[DUMTAB1SIZE], rsqrt_mantissatab[DUMTAB2SIZE];
  10
  11 static void init_rsqrt(void)
  12 {
  13     int i;
  14     for (i = 0; i < DUMTAB1SIZE; i++)
  15     {
  16         float f;
  17         long l = (i ? (i == DUMTAB1SIZE-1 ? DUMTAB1SIZE-2 : i) : 1)<< 23;
  18         *(long *)(&f) = l;
  19         rsqrt_exptab[i] = 1./sqrt(f);
  20     }
  21     for (i = 0; i < DUMTAB2SIZE; i++)
  22     {
  23         float f = 1 + (1./DUMTAB2SIZE) * i;
  24         rsqrt_mantissatab[i] = 1./sqrt(f);
  25     }
  26 }
  27
  28     /* these are used in externs like "bonk" */
  29
  30 float q8_rsqrt(float f)
  31 {
  32     long l = *(long *)(&f);
  33     if (f < 0) return (0);
  34     else return (rsqrt_exptab[(l >> 23) & 0xff] *
  35             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
  36 }
  37
  38 float q8_sqrt(float f)
  39 {
  40     long l = *(long *)(&f);
  41     if (f < 0) return (0);
  42     else return (f * rsqrt_exptab[(l >> 23) & 0xff] *
  43             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
  44 }
  45
  46     /* the old names are OK unless we're in IRIX N32 */
  47
  48 #ifndef N32
  49 float qsqrt(float f) {return (q8_sqrt(f)); }
  50 float qrsqrt(float f) {return (q8_rsqrt(f)); }
  51 #endif
  52
  53
  54
  55 typedef struct sigrsqrt
  56 {
  57     t_object x_obj;
  58     float x_f;
  59 } t_sigrsqrt;
  60
  61 static t_class *sigrsqrt_class;
  62
  63 static void *sigrsqrt_new(void)
  64 {
  65     t_sigrsqrt *x = (t_sigrsqrt *)pd_new(sigrsqrt_class);
  66     outlet_new(&x->x_obj, gensym("signal"));
  67     x->x_f = 0;
  68     return (x);
  69 }
  70
  71 static t_int *sigrsqrt_perform(t_int *w)
  72 {
  73     float *in = *(t_float **)(w+1), *out = *(t_float **)(w+2);
  74     t_int n = *(t_int *)(w+3);
  75     while (n--)
  76     {
  77         float f = *in;
  78         long l = *(long *)(in++);
  79         if (f < 0) *out++ = 0;
  80         else
  81         {
  82             float g = rsqrt_exptab[(l >> 23) & 0xff] *
  83                 rsqrt_mantissatab[(l >> 13) & 0x3ff];
  84             *out++ = 1.5 * g - 0.5 * g * g * g * f;
  85         }
  86     }
  87     return (w + 4);
  88 }
  89
  90 static void sigrsqrt_dsp(t_sigrsqrt *x, t_signal **sp)
  91 {
  92     dsp_add(sigrsqrt_perform, 3, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
  93 }
  94
  95 void rsqrt_tilde_setup(void)
  96 {
  97     init_rsqrt();
  98     sigrsqrt_class = class_new(gensym("rsqrt~"), (t_newmethod)sigrsqrt_new, 0,
  99         sizeof(t_sigrsqrt), 0, 0);
 100             /* an old name for it: */
 101     class_addcreator(sigrsqrt_new, gensym("q8_rsqrt~"), 0);
 102     CLASS_MAINSIGNALIN(sigrsqrt_class, t_sigrsqrt, x_f);
 103     class_addmethod(sigrsqrt_class, (t_method)sigrsqrt_dsp, gensym("dsp"), 0);
 104 }
 105
 106 #include <m_pd.h>
 107 #include <m_fixed.h>
 108
 109 /* sigrsqrt - reciprocal square root good to 8 mantissa bits  */
 110
 111 #define DUMTAB1SIZE 256
 112 #define DUMTAB2SIZE 1024
 113
 114 static float rsqrt_exptab[DUMTAB1SIZE], rsqrt_mantissatab[DUMTAB2SIZE];
 115
 116 static void init_rsqrt(void)
 117 {
 118     int i;
 119     for (i = 0; i < DUMTAB1SIZE; i++)
 120     {
 121         float f;
 122         long l = (i ? (i == DUMTAB1SIZE-1 ? DUMTAB1SIZE-2 : i) : 1)<< 23;
 123         *(long *)(&f) = l;
 124         rsqrt_exptab[i] = 1./sqrt(f);
 125     }
 126     for (i = 0; i < DUMTAB2SIZE; i++)
 127     {
 128         float f = 1 + (1./DUMTAB2SIZE) * i;
 129         rsqrt_mantissatab[i] = 1./sqrt(f);
 130     }
 131 }
 132
 133     /* these are used in externs like "bonk" */
 134
 135 float q8_rsqrt(float f)
 136 {
 137     long l = *(long *)(&f);
 138     if (f < 0) return (0);
 139     else return (rsqrt_exptab[(l >> 23) & 0xff] *
 140             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
 141 }
 142
 143 float q8_sqrt(float f)
 144 {
 145     long l = *(long *)(&f);
 146     if (f < 0) return (0);
 147     else return (f * rsqrt_exptab[(l >> 23) & 0xff] *
 148             rsqrt_mantissatab[(l >> 13) & 0x3ff]);
 149 }
 150
 151     /* the old names are OK unless we're in IRIX N32 */
 152
 153 #ifndef N32
 154 float qsqrt(float f) {return (q8_sqrt(f)); }
 155 float qrsqrt(float f) {return (q8_rsqrt(f)); }
 156 #endif
 157
 158
 159
 160 typedef struct sigrsqrt
 161 {
 162     t_object x_obj;
 163     float x_f;
 164 } t_sigrsqrt;
 165
 166 static t_class *sigrsqrt_class;
 167
 168 static void *sigrsqrt_new(void)
 169 {
 170     t_sigrsqrt *x = (t_sigrsqrt *)pd_new(sigrsqrt_class);
 171     outlet_new(&x->x_obj, gensym("signal"));
 172     x->x_f = 0;
 173     return (x);
 174 }
 175
 176 static t_int *sigrsqrt_perform(t_int *w)
 177 {
 178     float *in = *(t_float **)(w+1), *out = *(t_float **)(w+2);
 179     t_int n = *(t_int *)(w+3);
 180     while (n--)
 181     {
 182         float f = *in;
 183         long l = *(long *)(in++);
 184         if (f < 0) *out++ = 0;
 185         else
 186         {
 187             float g = rsqrt_exptab[(l >> 23) & 0xff] *
 188                 rsqrt_mantissatab[(l >> 13) & 0x3ff];
 189             *out++ = 1.5 * g - 0.5 * g * g * g * f;
 190         }
 191     }
 192     return (w + 4);
 193 }
 194
 195 static void sigrsqrt_dsp(t_sigrsqrt *x, t_signal **sp)
 196 {
 197     dsp_add(sigrsqrt_perform, 3, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
 198 }
 199
 200 void rsqrt_tilde_setup(void)
 201 {
 202     init_rsqrt();
 203     sigrsqrt_class = class_new(gensym("rsqrt~"), (t_newmethod)sigrsqrt_new, 0,
 204         sizeof(t_sigrsqrt), 0, 0);
 205             /* an old name for it: */
 206     class_addcreator(sigrsqrt_new, gensym("q8_rsqrt~"), 0);
 207     CLASS_MAINSIGNALIN(sigrsqrt_class, t_sigrsqrt, x_f);
 208     class_addmethod(sigrsqrt_class, (t_method)sigrsqrt_dsp, gensym("dsp"), 0);
 209 }
 210