#!/bin/bash

# An implementation of IEEE Single Precision Floating Point 
# entirely in Bash. Working, but probably not conforming :P
#
# - Vidar 'koala_man' Holen
#   2006-04-16



#Sample output (takes about 1.5min on 1ghz):
#
#Demo app for koala_man's IEEE Single Precision Floating Point bash script
#
#Basic numbers:
#Parsing  1.0: 1065353216 (1.0)
#Parsing -1337.666: 3299292494 (-1337.66577148)
#
#Flops:
#Doing 1000 multiply operations: 1000 ops per 15.820 secs, 63.2111 flops
#
#Scientific notation (based on exp/ln):
#Parsing   1E6: 1232348134 (999998.375)
#Parsing   1E30: 1900671618 (9.99980354E29)
#Parsing  -1E-30: 2376221322 (-1.00001072E-30)
#
#Various silly calculations: 
#E=2.71828031
#ln(2)=0.6931467
#Leibnitz' formula with 50 iterations: pi=3.18157386
#Wallis' Product with 50 iterations: pi=3.11091852


debug() {
	echo "$@" >&2
}

fp_get_sign() {
	echo $(($1>>31))
}

##set the sign -- num newsign
fp_set_sign() {
	echo $((($1&2147483647)|($2<<31)))
}

fp_get_exp() {
	echo $((($1>>23)&255))
}

fp_get_frac() {
	echo $((($1&8388607) | 8388608))
}

fp_is_nan() {
	[[ $(fp_get_exp $1) == 255 && $(fp_get_frac $1) != 8388608 ]]
}

fp_is_positive_infinity() {
	[[ $1 == $fp_CONST_PINF ]] 
}
fp_is_negative_infinity() {
	[[ $1 == $fp_CONST_NINF ]] 
}
fp_classify() {
	if fp_is_nan $1
	then
		echo 1
	elif fp_is_negative_infinity $1
	then
		echo 2
	elif fp_is_positive_infinity $1
	then
		echo 3
	else
		echo 0
	fi	
}

##is a less than b (must be normalized) -- a b
fp_lt() {
	if [[ $(fp_get_sign $1) == 1 ]] 
	then
		if [[ $(fp_get_sign $2) == 1 ]] 
		then
			fp_lt $(fp_set_sign $2 0) $(fp_set_sign $1 0)
		else
			true
		fi
	else
		if [[ $(fp_get_sign $2) == 1 ]]
		then
			false
		else
			if (( $(fp_get_exp $1) < $(fp_get_exp $2) ))
			then
				true
			elif (( $(fp_get_exp $1) > $(fp_get_exp $2) ))
			then
				false
			else 
				if (( $(fp_get_frac $1) < $(fp_get_frac $2) ))
				then
					true
				else
					false
				fi
			fi
		fi
	fi
}

## normalize and build fp -- sign exp frac
fp_normalize_build() {
	local fp_normalize_build_SIGN=$1
	local fp_normalize_build_EXP=$2
	local fp_normalize_build_FRAC=$3

	if (( fp_normalize_build_FRAC == 0 ))
	then
		fp_set_sign $fp_CONST_ZERO $fp_normalize_build_SIGN
		return
	fi

	while (( fp_normalize_build_FRAC > 8388608 )) 
	do
		(( fp_normalize_build_EXP++ ))
		(( fp_normalize_build_FRAC >>= 1 ))
	done
	if (( fp_normalize_build_EXP >= 255 )) 
	then
		echo $fp_CONST_PINF 
		return 
	fi
	while (( fp_normalize_build_FRAC < 8388608 ))
	do
		(( fp_normalize_build_EXP-- ))
		(( fp_normalize_build_FRAC <<= 1 ))
	done
	if (( fp_normalize_build_EXP <= 0 )) 
	then
		echo $fp_CONST_ZERO
		return 
	fi
	(( fp_normalize_build_FRAC &= 8388607 ))
	fp_build $fp_normalize_build_SIGN $fp_normalize_build_EXP $fp_normalize_build_FRAC
}

fp_build() {
	echo $((($1<<31)|($2<<23)|($3&8388607)))
}

fp_normalize() {
	fp_normalize_build $(fp_get_sign $1) $(fp_get_exp $1) $(fp_get_frac $1)
}


