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

   1 #define DUMTAB1SIZE 256
   2 #define DUMTAB2SIZE 1024
   3
   4 #include<m_pd.h>
   5 #include<m_fixed.h>
   6
   7 /* sigsqrt -  square root good to 8 mantissa bits  */
   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 typedef struct sigsqrt
  29 {
  30     t_object x_obj;
  31     float x_f;
  32 } t_sigsqrt;
  33
  34 static t_class *sigsqrt_class;
  35
  36 static void *sigsqrt_new(void)
  37 {
  38     t_sigsqrt *x = (t_sigsqrt *)pd_new(sigsqrt_class);
  39     outlet_new(&x->x_obj, gensym("signal"));
  40     x->x_f = 0;
  41     return (x);
  42 }
  43
  44 t_int *sigsqrt_perform(t_int *w)    /* not static; also used in d_fft.c */
  45 {
  46     float *in = *(t_float **)(w+1), *out = *(t_float **)(w+2);
  47     t_int n = *(t_int *)(w+3);
  48     while (n--)
  49     {
  50         float f = *in;
  51         long l = *(long *)(in++);
  52         if (f < 0) *out++ = 0;
  53         else
  54         {
  55             float g = rsqrt_exptab[(l >> 23) & 0xff] *
  56                 rsqrt_mantissatab[(l >> 13) & 0x3ff];
  57             *out++ = f * (1.5 * g - 0.5 * g * g * g * f);
  58         }
  59     }
  60     return (w + 4);
  61 }
  62
  63 static void sigsqrt_dsp(t_sigsqrt *x, t_signal **sp)
  64 {
  65     dsp_add(sigsqrt_perform, 3, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
  66 }
  67
  68 void sqrt_tilde_setup(void)
  69 {
  70     init_rsqrt();
  71     sigsqrt_class = class_new(gensym("sqrt~"), (t_newmethod)sigsqrt_new, 0,
  72         sizeof(t_sigsqrt), 0, 0);
  73     class_addcreator(sigsqrt_new, gensym("q8_sqrt~"), 0);   /* old name */
  74     CLASS_MAINSIGNALIN(sigsqrt_class, t_sigsqrt, x_f);
  75     class_addmethod(sigsqrt_class, (t_method)sigsqrt_dsp, gensym("dsp"), 0);
  76 }
  77
  78 #define DUMTAB1SIZE 256
  79 #define DUMTAB2SIZE 1024
  80
  81 #include<m_pd.h>
  82 #include<m_fixed.h>
  83
  84 /* sigsqrt -  square root good to 8 mantissa bits  */
  85
  86 static float rsqrt_exptab[DUMTAB1SIZE], rsqrt_mantissatab[DUMTAB2SIZE];
  87
  88 static void init_rsqrt(void)
  89 {
  90     int i;
  91     for (i = 0; i < DUMTAB1SIZE; i++)
  92     {
  93         float f;
  94         long l = (i ? (i == DUMTAB1SIZE-1 ? DUMTAB1SIZE-2 : i) : 1)<< 23;
  95         *(long *)(&f) = l;
  96         rsqrt_exptab[i] = 1./sqrt(f);
  97     }
  98     for (i = 0; i < DUMTAB2SIZE; i++)
  99     {
 100         float f = 1 + (1./DUMTAB2SIZE) * i;
 101         rsqrt_mantissatab[i] = 1./sqrt(f);
 102     }
 103 }
 104
 105 typedef struct sigsqrt
 106 {
 107     t_object x_obj;
 108     float x_f;
 109 } t_sigsqrt;
 110
 111 static t_class *sigsqrt_class;
 112
 113 static void *sigsqrt_new(void)
 114 {
 115     t_sigsqrt *x = (t_sigsqrt *)pd_new(sigsqrt_class);
 116     outlet_new(&x->x_obj, gensym("signal"));
 117     x->x_f = 0;
 118     return (x);
 119 }
 120
 121 t_int *sigsqrt_perform(t_int *w)    /* not static; also used in d_fft.c */
 122 {
 123     float *in = *(t_float **)(w+1), *out = *(t_float **)(w+2);
 124     t_int n = *(t_int *)(w+3);
 125     while (n--)
 126     {
 127         float f = *in;
 128         long l = *(long *)(in++);
 129         if (f < 0) *out++ = 0;
 130         else
 131         {
 132             float g = rsqrt_exptab[(l >> 23) & 0xff] *
 133                 rsqrt_mantissatab[(l >> 13) & 0x3ff];
 134             *out++ = f * (1.5 * g - 0.5 * g * g * g * f);
 135         }
 136     }
 137     return (w + 4);
 138 }
 139
 140 static void sigsqrt_dsp(t_sigsqrt *x, t_signal **sp)
 141 {
 142     dsp_add(sigsqrt_perform, 3, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
 143 }
 144
 145 void sqrt_tilde_setup(void)
 146 {
 147     init_rsqrt();
 148     sigsqrt_class = class_new(gensym("sqrt~"), (t_newmethod)sigsqrt_new, 0,
 149         sizeof(t_sigsqrt), 0, 0);
 150     class_addcreator(sigsqrt_new, gensym("q8_sqrt~"), 0);   /* old name */
 151     CLASS_MAINSIGNALIN(sigsqrt_class, t_sigsqrt, x_f);
 152     class_addmethod(sigsqrt_class, (t_method)sigsqrt_dsp, gensym("dsp"), 0);
 153 }
 154