cfloat.c

   1 /*
   2  * Floating point functions that are difficult or impossible to implement
   3  * in pure Haskell.
   4  *
   5  * Copyright (C) 2009-2010 Nick Bowler.
   6  *
   7  * License BSD2:  2-clause BSD license.  See LICENSE for full terms.
   8  * This is free software: you are free to change and redistribute it.
   9  * There is NO WARRANTY, to the extent permitted by law.
  10  */
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <errno.h>
  14 #include <math.h>
  15 #include <fenv.h>
  16
  17 #include "cfloat.h"
  18
  19 #pragma STDC FENV_ACCESS ON
  20
  21 int double_format(char *buf, char spec, int precision, double val)
  22 {
  23         char fmt[] = "%.*f";
  24         fmt[3] = spec;
  25
  26         if (buf == NULL)
  27                 return snprintf(NULL, 0, fmt, precision, val);
  28         return sprintf(buf, fmt, precision, val);
  29 }
  30
  31 double double_signum(double val)
  32 {
  33         if (signbit(val))
  34                 return -1;
  35         return 1;
  36 }
  37
  38 float float_signum(float val)
  39 {
  40         if (signbit(val))
  41                 return -1;
  42         return 1;
  43 }
  44
  45 int double_classify(double val)
  46 {
  47         switch (fpclassify(val)) {
  48         case FP_INFINITE:
  49                 return 0;
  50         case FP_NAN:
  51                 return 1;
  52         case FP_NORMAL:
  53                 return 2;
  54         case FP_SUBNORMAL:
  55                 return 3;
  56         case FP_ZERO:
  57                 return 4;
  58         }
  59
  60         return -1;
  61 }
  62
  63 int float_classify(float val)
  64 {
  65         switch (fpclassify(val)) {
  66         case FP_INFINITE:
  67                 return 0;
  68         case FP_NAN:
  69                 return 1;
  70         case FP_NORMAL:
  71                 return 2;
  72         case FP_SUBNORMAL:
  73                 return 3;
  74         case FP_ZERO:
  75                 return 4;
  76         }
  77
  78         return -1;
  79 }
  80
  81 int double_compare(double a, double b)
  82 {
  83         if (isless(a, b))
  84                 return 0;
  85         if (a == b)
  86                 return 1;
  87         if (isgreater(a, b))
  88                 return 2;
  89         if (isunordered(a, b))
  90                 return 3;
  91         return -1;
  92 }
  93
  94 int float_compare(float a, float b)
  95 {
  96         if (isless(a, b))
  97                 return 0;
  98         if (a == b)
  99                 return 1;
 100         if (isgreater(a, b))
 101                 return 2;
 102         if (isunordered(a, b))
 103                 return 3;
 104         return -1;
 105 }
 106
 107 int fenv_set_roundmode(int mode)
 108 {
 109         int cmode;
 110
 111         switch (mode) {
 112 #ifdef FE_TONEAREST
 113         case 0:
 114                 cmode = FE_TONEAREST;
 115                 break;
 116 #endif
 117 #ifdef FE_UPWARD
 118         case 1:
 119                 cmode = FE_UPWARD;
 120                 break;
 121 #endif
 122 #ifdef FE_DOWNWARD
 123         case 2:
 124                 cmode = FE_DOWNWARD;
 125                 break;
 126 #endif
 127 #ifdef FE_TOWARDZERO
 128         case 3:
 129                 cmode = FE_TOWARDZERO;
 130                 break;
 131 #endif
 132         default:
 133                 return -1;
 134         }
 135
 136         return fesetround(cmode);
 137 }
 138
 139 int fenv_get_roundmode(void)
 140 {
 141         int cmode = fegetround();
 142
 143         switch (cmode) {
 144 #ifdef FE_TONEAREST
 145         case FE_TONEAREST:
 146                 return 0;
 147 #endif
 148 #ifdef FE_UPWARD
 149         case FE_UPWARD:
 150                 return 1;
 151 #endif
 152 #ifdef FE_DOWNWARD
 153         case FE_DOWNWARD:
 154                 return 2;
 155 #endif
 156 #ifdef FE_TOWARDZERO
 157         case FE_TOWARDZERO:
 158                 return 3;
 159 #endif
 160         default:
 161                 return -1;
 162         }
 163 }
 164
 165 unsigned fenv_test_excepts(void)
 166 {
 167         int raw_excepts = fetestexcept(FE_ALL_EXCEPT);
 168         unsigned excepts = 0;
 169
 170 #ifdef FE_DIVBYZERO
 171         if (raw_excepts & FE_DIVBYZERO) excepts |= 0x01;
 172 #endif
 173 #ifdef FE_INEXACT
 174         if (raw_excepts & FE_INEXACT)   excepts |= 0x02;
 175 #endif
 176 #ifdef FE_INVALID
 177         if (raw_excepts & FE_INVALID)   excepts |= 0x04;
 178 #endif
 179 #ifdef FE_OVERFLOW
 180         if (raw_excepts & FE_OVERFLOW)  excepts |= 0x08;
 181 #endif
 182 #ifdef FE_UNDERFLOW
 183         if (raw_excepts & FE_UNDERFLOW) excepts |= 0x10;
 184 #endif
 185
 186         return excepts;
 187 }
 188
 189 int fenv_raise_excepts(unsigned excepts)
 190 {
 191         int raw_excepts = 0;
 192
 193 #ifdef FE_DIVBYZERO
 194         if (excepts & 0x01) raw_excepts |= FE_DIVBYZERO;
 195 #endif
 196 #ifdef FE_INEXACT
 197         if (excepts & 0x02) raw_excepts |= FE_INEXACT;
 198 #endif
 199 #ifdef FE_INVALID
 200         if (excepts & 0x04) raw_excepts |= FE_INVALID;
 201 #endif
 202 #ifdef FE_OVERFLOW
 203         if (excepts & 0x08) raw_excepts |= FE_OVERFLOW;
 204 #endif
 205 #ifdef FE_UNDERFLOW
 206         if (excepts & 0x10) raw_excepts |= FE_UNDERFLOW;
 207 #endif
 208         return feraiseexcept(raw_excepts);
 209 }