## multiply a and b -- a b
fp_mul() {
	case "$(fp_classify $1)$(fp_classify $2)" in
		01|10|11   ) echo $fp_CONST_NAN; return;; #NaN*x = NaN
		?3|3?|?2|2?) echo $(fp_set_sign $fp_CONST_PINF $(($(fp_get_sign $1)^$(fp_get_sign $2))))
			return;;
	esac
	
	local fp_mul_SIGN=$(($(fp_get_sign $1) ^ $(fp_get_sign $2)))
	local fp_mul_EXP=$(fp_get_exp  $1)
	local fp_mul_FRAC=$(fp_get_frac $1)
	
	
	(( fp_mul_EXP += $(fp_get_exp  $2) - 127 - 23 - 0 ))
	(( fp_mul_FRAC *= $(fp_get_frac $2) ))

	fp_normalize_build $fp_mul_SIGN $fp_mul_EXP $fp_mul_FRAC
}

## divide a and b -- a b
fp_div() {
	case "$(fp_classify $1)$(fp_classify $2)" in
		22|23|32|33) echo $fp_CONST_NAN; return;; #NaN*x = NaN
		2?|3?) echo $(fp_set_sign $fp_CONST_PINF $(($(fp_get_sign $1)^$(fp_get_sign $2)))) 
			return;; #+/- inf/x
	esac

#	fp_dbg $1 $2
	
	local fp_div_SIGN=$(($(fp_get_sign $1) ^ $(fp_get_sign $2)))
	local fp_div_EXP=$(fp_get_exp $1)
	local fp_div_FRAC=$(fp_get_frac $1)
	local fp_div_EXP2=$(fp_get_exp $2)

	if (( fp_div_EXP == 0 ))
	then
		if (( fp_div_EXP2 == 0 ))
		then
			echo $fp_const_NAN
			return
		fi
		fp_set_sign $fp_CONST_ZERO $fp_div_SIGN
		return
	fi
	if (( fp_div_EXP2 == 0 ))
	then
		fp_set_sign $fp_CONST_PINF $fp_div_SIGN
		return
	fi


	(( fp_div_FRAC <<= 23 )) 

	(( fp_div_EXP -= $(fp_get_exp $2) - 127 ))
	(( fp_div_FRAC /= ($(fp_get_frac $2)) ))

 
	fp_normalize_build $fp_div_SIGN $fp_div_EXP $fp_div_FRAC
}

fp_dbg() {
	for fp_dbg_I 
	do
		debug fp_dbg: $fp_dbg_I = $(fp_get_sign $fp_dbg_I) $(fp_get_exp $fp_dbg_I) $(fp_get_frac $fp_dbg_I) = $(fp_format $fp_dbg_I)
	done
}

## add a and b -- a b
fp_add() {
	case "$(fp_classify $1)$(fp_classify $2)" in
		01|10|11 ) echo $fp_CONST_NAN; return;; #NaN+x = NaN
		23|32    ) echo $fp_CONST_NAN; return;; #-inf+inf = NaN
		?3|3?    ) echo $fp_CONST_PINF; return;; #inf+x = inf
		?2|2?    ) echo $fp_CONST_NINF; return;; #-inf+x = -inf
	esac

	if fp_lt $(fp_set_sign $2 0) $(fp_set_sign $1 0)
	then
		set -- $2 $1
	fi

	fp_add_EXP1=$(fp_get_exp $1)
	fp_add_EXP2=$(fp_get_exp $2)

	if (( fp_add_EXP2 - fp_add_EXP1 >= 21 )) #$2 is too small
	then
		echo $2
		return
	fi
	
	fp_add_FRAC1=$(fp_get_frac $1)
	fp_add_FRAC2=$(fp_get_frac $2)
	fp_add_SIGN1=$(fp_get_sign $1)
	fp_add_SIGN2=$(fp_get_sign $2)

	while (( fp_add_EXP1 < fp_add_EXP2 )) 
	do
		(( fp_add_EXP1++ ))
		(( fp_add_FRAC1 >>= 1 )) 
	done

	if (( fp_add_SIGN1 ^ fp_add_SIGN2 == 0 ))
	then
		(( fp_add_FRAC2 += fp_add_FRAC1 ))
	else
		(( fp_add_FRAC2 -= fp_add_FRAC1 ))

		if (( fp_add_FRAC2 < 0 )) 
		then
			fp_add_SIGN2=$((fp_add_SIGN2^1))
		fi
	fi
	
	fp_normalize_build $fp_add_SIGN2 $fp_add_EXP2 $fp_add_FRAC2
}

