plots, code & comparison of methods

utils/plots/Makefile

@@ -0,0 +1,13 @@
+CFLAGS= -g -std=c99 -Wall -Wconversion -Wtype-limits -Wshadow -Wpointer-arith -Wunreachable-code -Werror -MMD -O0
+EXEC=CSV
+CC=gcc
+
+all: $(EXEC)
+
+$(EXEC):
+	$(CC) sr_stability.c $(CFLAGS) -o $(EXEC) -lm
+
+.PHONY: clean
+
+clean:
+	rm $(EXEC)

utils/plots/graph.R

@@ -0,0 +1,29 @@
+library("ggplot2")
+options(digits=3)
+setwd("/home/joaquintz/Desktop/projects/utils/utils/plots/")
+
+comparison <- read.table("benchmark.tsv",sep="\t",header=TRUE)
+
+comparison$err_rel <- (comparison$fsqrt-comparison$sqrt)/comparison$sqrt
+
+max_error <- max(comparison$err_rel,na.rm=TRUE)
+min_error <- min(comparison$err_rel,na.rm=TRUE)
+step_error <- (max_error-min_error)/5
+max_x_exp <- ceiling(max(log10(comparison$x),na.rm=TRUE))
+
+exponent_breaks <- function(coeficientes,pot) 
+    as.expression(parse(text=paste(coeficientes,"10^",pot,sep="")))
+
+comparison.plot <- (ggplot(data=comparison,aes(x=x,y=err_rel))
+                    + geom_line()
+                    + scale_y_continuous(limits=c(min_error,max_error),
+                                         breaks=seq(min_error,max_error,by=step_error))
+                    + scale_x_log10(breaks=10^seq(0,max_x_exp),
+                                    labels=exponent_breaks(NULL,seq(0,max_x_exp)))
+                    + theme_bw() 
+                    + ggtitle("Error relativo: sqrt vs fsqrt") 
+                    + ylab("Error relativo (%)") + xlab("x [log 10]"))
+
+print(comparison.plot)
+
+ggsave("comparison.png")

utils/plots/sr_stability.c

@@ -0,0 +1,32 @@
+#include <stdio.h>
+#include <math.h>
+#include <time.h>
+#include <stdlib.h>
+#define NUM_STEPS 100000000
+
+float fsqrt(float n) {
+    unsigned i = *(unsigned*)&n;
+    i = (i + 0x3f800000) >> 1;
+    float y = *(float*)&i;
+    //y = y*0.5f + n/(2*y);
+    //y = y*0.5f + n/(2*y);  
+    return y;
+}
+void next(size_t *n) {
+	if(*n<100)*n=*n+1;
+	else if(*n<10000)*n=*n+10;
+	else if(*n<1000000) *n=*n+100;
+	else *n=*n+1000;
+}
+
+int main() {
+    float res1, res2, x;
+	printf("x\tsqrt\tfsqrt\n");
+    for(size_t i= 0; i < NUM_STEPS ; next(&i)) {
+       	x = (float) i;        
+        res2 = (float)sqrt(x);
+        res1 = fsqrt(x);
+        printf("%.6f\t%.6f\t%.6f\n",x,res2,res1);
+    }
+    return 0;
+}