Octave, bc, and gnuplot scripts supporting math dev

This is a collection of scripts that are useful in developing math
approximations.

Refs: gh-206

Signed-off-by: Matthias Kretz <[email protected]>

Author: mattkretz

math/cosf.m

@@ -0,0 +1,12 @@
+function y = cosf (x)
+  x2 = x.*x;
+  y = 0;
+  y = y .* x2 + to_float(0x573F9F, -45);
+  y = y .* x2 - to_float(0x49CBA5, -37);
+  y = y .* x2 + to_float(0x0F76C7, -29);
+  y = y .* x2 - to_float(0x13F27E, -22);
+  y = y .* x2 + to_float(0x500D01, -16);
+  y = y .* x2 - to_float(0x360B61, -10);
+  y = y .* x2 + to_float(0x2AAAAB,  -5);
+  y = y .* (x2 .* x2) - .5 * x2 + 1.;
+endfunction

math/floattool

@@ -0,0 +1,139 @@
+#!/bin/zsh
+
+function usage() {
+  cat <<EOF
+Usage: $0 [options] <input number>
+
+This tool splits an arbitrary number into IEEE 768 floats.
+Options:
+  --help|-h        this message
+  -m <bits>        mantissa bits, 23 for sp and 52 for dp (default 23)
+  -n <count>       number of constants to split the input into (default 2)
+  -z <bits>        number of trailing zero bits in the first count-1 constants (default 12)
+  -p               only output positive numbers (or rather numbers with the same sign as the input number)
+  -o <offset>      add this offset to the mantissa of the first value
+  -f <function>    use the argument as output function name for easier cut and paste (default depends on -m)
+
+Available functions:
+  fak(x) : factorial of x
+  floor(x)
+  round(x)
+  log2(x)
+  abs(x)
+  sign(x)
+  exponent(x)
+  l(x)
+  s(x)
+  c(x)
+
+EOF
+}
+
+bits=23
+count=2
+zerobits=12
+mantissaReduction=round
+offset=0
+const_fun=
+
+while (( $# > 0 )); do
+  case "$1" in
+    --help|-h) usage; exit ;;
+    -m)  bits=$2; shift ;;
+    -m*) bits=${1#-m} ;;
+    -n)  count=$2; shift ;;
+    -n*) count=${1#-n} ;;
+    -z)  zerobits=$2; shift ;;
+    -z*) zerobits=${1#-z} ;;
+    -o)  offset=$2; shift ;;
+    -o*) offset=${1#-o} ;;
+    -f)  const_fun=$2; shift ;;
+    -f*) const_fun=${1#-o} ;;
+    -p)  mantissaReduction=floor ;;
+    *)   input="$input $1" ;;
+  esac
+  shift
+done
+
+if (($bits == 23)) && [[ -z "$const_fun" ]]; then
+  const_fun="const auto c## = Vc::Detail::floatConstant"
+elif [[ -z "$const_fun" ]]; then
+  const_fun="const auto c## = Vc::Detail::doubleConstant"
+fi
+
+BC_LINE_LENGTH=0 bc -l <<EOF
+scale=400
+pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936073
+oneoverlog2=1/l(2)
+scale=100
+define fak(x) { r=x; while(x>1) { r *= --x; }; return r; }
+define floor(x) { tmp=scale; scale=0; n = x/1; scale=tmp; return n; }
+define round(x) { return floor(x + 0.5); }
+define log2(x)  { return l(x)*oneoverlog2; }
+define abs(x)   { if(x < 0) return -x; return x; }
+define sign(x)  { tmp=scale; scale=0; n = x/abs(x); scale=tmp; return n; }
+define exponent(x) { x = abs(x); e = floor(log2(x)); while(x*2^-e < 1) e -= 1; while(x*2^-e >= 2) e += 1; return e; }
+define setprecision(x, p) { return round(x * 2^(p-exponent(x))) * 2^(exponent(x)-p); }
+define float(x) { return setprecision(x,${bits}); }
+define void printfloat(s,m,e) {
+  print "${const_fun}<", s, ", 0x"
+  obase=16
+  print m
+  obase=10
+  print ", ", e, ">(); // ", s*(m*2^(e-${bits})+2^e), "\n"
+}
+
+/* unused, but can be used to print full binary representation of floats */
+define void printfloat_int(x) {
+  h=0
+  if (x < 0) {
+    h = 2^31
+    x = -x;
+  }
+  e=exponent(x)
+  h+=${mantissaReduction}(x*2^(${bits}-e)-2^${bits})
+  h+=(e+127)*2^23
+  obase=16
+  print "0x", h
+  obase=10
+}
+
+in=${input}
+print "// ", in, "\n"
+x=in
+shift1=${bits}-${zerobits}
+shift2=${bits}-shift1
+for(i = 1; i < ${count}; ++i) {
+  s=sign(x)
+  x=abs(x)
+  e=exponent(x)
+  m=(${mantissaReduction}(x*2^(shift1-e)-2^shift1) + ${offset})*(2^shift2)
+  printfloat(s,m,e)
+  x=s*(x-(m*(2^-${bits})+1)*(2^e))
+  if(x == 0) {
+    break
+  }
+}
+if(x == 0) {
+  print "no remainder.\n"
+} else {
+  s=sign(x)
+  x=abs(x)
+  e=exponent(x)
+  m=${mantissaReduction}(x*2^(${bits}-e)-2^${bits})
+  printfloat(s,m,e)
+  x=s*(x-(m*(2^-${bits})+1)*(2^e))
+  if(x == 0) {
+    print "no remainder.\n"
+  } else {
+    err=abs(x/in)
+    e=floor(l(err)/l(10))-1
+    ulp=x * 2^(23 - exponent(in))
+    scale=1
+    print "relative error: ", (err*10^-e)/1, "e", e, "\n"
+    print (ulp*1000+0.5)/1*0.001, " ulp\n"
+  }
+}
+EOF
+
+# vim: sw=2 et

