summaryrefslogtreecommitdiffstats
path: root/baseline/source/fmref
diff options
context:
space:
mode:
authorJoshua Bakita <bakitajoshua@gmail.com>2019-10-07 19:13:39 -0400
committerJoshua Bakita <bakitajoshua@gmail.com>2019-10-07 19:13:39 -0400
commit386b7d3366f1359a265da207a9cafa3edf553b64 (patch)
treec76120c2c138faed822e4ae386be6ef22a738a78 /baseline/source/fmref
parent54a3f7091a2146b29c73a6fdc4b62a5c4ad7a3d8 (diff)
Reorganize and commit all the modified TACLeBench code and run scripts
Diffstat (limited to 'baseline/source/fmref')
-rw-r--r--baseline/source/fmref/Changelog.txt5
-rw-r--r--baseline/source/fmref/fmref.c283
-rw-r--r--baseline/source/fmref/license.txt21
-rw-r--r--baseline/source/fmref/math_private.h178
-rw-r--r--baseline/source/fmref/wcclibm.c522
-rw-r--r--baseline/source/fmref/wcclibm.h54
6 files changed, 1063 insertions, 0 deletions
diff --git a/baseline/source/fmref/Changelog.txt b/baseline/source/fmref/Changelog.txt
new file mode 100644
index 0000000..e9e7055
--- /dev/null
+++ b/baseline/source/fmref/Changelog.txt
@@ -0,0 +1,5 @@
12017-06-27
2- Remove unused variables
3- Introduce fmref_init and fmref_return functions.
4- Add prefix fmref_ to global variables.
5- Introduce dummy initialization in ieee754_rem_pio2f to please linter.
diff --git a/baseline/source/fmref/fmref.c b/baseline/source/fmref/fmref.c
new file mode 100644
index 0000000..6be1436
--- /dev/null
+++ b/baseline/source/fmref/fmref.c
@@ -0,0 +1,283 @@
1/*
2* fmref.c: C reference implementation of FM Radio
3* David Maze <dmaze@cag.lcs.mit.edu>
4* $Id: fmref.c,v 1.2 2010-10-04 21:21:26 garus Exp $
5*/
6
7
8#include "../extra.h"
9#include "wcclibm.h"
10#ifndef M_PI
11#define M_PI 3.1415926535897932384626433832795
12#endif
13
14// Defines
15#define SAMPLING_RATE 250000000
16#define CUTOFF_FREQUENCY 108000000
17#define NUM_TAPS 64
18#define MAX_AMPLITUDE 27000.0
19#define BANDWIDTH 10000
20#define DECIMATION 4
21/* Must be at least NUM_TAPS+1: */
22#define IN_BUFFER_LEN 200
23#define EQUALIZER_BANDS 10
24
25
26// Type declarations
27typedef struct FloatBuffer
28{
29 float buff[IN_BUFFER_LEN];
30 int rpos, rlen;
31}
32FloatBuffer;
33/* Low pass filter: */
34typedef struct LPFData
35{
36 float coeff[NUM_TAPS];
37 float freq;
38 int taps, decimation;
39}
40LPFData;
41typedef struct EqualizerData
42{
43 LPFData lpf[EQUALIZER_BANDS + 1];
44 FloatBuffer fb[EQUALIZER_BANDS + 1];
45 float gain[EQUALIZER_BANDS];
46}
47EqualizerData;
48
49// Global vars
50float fmref_lpf_coeff[NUM_TAPS];
51float fmref_eq_cutoffs[EQUALIZER_BANDS + 1] =
52 { 55.000004f, 77.78174f, 110.00001f, 155.56354f, 220.00002f, 311.12695f,
53 440.00003f, 622.25415f, 880.00006f, 1244.5078f, 1760.0001f };
54static int fmref_numiters = 2;
55
56// Forward declarations
57void fmref_fb_compact(FloatBuffer *fb);
58int fmref_fb_ensure_writable(FloatBuffer *fb, int amount);
59void fmref_get_floats(FloatBuffer *fb);
60void fmref_init_lpf_data(LPFData *data, float freq, int taps, int decimation);
61void fmref_run_lpf(FloatBuffer *fbin, FloatBuffer *fbout, LPFData *data);
62void fmref_run_demod(FloatBuffer *fbin, FloatBuffer *fbout);
63void fmref_init_equalizer(EqualizerData *data);
64void fmref_run_equalizer(FloatBuffer *fbin, FloatBuffer *fbout, EqualizerData *data);
65void fmref_main(void);
66
67void fmref_init(void)
68{
69 // dummy init function
70}
71
72int fmref_return(void)
73{
74 // dummy return value
75 return 0;
76}
77
78int main(int argc, char **argv){
79
80 SET_UP
81 for(jobsComplete=-1; jobsComplete<maxJobs; jobsComplete++){
82 START_LOOP
83 fmref_init();
84 fmref_main();
85 STOP_LOOP
86 }
87 WRITE_TO_FILE
88 return fmref_return();
89}
90
91FloatBuffer fmref_fb1, fmref_fb2, fmref_fb3, fmref_fb4;
92LPFData fmref_lpf_data;
93
94void fmref_main(void)
95{
96 int i;
97 EqualizerData eq_data;
98
99 fmref_fb1.rpos = fmref_fb1.rlen = 0;
100 fmref_fb2.rpos = fmref_fb2.rlen = 0;
101 fmref_fb3.rpos = fmref_fb3.rlen = 0;
102 fmref_fb4.rpos = fmref_fb4.rlen = 0;
103
104 fmref_init_lpf_data(&fmref_lpf_data, CUTOFF_FREQUENCY, NUM_TAPS, DECIMATION);
105 fmref_init_equalizer(&eq_data);
106
107 /* Startup: */
108 fmref_get_floats(&fmref_fb1);
109 /* LPF needs at least NUM_TAPS+1 inputs; fmref_get_floats is fine. */
110 fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
111 /* run_demod needs 1 input, OK here. */
112 /* run_equalizer needs 51 inputs (same reason as for LPF). This means
113 * running the pipeline up to demod 50 times in advance: */
114 _Pragma( "loopbound min 64 max 64" )
115 for (i = 0; i < 64; i++) {
116 if (fmref_fb1.rlen - fmref_fb1.rpos < NUM_TAPS + 1)
117 fmref_get_floats(&fmref_fb1);
118 fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
119 fmref_run_demod(&fmref_fb2, &fmref_fb3);
120 }
121
122 /* Main loop: */
123 _Pragma( "loopbound min 2 max 2" )
124 while (fmref_numiters-- > 0) {
125 /* The low-pass filter will need NUM_TAPS+1 items; read them if we
126 * need to. */
127 if (fmref_fb1.rlen - fmref_fb1.rpos < NUM_TAPS + 1)
128 fmref_get_floats(&fmref_fb1);
129 fmref_run_lpf(&fmref_fb1, &fmref_fb2, &fmref_lpf_data);
130 fmref_run_demod(&fmref_fb2, &fmref_fb3);
131 fmref_run_equalizer(&fmref_fb3, &fmref_fb4, &eq_data);
132
133 }
134}
135
136void fmref_fb_compact(FloatBuffer *fb)
137{
138
139 int i;
140 char *source;
141 char *target;
142 target = (char*)(fb->buff);
143 source = (char*)(fb->buff + fb->rpos);
144 _Pragma( "loopbound min 0 max 60" )
145 for (i = 0; i < fb->rlen - fb->rpos; i++) {
146 target[i] = source[i];
147 }
148 fb->rlen -= fb->rpos;
149 fb->rpos = 0;
150}
151
152int fmref_fb_ensure_writable(FloatBuffer *fb, int amount)
153{
154 int available = IN_BUFFER_LEN - fb->rlen;
155 if (available >= amount)
156 return 1;
157
158 /* Nope, not enough room, move current contents back to the beginning. */
159 fmref_fb_compact(fb);
160
161 available = IN_BUFFER_LEN - fb->rlen;
162 if (available >= amount)
163 return 1;
164
165 return 0;
166}
167
168void fmref_get_floats(FloatBuffer *fb)
169{
170 static int x = 0;
171 fmref_fb_compact(fb);
172
173 /* Fill the remaining space in fb with 1.0. */
174 _Pragma( "loopbound min 200 max 200" )
175 while (fb->rlen < IN_BUFFER_LEN) {
176 fb->buff[fb->rlen++] = (float)x;
177 x++;
178 }
179}
180
181void fmref_init_lpf_data(LPFData *data, float freq, int taps, int decimation)
182{
183 /* Assume that CUTOFF_FREQUENCY is non-zero. See comments in
184 * StreamIt LowPassFilter.java for origin. */
185 float w = 2 * M_PI * freq / SAMPLING_RATE;
186 int i;
187 float m = taps - 1.0f;
188
189 data->freq = freq;
190 data->taps = taps;
191 data->decimation = decimation;
192
193 _Pragma( "loopbound min 64 max 64" )
194 for (i = 0; i < taps; i++) {
195 if (i - m / 2 == 0.0f)
196 data->coeff[i] = w / M_PI;
197 else
198 data->coeff[i] = sin(w * (i - m / 2)) / M_PI / (i - m / 2) *
199 (0.54f - 0.46f * cos(2 * M_PI * i / m));
200 }
201}
202
203void fmref_run_lpf(FloatBuffer *fbin, FloatBuffer *fbout, LPFData *data)
204{
205 float sum = 0.0f;
206 int i = 0;
207
208 _Pragma( "loopbound min 64 max 64" )
209 for (i = 0; i < data->taps; i++) {
210 sum += fbin->buff[fbin->rpos + i] * data->coeff[i];
211 }
212
213 fbin->rpos += data->decimation + 1;
214
215 /* Check that there's room in the output buffer; move data if necessary. */
216 fmref_fb_ensure_writable(fbout, 1);
217 fbout->buff[fbout->rlen++] = sum;
218}
219
220void fmref_run_demod(FloatBuffer *fbin, FloatBuffer *fbout)
221{
222 float temp, gain;
223 gain = MAX_AMPLITUDE * SAMPLING_RATE / (BANDWIDTH * M_PI);
224 temp = fbin->buff[fbin->rpos] * fbin->buff[fbin->rpos + 1];
225 temp = gain * atan(temp);
226 fbin->rpos++;
227 fmref_fb_ensure_writable(fbout, 1);
228 fbout->buff[fbout->rlen++] = temp;
229}
230
231void fmref_init_equalizer(EqualizerData *data)
232{
233 int i;
234
235 /* Equalizer structure: there are ten band-pass filters, with
236 * cutoffs as shown below. The outputs of these filters get added
237 * together. Each band-pass filter is LPF(high)-LPF(low). */
238 _Pragma( "loopbound min 11 max 11" )
239 for (i = 0; i < EQUALIZER_BANDS + 1; i++)
240 fmref_init_lpf_data(&data->lpf[i], fmref_eq_cutoffs[i], 64, 0);
241
242 /* Also initialize member buffers. */
243 _Pragma( "loopbound min 11 max 11" )
244 for (i = 0; i < EQUALIZER_BANDS + 1; i++)
245 data->fb[i].rpos = data->fb[i].rlen = 0;
246
247 _Pragma( "loopbound min 10 max 10" )
248 for (i = 0; i < EQUALIZER_BANDS; i++) {
249 // the gain amplifies the middle bands the most
250 float val = (((float)i) - (((float)(EQUALIZER_BANDS - 1)) / 2.0f)) / 5.0f;
251 data->gain[i] = val > 0 ? 2.0f - val : 2.0f + val;
252 }
253}
254
255void fmref_run_equalizer(FloatBuffer *fbin, FloatBuffer *fbout, EqualizerData *data)
256{
257 int i, rpos;
258 float lpf_out[EQUALIZER_BANDS + 1];
259 float sum = 0.0;
260
261 /* Save the input read location; we can reuse the same input data on all
262 * of the LPFs. */
263 rpos = fbin->rpos;
264
265 /* Run the child filters. */
266 _Pragma( "loopbound min 11 max 11" )
267 for (i = 0; i < EQUALIZER_BANDS + 1; i++) {
268 fbin->rpos = rpos;
269 fmref_run_lpf(fbin, &data->fb[i], &data->lpf[i]);
270 lpf_out[i] = data->fb[i].buff[data->fb[i].rpos++];
271 }
272
273 /* Now process the results of the filters. Remember that each band is
274 * output(hi)-output(lo). */
275 _Pragma( "loopbound min 10 max 10" )
276 for (i = 0; i < EQUALIZER_BANDS; i++)
277 sum += (lpf_out[i + 1] - lpf_out[i]) * data->gain[i];
278
279 /* Write that result. */
280 fmref_fb_ensure_writable(fbout, 1);
281 fbout->buff[fbout->rlen++] = sum;
282}
283
diff --git a/baseline/source/fmref/license.txt b/baseline/source/fmref/license.txt
new file mode 100644
index 0000000..7029925
--- /dev/null
+++ b/baseline/source/fmref/license.txt
@@ -0,0 +1,21 @@
1The MIT License
2
3Copyright (c) <year> <copyright holders>
4
5Permission is hereby granted, free of charge, to any person obtaining a copy
6of this software and associated documentation files (the "Software"), to deal
7in the Software without restriction, including without limitation the rights
8to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9copies of the Software, and to permit persons to whom the Software is
10furnished to do so, subject to the following conditions:
11
12The above copyright notice and this permission notice shall be included in
13all copies or substantial portions of the Software.
14
15THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21THE SOFTWARE.
diff --git a/baseline/source/fmref/math_private.h b/baseline/source/fmref/math_private.h
new file mode 100644
index 0000000..58d4354
--- /dev/null
+++ b/baseline/source/fmref/math_private.h
@@ -0,0 +1,178 @@
1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/*
13 * from: @(#)fdlibm.h 5.1 93/09/24
14 */
15
16#ifndef _MATH_PRIVATE_H_
17#define _MATH_PRIVATE_H_
18
19#include "wcclibm.h"
20
21//#include <endian.h>
22//#include <sys/types.h>
23
24/* A representation of a double as a union. */
25union ieee754_double
26{
27 double d;
28
29 /* This is the IEEE 754 double-precision format. */
30 struct {
31 /* Together these comprise the mantissa. */
32 unsigned int mantissa1:32;
33 unsigned int mantissa0:20;
34 unsigned int exponent:11;
35 unsigned int negative:1;
36 } ieee;
37
38 /* This format makes it easier to see if a NaN is a signalling NaN. */
39 struct {
40 /* Together these comprise the mantissa. */
41 unsigned int mantissa1:32;
42 unsigned int mantissa0:19;
43 unsigned int quiet_nan:1;
44 unsigned int exponent:11;
45 unsigned int negative:1;
46 } ieee_nan;
47};
48
49/* The original fdlibm code used statements like:
50 n0 = ((*(int*)&one)>>29)^1; * index of high word *
51 ix0 = *(n0+(int*)&x); * high word of x *
52 ix1 = *((1-n0)+(int*)&x); * low word of x *
53 to dig two 32 bit words out of the 64 bit IEEE floating point
54 value. That is non-ANSI, and, moreover, the gcc instruction
55 scheduler gets it wrong. We instead use the following macros.
56 Unlike the original code, we determine the endianness at compile
57 time, not at run time; I don't see much benefit to selecting
58 endianness at run time. */
59
60/* A union which permits us to convert between a double and two 32 bit
61 ints. */
62
63/* #if __FLOAT_WORD_ORDER == BIG_ENDIAN */
64/* #warning USING Big Endian float word order */
65/* typedef union */
66/* { */
67/* double value; */
68/* struct */
69/* { */
70/* u_int16_t msw; */
71/* u_int16_t lsw; */
72/* } parts; */
73/* } ieeeDoubleShapeType; */
74
75/* #endif */
76
77/* #if __FLOAT_WORD_ORDER == LITTLE_ENDIAN */
78/* #warning USING Little Endian float word order */
79
80typedef union
81{
82 double value;
83 struct
84 {
85 u_int16_t lsw;
86 u_int16_t msw;
87 } parts;
88} ieeeDoubleShapeType;
89
90/* #endif */
91
92/* Get two 32 bit ints from a double. */
93
94#define EXTRACT_WORDS(ix0,ix1,d) \
95{ \
96 ieeeDoubleShapeType ew_u; \
97 ew_u.value = (d); \
98 (ix0) = ew_u.parts.msw; \
99 (ix1) = ew_u.parts.lsw; \
100}
101
102/* Get the more significant 32 bit int from a double. */
103
104#define GET_HIGH_WORD(i,d) \
105{ \
106 ieeeDoubleShapeType gh_u; \
107 gh_u.value = (d); \
108 (i) = gh_u.parts.msw; \
109}
110
111/* Get the less significant 32 bit int from a double. */
112
113#define GET_LOW_WORD(i,d) \
114{ \
115 ieeeDoubleShapeType gl_u; \
116 gl_u.value = (d); \
117 (i) = gl_u.parts.lsw; \
118}
119
120/* Set a double from two 32 bit ints. */
121
122#define INSERT_WORDS(d,ix0,ix1) \
123{ \
124 ieeeDoubleShapeType iw_u; \
125 iw_u.parts.msw = (ix0); \
126 iw_u.parts.lsw = (ix1); \
127 (d) = iw_u.value; \
128}
129
130/* Set the more significant 32 bits of a double from an int. */
131
132#define SET_HIGH_WORD(d,v) \
133{ \
134 ieeeDoubleShapeType sh_u; \
135 sh_u.value = (d); \
136 sh_u.parts.msw = (v); \
137 (d) = sh_u.value; \
138}
139
140/* Set the less significant 32 bits of a double from an int. */
141
142#define SET_LOW_WORD(d,v) \
143{ \
144 ieeeDoubleShapeType sl_u; \
145 sl_u.value = (d); \
146 sl_u.parts.lsw = (v); \
147 (d) = sl_u.value; \
148}
149
150/* A union which permits us to convert between a float and a 32 bit
151 int. */
152
153typedef union
154{
155 float value;
156 u_int32_t word;
157} ieee_float_shape_type;
158
159/* Get a 32 bit int from a float. */
160
161#define GET_FLOAT_WORD(i,d) \
162{ \
163 ieee_float_shape_type gf_u; \
164 gf_u.value = (d); \
165 (i) = gf_u.word; \
166}
167
168/* Set a float from a 32 bit int. */
169
170#define SET_FLOAT_WORD(d,i) \
171{ \
172 ieee_float_shape_type sf_u; \
173 sf_u.word = (i); \
174 (d) = sf_u.value; \
175}
176
177
178#endif /* _MATH_PRIVATE_H_ */
diff --git a/baseline/source/fmref/wcclibm.c b/baseline/source/fmref/wcclibm.c
new file mode 100644
index 0000000..39b8cbe
--- /dev/null
+++ b/baseline/source/fmref/wcclibm.c
@@ -0,0 +1,522 @@
1#include "math_private.h"
2#include "wcclibm.h"
3
4
5
6/* e_rem_pio2f.c -- float version of e_rem_pio2.c
7 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
8 */
9
10/*
11 * ====================================================
12 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
13 *
14 * Developed at SunPro, a Sun Microsystems, Inc. business.
15 * Permission to use, copy, modify, and distribute this
16 * software is freely granted, provided that this notice
17 * is preserved.
18 * ====================================================
19 */
20
21#if defined(LIBM_SCCS) && !defined(lint)
22static char rcsid[] = "$NetBSD: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
23#endif
24
25/* __ieee754_rem_pio2f(x,y)
26 *
27 * return the remainder of x rem pi/2 in y[0]+y[1]
28 * use __kernel_rem_pio2f()
29 */
30
31/* This array is like the one in e_rem_pio2.c, but the numbers are
32 single precision and the last 8 bits are forced to 0. */
33#ifdef __STDC__
34static const int32_t fmref_npio2_hw[] = {
35#else
36static int32_t fmref_npio2_hw[] = {
37#endif
380x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
390x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
400x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
410x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
420x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
430x4242c700, 0x42490f00
44};
45
46/*
47 * invpio2: 24 bits of 2/pi
48 * pio2_1: first 17 bit of pi/2
49 * pio2_1t: pi/2 - pio2_1
50 * pio2_2: second 17 bit of pi/2
51 * pio2_2t: pi/2 - (pio2_1+pio2_2)
52 * pio2_3: third 17 bit of pi/2
53 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
54 */
55
56#ifdef __STDC__
57static const float
58#else
59static float
60#endif
61/* zero = 0.0000000000e+00f, /\* 0x00000000 *\/ */
62/* half = 5.0000000000e-01f, /\* 0x3f000000 *\/ */
63/* two8 = 2.5600000000e+02f, /\* 0x43800000 *\/ */
64fmref_invpio2 = 6.3661980629e-01f, /* 0x3f22f984 */
65fmref_pio2_1 = 1.5707855225e+00f, /* 0x3fc90f80 */
66fmref_pio2_1t = 1.0804334124e-05f, /* 0x37354443 */
67fmref_pio2_2 = 1.0804273188e-05f, /* 0x37354400 */
68fmref_pio2_2t = 6.0770999344e-11f, /* 0x2e85a308 */
69fmref_pio2_3 = 6.0770943833e-11f, /* 0x2e85a300 */
70fmref_pio2_3t = 6.1232342629e-17f; /* 0x248d3132 */
71
72#ifdef __STDC__
73 int32_t fmref___ieee754_rem_pio2f(float x, float *y)
74#else
75 int32_t fmref___ieee754_rem_pio2f(x,y)
76 float x,y[];
77#endif
78{
79 float z,w,t,r,fn;
80 int32_t i,j,n,ix,hx;
81
82 GET_FLOAT_WORD(hx,x);
83 ix = hx&0x7fffffff;
84 if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
85 {y[0] = x; y[1] = 0; return 0;}
86 if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
87 if(hx>0) {
88 z = x - fmref_pio2_1;
89 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
90 y[0] = z - fmref_pio2_1t;
91 y[1] = (z-y[0])-fmref_pio2_1t;
92 } else { /* near pi/2, use 24+24+24 bit pi */
93 z -= fmref_pio2_2;
94 y[0] = z - fmref_pio2_2t;
95 y[1] = (z-y[0])-fmref_pio2_2t;
96 }
97 return 1;
98 } else { /* negative x */
99 z = x + fmref_pio2_1;
100 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
101 y[0] = z + fmref_pio2_1t;
102 y[1] = (z-y[0])+fmref_pio2_1t;
103 } else { /* near pi/2, use 24+24+24 bit pi */
104 z += fmref_pio2_2;
105 y[0] = z + fmref_pio2_2t;
106 y[1] = (z-y[0])+fmref_pio2_2t;
107 }
108 return -1;
109 }
110 }
111 if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
112 t = fabsf(x);
113 n = (int32_t) (t*fmref_invpio2+fmref_half);
114 fn = (float)n;
115 r = t-fn*fmref_pio2_1;
116 w = fn*fmref_pio2_1t; /* 1st round good to 40 bit */
117 if(n<32&&(int32_t)(ix&0xffffff00)!=fmref_npio2_hw[n-1]) {
118 y[0] = r-w; /* quick check no cancellation */
119 } else {
120 u_int32_t high;
121 j = ix>>23;
122 y[0] = r-w;
123 GET_FLOAT_WORD(high,y[0]);
124 i = j-((high>>23)&0xff);
125 if(i>8) { /* 2nd iteration needed, good to 57 */
126 t = r;
127 w = fn*fmref_pio2_2;
128 r = t-w;
129 w = fn*fmref_pio2_2t-((t-r)-w);
130 y[0] = r-w;
131 GET_FLOAT_WORD(high,y[0]);
132 i = j-((high>>23)&0xff);
133 if(i>25) { /* 3rd iteration need, 74 bits acc */
134 t = r; /* will cover all possible cases */
135 w = fn*fmref_pio2_3;
136 r = t-w;
137 w = fn*fmref_pio2_3t-((t-r)-w);
138 y[0] = r-w;
139 }
140 }
141 }
142 y[1] = (r-y[0])-w;
143 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
144 else return n;
145 }
146 /*
147 * all other (large) arguments
148 */
149 if(ix>=0x7f800000) { /* x is inf or NaN */
150 y[0]=y[1]=x-x; return 0;
151 }
152
153 y[0]=y[1]=x-x; /* dummy initialization */
154 return 0; /* doesn't happen for our input */
155}
156
157/* k_cosf.c -- float version of k_cos.c
158 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
159 */
160
161/*
162 * ====================================================
163 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
164 *
165 * Developed at SunPro, a Sun Microsystems, Inc. business.
166 * Permission to use, copy, modify, and distribute this
167 * software is freely granted, provided that this notice
168 * is preserved.
169 * ====================================================
170 */
171
172#if defined(LIBM_SCCS) && !defined(lint)
173static char rcsid[] = "$NetBSD: k_cosf.c,v 1.4 1995/05/10 20:46:23 jtc Exp $";
174#endif
175
176
177#ifdef __STDC__
178static const float
179#else
180static float
181#endif
182/* one = 1.0000000000e+00, /\* 0x3f800000 *\/ */
183fmref_C1 = 4.1666667908e-02f, /* 0x3d2aaaab */
184fmref_C2 = -1.3888889225e-03f, /* 0xbab60b61 */
185fmref_C3 = 2.4801587642e-05f, /* 0x37d00d01 */
186fmref_C4 = -2.7557314297e-07f, /* 0xb493f27c */
187fmref_C5 = 2.0875723372e-09f, /* 0x310f74f6 */
188fmref_C6 = -1.1359647598e-11f; /* 0xad47d74e */
189
190#ifdef __STDC__
191 float fmref___kernel_cosf(float x, float y)
192#else
193 float fmref___kernel_cosf(x, y)
194 float x,y;
195#endif
196{
197 float a,hz,z,r,qx;
198 int32_t ix;
199 GET_FLOAT_WORD(ix,x);
200 ix &= 0x7fffffff; /* ix = |x|'s high word*/
201 if(ix<0x32000000) { /* if x < 2**27 */
202 if(((int)x)==0) return fmref_one; /* generate inexact */
203 }
204 z = x*x;
205 r = z*(fmref_C1+z*(fmref_C2+z*(fmref_C3+z*(fmref_C4+z*(fmref_C5+z*fmref_C6)))));
206 if(ix < 0x3e99999a) /* if |x| < 0.3 */
207 return fmref_one - ((float)0.5f*z - (z*r - x*y));
208 else {
209 if(ix > 0x3f480000) { /* x > 0.78125 */
210 qx = (float)0.28125f;
211 } else {
212 SET_FLOAT_WORD(qx,ix-0x01000000); /* x/4 */
213 }
214 hz = (float)0.5f*z-qx;
215 a = fmref_one-qx;
216 return a - (hz - (z*r-x*y));
217 }
218}
219
220/* k_sinf.c -- float version of k_sin.c
221 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
222 */
223
224/*
225 * ====================================================
226 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
227 *
228 * Developed at SunPro, a Sun Microsystems, Inc. business.
229 * Permission to use, copy, modify, and distribute this
230 * software is freely granted, provided that this notice
231 * is preserved.
232 * ====================================================
233 */
234
235#if defined(LIBM_SCCS) && !defined(lint)
236static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
237#endif
238
239
240#ifdef __STDC__
241static const float
242#else
243static float
244#endif
245/* half = 5.0000000000e-01f,/\* 0x3f000000 *\/ */
246fmref_S1 = -1.6666667163e-01f, /* 0xbe2aaaab */
247fmref_S2 = 8.3333337680e-03f, /* 0x3c088889 */
248fmref_S3 = -1.9841270114e-04f, /* 0xb9500d01 */
249fmref_S4 = 2.7557314297e-06f, /* 0x3638ef1b */
250fmref_S5 = -2.5050759689e-08f, /* 0xb2d72f34 */
251fmref_S6 = 1.5896910177e-10f; /* 0x2f2ec9d3 */
252
253#ifdef __STDC__
254 float fmref___kernel_sinf(float x, float y, int iy)
255#else
256 float fmref___kernel_sinf(x, y, iy)
257 float x,y; int iy; /* iy=0 if y is zero */
258#endif
259{
260 float z,r,v;
261 int32_t ix;
262 GET_FLOAT_WORD(ix,x);
263 ix &= 0x7fffffff; /* high word of x */
264 if(ix<0x32000000) /* |x| < 2**-27 */
265 {if((int)x==0) return x;} /* generate inexact */
266 z = x*x;
267 v = z*x;
268 r = fmref_S2+z*(fmref_S3+z*(fmref_S4+z*(fmref_S5+z*fmref_S6)));
269 if(iy==0) return x+v*(fmref_S1+z*r);
270 else return x-((z*(fmref_half*y-v*r)-y)-v*fmref_S1);
271}
272/* s_atanf.c -- float version of s_atan.c.
273 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
274 */
275
276/*
277 * ====================================================
278 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
279 *
280 * Developed at SunPro, a Sun Microsystems, Inc. business.
281 * Permission to use, copy, modify, and distribute this
282 * software is freely granted, provided that this notice
283 * is preserved.
284 * ====================================================
285 */
286
287#if defined(LIBM_SCCS) && !defined(lint)
288static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $";
289#endif
290
291
292#ifdef __STDC__
293static const float fmref_atanhi[] = {
294#else
295static float fmref_atanhi[] = {
296#endif
297 4.6364760399e-01f, /* atan(0.5)hi 0x3eed6338 */
298 7.8539812565e-01f, /* atan(1.0)hi 0x3f490fda */
299 9.8279368877e-01f, /* atan(1.5)hi 0x3f7b985e */
300 1.5707962513e+00f, /* atan(inf)hi 0x3fc90fda */
301};
302
303#ifdef __STDC__
304static const float fmref_atanlo[] = {
305#else
306static float fmref_atanlo[] = {
307#endif
308 5.0121582440e-09f, /* atan(0.5)lo 0x31ac3769 */
309 3.7748947079e-08f, /* atan(1.0)lo 0x33222168 */
310 3.4473217170e-08f, /* atan(1.5)lo 0x33140fb4 */
311 7.5497894159e-08f, /* atan(inf)lo 0x33a22168 */
312};
313
314#ifdef __STDC__
315static const float fmref_aT[] = {
316#else
317static float fmref_aT[] = {
318#endif
319 3.3333334327e-01f, /* 0x3eaaaaaa */
320 -2.0000000298e-01f, /* 0xbe4ccccd */
321 1.4285714924e-01f, /* 0x3e124925 */
322 -1.1111110449e-01f, /* 0xbde38e38 */
323 9.0908870101e-02f, /* 0x3dba2e6e */
324 -7.6918758452e-02f, /* 0xbd9d8795 */
325 6.6610731184e-02f, /* 0x3d886b35 */
326 -5.8335702866e-02f, /* 0xbd6ef16b */
327 4.9768779427e-02f, /* 0x3d4bda59 */
328 -3.6531571299e-02f, /* 0xbd15a221 */
329 1.6285819933e-02f, /* 0x3c8569d7 */
330};
331
332/* #ifdef __STDC__ */
333/* static const float */
334/* #else */
335/* static float */
336/* #endif */
337/* one = 1.0, */
338/* huge = 1.0e30; */
339
340#ifdef __STDC__
341 float fmref___atanf(float x)
342#else
343 float fmref___atanf(x)
344 float x;
345#endif
346{
347 float w,s1,s2,z;
348 int32_t ix,hx,id;
349
350 GET_FLOAT_WORD(hx,x);
351 ix = hx&0x7fffffff;
352 if(ix>=0x50800000) { /* if |x| >= 2^34 */
353 if(ix>0x7f800000)
354 return x+x; /* NaN */
355 if(hx>0) return fmref_atanhi[3]+fmref_atanlo[3];
356 else return -fmref_atanhi[3]-fmref_atanlo[3];
357 } if (ix < 0x3ee00000) { /* |x| < 0.4375 */
358 if (ix < 0x31000000) { /* |x| < 2^-29 */
359 if(fmref_huge+x>fmref_one) return x; /* raise inexact */
360 }
361 id = -1;
362 } else {
363 x = fabsf(x);
364 if (ix < 0x3f980000) { /* |x| < 1.1875 */
365 if (ix < 0x3f300000) { /* 7/16 <=|x|<11/16 */
366 id = 0; x = ((float)2.0f*x-fmref_one)/((float)2.0f+x);
367 } else { /* 11/16<=|x|< 19/16 */
368 id = 1; x = (x-fmref_one)/(x+fmref_one);
369 }
370 } else {
371 if (ix < 0x401c0000) { /* |x| < 2.4375 */
372 id = 2; x = (x-(float)1.5f)/(fmref_one+(float)1.5f*x);
373 } else { /* 2.4375 <= |x| < 2^66 */
374 id = 3; x = -(float)1.0f/x;
375 }
376 }}
377 /* end of argument reduction */
378 z = x*x;
379 w = z*z;
380 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
381 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])))));
382 s2 = w*(fmref_aT[1]+w*(fmref_aT[3]+w*(fmref_aT[5]+w*(fmref_aT[7]+w*fmref_aT[9]))));
383 if (id<0) return x - x*(s1+s2);
384 else {
385 z = fmref_atanhi[id] - ((x*(s1+s2) - fmref_atanlo[id]) - x);
386 return (hx<0)? -z:z;
387 }
388}
389//weak_alias (__atanf, atanf)
390
391/* s_cosf.c -- float version of s_cos.c.
392 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
393 */
394
395/*
396 * ====================================================
397 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
398 *
399 * Developed at SunPro, a Sun Microsystems, Inc. business.
400 * Permission to use, copy, modify, and distribute this
401 * software is freely granted, provided that this notice
402 * is preserved.
403 * ====================================================
404 */
405
406/* #ifdef __STDC__ */
407/* static const float one=1.0; */
408/* #else */
409/* static float one=1.0; */
410/* #endif */
411
412#ifdef __STDC__
413 float fmref___cosf(float x)
414#else
415 float fmref___cosf(x)
416 float x;
417#endif
418{
419 float y[2],z=0.0f;
420 int32_t n,ix;
421
422 GET_FLOAT_WORD(ix,x);
423
424 /* |x| ~< pi/4 */
425 ix &= 0x7fffffff;
426 if(ix <= 0x3f490fd8) return fmref___kernel_cosf(x,z);
427
428 /* cos(Inf or NaN) is NaN */
429 else if (ix>=0x7f800000) return x-x;
430
431 /* argument reduction needed */
432 else {
433 n = fmref___ieee754_rem_pio2f(x,y);
434 switch(n&3) {
435 case 0: return fmref___kernel_cosf(y[0],y[1]);
436 case 1: return -fmref___kernel_sinf(y[0],y[1],1);
437 case 2: return -fmref___kernel_cosf(y[0],y[1]);
438 default:
439 return fmref___kernel_sinf(y[0],y[1],1);
440 }
441 }
442}
443
444/* s_sinf.c -- float version of s_sin.c.
445 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
446 */
447
448/*
449 * ====================================================
450 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
451 *
452 * Developed at SunPro, a Sun Microsystems, Inc. business.
453 * Permission to use, copy, modify, and distribute this
454 * software is freely granted, provided that this notice
455 * is preserved.
456 * ====================================================
457 */
458
459 #ifdef __STDC__
460 float fmref___sinf(float x)
461#else
462 float fmref___sinf(x)
463 float x;
464#endif
465{
466 float y[2],z=0.0;
467 int32_t n, ix;
468
469 GET_FLOAT_WORD(ix,x);
470
471 /* |x| ~< pi/4 */
472 ix &= 0x7fffffff;
473 if(ix <= 0x3f490fd8) return fmref___kernel_sinf(x,z,0);
474
475 /* sin(Inf or NaN) is NaN */
476 else if (ix>=0x7f800000) return x-x;
477
478 /* argument reduction needed */
479 else {
480 n = fmref___ieee754_rem_pio2f(x,y);
481 switch(n&3) {
482 case 0: return fmref___kernel_sinf(y[0],y[1],1);
483 case 1: return fmref___kernel_cosf(y[0],y[1]);
484 case 2: return -fmref___kernel_sinf(y[0],y[1],1);
485 default:
486 return -fmref___kernel_cosf(y[0],y[1]);
487 }
488 }
489}
490
491/* s_fabsf.c -- float version of s_fabs.c.
492 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
493 */
494
495/*
496 * ====================================================
497 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
498 *
499 * Developed at SunPro, a Sun Microsystems, Inc. business.
500 * Permission to use, copy, modify, and distribute this
501 * software is freely granted, provided that this notice
502 * is preserved.
503 * ====================================================
504 */
505
506/*
507 * fabsf(x) returns the absolute value of x.
508 */
509
510
511#ifdef __STDC__
512 float fmref___fabsf(float x)
513#else
514 float fmref___fabsf(x)
515 float x;
516#endif
517{
518 u_int32_t ix;
519 GET_FLOAT_WORD(ix,x);
520 SET_FLOAT_WORD(x,ix&0x7fffffff);
521 return x;
522}
diff --git a/baseline/source/fmref/wcclibm.h b/baseline/source/fmref/wcclibm.h
new file mode 100644
index 0000000..c00738b
--- /dev/null
+++ b/baseline/source/fmref/wcclibm.h
@@ -0,0 +1,54 @@
1#ifndef _WCCLIBM
2#define _WCCLIBM
3
4#define size_x unsigned long
5#define int32_t int
6#define uint32_t unsigned int
7#define u_int16_t unsigned short
8#define u_int32_t unsigned int
9
10// Often used variables/consts
11#ifdef __STDC__
12static const float
13#else
14static float
15#endif
16fmref_one = 1.0f,
17fmref_half = 5.0000000000e-01f, /* 0x3f000000 */
18fmref_zero = 0.0f,
19fmref_huge = 1.0e30,
20fmref_two8 = 2.5600000000e+02f, /* 0x43800000 */
21fmref_twon8 = 3.9062500000e-03f; /* 0x3b800000 */
22
23// The following defines map the math functions to specialized calls
24#define acos fmref___ieee754_acosf
25#define atan fmref___atanf
26#define cos fmref___cosf
27#define fabs fmref___fabsf
28#define fabsf fmref___fabsf
29#define isinf fmref___isinff
30#define pow fmref___ieee754_powf
31#define sqrt fmref___ieee754_sqrtf
32#define log10 fmref___ieee754_log10f
33#define log fmref___ieee754_logf
34#define sin fmref___sinf
35
36float fmref___atanf(float x);
37float fmref___copysignf(float x, float y);
38float fmref___cosf(float x);
39float fmref___fabsf(float x);
40float fmref___floorf(float x);
41float fmref___ieee754_acosf(float x);
42float fmref___ieee754_powf(float x, float y);
43int32_t fmref___ieee754_rem_pio2f(float x, float *y);
44float fmref___ieee754_sqrtf(float x);
45int fmref___isinff (float x);
46float fmref___kernel_cosf(float x, float y);
47float fmref___kernel_sinf(float x, float y, int iy);
48int fmref___kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2);
49float fmref___scalbnf (float x, int n);
50float fmref___ieee754_logf(float x);
51float fmref___ieee754_log10f(float x);
52float fmref___sinf(float x);
53
54#endif // _WCCLIBM