Jump to content


Photo

Some Float Maths


13 replies to this topic

#1 Neil

Neil

    Regular

  • EstablishedMember
  • Pip
  • 26 posts

Posted 24 November 2008 - 09:43 PM

Hi all

I dont know if this is where I should do this, so sorry if Im in error.

I have seen some posts regarding float maths functions. Here is some code that I have writen to provide some functions. I dont clame it to be fast or even bug free, but it works for me. If you find it usfull all well and good.

Ive added an sprintf32 function. The code is derived from BoostC stdio.c.

e.g.
sprintf32f(myStr, "%8.2f", myFloat);

The format string is:- "%[optional '-'][optional ' ' or '+'][optional width][optional .precision]f"
were '-' gives left justification
were '+' always shows the sign (except if it is zero)
were ' ' pads a positive sign value with one space were the - or + would be
were width is the total width of the value, all chars including sign, decimal point and padding MAY TRUNCATE
were . (period) precision is the number of chars after the decimal point. a value of 0 supresses the decimal point

Also added atof

Added bug fixes in both sprintf32f and atof
Added yet another bug fix to sprintf32f


Neil

#ifndef _mcuMath_h
#define _mcuMath_h

#include <float.h>

unsigned char sprintf32f&#40; char* buffer, const char *format, float val &#41;;

float atof&#40;char* str&#41;;


float
mcu_sin&#40;float x&#41;;

float
mcu_cos&#40;float x&#41;;

float
mcu_tan&#40;float x&#41;;

float
mcu_asin&#40;float x&#41;;

float
mcu_acos&#40;float x&#41;;

float
mcu_atan&#40;float x&#41;;

float
mcu_epowx&#40;float x&#41;;

float
mcu_ln&#40;float x&#41;;

float
mcu_log10&#40;float x&#41;;

float
mcu_xpown&#40;float x, float n&#41;;


#endif // _mcuMath_h

The code

#include <float.h>
#include &#34;mcuMath.h&#34;

const float mcu_PI  = 3.14159265359;
const float mcu_hPI = 1.57079632679;
const float mcu_E   = 2.71828182846;


float
mcu_sin&#40;float x&#41;
{
 	float_rounding_mode = float_round_down;

	float ac = 0.0;
	float nv = 1.0;
	float xv = x;
	
	for&#40;unsigned char cnt = 2; cnt < 50; cnt+=2&#41;
	{
		ac = float32_add&#40;ac, float32_div&#40;xv, nv&#41;&#41;;
		nv = float32_mul&#40;nv, float32_mul&#40;float32_from_int32&#40;&#40;int32&#41;cnt&#41;, float32_from_int32&#40;&#40;int32&#41;&#40;cnt + 1&#41;&#41;&#41;&#41;;
		xv = float32_mul&#40;xv, float32_mul&#40;x, x&#41;&#41;;
		xv = float32_mul&#40;xv, -1.0&#41;;
	}
	
	return ac;
}


float
mcu_cos&#40;float x&#41;
{
 	float_rounding_mode = float_round_down;

	float ac = 0.0;
	float nv = 1.0;
	float xv = 1.0;
	
	for&#40;unsigned char cnt = 1; cnt < 50; cnt+=2&#41;
	{
		ac = float32_add&#40;ac, float32_div&#40;xv, nv&#41;&#41;;
		nv = float32_mul&#40;nv, float32_mul&#40;float32_from_int32&#40;&#40;int32&#41;cnt&#41;, float32_from_int32&#40;&#40;int32&#41;&#40;cnt + 1&#41;&#41;&#41;&#41;;
		xv = float32_mul&#40;xv, float32_mul&#40;x, x&#41;&#41;;
		xv = float32_mul&#40;xv, -1.0&#41;;
	}
	
	return ac;
}


float
mcu_tan&#40;float x&#41;
{
	return float32_div&#40;mcu_sin&#40;x&#41;, mcu_cos&#40;x&#41;&#41;;
}


float
mcu_asin&#40;float x&#41;
{
	return mcu_atan&#40;float32_div&#40;x,  float32_sqrt&#40;float32_sub&#40;1.0, float32_mul&#40;x, x&#41;&#41;&#41;&#41;&#41;;
}


float
mcu_acos&#40;float x&#41;
{
 	float_rounding_mode = float_round_down;

	float ac = 0;
	ac = mcu_atan&#40;float32_div&#40;float32_sqrt&#40;float32_sub&#40;1.0, float32_mul&#40;x, x&#41;&#41;&#41;,  x&#41;&#41;;
	if&#40;float32_lt&#40;x, 0&#41; != 0&#41; ac = float32_add&#40;mcu_PI, ac&#41;;
	return ac;
}