math/fold_pi4_single.m

@@ -0,0 +1,17 @@
+function x = fold_pi4_single(x)
+  x = round_to_float(x);
+  expected = round_to_float(x - floor(x / (pi/2) + .5)*(pi/2));
+  y = round_to_float(floor(x * round_to_float(2/pi) + .5)) * 2;
+  printf("y = %i = %s\n", y, hexfloat(y));
+  while(y > 0)
+    p2 = 2^floor(log2(y));
+    y -= p2;
+    diff = p2 * round_to_float(pi/4);
+    x -= diff;
+    printf("(x -= %s) -> %s\n", hexfloat(diff), hexfloat(x));
+  endwhile
+  diff = p2 * round_to_float(pi/4 - round_to_float(pi/4));
+  x -= diff;
+  printf("(x -= %s) -> %s\n", hexfloat(diff), hexfloat(x));
+  printf("expected result:         %s\n", hexfloat(expected));
+endfunction

math/generate_sindat.sh

@@ -0,0 +1,36 @@
+#!/bin/zsh
+
+cd "${0%/*}"
+
+function usage() {
+  cat <<EOF
+Usage: $0
+
+Options:
+  --help|-h        this message
+
+EOF
+}
+
+while (( $# > 0 )); do
+  case "$1" in
+    --help|-h) usage; exit ;;
+  esac
+  shift
+done
+
+octave <<EOF
+ff=fopen("sincosf.dat", "w");
+fd=fopen("sincosd.dat", "w");
+
+printf("Generating 100000 dp & sp sincos values [0, π/4]     ");
+for i = 1:100000
+  x=i * pi / 400000;
+  fprintf(fd, "%.60f\t%.60f\t%.60f\n", x, sin(x), cos(x));
+  x=round_to_float(x);
+  fprintf(ff, "%.40f\t%.50f\t%.50f\n", x, sin(x), cos(x));
+  printf("\033[4D%3i%%", floor(i/1000));
+endfor
+EOF
+
+# vim: sw=2 et

math/hexdouble.m

@@ -0,0 +1,11 @@
+function r = hexdouble(x)
+  if (signbit(x))
+    s = "-";
+  else
+    s = "";
+  endif
+  x = abs(x);
+  e = floor(log2(x));
+  m = round((x*2^-e - 1) * 2^52);
+  r = sprintf("%s0x1.%013xp%i", s, m, e);
+endfunction

math/hexfloat.m

@@ -0,0 +1,15 @@
+function r = hexfloat(x)
+  if (x == 0)
+    r = "0x0.000000p0";
+    return;
+  endif
+  if (signbit(x))
+    s = "-";
+  else
+    s = "";
+  endif
+  x = abs(x);
+  e = floor(log2(x));
+  m = roundb((x*2^-e - 1) * 2^23) * 2;
+  r = sprintf("%s0x1.%06xp%i", s, m, e);
+endfunction

math/mycos.m

@@ -0,0 +1,12 @@
+function y = mycos (x)
+  x2 = x.*x;
+  y = 0;
+  y = y .* x2 + to_double(0xAC00000000000, -45);
+  y = y .* x2 - to_double(0x9394000000000, -37);
+  y = y .* x2 + to_double(0x1EED8C0000000, -29);
+  y = y .* x2 - to_double(0x27E4FB7400000, -22);
+  y = y .* x2 + to_double(0xA01A01A018000, -16);
+  y = y .* x2 - to_double(0x6C16C16C16C00, -10);
+  y = y .* x2 + to_double(0x5555555555554,  -5);
+  y = y .* (x2 .* x2) - .5 * x2 + 1.;
+endfunction

math/mysin.m

@@ -0,0 +1,12 @@
+function y = mysin (x)
+  x2 = x.*x;
+  y = 0;
+  y = (y - to_double(0xACF0000000000 + 0, -41)) .* x2;
+  y = (y + to_double(0x6124400000000 + 0, -33)) .* x2;
+  y = (y - to_double(0xAE64567000000 + 0, -26)) .* x2;
+  y = (y + to_double(0x71DE3A5540000 + 0, -19)) .* x2;
+  y = (y - to_double(0xa01a01a01a000 + 0, -13)) .* x2;
+  y = (y + to_double(0x1111111111110 + 0, -7)) .* x2;
+  y = (y - to_double(0x5555555555555 + 0, -3)) .* x2;
+  y = (y + 1) .* x;
+endfunction

math/print_hexfloat.m

@@ -0,0 +1,11 @@
+function r = hexfloat(x)
+  if (signbit(x))
+    s = "-";
+  else
+    s = "";
+  endif
+  x = abs(x);
+  e = floor(log2(x));
+  m = round((x*2^-e - 1) * 2^23) * 2;
+  r = sprintf("%s0x1.%xp%i", s, m, e);
+endfunction

math/round_to_float.m

@@ -0,0 +1,10 @@
+function x = round_to_float(x)
+  if(x == 0)
+    return;
+  endif
+  s = sign(x);
+  x = abs(x);
+  e = floor(log2(x));
+  m = roundb((x.*2.^-e - 1) .* 2.^23);
+  x = s .* (m .* 2.^-23 + 1) .* 2.^e;
+endfunction