fp_sub() {
#	debug fp_add $1 $(fp_set_sign $2 $(($(fp_get_sign $2)^1)))
	fp_add $1 $(fp_set_sign $2 $(($(fp_get_sign $2)^1)))
}

fp_parse() {
	local fp_parse_NUM=$1
	local fp_parse_SIGN=0
	local fp_parse_EXP10=0
	local fp_parse_EXP2=127
	local fp_parse_FRAC=0
	local fp_parse_RES=0
	
	if [[ $fp_parse_NUM == -* ]] 
	then
		fp_parse_SIGN=1
		fp_parse_NUM=${fp_parse_NUM#-}
	fi	
	if [[ $fp_parse_NUM == *E* ]] 
	then
		fp_parse_EXP10=${fp_parse_NUM##*E}
		fp_parse_NUM=${fp_parse_NUM%%E*}
	fi
	
	fp_parse_RES=${fp_parse_NUM%%.*}
	fp_parse_RES=$(fp_parse_int $fp_parse_RES)

	if [[ $fp_parse_NUM == *.* ]] 
	then
		fp_parse_DEC=${fp_parse_NUM##*.}
	else
		fp_parse_DEC=0
	fi
	
	fp_parse_TMP=${#fp_parse_DEC}
	fp_parse_TMP2=$(fp_div $(fp_parse_int $fp_parse_DEC) $(fp_parse_int $((10**$fp_parse_TMP))))
	fp_parse_RES=$(fp_add $fp_parse_RES $fp_parse_TMP2)

	fp_parse_RES=$(fp_set_sign $fp_parse_RES $fp_parse_SIGN)

	if [[ $fp_parse_EXP10 == 0 ]]
	then
		echo $fp_parse_RES
	else
		echo $(fp_mulp10 $fp_parse_RES $fp_parse_EXP10)
	fi
}

##multiply a by 10^b (cleverly) -- a b
fp_mulp10() {
	fp_mulp10_NUM=$1
	fp_mulp10_EXP=$2

	fp_mulp10_NUM_EXP=$(fp_get_exp $fp_mulp10_NUM)
	fp_mulp10_NUM_FRAC=$(fp_get_frac $fp_mulp10_NUM)
	fp_mulp10_EXP2=$(fp_mul $(fp_parse $fp_mulp10_EXP) $(fp_parse 3.321928094))

	(( fp_mulp10_NUM_EXP += $(fp_int $fp_mulp10_EXP2) ))
	
	fp_mulp10_EXP2=$(fp_sub $fp_mulp10_EXP2  $(fp_parse $(fp_int $fp_mulp10_EXP2)))
	fp_mulp10_TMP=$(fp_pow $(fp_parse 2) $fp_mulp10_EXP2)

	fp_mulp10_RES=$(fp_normalize_build $(fp_get_sign $fp_mulp10_NUM) $fp_mulp10_NUM_EXP $fp_mulp10_NUM_FRAC)
	fp_mulp10_RES=$(fp_mul $fp_mulp10_RES $fp_mulp10_TMP)

	echo $fp_mulp10_RES
}

## find the integer value of x, rounds towards zero -- x
## the answer is returned as a bash int, not a bash float
fp_int() {
	fp_int_N=$1
	fp_int_EXP=$(fp_get_exp $1)

	if (( fp_int_EXP < 127 )) 
	then
		echo 0
		return
	fi
	
	fp_int_FRAC=$(fp_get_frac $1)
	fp_int_SIGN=$(fp_get_sign $1)

	(( fp_int_FRAC <<= (fp_int_EXP-127) ))

	echo $(((fp_int_FRAC >> 23)*(fp_int_SIGN?-1:1)))
}

fp_log2() {
	fp_log2_N=$1
	fp_log2_L=0
	while (( $fp_log2_N > 0 ))
	do
		(( fp_log2_L++ ))
		(( fp_log2_N>>=1 ))
	done
	echo $fp_log2_L
}

#parse an int to a float
fp_parse_int() {
	fp_parse_int_N=$1
	fp_parse_int_SIGN=0
	if [[ $fp_parse_int_N == -* ]]
	then
		fp_parse_int_SIGN=1
		fp_parse_int_N=${fp_parse_int_N#-}
	fi
	fp_normalize_build $fp_parse_int_SIGN 150 $fp_parse_int_N
}

##Format the float into a string
fp_format() {
	case "$(fp_classify $1)" in
		1 ) echo "NaN"; return;;
		2 ) echo "-Infinity"; return;;
		3 ) echo "Infinity"; return;;
	esac
	local fp_format_DIGITS=$2
	[[ $fp_format_DIGITS == "" ]] && fp_format_DIGITS=8
	
	fp_format_SIGN=$(fp_get_sign $1)
	fp_format_EXP=$(fp_get_exp $1)
	fp_format_FRAC=$(fp_get_frac $1)
	fp_format_EXP10=$fp_CONST_ZERO
	
	if (( fp_format_EXP-127 > 50 )) || (( fp_format_EXP-127 < -50 ))
	then
		fp_format_LO=$(fp_mul $(fp_parse $((fp_format_EXP-127))) $(fp_parse 0.301029995664))
		fp_format_EXP10=$(fp_int $fp_format_LO)
		fp_format_EXP=127
		fp_format_LO=$(fp_sub $fp_format_LO $(fp_parse $fp_format_EXP10))

		fp_format_NUM=$(fp_build 0 $fp_format_EXP $fp_format_FRAC)
		fp_format_NUM=$(fp_mul $fp_format_NUM $(fp_pow $(fp_parse 10) $fp_format_LO))
	else
		fp_format_NUM=$(fp_set_sign $1 0)
	fi

	fp_format_ADJ="0."
	for (( fp_format_I=0; fp_format_I<fp_format_DIGITS; fp_format_I++))
	do
		fp_format_ADJ="${fp_format_ADJ}0"
	done
	fp_format_ADJ="${fp_format_ADJ}5"
	fp_format_NUM=$(fp_add $fp_format_NUM $(fp_set_sign $(fp_parse $fp_format_ADJ) $fp_format_SIGN))
	

	fp_format_STR=$(fp_int $fp_format_NUM)
	(( fp_format_DIGITS > 0 )) && fp_format_STR="$fp_format_STR."
	

	for((fp_format_I=0; fp_format_I < fp_format_DIGITS; fp_format_I++))
	do
		fp_format_NUM=$(fp_mul $(fp_sub $fp_format_NUM $(fp_parse $(fp_int $fp_format_NUM))) $(fp_parse 10))
		fp_format_STR="$fp_format_STR$(fp_int $fp_format_NUM)"
		(( fp_format_NUM == $fp_CONST_ZERO )) && break
	done

	while [[ ${fp_format_STR%0} != $fp_format_STR && $fp_format_STR != *.0 ]]
	do
		fp_format_STR=${fp_format_STR%0}
	done


	[[ $(fp_get_sign $1) == 1 ]] && fp_format_STR="-$fp_format_STR"
	(( fp_format_EXP10 != 0 )) && fp_format_STR="${fp_format_STR}E$fp_format_EXP10"
	echo $fp_format_STR

#	echo "scale=$fp_format_DIGITS; (-1)^$fp_format_SIGN * $fp_format_FRAC/2^23 * 2^($fp_format_EXP-127)" | bc
}

fp_inc() {
	fp_add $1 1065353216;
}
fp_dec() {
	fp_add $1 3212836864;
}

##Important constants
fp_CONST_NAN=$(fp_build 0 255 1)
fp_CONST_PINF=$(fp_build 0 255 0)
fp_CONST_NINF=$(fp_build 1 255 0)
fp_CONST_ZERO=$(fp_build 0 0 0)

##Find e^x -- x iter
fp_exp() {
	fp_exp_EX=$(fp_parse 1)
	fp_exp_X=${1:-$(fp_parse 1)}
	fp_exp_XN=$fp_exp_X
	fp_exp_NF=$(fp_parse 1)
	fp_exp_ITER=${2:-50}

	for ((fp_exp_I=2; fp_exp_I < fp_exp_ITER; fp_exp_I++ ))
	do
#		fp_dbg $fp_exp_EX $fp_exp_NF $fp_exp_XN
		fp_exp_EX=$(fp_add $(fp_div $fp_exp_XN $fp_exp_NF) $fp_exp_EX)
		fp_exp_NF=$(fp_mul $fp_exp_NF $(fp_parse $fp_exp_I))
		fp_exp_XN=$(fp_mul $fp_exp_XN $fp_exp_X)
		fp_is_positive_infinity $fp_exp_XN && break
		fp_is_positive_infinity $fp_exp_NF && break
	done

	echo $fp_exp_EX
}

##Find ln x -- x iter
fp_ln() {
	fp_ln_X=$1
	fp_ln_LNX=$fp_CONST_ZERO
	fp_ln_ITER=${2:-50}
	
	fp_ln_XN=$(fp_div $(fp_dec $fp_ln_X) $(fp_inc $fp_ln_X)) 
	fp_ln_XPNS=$(fp_mul $fp_ln_XN $fp_ln_XN)

	for((fp_ln_I=0; fp_ln_I < fp_ln_ITER; fp_ln_I++))
	do
#		fp_dbg $fp_ln_XN
		fp_ln_LNX=$(fp_add $fp_ln_LNX $(fp_div $fp_ln_XN $(fp_parse $((fp_ln_I*2+1))) ) )
		fp_ln_XN=$(fp_mul $fp_ln_XN $fp_ln_XPNS)
	done
	fp_mul fp_ln_LNX $(fp_parse 2)
}

##Find x^y -- x y
fp_pow() {
	fp_exp $(fp_mul $2 $(fp_ln $1))
}

###################################################################################
#Examples:


##calculate pi using Leibnitz' Formula -- iterations
pi_leibnitz() {
	PI=0
	K=1
	ITER=${1:-100}
	for((I=1; I<$ITER; I=I+2))
	do
		C=$(fp_div $(fp_parse 4) $(fp_parse $((I*K))))
		PI=$(fp_add $PI $C)
		(( K*=-1 ))
	done
	fp_format $PI
}

##calculate pi using Wallis' Product -- iterations
##This function will converge from below
pi_wallis() {
	PI=$(fp_parse 2)
	
	K=$(fp_parse 2)
	I=$(fp_parse 1)
	LIM=$(fp_parse ${1:-100})
	
	while fp_lt $I $LIM
	do
		PI=$(fp_mul $PI $(fp_div $K $I))
		I=$(fp_add $I $(fp_parse 2))
		PI=$(fp_mul $PI $(fp_div $K $I))
		K=$(fp_add $K $(fp_parse 2))
	done

	fp_format $PI
}

##do 1000 multiply operations
flops_mul() {
	for((I=0; I<100; I++))
	do
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
		fp_mul $1 $2
	done
}

echo "Demo app for koala_man's IEEE Single Precision Floating Point bash script"
echo
echo "Basic numbers:"
echo "Parsing  1.0: $(fp_parse 1.0) ($(fp_format $(fp_parse 1.0)))"
echo "Parsing -1337.666: $(fp_parse -1337.666) ($(fp_format $(fp_parse -1337.666)))"
echo
echo "Flops:"
echo -n "Doing 1000 multiply operations: "; 
FLOPS=$(TIMEFORMAT=%R; { time flops_mul $(fp_parse 42) $(fp_parse 42); } 2>&1 > /dev/null)
echo "1000 ops per $FLOPS secs, $(awk '{ print 1000.0/$0 }' <<< $FLOPS) flops"

echo 
echo "Scientific notation (based on exp/ln):"
echo -n "Parsing   1E6: "; A=$(fp_parse 1E6); echo "$A ($(fp_format $A))"
echo -n "Parsing   1E30: "; A=$(fp_parse 1E30); echo "$A ($(fp_format $A))"
echo -n "Parsing  -1E-30: "; A=$(fp_parse -1E-30); echo "$A ($(fp_format $A))"
echo
echo "Various silly calculations: "
echo -n "E="; fp_format $(fp_exp $(fp_parse 1))
echo -n "ln(2)="; fp_format $(fp_ln $(fp_parse 2))
echo -n "Leibnitz' formula with 50 iterations: pi="; pi_leibnitz 50
echo -n "Wallis' Product with 50 iterations: pi="; pi_wallis 50