impl-mpfr/impl-mpfr.c

   1 #include <errno.h>
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <string.h>
   5
   6 #include <unistd.h>
   7
   8 #include <mpfr.h>
   9
  10 /* Whether we're looking at 32-bit or 64-bit floats */
  11 typedef enum { P_SINGLE, P_DOUBLE } precision;
  12
  13 /* What type of arguments we expect, and what we'll give back. */
  14 typedef enum {
  15         /* */
  16         A_UNKNOWN,
  17         A__FLT__FLT,
  18         A__FLT_FLT_FLT_RND__FLT,
  19         A__FLT_RND__FLT,
  20 } argtype;
  21
  22 /* Types of functions we could call */
  23 typedef int (*f__flt__flt)(mpfr_ptr, mpfr_srcptr);
  24 typedef int (*f__flt_flt_flt_rnd__flt)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
  25                                        mpfr_srcptr, mpfr_rnd_t);
  26 typedef int (*f__flt_rnd__flt)(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
  27
  28 /* Wrapper around a function pointer */
  29 typedef struct {
  30         /* */
  31         precision p;
  32         argtype a;
  33
  34         union {
  35                 /* */
  36                 f__flt__flt flt__flt;
  37                 f__flt_flt_flt_rnd__flt flt_flt_flt_rnd__flt;
  38                 f__flt_rnd__flt flt_rnd__flt;
  39         } f;
  40
  41 } action;
  42
  43 void usage(void)
  44 {
  45         fprintf(stderr,
  46                 "usage: impl-mpfr [-s|-d] -f <function_name> -n <num_inputs>\n");
  47         _exit(1);
  48 }
  49
  50 void determine_function(const char *f, action *a)
  51 {
  52         if (!strcmp(f, "zzzzzzzzz")) {
  53                 a->a = A_UNKNOWN;
  54         } else if (!strcmp(f, "id")) {
  55                 a->a = A__FLT_RND__FLT;
  56                 a->f.flt_rnd__flt = mpfr_set;
  57         } else if (!strcmp(f, "ceil")) {
  58                 a->a = A__FLT__FLT;
  59                 a->f.flt__flt = mpfr_ceil;
  60         } else if (!strcmp(f, "cos")) {
  61                 a->a = A__FLT_RND__FLT;
  62                 a->f.flt_rnd__flt = mpfr_cos;
  63         } else if (!strcmp(f, "floor")) {
  64                 a->a = A__FLT__FLT;
  65                 a->f.flt__flt = mpfr_floor;
  66         } else if (!strcmp(f, "exp")) {
  67                 a->a = A__FLT_RND__FLT;
  68                 a->f.flt_rnd__flt = mpfr_exp;
  69         } else if (!strcmp(f, "expm1")) {
  70                 a->a = A__FLT_RND__FLT;
  71                 a->f.flt_rnd__flt = mpfr_expm1;
  72         } else if (!strcmp(f, "fma")) {
  73                 a->a = A__FLT_FLT_FLT_RND__FLT;
  74                 a->f.flt_flt_flt_rnd__flt = mpfr_fma;
  75         } else if (!strcmp(f, "log")) {
  76                 a->a = A__FLT_RND__FLT;
  77                 a->f.flt_rnd__flt = mpfr_log;
  78         } else if (!strcmp(f, "log1p")) {
  79                 a->a = A__FLT_RND__FLT;
  80                 a->f.flt_rnd__flt = mpfr_log1p;
  81         } else if (!strcmp(f, "sin")) {
  82                 a->a = A__FLT_RND__FLT;
  83                 a->f.flt_rnd__flt = mpfr_sin;
  84         } else if (!strcmp(f, "sqrt")) {
  85                 a->a = A__FLT_RND__FLT;
  86                 a->f.flt_rnd__flt = mpfr_sqrt;
  87         } else if (!strcmp(f, "trunc")) {
  88                 a->a = A__FLT__FLT;
  89                 a->f.flt__flt = mpfr_trunc;
  90         } else {
  91                 fprintf(stderr, "impl-mpfr: unknown function \"%s\"\n", f);
  92                 _exit(1);
  93         }
  94 }
  95
  96 void read_buf(char *b, ssize_t len)
  97 {
  98         ssize_t r;
  99         ssize_t total = 0;
 100
 101         while (total < len) {
 102                 r = read(0, (b + total), (len - total));
 103
 104                 if (!r) {
 105                         /* EOF */
 106                         _exit(0);
 107                 } else if (r == -1) {
 108                         perror("impl-mpfr: read");
 109                         _exit(1);
 110                 } else {
 111                         total += r;
 112                 }
 113         }
 114 }
 115
 116 void write_buf(const char *b, ssize_t len)
 117 {
 118         ssize_t r;
 119         ssize_t total = 0;
 120
 121         while (total < len) {
 122                 r = write(1, (b + total), (len - total));
 123
 124                 if (r == -1) {
 125                         perror("impl-mpfr: write");
 126                         _exit(1);
 127                 } else {
 128                         total += r;
 129                 }
 130         }
 131 }
 132
 133 size_t input_width(argtype a, precision p)
 134 {
 135         size_t w = (p == P_SINGLE) ? 4 : 8;
 136
 137         switch (a) {
 138         case A_UNKNOWN:
 139                 break;
 140         case A__FLT__FLT:
 141
 142                 return 1 * w;
 143         case A__FLT_FLT_FLT_RND__FLT:
 144
 145                 return 3 * w;
 146         case A__FLT_RND__FLT:
 147
 148                 return 1 * w;
 149         }
 150
 151         return (size_t) -1;
 152 }
 153
 154 size_t output_width(argtype a, precision p)
 155 {
 156         size_t w = (p == P_SINGLE) ? 4 : 8;
 157
 158         switch (a) {
 159         case A_UNKNOWN:
 160                 break;
 161         case A__FLT__FLT:
 162
 163                 return 1 * w;
 164         case A__FLT_FLT_FLT_RND__FLT:
 165
 166                 return 1 * w;
 167         case A__FLT_RND__FLT:
 168
 169                 return 1 * w;
 170         }
 171
 172         return (size_t) -1;
 173 }
 174
 175 void io_loop(action a, size_t n)
 176 {
 177         char *in_buf = 0;
 178         char *out_buf = 0;
 179         size_t in_sz = input_width(a.a, a.p);
 180         size_t out_sz = output_width(a.a, a.p);
 181         mpfr_t x1;
 182         mpfr_t x2;
 183         mpfr_t x3;
 184         mpfr_t y;
 185
 186         if ((in_sz * n) / n != in_sz) {
 187                 fprintf(stderr, "impl-libc: input length overflow\n");
 188                 _exit(1);
 189         }
 190
 191         if ((out_sz * n) / n != out_sz) {
 192                 fprintf(stderr, "impl-libc: output length overflow\n");
 193                 _exit(1);
 194         }
 195
 196         if (!(in_buf = malloc(in_sz * n))) {
 197                 perror("impl-libc: malloc");
 198                 _exit(1);
 199         }
 200
 201         if (!(out_buf = malloc(out_sz * n))) {
 202                 perror("impl-libc: malloc");
 203                 _exit(1);
 204         }
 205
 206         /* I'm pretty sure 53 precision would be enough */
 207         mpfr_init2(x1, 75);
 208         mpfr_init2(x2, 75);
 209         mpfr_init2(x3, 75);
 210         mpfr_init2(y, 75);
 211
 212         while (1) {
 213                 read_buf(in_buf, in_sz * n);
 214
 215                 switch (a.a) {
 216                 case A_UNKNOWN:
 217                         fprintf(stderr, "impl-libc: impossible\n");
 218                         _exit(1);
 219                         break;
 220                 case A__FLT__FLT:
 221
 222                         switch (a.p) {
 223                         case P_SINGLE:
 224
 225                                 for (size_t j = 0; j < n; ++j) {
 226                                         float *xf1 = (float *) (in_buf +
 227                                                                 (in_sz * j));
 228                                         float *yf = (float *) (out_buf +
 229                                                                (out_sz * j));
 230
 231                                         mpfr_set_flt(x1, *xf1, MPFR_RNDN);
 232                                         a.f.flt__flt(y, x1);
 233                                         *yf = mpfr_get_flt(y, MPFR_RNDN);
 234                                 }
 235
 236                                 break;
 237                         case P_DOUBLE:
 238
 239                                 for (size_t j = 0; j < n; ++j) {
 240                                         double *xf1 = (double *) (in_buf +
 241                                                                   (in_sz * j));
 242                                         double *yf = (double *) (out_buf +
 243                                                                  (out_sz * j));
 244
 245                                         mpfr_set_d(x1, *xf1, MPFR_RNDN);
 246                                         a.f.flt__flt(y, x1);
 247                                         *yf = mpfr_get_d(y, MPFR_RNDN);
 248                                 }
 249
 250                                 break;
 251                         }
 252
 253                         break;
 254                 case A__FLT_FLT_FLT_RND__FLT:
 255
 256                         switch (a.p) {
 257                         case P_SINGLE:
 258
 259                                 for (size_t j = 0; j < n; ++j) {
 260                                         float *xf1 = (float *) (in_buf +
 261                                                                 (in_sz * j));
 262                                         float *xf2 = (float *) (in_buf +
 263                                                                 (in_sz * j) +
 264                                                                 4);
 265                                         float *xf3 = (float *) (in_buf +
 266                                                                 (in_sz * j) +
 267                                                                 8);
 268                                         float *yf = (float *) (out_buf +
 269                                                                (out_sz * j));
 270
 271                                         mpfr_set_flt(x1, *xf1, MPFR_RNDN);
 272                                         mpfr_set_flt(x2, *xf2, MPFR_RNDN);
 273                                         mpfr_set_flt(x3, *xf3, MPFR_RNDN);
 274                                         a.f.flt_flt_flt_rnd__flt(y, x1, x2, x3,
 275                                                                  MPFR_RNDN);
 276                                         *yf = mpfr_get_flt(y, MPFR_RNDN);
 277                                 }
 278
 279                                 break;
 280                         case P_DOUBLE:
 281
 282                                 for (size_t j = 0; j < n; ++j) {
 283                                         double *xf1 = (double *) (in_buf +
 284                                                                   (in_sz * j));
 285                                         double *xf2 = (double *) (in_buf +
 286                                                                   (in_sz * j) +
 287                                                                   8);
 288                                         double *xf3 = (double *) (in_buf +
 289                                                                   (in_sz * j) +
 290                                                                   16);
 291                                         double *yf = (double *) (out_buf +
 292                                                                  (out_sz * j));
 293
 294                                         mpfr_set_d(x1, *xf1, MPFR_RNDN);
 295                                         mpfr_set_d(x2, *xf2, MPFR_RNDN);
 296                                         mpfr_set_d(x3, *xf3, MPFR_RNDN);
 297                                         a.f.flt_flt_flt_rnd__flt(y, x1, x2, x3,
 298                                                                  MPFR_RNDN);
 299                                         *yf = mpfr_get_d(y, MPFR_RNDN);
 300                                 }
 301
 302                                 break;
 303                         }
 304
 305                         break;
 306                 case A__FLT_RND__FLT:
 307
 308                         switch (a.p) {
 309                         case P_SINGLE:
 310
 311                                 for (size_t j = 0; j < n; ++j) {
 312                                         float *xf1 = (float *) (in_buf +
 313                                                                 (in_sz * j));
 314                                         float *yf = (float *) (out_buf +
 315                                                                (out_sz * j));
 316
 317                                         mpfr_set_flt(x1, *xf1, MPFR_RNDN);
 318                                         a.f.flt_rnd__flt(y, x1, MPFR_RNDN);
 319                                         *yf = mpfr_get_flt(y, MPFR_RNDN);
 320                                 }
 321
 322                                 break;
 323                         case P_DOUBLE:
 324
 325                                 for (size_t j = 0; j < n; ++j) {
 326                                         double *xf1 = (double *) (in_buf +
 327                                                                   (in_sz * j));
 328                                         double *yf = (double *) (out_buf +
 329                                                                  (out_sz * j));
 330
 331                                         mpfr_set_d(x1, *xf1, MPFR_RNDN);
 332                                         a.f.flt_rnd__flt(y, x1, MPFR_RNDN);
 333                                         *yf = mpfr_get_d(y, MPFR_RNDN);
 334                                 }
 335
 336                                 break;
 337                         }
 338
 339                         break;
 340                 }
 341
 342                 write_buf(out_buf, out_sz * n);
 343         }
 344 }
 345
 346 int main(int argc, char **argv)
 347 {
 348         int c = 0;
 349         action a = { .p = P_SINGLE };
 350         long long n = 0;
 351
 352         while ((c = getopt(argc, argv, "sdf:n:")) != -1) {
 353                 switch (c) {
 354                 case 's':
 355                         a.p = P_SINGLE;
 356                         break;
 357                 case 'd':
 358                         a.p = P_DOUBLE;
 359                         break;
 360                 case 'f':
 361                         determine_function(optarg, &a);
 362                         break;
 363                 case 'n':
 364                         errno = 0;
 365                         n = strtoll(optarg, 0, 0);
 366
 367                         if (errno) {
 368                                 perror("impl-libc: unparsable");
 369
 370                                 return 1;
 371                         }
 372
 373                         break;
 374                 default:
 375                         usage();
 376                         break;
 377                 }
 378         }
 379
 380         if (a.a == A_UNKNOWN) {
 381                 usage();
 382         }
 383
 384         io_loop(a, n);
 385 }