float
mcu_atan&#40;float x&#41;
{
 	float_rounding_mode = float_round_down;

	unsigned char inv = 0;
	if&#40;float32_gt&#40;x, 1.0&#41; != 0&#41;
	{
		x = float32_div&#40;1.0, x&#41;;
		inv = 1;
	}
	
	if&#40;float32_lt&#40;x, 1.0&#41; != 0&#41;
	{
		x = float32_div&#40;1.0, x&#41;;
		inv = 2;
	}

	float ac = 0.0;
	float nv = 1.0;
	float xv = x;
	
	for&#40;unsigned int cnt = 1; cnt < 500; cnt++&#41;
	{
		ac = float32_add&#40;ac, float32_div&#40;xv, nv&#41;&#41;;
		nv = float32_add&#40;nv, 2.0&#41;;
		xv = float32_mul&#40;xv, float32_mul&#40;x, x&#41;&#41;;
		xv = float32_mul&#40;xv, -1.0&#41;;
	}
	
	if&#40;inv == 1&#41;
	{
		ac = float32_sub&#40;mcu_hPI, ac&#41;;
	}
	if&#40;inv == 2&#41;
	{
		ac = float32_add&#40;mcu_hPI, ac&#41;;
		ac = float32_mul&#40;ac, -1.0&#41;;
	}
	
	return ac;
}


float
mcu_epowx&#40;float x&#41;
{
	float_rounding_mode = float_round_down;

	int val = &#40;int&#41;float32_to_int32&#40;x&#41;;
	x = float32_sub&#40;x, float32_from_int32&#40;&#40;int32&#41;val&#41;&#41;;
	
	float ac = 1.0;
	float nv = 1.0;
	float xv = x;
	
	for&#40;unsigned char cnt = 1; cnt < 8; ++cnt&#41;
	{
		nv = float32_mul&#40;nv, float32_from_int32&#40;&#40;int32&#41;cnt&#41;&#41;;
		ac = float32_add&#40;ac, float32_div&#40;xv, nv&#41;&#41;;
		xv = float32_mul&#40;xv, x&#41;;
	}
	
	float esum = 1.0;
	unsigned int cntL;
	if&#40;val < 0&#41; cntL = &#40;unsigned int&#41;&#40;val * -1&#41;;
	else		cntL = &#40;unsigned int&#41; val;
	
	for&#40;unsigned int cnt2 = 0; cnt2 < cntL; ++cnt2&#41;
	{
		if&#40;val > 0&#41; esum = float32_mul&#40;esum, mcu_E&#41;;
		else		esum = float32_div&#40;esum, mcu_E&#41;;
	}
		
	return float32_mul&#40;ac, esum&#41;;
}


float
mcu_ln&#40;float x&#41;
{
	float_rounding_mode = float_round_down;
	
	int cnt = 0;
	while&#40;float32_gt&#40;x, 1.0&#41; != 0&#41;
	{
		x = float32_div&#40;x, 2.71828182846&#41;;
		cnt++;
	}
	
	x = float32_sub&#40;x, 1.0&#41;;
	
	float ac = 0;
	float xv = x;
	
	for&#40;unsigned char cnt2 = 1; cnt2 < 100; ++cnt2&#41;
	{
		ac = float32_add&#40;ac, float32_div&#40;xv, float32_from_int32&#40;&#40;int32&#41;cnt2&#41;&#41;&#41;;
		xv = float32_mul&#40;xv, x&#41;;
		xv = float32_mul&#40;xv, -1.0&#41;;
	 }
	
	ac = float32_add&#40;ac, float32_from_int32&#40;&#40;int32&#41;cnt&#41;&#41;;
	
	return ac;
}


float
mcu_log10&#40;float x&#41;
{
	return float32_div&#40;mcu_ln&#40;x&#41;, 2.30258509299&#41;;   // = ln&#40;10.0&#41;;
}


