1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <stdint.h>
5
6 #define BUFSZ (64 * 64 * sizeof(uint16_t))
7 #define NUM_COEFF_BUCKETS (4)
8 #define NUM_OTHER_BUCKETS (0)
9 #define NUM_TOTAL_BUCKETS ((NUM_COEFF_BUCKETS) + (NUM_OTHER_BUCKETS))
10 #ifdef ERR_SQUARED
11 #define STEPSIZE (0.00000001f * 0.000001f)
12 #else
13 #define STEPSIZE (0.00000001f)
14 #endif
15
16 #define clz(x) __builtin_clz(x)
17 #define ilog2(x) (sizeof(x) * 8 - clz(x) - 1)
18 #define coord(x,y,w) ((x)+((y)*(w)))
19
update_result(const uint64_t * buckets,uint64_t ccc,const double * mat,double * res)20 void update_result(const uint64_t *buckets, uint64_t ccc, const double *mat, double *res)
21 {
22 for (int y = 0; y < NUM_TOTAL_BUCKETS; y++) {
23 double addend = 0.0;
24 for (int x = 0; x < NUM_TOTAL_BUCKETS; x++) {
25 addend += mat[coord(x, y, NUM_TOTAL_BUCKETS)] * (double)buckets[x];
26 }
27 addend *= (double)ccc;
28 res[y] += addend;
29 }
30 }
31
read_matrix(const char * fn,double * mat)32 void read_matrix(const char *fn, double *mat)
33 {
34 FILE *f = fopen(fn, "r");
35 for (int y = 0; y < NUM_TOTAL_BUCKETS; y++) {
36 for (int x = 0; x < NUM_TOTAL_BUCKETS; x++) {
37 float curr;
38 fscanf(f, "%f", &curr);
39 mat[x + y * NUM_TOTAL_BUCKETS] = curr;
40 }
41 }
42 fclose(f);
43 }
44
count_coeffs(const int16_t * buf,uint32_t size,uint64_t * buckets,uint64_t * num_signs)45 void count_coeffs(const int16_t *buf, uint32_t size, uint64_t *buckets, uint64_t *num_signs)
46 {
47 uint32_t i;
48 for (i = 0; i < size; i++) {
49 int16_t curr = buf[i];
50 int16_t is_signed = curr >> 15;
51 *num_signs += (is_signed & 1);
52
53 uint16_t abs = (curr ^ is_signed) - is_signed;
54 if (abs >= NUM_COEFF_BUCKETS)
55 abs = NUM_COEFF_BUCKETS - 1;
56
57 buckets[abs]++;
58 }
59 }
60
is_power_of_two(uint32_t u)61 static inline int is_power_of_two(uint32_t u)
62 {
63 return (u & (u - 1)) == 0;
64 }
65
process_rdcosts(FILE * in,FILE * out,const double * mat)66 int process_rdcosts(FILE *in, FILE *out, const double *mat)
67 {
68 void *buf = malloc(BUFSZ);
69 uint32_t *u32buf = (uint32_t *)buf;
70 int16_t *i16buf = (int16_t *)buf;
71 int rv = 0;
72
73 double res[NUM_TOTAL_BUCKETS] = {0.0};
74
75 while (!feof(in)) {
76 uint32_t size, ccc, size_sqrt;
77 uint64_t cg_buckets[NUM_TOTAL_BUCKETS] = {0};
78 uint64_t cg_num_signs = 0;
79 size_t n_read;
80
81 n_read = fread(buf, sizeof(uint32_t), 2, in);
82 size = u32buf[0];
83 ccc = u32buf[1];
84
85 // Can't rely on feof() alone when reading from a pipe that might only get
86 // closed long after the last data has been poured in
87 if (n_read == 0) {
88 break;
89 }
90 if (feof(in) || n_read < 2) {
91 fprintf(stderr, "Unexpected EOF when reading header, managed still to read %u u32's\n", n_read);
92 rv = 1;
93 goto out;
94 }
95 if (!is_power_of_two(size)) {
96 fprintf(stderr, "Errorneous block size %u\n", size);
97 rv = 1;
98 goto out;
99 }
100
101 size_sqrt = 1 << (ilog2(size) >> 1);
102
103 n_read = fread(buf, sizeof(int16_t), size, in);
104 if (n_read != size) {
105 fprintf(stderr, "Unexpected EOF when reading block, managed still to read %u i16's\n", n_read);
106 rv = 1;
107 goto out;
108 }
109
110 count_coeffs(i16buf, size, cg_buckets, &cg_num_signs);
111 update_result(cg_buckets, ccc, mat, res);
112 }
113
114 for (int y = 0; y < NUM_TOTAL_BUCKETS; y++)
115 fprintf(out, "%g\n", (float)(res[y]));
116
117 out:
118 free(buf);
119 return rv;
120 }
121
main(int ar,char ** av)122 int main(int ar, char **av)
123 {
124 double mat[NUM_TOTAL_BUCKETS * NUM_TOTAL_BUCKETS] = {0.0};
125 if (ar != 2) {
126 fprintf(stderr, "gib matrix plz\n");
127 return 1;
128 }
129 read_matrix(av[1], mat);
130 return process_rdcosts(stdin, stdout, mat);
131 }
132
133