example of new mathfunc (fast) "invsqrt" (with given accuracy) as c-module + tcl examples of usage

/*
 * tclfastinvsqrt.c --
 * 
 * Test module to get native handle of channel.
 *
 * Compile:
 *   mingw: gcc -O2 -DUSE_TCL_STUBS=1 -I$tcl/win  -I$tcl/generic tclfastinvsqrt.c -shared -o tclfastinvsqrt.dll libtclstub86.a
 *   *nix:  gcc -O2 -DUSE_TCL_STUBS=1 -I$tcl/unix -I$tcl/generic tclfastinvsqrt.c -shared -o tclfastinvsqrt.so  libtclstub86.a
 * 
 * Usage:
 *   $ tclsh86
 *   % load tclfastinvsqrt
 *   % expr {invsqrt(10)}
 */

#include "tcl.h"

static inline double
FastInvSqrt(double x, int accuracy) {
  char *buf = (char *)&x;
  double xhalf = 0.5f * x;
  Tcl_WideInt i = *(Tcl_WideInt*)buf;
  i = 0x5fe6eb50c7b537a9 - (i >> 1);
  buf = (char *)&i;
  x = *(double*)buf;
  while (accuracy--) {
    x = x*(1.5-(xhalf*x*x));
  }
  return x;
}

int
FastInvSqrtCmd(
  ClientData dummy,
  Tcl_Interp* interp,
  int objc,
  Tcl_Obj * const objv[]
) {
  double dbl;
  int accuracy = 4;

  if (objc < 2 || objc > 3) {
    Tcl_WrongNumArgs(interp, 1, objv, "value ?accuracy?");
    return TCL_ERROR;
  }
  
  if (Tcl_GetDoubleFromObj(interp, objv[1], &dbl) != TCL_OK) {
    return TCL_ERROR;
  }
  if (objc > 2 && Tcl_GetIntFromObj(interp, objv[2], &accuracy) != TCL_OK) {
    return TCL_ERROR;
  }
  
  dbl = FastInvSqrt(dbl, accuracy);

  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(dbl));
  return TCL_OK;
}

int
Tclfastinvsqrt_Init(Tcl_Interp *interp) {
  if (!Tcl_InitStubs(interp, "8.5", 0)) {
    return TCL_ERROR;
  }
  Tcl_CreateObjCommand(interp, "::tcl::mathfunc::invsqrt", FastInvSqrtCmd, NULL, NULL);
  return TCL_OK;
}

/*

load tclfastinvsqrt

# speed vs accuracy:

set tm 500
for {set x 2} {$x < 100000} {set x [expr {$x * 10}]} {
  puts $x:\t[timerate { expr { 1 / sqrt($x) } } $tm]
  puts \t[timerate { expr { invsqrt($x,1) } } $tm]
  puts \t[timerate { expr { invsqrt($x) } } $tm]
}

# accuracy :

for {set x 2} {$x < 100000} {set x [expr {$x * 10}]} {
  puts $x:\t[expr { 1 / sqrt($x) }]
  puts \t[expr { invsqrt($x) }]
}
for {set x 2} {$x < 100000} {set x [expr {$x * 10}]} {
  puts $x:\t[expr { 1 / sqrt($x) }]
  puts \t[expr { invsqrt($x,1) }]
}

# Herbie example:

for {set x 2} {$x < 100000} {set x [expr {$x * 10}]} {
  puts $x:\t[expr { sqrt($x+1) - sqrt($x) }]
  puts \t[expr { 1 / (sqrt($x+1) + sqrt($x)) }]
  puts \t[expr { invsqrt( (2*$x + 1 + 2/invsqrt(($x+1)*$x)) ) }]
}

*/