float
mcu_xpown&#40;float x, float n&#41;
{
 	float_rounding_mode = float_round_down;

	bool inv = false;
	if&#40;float32_lt&#40;n, 0.0&#41; != 0&#41;
	{
		n = float32_mul&#40;n, -1.0&#41;;
		inv = true;
	}

	int val = &#40;int&#41;float32_to_int32&#40;n&#41;;
	n = float32_sub&#40;n, float32_from_int32&#40;&#40;int32&#41;val&#41;&#41;;
	
	x = float32_sub&#40;1.0, x&#41;;

	float ac = 0.0;
	float xv = 1.0;
	float nv = float32_mul&#40;-1.0, n&#41;;
	float cv = 1.0;
	
	ac = float32_add&#40;ac, xv&#41;;
	for&#40;unsigned char cnt = 1; cnt < 25; ++cnt&#41;
	{
		cv = float32_mul&#40;cv, float32_from_int32&#40;&#40;int32&#41;cnt&#41;&#41;;
		xv = float32_mul&#40;xv, x&#41;;
		ac = float32_add&#40;ac, float32_div&#40;float32_mul&#40;nv, xv&#41;, cv&#41;&#41;;
		nv = float32_mul&#40;nv, float32_sub&#40;n, float32_from_int32&#40;&#40;int32&#41;cnt&#41;&#41;&#41;;
		nv = float32_mul&#40;nv, -1.0&#41;;
	}
	x = float32_sub&#40;1.0, x&#41;;
	
	float xsum = 1.0;
	unsigned int cntL;
	if&#40;val < 0&#41; cntL = &#40;unsigned int&#41;&#40;val * -1&#41;;
	else		cntL = &#40;unsigned int&#41; val;
	
	for&#40;unsigned int cnt2 = 1; cnt2 < cntL; ++cnt2&#41;
	{
		if&#40;val > 0&#41; xsum = float32_mul&#40;xsum, mcu_E&#41;;
		else		xsum = float32_div&#40;xsum, mcu_E&#41;;
	}
	ac = float32_mul&#40;ac, xsum&#41;;
	
	if&#40;inv == true&#41; ac = float32_div&#40;1.0, ac&#41;;
 
	return ac;   
}

///////////////////////////////////////////////////////
// stdio.c
///////////////////////////////////////////////////////
// This file is a part of the system libraries
// supplied with BoostC C compiler distribution. 
// More information can be found on&#58;
// http&#58;//www.sourceboost.com
//
// You are allowed to use or modify this code as
// long as the following conditions are met&#58;
// 1. this entire text from this comment block 
//	remains intact;
// 2. this code and derivative work are compiled and
//	linked with products from SourceBoost Technologies;
// 3. the source code for the derivative work includes 
//	prominent notice that the work is derivative.
//
// THIS SOFTWARE IS PROVIDED IN AN &#34;AS IS&#34; CONDITION. 
// NO WARRANTIES, WHETHER EXPRESS, IMPLIED OR STATUTORY, 
// INCLUDING, BUT NOT LIMITED TO, IMPLIED WARRANTIES OF
// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
// APPLY TO THIS SOFTWARE. SOURCEBOOST TECHNOLOGIES SHALL 
// NOT, IN ANY CIRCUMSTANCES, BE LIABLE FOR SPECIAL, 
// INCIDENTAL OR CONSEQUENTIAL DAMAGES, FOR ANY REASON 
// WHATSOEVER.
//
// Copyright&#40;C&#41; 2004-2008 Pavel Baranov
// Copyright&#40;C&#41; 2004-2008 David Hobday
//
/////////////////////////////////////////////////////////////
//
// Author&#58; David Hobday
// Date&#58; 11/11/2006
// V1.00 11/11/2006
// Initial version.
//
// V1.01 18/02/2007 
// Added sprintf32 function for 32bit values.
//
//
// Auther&#58; Neil Robinson 25-11-2008
// Added sprintf32f function for floats
//

