Reorganize and commit all the modified TACLeBench code and run scripts

6 files changed, 1063 insertions, 0 deletions

diff --git a/all_pairs/source/fmref/Changelog.txt b/all_pairs/source/fmref/Changelog.txt
new file mode 100644
index 0000000..e9e7055
--- /dev/null
+++ b/all_pairs/source/fmref/Changelog.txt
@@ -0,0 +1,5 @@
+2017-06-27
+- Remove unused variables
+- Introduce fmref_init and fmref_return functions.
+- Add prefix fmref_ to global variables.
+- Introduce dummy initialization in ieee754_rem_pio2f to please linter.
diff --git a/all_pairs/source/fmref/fmref.c b/all_pairs/source/fmref/fmref.c
new file mode 100644
index 0000000..6be1436
--- /dev/null
+++ b/all_pairs/source/fmref/fmref.c
@@ -0,0 +1,283 @@
+/*
+* fmref.c: C reference implementation of FM Radio
+* David Maze <dmaze@cag.lcs.mit.edu>
+* $Id: fmref.c,v 1.2 2010-10-04 21:21:26 garus Exp $
+*/
+
+
+#include "../extra.h"
+#include "wcclibm.h"
+#ifndef M_PI
+#define M_PI        3.1415926535897932384626433832795
+#endif
+
+// Defines
+#define SAMPLING_RATE 250000000
+#define CUTOFF_FREQUENCY 108000000
+#define NUM_TAPS 64
+#define MAX_AMPLITUDE 27000.0
+#define BANDWIDTH 10000
+#define DECIMATION 4 
+/* Must be at least NUM_TAPS+1: */
+#define IN_BUFFER_LEN 200
+#define EQUALIZER_BANDS 10
+
+
+// Type declarations
+typedef struct FloatBuffer
+{
+  float buff[IN_BUFFER_LEN];
+  int rpos, rlen;
+}
+FloatBuffer;
+/* Low pass filter: */
+typedef struct LPFData
+{
+  float coeff[NUM_TAPS];
+  float freq;
+  int taps, decimation;
+}
+LPFData;
+typedef struct EqualizerData
+{
+  LPFData lpf[EQUALIZER_BANDS + 1];
+  FloatBuffer fb[EQUALIZER_BANDS + 1];
+  float gain[EQUALIZER_BANDS];
+}
+EqualizerData;
+
+// Global vars
+float fmref_lpf_coeff[NUM_TAPS];
+float fmref_eq_cutoffs[EQUALIZER_BANDS + 1] =
+  { 55.000004f, 77.78174f, 110.00001f, 155.56354f, 220.00002f, 311.12695f,
+    440.00003f, 622.25415f, 880.00006f, 1244.5078f, 1760.0001f };
+static int fmref_numiters = 2;
+
+// Forward declarations
+void fmref_fb_compact(FloatBuffer *fb);
+int fmref_fb_ensure_writable(FloatBuffer *fb, int amount);
+void fmref_get_floats(FloatBuffer *fb);
+void fmref_init_lpf_data(LPFData *data, float freq, int taps, int decimation);
+void fmref_run_lpf(FloatBuffer *fbin, FloatBuffer *fbout, LPFData *data);
+void fmref_run_demod(FloatBuffer *fbin, FloatBuffer *fbout);
+void fmref_init_equalizer(EqualizerData *data);
+void fmref_run_equalizer(FloatBuffer *fbin, FloatBuffer *fbout, EqualizerData *data);
+void fmref_main(void);
+
+void fmref_init(void)
+{
+  // dummy init function
+}
+
+int fmref_return(void)
+{
+  // dummy return value
+  return 0;
+}
+
+int main(int argc, char **argv){
+
+        SET_UP
+        for(jobsComplete=-1; jobsComplete<maxJobs; jobsComplete++){
+                START_LOOP
+                fmref_init();
+                fmref_main();
+                STOP_LOOP
+        }
+        WRITE_TO_FILE
+  return fmref_return();
+}
+
+FloatBuffer fmref_fb1, fmref_fb2, fmref_fb3, fmref_fb4;
+LPFData fmref_lpf_data;
+
+void fmref_main(void)
+{
+  int i;
+  EqualizerData eq_data;
+
+  fmref_fb1.rpos = fmref_fb1.rlen = 0;
+  fmref_fb2.rpos = fmref_fb2.rlen = 0;
+  fmref_fb3.rpos = fmref_fb3.rlen = 0;
+  fmref_fb4.rpos = fmref_fb4.rlen = 0;
+
+  fmref_init_lpf_data(&fmref_lpf_data, CUTOFF_FREQUENCY, NUM_TAPS, DECIMATION);
+  fmref_init_equalizer(&eq_data);
+
+  /* Startup: */
+  fmref_get_floats(&fmref_fb1);
+  /* LPF needs at least NUM_TAPS+1 inputs; fmref_get_floats is fine. */
+  fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
+  /* run_demod needs 1 input, OK here. */
+  /* run_equalizer needs 51 inputs (same reason as for LPF).  This means
+   * running the pipeline up to demod 50 times in advance: */
+  _Pragma( "loopbound min 64 max 64" )
+  for (i = 0; i < 64; i++) {
+    if (fmref_fb1.rlen - fmref_fb1.rpos < NUM_TAPS + 1)
+      fmref_get_floats(&fmref_fb1);
+    fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
+    fmref_run_demod(&fmref_fb2, &fmref_fb3);
+  }
+
+  /* Main loop: */
+  _Pragma( "loopbound min 2 max 2" )
+  while (fmref_numiters-- > 0) {
+    /* The low-pass filter will need NUM_TAPS+1 items; read them if we
+     * need to. */
+    if (fmref_fb1.rlen - fmref_fb1.rpos < NUM_TAPS + 1)
+      fmref_get_floats(&fmref_fb1);
+    fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
+    fmref_run_demod(&fmref_fb2, &fmref_fb3);
+    fmref_run_equalizer(&fmref_fb3, &fmref_fb4, &eq_data);
+    
+  }
+}
+
+void fmref_fb_compact(FloatBuffer *fb)
+{
+  
+  int i;
+  char *source;
+  char *target;
+  target = (char*)(fb->buff);
+  source = (char*)(fb->buff + fb->rpos);
+  _Pragma( "loopbound min 0 max 60" )
+  for (i = 0; i < fb->rlen - fb->rpos; i++) {
+    target[i] = source[i];
+  }
+  fb->rlen -= fb->rpos;
+  fb->rpos = 0;
+}
+
+int fmref_fb_ensure_writable(FloatBuffer *fb, int amount)
+{
+  int available = IN_BUFFER_LEN - fb->rlen;
+  if (available >= amount)
+    return 1;
+
+  /* Nope, not enough room, move current contents back to the beginning. */
+  fmref_fb_compact(fb);
+
+  available = IN_BUFFER_LEN - fb->rlen;
+  if (available >= amount)
+    return 1;
+
+  return 0;
+}
+
+void fmref_get_floats(FloatBuffer *fb)
+{
+  static int x = 0;
+  fmref_fb_compact(fb);
+
+  /* Fill the remaining space in fb with 1.0. */
+  _Pragma( "loopbound min 200 max 200" )
+  while (fb->rlen < IN_BUFFER_LEN) {
+    fb->buff[fb->rlen++] = (float)x;
+    x++;
+  }
+}
+
+void fmref_init_lpf_data(LPFData *data, float freq, int taps, int decimation)
+{
+  /* Assume that CUTOFF_FREQUENCY is non-zero.  See comments in
+   * StreamIt LowPassFilter.java for origin. */
+  float w = 2 * M_PI * freq / SAMPLING_RATE;
+  int i;
+  float m = taps - 1.0f;
+
+  data->freq = freq;
+  data->taps = taps;
+  data->decimation = decimation;
+
+  _Pragma( "loopbound min 64 max 64" )
+  for (i = 0; i < taps; i++) {
+    if (i - m / 2 == 0.0f)
+      data->coeff[i] = w / M_PI;
+    else
+      data->coeff[i] = sin(w * (i - m / 2)) / M_PI / (i - m / 2) *
+                       (0.54f - 0.46f * cos(2 * M_PI * i / m));
+  }
+}
+
+void fmref_run_lpf(FloatBuffer *fbin, FloatBuffer *fbout, LPFData *data)
+{
+  float sum = 0.0f;
+  int i = 0;
+
+  _Pragma( "loopbound min 64 max 64" )
+  for (i = 0; i < data->taps; i++) {
+    sum += fbin->buff[fbin->rpos + i] * data->coeff[i];
+  }
+
+  fbin->rpos += data->decimation + 1;
+
+  /* Check that there's room in the output buffer; move data if necessary. */
+  fmref_fb_ensure_writable(fbout, 1);
+  fbout->buff[fbout->rlen++] = sum;
+}
+
+void fmref_run_demod(FloatBuffer *fbin, FloatBuffer *fbout)
+{
+  float temp, gain;
+  gain = MAX_AMPLITUDE * SAMPLING_RATE / (BANDWIDTH * M_PI);
+  temp = fbin->buff[fbin->rpos] * fbin->buff[fbin->rpos + 1];
+  temp = gain * atan(temp);
+  fbin->rpos++;
+  fmref_fb_ensure_writable(fbout, 1);
+  fbout->buff[fbout->rlen++] = temp;
+}
+
+void fmref_init_equalizer(EqualizerData *data)
+{
+  int i;
+
+  /* Equalizer structure: there are ten band-pass filters, with
+   * cutoffs as shown below.  The outputs of these filters get added
+   * together.  Each band-pass filter is LPF(high)-LPF(low). */
+  _Pragma( "loopbound min 11 max 11" )
+  for (i = 0; i < EQUALIZER_BANDS + 1; i++)
+    fmref_init_lpf_data(&data->lpf[i], fmref_eq_cutoffs[i], 64, 0);
+
+  /* Also initialize member buffers. */
+  _Pragma( "loopbound min 11 max 11" )
+  for (i = 0; i < EQUALIZER_BANDS + 1; i++)
+    data->fb[i].rpos = data->fb[i].rlen = 0;
+
+  _Pragma( "loopbound min 10 max 10" )
+  for (i = 0; i < EQUALIZER_BANDS; i++) {
+    // the gain amplifies the middle bands the most
+    float val = (((float)i) - (((float)(EQUALIZER_BANDS - 1)) / 2.0f)) / 5.0f;
+    data->gain[i] = val > 0 ? 2.0f - val : 2.0f + val;
+  }
+}
+
+void fmref_run_equalizer(FloatBuffer *fbin, FloatBuffer *fbout, EqualizerData *data)
+{
+  int i, rpos;
+  float lpf_out[EQUALIZER_BANDS + 1];
+  float sum = 0.0;
+
+  /* Save the input read location; we can reuse the same input data on all
+   * of the LPFs. */
+  rpos = fbin->rpos;
+
+  /* Run the child filters. */
+  _Pragma( "loopbound min 11 max 11" )
+  for (i = 0; i < EQUALIZER_BANDS + 1; i++) {
+    fbin->rpos = rpos;
+    fmref_run_lpf(fbin, &data->fb[i], &data->lpf[i]);
+    lpf_out[i] = data->fb[i].buff[data->fb[i].rpos++];
+  }
+
+  /* Now process the results of the filters.  Remember that each band is
+   * output(hi)-output(lo). */
+   _Pragma( "loopbound min 10 max 10" )
+  for (i = 0; i < EQUALIZER_BANDS; i++)
+    sum += (lpf_out[i + 1] - lpf_out[i]) * data->gain[i];
+
+  /* Write that result.  */
+  fmref_fb_ensure_writable(fbout, 1);
+  fbout->buff[fbout->rlen++] = sum;
+}
+
diff --git a/all_pairs/source/fmref/license.txt b/all_pairs/source/fmref/license.txt
new file mode 100644
index 0000000..7029925
--- /dev/null
+++ b/all_pairs/source/fmref/license.txt
@@ -0,0 +1,21 @@
+The MIT License
+
+Copyright (c) <year> <copyright holders>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/all_pairs/source/fmref/math_private.h b/all_pairs/source/fmref/math_private.h
new file mode 100644
index 0000000..58d4354
--- /dev/null
+++ b/all_pairs/source/fmref/math_private.h
@@ -0,0 +1,178 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ */
+
+#ifndef _MATH_PRIVATE_H_
+#define _MATH_PRIVATE_H_
+
+#include "wcclibm.h"
+
+//#include <endian.h>
+//#include <sys/types.h>
+
+/* A representation of a double as a union. */
+union ieee754_double
+{
+  double d;
+
+  /* This is the IEEE 754 double-precision format.  */
+  struct {
+    /* Together these comprise the mantissa.  */
+    unsigned int mantissa1:32;
+    unsigned int mantissa0:20;
+    unsigned int exponent:11;
+    unsigned int negative:1;
+  } ieee;
+
+  /* This format makes it easier to see if a NaN is a signalling NaN.  */
+  struct {
+    /* Together these comprise the mantissa.  */
+    unsigned int mantissa1:32;
+    unsigned int mantissa0:19;
+    unsigned int quiet_nan:1;
+    unsigned int exponent:11;
+    unsigned int negative:1;
+  } ieee_nan;
+};
+
+/* The original fdlibm code used statements like:
+        n0 = ((*(int*)&one)>>29)^1;             * index of high word *
+        ix0 = *(n0+(int*)&x);                   * high word of x *
+        ix1 = *((1-n0)+(int*)&x);               * low word of x *
+   to dig two 32 bit words out of the 64 bit IEEE floating point
+   value.  That is non-ANSI, and, moreover, the gcc instruction
+   scheduler gets it wrong.  We instead use the following macros.
+   Unlike the original code, we determine the endianness at compile
+   time, not at run time; I don't see much benefit to selecting
+   endianness at run time.  */
+
+/* A union which permits us to convert between a double and two 32 bit
+   ints.  */
+
+/* #if __FLOAT_WORD_ORDER == BIG_ENDIAN */
+/* #warning USING Big Endian float word order */
+/* typedef union */
+/* { */
+/*   double value; */
+/*   struct */
+/*   { */
+/*     u_int16_t msw; */
+/*     u_int16_t lsw; */
+/*   } parts; */
+/* } ieeeDoubleShapeType; */
+
+/* #endif */
+
+/* #if __FLOAT_WORD_ORDER == LITTLE_ENDIAN */
+/* #warning USING Little Endian float word order */
+
+typedef union
+{
+  double value;
+  struct
+  {
+    u_int16_t lsw;
+    u_int16_t msw;
+  } parts;
+} ieeeDoubleShapeType;
+
+/* #endif */
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)                                \
+{                                                               \
+  ieeeDoubleShapeType ew_u;                                     \
+  ew_u.value = (d);                                             \
+  (ix0) = ew_u.parts.msw;                                       \
+  (ix1) = ew_u.parts.lsw;                                       \
+}
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)                                      \
+{                                                               \
+  ieeeDoubleShapeType gh_u;                                     \
+  gh_u.value = (d);                                             \
+  (i) = gh_u.parts.msw;                                         \
+}
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i,d)                                       \
+{                                                               \
+  ieeeDoubleShapeType gl_u;                                     \
+  gl_u.value = (d);                                             \
+  (i) = gl_u.parts.lsw;                                         \
+}
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d,ix0,ix1)                                 \
+{                                                               \
+  ieeeDoubleShapeType iw_u;                                     \
+  iw_u.parts.msw = (ix0);                                       \
+  iw_u.parts.lsw = (ix1);                                       \
+  (d) = iw_u.value;                                             \
+}
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)                                      \
+{                                                               \
+  ieeeDoubleShapeType sh_u;                                     \
+  sh_u.value = (d);                                             \
+  sh_u.parts.msw = (v);                                         \
+  (d) = sh_u.value;                                             \
+}
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d,v)                                       \
+{                                                               \
+  ieeeDoubleShapeType sl_u;                                     \
+  sl_u.value = (d);                                             \
+  sl_u.parts.lsw = (v);                                         \
+  (d) = sl_u.value;                                             \
+}
+
+/* A union which permits us to convert between a float and a 32 bit
+   int.  */
+
+typedef union
+{
+  float value;
+  u_int32_t word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)                                     \
+{                                                               \
+  ieee_float_shape_type gf_u;                                   \
+  gf_u.value = (d);                                             \
+  (i) = gf_u.word;                                              \
+}
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)                                     \
+{                                                               \
+  ieee_float_shape_type sf_u;                                   \
+  sf_u.word = (i);                                              \
+  (d) = sf_u.value;                                             \
+}
+
+
+#endif /* _MATH_PRIVATE_H_ */
diff --git a/all_pairs/source/fmref/wcclibm.c b/all_pairs/source/fmref/wcclibm.c
new file mode 100644
index 0000000..39b8cbe
--- /dev/null
+++ b/all_pairs/source/fmref/wcclibm.c
@@ -0,0 +1,522 @@
+#include "math_private.h"
+#include "wcclibm.h"
+
+
+
+/* e_rem_pio2f.c -- float version of e_rem_pio2.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
+#endif
+
+/* __ieee754_rem_pio2f(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2f()
+ */
+
+/* This array is like the one in e_rem_pio2.c, but the numbers are
+   single precision and the last 8 bits are forced to 0.  */
+#ifdef __STDC__
+static const int32_t fmref_npio2_hw[] = {
+#else
+static int32_t fmref_npio2_hw[] = {
+#endif
+0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
+0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
+0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
+0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
+0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
+0x4242c700, 0x42490f00
+};
+
+/*
+ * invpio2:  24 bits of 2/pi
+ * pio2_1:   first  17 bit of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 17 bit of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  17 bit of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+/* zero =  0.0000000000e+00f, /\* 0x00000000 *\/ */
+/* half =  5.0000000000e-01f, /\* 0x3f000000 *\/ */
+/* two8 =  2.5600000000e+02f, /\* 0x43800000 *\/ */
+fmref_invpio2 =  6.3661980629e-01f, /* 0x3f22f984 */
+fmref_pio2_1  =  1.5707855225e+00f, /* 0x3fc90f80 */
+fmref_pio2_1t =  1.0804334124e-05f, /* 0x37354443 */
+fmref_pio2_2  =  1.0804273188e-05f, /* 0x37354400 */
+fmref_pio2_2t =  6.0770999344e-11f, /* 0x2e85a308 */
+fmref_pio2_3  =  6.0770943833e-11f, /* 0x2e85a300 */
+fmref_pio2_3t =  6.1232342629e-17f; /* 0x248d3132 */
+
+#ifdef __STDC__
+        int32_t fmref___ieee754_rem_pio2f(float x, float *y)
+#else
+        int32_t fmref___ieee754_rem_pio2f(x,y)
+        float x,y[];
+#endif
+{
+        float z,w,t,r,fn;
+        int32_t i,j,n,ix,hx;
+
+        GET_FLOAT_WORD(hx,x);
+        ix = hx&0x7fffffff;
+        if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
+            {y[0] = x; y[1] = 0; return 0;}
+        if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
+            if(hx>0) {
+                z = x - fmref_pio2_1;
+                if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+                    y[0] = z - fmref_pio2_1t;
+                    y[1] = (z-y[0])-fmref_pio2_1t;
+                } else {                /* near pi/2, use 24+24+24 bit pi */
+                    z -= fmref_pio2_2;
+                    y[0] = z - fmref_pio2_2t;
+                    y[1] = (z-y[0])-fmref_pio2_2t;
+                }
+                return 1;
+            } else {    /* negative x */
+                z = x + fmref_pio2_1;
+                if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+                    y[0] = z + fmref_pio2_1t;
+                    y[1] = (z-y[0])+fmref_pio2_1t;
+                } else {                /* near pi/2, use 24+24+24 bit pi */
+                    z += fmref_pio2_2;
+                    y[0] = z + fmref_pio2_2t;
+                    y[1] = (z-y[0])+fmref_pio2_2t;
+                }
+                return -1;
+            }
+        }
+        if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
+            t  = fabsf(x);
+            n  = (int32_t) (t*fmref_invpio2+fmref_half);
+            fn = (float)n;
+            r  = t-fn*fmref_pio2_1;
+            w  = fn*fmref_pio2_1t;      /* 1st round good to 40 bit */
+            if(n<32&&(int32_t)(ix&0xffffff00)!=fmref_npio2_hw[n-1]) {
+                y[0] = r-w;     /* quick check no cancellation */
+            } else {
+                u_int32_t high;
+                j  = ix>>23;
+                y[0] = r-w;
+                GET_FLOAT_WORD(high,y[0]);
+                i = j-((high>>23)&0xff);
+                if(i>8) {  /* 2nd iteration needed, good to 57 */
+                    t  = r;
+                    w  = fn*fmref_pio2_2;
+                    r  = t-w;
+                    w  = fn*fmref_pio2_2t-((t-r)-w);
+                    y[0] = r-w;
+                    GET_FLOAT_WORD(high,y[0]);
+                    i = j-((high>>23)&0xff);
+                    if(i>25)  { /* 3rd iteration need, 74 bits acc */
+                        t  = r; /* will cover all possible cases */
+                        w  = fn*fmref_pio2_3;
+                        r  = t-w;
+                        w  = fn*fmref_pio2_3t-((t-r)-w);
+                        y[0] = r-w;
+                    }
+                }
+            }
+            y[1] = (r-y[0])-w;
+            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
+            else         return n;
+        }
+    /*
+     * all other (large) arguments
+     */
+        if(ix>=0x7f800000) {            /* x is inf or NaN */
+            y[0]=y[1]=x-x; return 0;
+        }
+
+        y[0]=y[1]=x-x; /* dummy initialization */
+        return 0; /* doesn't happen for our input */
+}
+
+/* k_cosf.c -- float version of k_cos.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: k_cosf.c,v 1.4 1995/05/10 20:46:23 jtc Exp $";
+#endif
+
+
+#ifdef __STDC__
+static const float 
+#else
+static float 
+#endif
+/* one =  1.0000000000e+00, /\* 0x3f800000 *\/ */
+fmref_C1  =  4.1666667908e-02f, /* 0x3d2aaaab */
+fmref_C2  = -1.3888889225e-03f, /* 0xbab60b61 */
+fmref_C3  =  2.4801587642e-05f, /* 0x37d00d01 */
+fmref_C4  = -2.7557314297e-07f, /* 0xb493f27c */
+fmref_C5  =  2.0875723372e-09f, /* 0x310f74f6 */
+fmref_C6  = -1.1359647598e-11f; /* 0xad47d74e */
+
+#ifdef __STDC__
+        float fmref___kernel_cosf(float x, float y)
+#else
+        float fmref___kernel_cosf(x, y)
+        float x,y;
+#endif
+{
+        float a,hz,z,r,qx;
+        int32_t ix;
+        GET_FLOAT_WORD(ix,x);
+        ix &= 0x7fffffff;                       /* ix = |x|'s high word*/
+        if(ix<0x32000000) {                     /* if x < 2**27 */
+            if(((int)x)==0) return fmref_one;           /* generate inexact */
+        }
+        z  = x*x;
+        r  = z*(fmref_C1+z*(fmref_C2+z*(fmref_C3+z*(fmref_C4+z*(fmref_C5+z*fmref_C6)))));
+        if(ix < 0x3e99999a)                     /* if |x| < 0.3 */ 
+            return fmref_one - ((float)0.5f*z - (z*r - x*y));
+        else {
+            if(ix > 0x3f480000) {               /* x > 0.78125 */
+                qx = (float)0.28125f;
+            } else {
+                SET_FLOAT_WORD(qx,ix-0x01000000);       /* x/4 */
+            }
+            hz = (float)0.5f*z-qx;
+            a  = fmref_one-qx;
+            return a - (hz - (z*r-x*y));
+        }
+}
+
+/* k_sinf.c -- float version of k_sin.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
+#endif
+
+
+#ifdef __STDC__
+static const float 
+#else
+static float 
+#endif
+/* half =  5.0000000000e-01f,/\* 0x3f000000 *\/ */
+fmref_S1  = -1.6666667163e-01f, /* 0xbe2aaaab */
+fmref_S2  =  8.3333337680e-03f, /* 0x3c088889 */
+fmref_S3  = -1.9841270114e-04f, /* 0xb9500d01 */
+fmref_S4  =  2.7557314297e-06f, /* 0x3638ef1b */
+fmref_S5  = -2.5050759689e-08f, /* 0xb2d72f34 */
+fmref_S6  =  1.5896910177e-10f; /* 0x2f2ec9d3 */
+
+#ifdef __STDC__
+        float fmref___kernel_sinf(float x, float y, int iy)
+#else
+        float fmref___kernel_sinf(x, y, iy)
+        float x,y; int iy;              /* iy=0 if y is zero */
+#endif
+{
+        float z,r,v;
+        int32_t ix;
+        GET_FLOAT_WORD(ix,x);
+        ix &= 0x7fffffff;                       /* high word of x */
+        if(ix<0x32000000)                       /* |x| < 2**-27 */
+           {if((int)x==0) return x;}            /* generate inexact */
+        z       =  x*x;
+        v       =  z*x;
+        r       =  fmref_S2+z*(fmref_S3+z*(fmref_S4+z*(fmref_S5+z*fmref_S6)));
+        if(iy==0) return x+v*(fmref_S1+z*r);
+        else      return x-((z*(fmref_half*y-v*r)-y)-v*fmref_S1);
+}
+/* s_atanf.c -- float version of s_atan.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $";
+#endif
+
+
+#ifdef __STDC__
+static const float fmref_atanhi[] = {
+#else
+static float fmref_atanhi[] = {
+#endif
+  4.6364760399e-01f, /* atan(0.5)hi 0x3eed6338 */
+  7.8539812565e-01f, /* atan(1.0)hi 0x3f490fda */
+  9.8279368877e-01f, /* atan(1.5)hi 0x3f7b985e */
+  1.5707962513e+00f, /* atan(inf)hi 0x3fc90fda */
+};
+
+#ifdef __STDC__
+static const float fmref_atanlo[] = {
+#else
+static float fmref_atanlo[] = {
+#endif
+  5.0121582440e-09f, /* atan(0.5)lo 0x31ac3769 */
+  3.7748947079e-08f, /* atan(1.0)lo 0x33222168 */
+  3.4473217170e-08f, /* atan(1.5)lo 0x33140fb4 */
+  7.5497894159e-08f, /* atan(inf)lo 0x33a22168 */
+};
+
+#ifdef __STDC__
+static const float fmref_aT[] = {
+#else
+static float fmref_aT[] = {
+#endif
+  3.3333334327e-01f, /* 0x3eaaaaaa */
+ -2.0000000298e-01f, /* 0xbe4ccccd */
+  1.4285714924e-01f, /* 0x3e124925 */
+ -1.1111110449e-01f, /* 0xbde38e38 */
+  9.0908870101e-02f, /* 0x3dba2e6e */
+ -7.6918758452e-02f, /* 0xbd9d8795 */
+  6.6610731184e-02f, /* 0x3d886b35 */
+ -5.8335702866e-02f, /* 0xbd6ef16b */
+  4.9768779427e-02f, /* 0x3d4bda59 */
+ -3.6531571299e-02f, /* 0xbd15a221 */
+  1.6285819933e-02f, /* 0x3c8569d7 */
+};
+
+/* #ifdef __STDC__ */
+/*      static const float  */
+/* #else */
+/*      static float  */
+/* #endif */
+/* one   = 1.0, */
+/* huge   = 1.0e30; */
+
+#ifdef __STDC__
+        float fmref___atanf(float x)
+#else
+        float fmref___atanf(x)
+        float x;
+#endif
+{
+        float w,s1,s2,z;
+        int32_t ix,hx,id;
+
+        GET_FLOAT_WORD(hx,x);
+        ix = hx&0x7fffffff;
+        if(ix>=0x50800000) {    /* if |x| >= 2^34 */
+            if(ix>0x7f800000)
+                return x+x;             /* NaN */
+            if(hx>0) return  fmref_atanhi[3]+fmref_atanlo[3];
+            else     return -fmref_atanhi[3]-fmref_atanlo[3];
+        } if (ix < 0x3ee00000) {        /* |x| < 0.4375 */
+            if (ix < 0x31000000) {      /* |x| < 2^-29 */
+                if(fmref_huge+x>fmref_one) return x;    /* raise inexact */
+            }
+            id = -1;
+        } else {
+        x = fabsf(x);
+        if (ix < 0x3f980000) {          /* |x| < 1.1875 */
+            if (ix < 0x3f300000) {      /* 7/16 <=|x|<11/16 */
+                id = 0; x = ((float)2.0f*x-fmref_one)/((float)2.0f+x); 
+            } else {                    /* 11/16<=|x|< 19/16 */
+                id = 1; x  = (x-fmref_one)/(x+fmref_one); 
+            }
+        } else {
+            if (ix < 0x401c0000) {      /* |x| < 2.4375 */
+                id = 2; x  = (x-(float)1.5f)/(fmref_one+(float)1.5f*x);
+            } else {                    /* 2.4375 <= |x| < 2^66 */
+                id = 3; x  = -(float)1.0f/x;
+            }
+        }}
+    /* end of argument reduction */
+        z = x*x;
+        w = z*z;
+    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+        s1 = z*(fmref_aT[0]+w*(fmref_aT[2]+w*(fmref_aT[4]+w*(fmref_aT[6]+w*(fmref_aT[8]+w*fmref_aT[10])))));
+        s2 = w*(fmref_aT[1]+w*(fmref_aT[3]+w*(fmref_aT[5]+w*(fmref_aT[7]+w*fmref_aT[9]))));
+        if (id<0) return x - x*(s1+s2);
+        else {
+            z = fmref_atanhi[id] - ((x*(s1+s2) - fmref_atanlo[id]) - x);
+            return (hx<0)? -z:z;
+        }
+}
+//weak_alias (__atanf, atanf)
+
+/* s_cosf.c -- float version of s_cos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* #ifdef __STDC__ */
+/* static const float one=1.0; */
+/* #else */
+/* static float one=1.0; */
+/* #endif */
+
+#ifdef __STDC__
+        float fmref___cosf(float x)
+#else
+        float fmref___cosf(x)
+        float x;
+#endif
+{
+        float y[2],z=0.0f;
+        int32_t n,ix;
+
+        GET_FLOAT_WORD(ix,x);
+
+    /* |x| ~< pi/4 */
+        ix &= 0x7fffffff;
+        if(ix <= 0x3f490fd8) return fmref___kernel_cosf(x,z);
+
+    /* cos(Inf or NaN) is NaN */
+        else if (ix>=0x7f800000) return x-x;
+
+    /* argument reduction needed */
+        else {
+            n = fmref___ieee754_rem_pio2f(x,y);
+            switch(n&3) {
+                case 0: return  fmref___kernel_cosf(y[0],y[1]);
+                case 1: return -fmref___kernel_sinf(y[0],y[1],1);
+                case 2: return -fmref___kernel_cosf(y[0],y[1]);
+                default:
+                        return  fmref___kernel_sinf(y[0],y[1],1);
+            }
+        }
+}
+
+/* s_sinf.c -- float version of s_sin.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+ #ifdef __STDC__
+        float fmref___sinf(float x)
+#else
+        float fmref___sinf(x)
+        float x;
+#endif
+{
+        float y[2],z=0.0;
+        int32_t n, ix;
+
+        GET_FLOAT_WORD(ix,x);
+
+    /* |x| ~< pi/4 */
+        ix &= 0x7fffffff;
+        if(ix <= 0x3f490fd8) return fmref___kernel_sinf(x,z,0);
+
+    /* sin(Inf or NaN) is NaN */
+        else if (ix>=0x7f800000) return x-x;
+
+    /* argument reduction needed */
+        else {
+            n = fmref___ieee754_rem_pio2f(x,y);
+            switch(n&3) {
+                case 0: return  fmref___kernel_sinf(y[0],y[1],1);
+                case 1: return  fmref___kernel_cosf(y[0],y[1]);
+                case 2: return -fmref___kernel_sinf(y[0],y[1],1);
+                default:
+                        return -fmref___kernel_cosf(y[0],y[1]);
+            }
+        }
+}
+
+/* s_fabsf.c -- float version of s_fabs.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * fabsf(x) returns the absolute value of x.
+ */
+
+
+#ifdef __STDC__
+        float fmref___fabsf(float x)
+#else
+        float fmref___fabsf(x)
+        float x;
+#endif
+{
+        u_int32_t ix;
+        GET_FLOAT_WORD(ix,x);
+        SET_FLOAT_WORD(x,ix&0x7fffffff);
+        return x;
+}
diff --git a/all_pairs/source/fmref/wcclibm.h b/all_pairs/source/fmref/wcclibm.h
new file mode 100644
index 0000000..c00738b
--- /dev/null
+++ b/all_pairs/source/fmref/wcclibm.h
@@ -0,0 +1,54 @@
+#ifndef _WCCLIBM
+#define _WCCLIBM
+
+#define size_x       unsigned long
+#define int32_t      int
+#define uint32_t     unsigned int
+#define u_int16_t    unsigned short
+#define u_int32_t    unsigned int
+
+// Often used variables/consts
+#ifdef __STDC__
+static  const float     
+#else
+static  float   
+#endif
+fmref_one          =  1.0f,
+fmref_half   =  5.0000000000e-01f, /* 0x3f000000 */
+fmref_zero   =  0.0f,
+fmref_huge   = 1.0e30,
+fmref_two8   =  2.5600000000e+02f, /* 0x43800000 */
+fmref_twon8  =  3.9062500000e-03f; /* 0x3b800000 */
+
+// The following defines map the math functions to specialized calls
+#define acos     fmref___ieee754_acosf
+#define atan     fmref___atanf
+#define cos      fmref___cosf
+#define fabs     fmref___fabsf
+#define fabsf    fmref___fabsf
+#define isinf    fmref___isinff
+#define pow      fmref___ieee754_powf
+#define sqrt     fmref___ieee754_sqrtf
+#define log10    fmref___ieee754_log10f
+#define log      fmref___ieee754_logf
+#define sin      fmref___sinf
+
+float   fmref___atanf(float x);
+float   fmref___copysignf(float x, float y);
+float   fmref___cosf(float x);
+float   fmref___fabsf(float x);
+float   fmref___floorf(float x);
+float   fmref___ieee754_acosf(float x);
+float   fmref___ieee754_powf(float x, float y);
+int32_t fmref___ieee754_rem_pio2f(float x, float *y);
+float   fmref___ieee754_sqrtf(float x);
+int     fmref___isinff (float x);
+float   fmref___kernel_cosf(float x, float y);
+float   fmref___kernel_sinf(float x, float y, int iy);
+int     fmref___kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2);
+float   fmref___scalbnf (float x, int n);
+float   fmref___ieee754_logf(float x);
+float   fmref___ieee754_log10f(float x);
+float   fmref___sinf(float x);
+
+#endif // _WCCLIBM