unsigned char
sprintf32f&#40;char* buffer, const char *format, float val&#41;
{
	float_rounding_mode = float_round_down;

	unsigned char bi = 0, fi = 0;
	char c;
	
	while&#40; 1 &#41;
	{
		c = format&#91;fi++&#93;; // read next string format character
		
		if &#40; !c &#41; goto end; // end of format string reached
		
		if&#40; 	c != &#39;%&#39;  // not format specifier
			|| &#40;c = format&#91;fi++&#93;&#41; == &#39;%&#39; &#41; // %% sequence		
		{
			buffer&#91; bi++ &#93; = c; // out character raw
		}
		else
		{			
			unsigned char fillChar, fmatExp, signChar = 0;
			bit alignLeft = false;
			
			fillChar = &#39; &#39;; // default fill character
			
			// process flags
			while&#40; 1 &#41;
			{				
				switch&#40; c &#41;
				{
					case 0&#58; goto end; // end of format string reached
					case &#39;0&#39;&#58; fillChar = &#39;0&#39;; break;
					case &#39; &#39;&#58; signChar = &#39; &#39;; break; // show space in front of signed numbers when not negative						
					case &#39;+&#39;&#58; signChar = &#39;+&#39;; break; // show sign on signed numbers						
					case &#39;-&#39;&#58; alignLeft = true; break;
					default&#58; goto processWidth;
				}
				
				c = format&#91;fi++&#93;; // read next string format character
			}

processWidth&#58;
			// Next character if valid digit indicates field width
			unsigned char fieldWidth = 0;
			while&#40; c >= &#39;0&#39; && c <= &#39;9&#39; &#41;
			{			
				fieldWidth *= 10;
				fieldWidth += c;
				fieldWidth -= &#39;0&#39;;				
				c = format&#91;fi++&#93;; // read next string format character
				if &#40; !c &#41; goto end; // end of format string reached
			}	
			
			if&#40;c == &#39;.&#39;&#41;
			{
				c = format&#91;fi++&#93;; // read next string format character
				if &#40; !c &#41; goto end; // end of format string reached
			}
			
			unsigned char precision = 0;
			signed char preTmp = -1;
			while&#40; c >= &#39;0&#39; && c <= &#39;9&#39; &#41;
			{
				if&#40;preTmp == -1&#41; preTmp = 0;
				preTmp *= 10;
				preTmp += c;
				preTmp -= &#39;0&#39;;				
				c = format&#91;fi++&#93;; // read next string format character
				if &#40; !c &#41; goto end; // end of format string reached
			}				
			
			if&#40;preTmp == -1&#41; precision = 6;
			else			 precision = &#40;unsigned char&#41;preTmp;
			
			unsigned char floatSignChar = 0;
			
			// Next character indicates the format
			switch&#40; c &#41;
			{
				case &#39;f&#39;&#58;
				{
					fmatExp = 0;		
					if&#40;float32_lt&#40;val, 0&#41; != 0&#41;		// Negative values must be adjusted to be positive
					{
						floatSignChar = &#39;-&#39;;
						val = float32_mul&#40;val, -1.0&#41;;
					}
				}
				break;
					// fall through to next case intended
				default&#58; goto end; // no radix specified - can&#39;t continue
			}
			
			long valFr = 0;
			long valIn = 0;
			
			long frMag = 1;
			
			for&#40;unsigned char ct = 0; ct < precision; ++ct&#41;
			{
				valFr *= 10;
				frMag *= 10;
			}
			
			valIn = &#40;long&#41;float32_to_int32&#40;val&#41;;
			valFr = &#40;long&#41;float32_to_int32&#40;float32_add&#40;float32_mul&#40;float32_sub&#40;val, float32_from_int32&#40;&#40;int32&#41;valIn&#41;&#41;, float32_from_int32&#40;frMag&#41;&#41;, 0.5&#41;&#41;;
			frMag--;
			
			if&#40;&#40;valIn == 0&#41; && &#40;valFr == 0&#41;&#41; floatSignChar = 0;
			
			if&#40;floatSignChar == &#39;-&#39;&#41; signChar = floatSignChar;
			
			unsigned char digitWidth = 0;
			{	// determine number of digits in the number
				unsigned long valx = &#40;unsigned long&#41;float32_to_int32&#40;val&#41;;
				do
				{
					digitWidth++;
					valx /= 10;
				}
				while&#40; valx != 0 &#41;;
			}
				
			// sign takes extra space
			if&#40; signChar != 0 &#41;
				digitWidth++;
				
			// precision takes extra space
			digitWidth += precision;
			
			// one for the decimal point
			if&#40;precision > 0&#41; digitWidth++;
			
			
			// limit digit width to prevent buffer overflow - not standard, but will prevent buffer overrun corruption
			unsigned char fillWidth = 0;
			if&#40; fieldWidth > 0 &#41;
			{
				if&#40; digitWidth > fieldWidth &#41; // limit number of digits
					digitWidth = fieldWidth;
					
				fillWidth = fieldWidth - digitWidth;
			}
			else
			{
				fieldWidth = digitWidth;
				fillWidth = 0;
			}
			
			if&#40; fillWidth == 0 &#41;
				alignLeft = false; // alignment serves no purpose if no fill
				
			// we reverse fill string, so start at the end
			bi += fieldWidth;
			
			bit doFill = alignLeft; // if left aligning do fill first.
			unsigned char fracFlg = precision;
			
			while&#40; digitWidth != 0 || fillWidth != 0  &#41;
			{
				bi--;
				if&#40; doFill &#41;
				{
					c = fillChar;
					if&#40; --fillWidth == 0 &#41;
					 doFill = false;
				}
				else
				{
					if&#40; --digitWidth == 0 && signChar &#41; // last digit and sign character
						c = signChar; // add sign character
					else if &#40;frMag > 0&#41;
					{
						c = &#40;char&#41;&#40;valFr % 10 &#41;;											
						c += &#39;0&#39;; // convert to digit character 0-9
						valFr /= 10&#59;
						frMag /= 10&#59;
						fracFlg = 0;
					}
					else if &#40;&#40;valFr == 0&#41; && &#40;precision > 0&#41;&#41;
					{
						if&#40;fracFlg > 0&#41;
						{
							c = &#39;0&#39;;											
							fracFlg--;
						}
						else
						{
							c = &#39;.&#39;;											
							valFr = -1;
						}
					}
					else
					{
						c = &#40;char&#41;&#40;valIn % 10 &#41;;											
						c += &#39;0&#39;; // convert to digit character 0-9
						valIn /= 10&#59;					
					}
					
					if&#40; digitWidth == 0 &#41;
						doFill = true;
				}
								
				buffer&#91; bi &#93; = c;
			}
			bi += fieldWidth;
		}	
	}
end&#58;	
	buffer&#91; bi &#93; = 0; // null terminate
	return bi;
}



float
atof&#40;char* str&#41;
{
	if&#40;&#40;str == 0&#41; || &#40;strlen&#40;str&#41; == 0&#41;&#41; return 0;
	
	int hl = 0;
	int fr = 0;
	float retVal = 0;
	
	char* pntPtr = strrchr&#40;str,&#39;.&#39;&#41;;
	for&#40;char ct = strlen&#40;str&#41; - 1; ct >=0; --ct&#41;
	{
		if&#40;str&#91;ct&#93; != &#39;0&#39;&#41; break;
		
		str&#91;ct&#93; = 0;
	}
	
	
	if&#40;pntPtr != 0&#41;
	{
		fr = atoi&#40;pntPtr + 1&#41;;
		*pntPtr = 0;
	}
	
	hl = atoi&#40;str&#41;;
	
	char foo&#91;32&#93;;
	if&#40;pntPtr != 0&#41;
	{
		*pntPtr = &#39;.&#39;;
		
		float pwr = 1.0;
		for&#40;unsigned char ct = 0; ct < strlen&#40;pntPtr + 1&#41;; ++ct&#41;
		{
			pwr = float32_mul&#40;pwr, 10.0&#41;;
		}
		
		retVal = float32_div&#40;float32_from_int32&#40;&#40;int32&#41; fr&#41;, pwr&#41;;
	}
	

	if&#40;hl < 0&#41;
	{
		retVal = float32_sub&#40;float32_from_int32&#40;&#40;int32&#41; hl&#41;, retVal&#41;;
	}
	else
	{
		retVal = float32_add&#40;float32_from_int32&#40;&#40;int32&#41; hl&#41;, retVal&#41;;
	}

	return retVal;
}

Edited by Neil, 27 November 2008 - 09:26 PM.


#2 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 27 November 2008 - 06:13 PM

Hi Neil,

Thanks for the float math and prinf functions.

I thought I would put it to a small test with the sin function but ran into a small problem.

nv = float32_mul&#40;nv, float32_mul&#40;float32_from_int32&#40;&#40;int32&#41;cnt&#41;, float32_from_int32&#40;&#40;int32&#41;&#40;cnt + 1&#41;&#41;&#41;&#41;;

This line does not work correctly. I think it may be due to recursion with the float32_mul function. The first time through the loop, nv has the value 12.00000 instead of 6.00000.

Splitting these two multiplications up into separate calls seems to work OK.

Cheers

Reynard

#3 Neil

Neil

    Regular

  • EstablishedMember
  • Pip
  • 26 posts

Posted 27 November 2008 - 09:30 PM

Hi Reynard

Thanks for the bug fix. I have not tested the functions to distruction, so I expect some bugs. As I use a function I test it more!!

I first wrote these under gcc/linux and compared the resulst over many values.

Porting to PIC float lib was not as simple as I expected, mainly due to rounding etc.

If there is a better place for this code Im happy to donate it.

Also, I just put a bug fix into sprintf32f. Id would print 1.002 as 001.2 !!! Fixed now!!!!

Neil

#4 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 27 November 2008 - 09:39 PM

Hi Neil,

Another little tweak would be to take the x squared outside of the for loop and you gain 12 percent in speed.

x = float32_mul&#40;x, x&#41;;
for&#40;unsigned char cnt = 2; cnt < 50; cnt+=2&#41;

Cheers

Reynard

#5 Dave

Dave

    Super Maniac

  • Administrators
  • PipPipPipPipPip
  • 2,091 posts
  • Gender:Male
  • Location:UK
  • Interests:How things work, Electronics, Software, Cycling.

Posted 28 November 2008 - 11:16 AM

Neil,

I first wrote these under gcc/linux and compared the resulst over many values.

Porting to PIC float lib was not as simple as I expected, mainly due to rounding etc.

If there is a better place for this code Im happy to donate it.

Also, I just put a bug fix into sprintf32f. Id would print 1.002 as 001.2 !!! Fixed now!!!!

This code looks like implementation of some missing features. We would like to consider them to be added to the Boostc compiler libraries included in the SourceBoost installation. Once your are happy that these routines provide accurate results and are fully tested, please let supply the source code to support@sourceboost.com for further evaluation.

Keep up the good work.

Regards
Dave

#6 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 28 November 2008 - 02:49 PM

Nice to see that Dave is interested and I hope others are too.

I have made a few more tweaks to the sine function which increased the speed by 42 percent from the original working source.
float
mcu_sin&#40;float x&#41;
{
	float_rounding_mode = float_round_down;

	float ac = 0.0;
	float nv = 1.0;
	float xv = x;

	x = float32_mul&#40;x, x&#41;;
	for&#40;unsigned char cnt = 2; cnt < 50; cnt+=2&#41;
	{
		ac = float32_add&#40;ac, float32_div&#40;xv, nv&#41;&#41;;
		nv = float32_mul&#40;nv, float32_from_int32&#40;&#40;int32&#41;&#40;&#40;unsigned int&#41;&#40;cnt * &#40;cnt + 1&#41;&#41;&#41;&#41;&#41;;
		xv = float32_mul&#40;xv, x&#41;;
		xv = xv ^ 0x80000000l;
	}
	
	return ac;
}

The last line of the loop is a bit dirty for the multiplication by -1.0 but it works (I hope). It just diddles with the sign bit.

Accuracy is OK to 5 places which is probably all you can expect with single precision floats I suppose.
My speed tests where with a PIC 18Fxxx using the simulator tick counter.
For those who don't require the accuracy can reduce the iteration count and go faster.

Cheers

Reynard

#7 Neil

Neil

    Regular

  • EstablishedMember
  • Pip
  • 26 posts

Posted 28 November 2008 - 07:07 PM

Hi Dave and Reynard


Thanks for the interest. Im more than happy for the code to go into BoostC.

From the tests I have doen, all the paths functions work correctly. The could be faster, but Reynard seems to be on the case.

I feel they really need some testing by others to shake out any other bugs.

Feel free to include my code into BoostC.

I dont know when I will next get time to work on this code. It just depends on what project Im working on.

Best regards

Neil

PS I do think sprintf32 is now working!!

#8 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 28 November 2008 - 09:04 PM

Hi Neil,

Done the 3 basic trig functions (sin,cos,tan) and they seem to be working OK.

Worked in the atan function but with 500 iterations it takes a long time and accuracy only good to 3 places. Nothing wrong with the formula, Taylor series or something, but the build up of floating point errors using short floats.

I found another bit of code that I have converted over to BoostC. It is very fast by a factor of 10 over the original code I have been working with and alot more accurate. The accuracy probably due to less arithmetic hence less floating point errors to build up.

It uses a bit more ram, but for a big PIC 18 it is worth it for the speed gained.

#define PI_DIV_BY_TWO	1.570796326794896

float
cos&#40;float x&#41;
{
	float p&#91;4&#93; = {
		-0.499999993585,
		 0.041666636258,
		-0.0013888361399,
		 0.00002476016134
	};
	float frac, t, y, t2 = 1.0;
	unsigned char quad, i;

	if &#40;float32_lt&#40;x,  0&#41;&#41; {
		x = x ^ 0x80000000l;			 // absolute value of input
	}

	quad = &#40;unsigned char&#41;float32_to_int32&#40;float32_div&#40;x , PI_DIV_BY_TWO&#41;&#41;;
	frac = float32_sub&#40;float32_div&#40;x , PI_DIV_BY_TWO&#41;, float32_from_int32&#40;&#40;int32&#41;quad&#41;&#41;;
	quad = quad % 4;					// quadrant &#40;0 to 3&#41;

	if &#40;quad == 0 || quad == 2&#41;
		t = float32_mul&#40;frac, PI_DIV_BY_TWO&#41;;
	else if &#40;quad == 1&#41;
		t = float32_mul&#40;float32_sub&#40;1.0, frac&#41;, PI_DIV_BY_TWO&#41;;
	else // should be 3
		t = float32_mul&#40;float32_sub&#40;frac, 1.0&#41;, PI_DIV_BY_TWO&#41;;
	y = 0.999999999781;
	t = float32_mul&#40;t, t&#41;;

	for &#40;i = 0; i <= 3; i++&#41;
	{
		t2 = float32_mul&#40;t2, t&#41;;
		y = float32_add&#40;y, float32_mul&#40;p&#91;i&#93;, t2&#41;&#41;;
	}

	if &#40;quad == 2 || quad == 1&#41;
	y = y ^ 0x80000000l;  // correct sign

	return y;
}

float
sin&#40;float x&#41;
{
	return cos&#40;float32_sub&#40;x, PI_DIV_BY_TWO&#41;&#41;;
}

Cheers

Reynard

#9 PeterK

PeterK

    Newbrie

  • Members
  • 2 posts

Posted 24 April 2010 - 11:52 PM

I know this is an old topic, however I would like to comment on the implementation of the math functions as shown in this topic. The method used is not the proper way to do it. You can write for every math function a routine which uses at average 7 multiplications and 7 additions to calculate the value with 6 digits accuracy! Its very simple to create a polynome which has great accuracy for a small range of input values. With simple math rules you can translate all possible input values to this small range. All math libraries are written in this manner !!!

I stumbled on this topic because I needed a natural log function. I was not happy with the supplied code, therefore I have written my own "log" function which uses only 7 multiplications and 7 additions with 6 digits accuracy for ALL possibly input values. If anyone is interested I will post my "log" version.

The same method can be used for all math functions, however you need some math knowledge to accomplish this.

#10 Dave

Dave

    Super Maniac

  • Administrators
  • PipPipPipPipPip
  • 2,091 posts
  • Gender:Male
  • Location:UK
  • Interests:How things work, Electronics, Software, Cycling.

Posted 25 April 2010 - 11:36 AM

PeterK,

I know this is an old topic, however I would like to comment on the implementation of the math functions as shown in this topic. The method used is not the proper way to do it. You can write for every math function a routine which uses at average 7 multiplications and 7 additions to calculate the value with 6 digits accuracy! Its very simple to create a polynome which has great accuracy for a small range of input values. With simple math rules you can translate all possible input values to this small range. All math libraries are written in this manner !!!

I stumbled on this topic because I needed a natural log function. I was not happy with the supplied code, therefore I have written my own "log" function which uses only 7 multiplications and 7 additions with 6 digits accuracy for ALL possibly input values. If anyone is interested I will post my "log" version.

The same method can be used for all math functions, however you need some math knowledge to accomplish this.

We are very interested in this, we want to improve on the supplied maths function.
Please post your log funtion.

Regards
Dave

#11 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 25 April 2010 - 05:16 PM

Most compilers will use a version of Henri Pade approximant.

As Peter says it can be done with 7 muls and adds or there abouts to produce an accurate value.

Finding the correct coefficients without stealing someone elses is the secret.

Once BoostC supports floats as a basic data type it should be vary fast for all trig functions.

Cheers

Reynard

#12 Reynard

Reynard

    Maniac

  • EstablishedMember
  • PipPipPipPip
  • 659 posts
  • Gender:Male
  • Location:Scotland
  • Interests:Archery - target and field

Posted 25 April 2010 - 05:21 PM

Here is a little something I was playing around with a year or so ago.

It may want some refining.

/************************************************************
*   Return the exp of the fp number x.					  *
************************************************************/
#define LN2 0.69314718056		// ln&#40;2&#41;

float exp&#40;float x&#41;
{
	float y, r, twopwrn;
	int n;
	unsigned char sign;

	float_rounding_mode = float_round_down;

	// check if result would exceed max fp number &#40;3.4028235e38 ish&#41;.
	if&#40;float32_gt&#40;x, 88.722839&#41;&#41;
	{
		float_raise &#40;float_flag_overflow&#41;
		return&#40;0.0&#41;;
	}

	// calc. raise power value &#40;2^n&#41;
	n = &#40;int&#41;float32_to_int32&#40;float32_div&#40;x, LN2&#41;&#41;;
	sign = 0;					// assume x is positive.
	y = x;						// take local copy of x.

	if &#40;float32_lt&#40;x, 0.0&#41;&#41;		// x negative ?
	{
		sign = 1;				// yes, set sign flag.
		n = -n;					// make power absolute value.
		asm {
			bcf _y+3,7			// make y absolute value.
		}
	}

	// Calc. 2^n.
	twopwrn = 0.0;					// set res to all zeroes.
	// play some tricks to compute the exponent value &#40;2^n&#41;.
	*&#40;&#40;long*&#41;&twopwrn&#41; = &#40;long&#41;&#40;&#40;n + 127&#41; << 23&#41;;

	// y = &#40;y/ln&#40;2&#41;&#41; - n
	y = float32_sub&#40;float32_div&#40;y, LN2&#41;, float32_from_int32&#40;n&#41;&#41;;

	// compute a series.
	// y = P0y^6 + P1y^5 + P2y^4 + P3y^3 + P4y2 + P5y
	r = float32_add&#40;float32_mul&#40;0.000207455774, y&#41;, 0.00127100575&#41;;
	r = float32_add&#40;float32_mul&#40;r, y&#41;, 0.00965065093&#41;;
	r = float32_add&#40;float32_mul&#40;r, y&#41;, 0.0554965651&#41;;
	r = float32_add&#40;float32_mul&#40;r, y&#41;, 0.240227138&#41;;
	r = float32_add&#40;float32_mul&#40;r, y&#41;, 0.693147172&#41;;
	y = float32_mul&#40;y, r&#41;;

	// Calc. res = 2^n&#40;1 + y&#41;
	y = float32_mul&#40;twopwrn, float32_add&#40;1.0, y&#41;&#41;;

	if &#40;sign&#41;		// was x negative ?
		y = float32_div&#40;1.0, y&#41;;	// yes, invert. y = 1/y

	return&#40;y&#41;;
}

Cheers

Reynard

#13 PeterK

PeterK

    Newbrie

  • Members
  • 2 posts

Posted 25 April 2010 - 06:29 PM

You are right about "Pade approximant", which is even better then "Chebyshev Approximation". I use the "Polynomial Approximation from Chebyshev Coefficients" which are easy to determine. If you google at "numerical recipes in c pdf", you will find the PDF version of the book I use. In this document you find full documented C-functions of "Chebyshev Approximation", "Polynomial Approximation from Chebyshev Coefficients", and of course "Pade approximant". You can copy and paste these functions directly from the PDF file, and you are in business in calculating your own polynomial coefficients. After calculating your coefficients, you will have a polynomial which is accurate for a small range of input values. The trick is to convert all input values to this range, for most functions this is easy to accomplish with simple math.

#14 Neil

Neil

    Regular

  • EstablishedMember
  • Pip
  • 26 posts

Posted 03 June 2010 - 09:16 AM

You are right about "Pade approximant", which is even better then "Chebyshev Approximation". I use the "Polynomial Approximation from Chebyshev Coefficients" which are easy to determine. If you google at "numerical recipes in c pdf", you will find the PDF version of the book I use. In this document you find full documented C-functions of "Chebyshev Approximation", "Polynomial Approximation from Chebyshev Coefficients", and of course "Pade approximant". You can copy and paste these functions directly from the PDF file, and you are in business in calculating your own polynomial coefficients. After calculating your coefficients, you will have a polynomial which is accurate for a small range of input values. The trick is to convert all input values to this range, for most functions this is easy to accomplish with simple math.


Hi PeterK

Its been a while since I have be on the boards, so sorry for the late reply.

While developing my own code, I had a need for float math. I soon found, as may had that BoostC++ does not provide any. So, as a service, I freely gave away my source code. No license, do with as you see fit. Kindly, some have helped to improve it. I am no compiler writer, so excude me for a poor implementaton. I will look into Henry Pade, may be I can write a better lib.

While I apreciate your comments, and your clear knowledge on this point, I do not see any source code to help the rest of the community. May be you could share your code as I have done.

regards

Neil



Reply to this topic



  


0 user(s) are reading this topic

0 members, 0 guests, 0 anonymous users