Jump to content
Neil

Some Float Maths

Recommended Posts

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( char* buffer, const char *format, float val );

float atof(char* str);


float
mcu_sin(float x);

float
mcu_cos(float x);

float
mcu_tan(float x);

float
mcu_asin(float x);

float
mcu_acos(float x);

float
mcu_atan(float x);

float
mcu_epowx(float x);

float
mcu_ln(float x);

float
mcu_log10(float x);

float
mcu_xpown(float x, float n);


#endif // _mcuMath_h

 

The code

 

#include <float.h>
#include "mcuMath.h"

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


float
mcu_sin(float x)
{
	float_rounding_mode = float_round_down;

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

for(unsigned char cnt = 2; cnt < 50; cnt+=2)
{
	ac = float32_add(ac, float32_div(xv, nv));
	nv = float32_mul(nv, float32_mul(float32_from_int32((int32)cnt), float32_from_int32((int32)(cnt + 1))));
	xv = float32_mul(xv, float32_mul(x, x));
	xv = float32_mul(xv, -1.0);
}

return ac;
}


float
mcu_cos(float x)
{
	float_rounding_mode = float_round_down;

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

for(unsigned char cnt = 1; cnt < 50; cnt+=2)
{
	ac = float32_add(ac, float32_div(xv, nv));
	nv = float32_mul(nv, float32_mul(float32_from_int32((int32)cnt), float32_from_int32((int32)(cnt + 1))));
	xv = float32_mul(xv, float32_mul(x, x));
	xv = float32_mul(xv, -1.0);
}

return ac;
}


float
mcu_tan(float x)
{
return float32_div(mcu_sin(x), mcu_cos(x));
}


float
mcu_asin(float x)
{
return mcu_atan(float32_div(x,  float32_sqrt(float32_sub(1.0, float32_mul(x, x)))));
}


float
mcu_acos(float x)
{
	float_rounding_mode = float_round_down;

float ac = 0;
ac = mcu_atan(float32_div(float32_sqrt(float32_sub(1.0, float32_mul(x, x))),  x));
if(float32_lt(x, 0) != 0) ac = float32_add(mcu_PI, ac);
return ac;
}


float
mcu_atan(float x)
{
	float_rounding_mode = float_round_down;

unsigned char inv = 0;
if(float32_gt(x, 1.0) != 0)
{
	x = float32_div(1.0, x);
	inv = 1;
}

if(float32_lt(x, 1.0) != 0)
{
	x = float32_div(1.0, x);
	inv = 2;
}

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

for(unsigned int cnt = 1; cnt < 500; cnt++)
{
	ac = float32_add(ac, float32_div(xv, nv));
	nv = float32_add(nv, 2.0);
	xv = float32_mul(xv, float32_mul(x, x));
	xv = float32_mul(xv, -1.0);
}

if(inv == 1)
{
	ac = float32_sub(mcu_hPI, ac);
}
if(inv == 2)
{
	ac = float32_add(mcu_hPI, ac);
	ac = float32_mul(ac, -1.0);
}

return ac;
}


float
mcu_epowx(float x)
{
float_rounding_mode = float_round_down;

int val = (int)float32_to_int32(x);
x = float32_sub(x, float32_from_int32((int32)val));

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

for(unsigned char cnt = 1; cnt < 8; ++cnt)
{
	nv = float32_mul(nv, float32_from_int32((int32)cnt));
	ac = float32_add(ac, float32_div(xv, nv));
	xv = float32_mul(xv, x);
}

float esum = 1.0;
unsigned int cntL;
if(val < 0) cntL = (unsigned int)(val * -1);
else		cntL = (unsigned int) val;

for(unsigned int cnt2 = 0; cnt2 < cntL; ++cnt2)
{
	if(val > 0) esum = float32_mul(esum, mcu_E);
	else		esum = float32_div(esum, mcu_E);
}

return float32_mul(ac, esum);
}


float
mcu_ln(float x)
{
float_rounding_mode = float_round_down;

int cnt = 0;
while(float32_gt(x, 1.0) != 0)
{
	x = float32_div(x, 2.71828182846);
	cnt++;
}

x = float32_sub(x, 1.0);

float ac = 0;
float xv = x;

for(unsigned char cnt2 = 1; cnt2 < 100; ++cnt2)
{
	ac = float32_add(ac, float32_div(xv, float32_from_int32((int32)cnt2)));
	xv = float32_mul(xv, x);
	xv = float32_mul(xv, -1.0);
 }

ac = float32_add(ac, float32_from_int32((int32)cnt));

return ac;
}


float
mcu_log10(float x)
{
return float32_div(mcu_ln(x), 2.30258509299);   // = ln(10.0);
}


float
mcu_xpown(float x, float n)
{
	float_rounding_mode = float_round_down;

bool inv = false;
if(float32_lt(n, 0.0) != 0)
{
	n = float32_mul(n, -1.0);
	inv = true;
}

int val = (int)float32_to_int32(n);
n = float32_sub(n, float32_from_int32((int32)val));

x = float32_sub(1.0, x);

float ac = 0.0;
float xv = 1.0;
float nv = float32_mul(-1.0, n);
float cv = 1.0;

ac = float32_add(ac, xv);
for(unsigned char cnt = 1; cnt < 25; ++cnt)
{
	cv = float32_mul(cv, float32_from_int32((int32)cnt));
	xv = float32_mul(xv, x);
	ac = float32_add(ac, float32_div(float32_mul(nv, xv), cv));
	nv = float32_mul(nv, float32_sub(n, float32_from_int32((int32)cnt)));
	nv = float32_mul(nv, -1.0);
}
x = float32_sub(1.0, x);

float xsum = 1.0;
unsigned int cntL;
if(val < 0) cntL = (unsigned int)(val * -1);
else		cntL = (unsigned int) val;

for(unsigned int cnt2 = 1; cnt2 < cntL; ++cnt2)
{
	if(val > 0) xsum = float32_mul(xsum, mcu_E);
	else		xsum = float32_div(xsum, mcu_E);
}
ac = float32_mul(ac, xsum);

if(inv == true) ac = float32_div(1.0, ac);

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:
// http://www.sourceboost.com
//
// You are allowed to use or modify this code as
// long as the following conditions are met:
// 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 "AS IS" 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(C) 2004-2008 Pavel Baranov
// Copyright(C) 2004-2008 David Hobday
//
/////////////////////////////////////////////////////////////
//
// Author: David Hobday
// Date: 11/11/2006
// V1.00 11/11/2006
// Initial version.
//
// V1.01 18/02/2007 
// Added sprintf32 function for 32bit values.
//
//
// Auther: Neil Robinson 25-11-2008
// Added sprintf32f function for floats
//

unsigned char
sprintf32f(char* buffer, const char *format, float val)
{
float_rounding_mode = float_round_down;

unsigned char bi = 0, fi = 0;
char c;

while( 1 )
{
	c = format[fi++]; // read next string format character

	if ( !c ) goto end; // end of format string reached

	if( 	c != '%'  // not format specifier
		|| (c = format[fi++]) == '%' ) // %% sequence		
	{
		buffer[ bi++ ] = c; // out character raw
	}
	else
	{			
		unsigned char fillChar, fmatExp, signChar = 0;
		bit alignLeft = false;

		fillChar = ' '; // default fill character

		// process flags
		while( 1 )
		{				
			switch( c )
			{
				case 0: goto end; // end of format string reached
				case '0': fillChar = '0'; break;
				case ' ': signChar = ' '; break; // show space in front of signed numbers when not negative						
				case '+': signChar = '+'; break; // show sign on signed numbers						
				case '-': alignLeft = true; break;
				default: goto processWidth;
			}

			c = format[fi++]; // read next string format character
		}

processWidth:
		// Next character if valid digit indicates field width
		unsigned char fieldWidth = 0;
		while( c >= '0' && c <= '9' )
		{			
			fieldWidth *= 10;
			fieldWidth += c;
			fieldWidth -= '0';				
			c = format[fi++]; // read next string format character
			if ( !c ) goto end; // end of format string reached
		}	

		if(c == '.')
		{
			c = format[fi++]; // read next string format character
			if ( !c ) goto end; // end of format string reached
		}

		unsigned char precision = 0;
		signed char preTmp = -1;
		while( c >= '0' && c <= '9' )
		{
			if(preTmp == -1) preTmp = 0;
			preTmp *= 10;
			preTmp += c;
			preTmp -= '0';				
			c = format[fi++]; // read next string format character
			if ( !c ) goto end; // end of format string reached
		}				

		if(preTmp == -1) precision = 6;
		else			 precision = (unsigned char)preTmp;

		unsigned char floatSignChar = 0;

		// Next character indicates the format
		switch( c )
		{
			case 'f':
			{
				fmatExp = 0;		
				if(float32_lt(val, 0) != 0)		// Negative values must be adjusted to be positive
				{
					floatSignChar = '-';
					val = float32_mul(val, -1.0);
				}
			}
			break;
				// fall through to next case intended
			default: goto end; // no radix specified - can't continue
		}

		long valFr = 0;
		long valIn = 0;

		long frMag = 1;

		for(unsigned char ct = 0; ct < precision; ++ct)
		{
			valFr *= 10;
			frMag *= 10;
		}

		valIn = (long)float32_to_int32(val);
		valFr = (long)float32_to_int32(float32_add(float32_mul(float32_sub(val, float32_from_int32((int32)valIn)), float32_from_int32(frMag)), 0.5));
		frMag--;

		if((valIn == 0) && (valFr == 0)) floatSignChar = 0;

		if(floatSignChar == '-') signChar = floatSignChar;

		unsigned char digitWidth = 0;
		{	// determine number of digits in the number
			unsigned long valx = (unsigned long)float32_to_int32(val);
			do
			{
				digitWidth++;
				valx /= 10;
			}
			while( valx != 0 );
		}

		// sign takes extra space
		if( signChar != 0 )
			digitWidth++;

		// precision takes extra space
		digitWidth += precision;

		// one for the decimal point
		if(precision > 0) digitWidth++;


		// limit digit width to prevent buffer overflow - not standard, but will prevent buffer overrun corruption
		unsigned char fillWidth = 0;
		if( fieldWidth > 0 )
		{
			if( digitWidth > fieldWidth ) // limit number of digits
				digitWidth = fieldWidth;

			fillWidth = fieldWidth - digitWidth;
		}
		else
		{
			fieldWidth = digitWidth;
			fillWidth = 0;
		}

		if( fillWidth == 0 )
			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( digitWidth != 0 || fillWidth != 0  )
		{
			bi--;
			if( doFill )
			{
				c = fillChar;
				if( --fillWidth == 0 )
				 doFill = false;
			}
			else
			{
				if( --digitWidth == 0 && signChar ) // last digit and sign character
					c = signChar; // add sign character
				else if (frMag > 0)
				{
					c = (char)(valFr % 10 );											
					c += '0'; // convert to digit character 0-9
					valFr /= 10;
					frMag /= 10;
					fracFlg = 0;
				}
				else if ((valFr == 0) && (precision > 0))
				{
					if(fracFlg > 0)
					{
						c = '0';											
						fracFlg--;
					}
					else
					{
						c = '.';											
						valFr = -1;
					}
				}
				else
				{
					c = (char)(valIn % 10 );											
					c += '0'; // convert to digit character 0-9
					valIn /= 10;					
				}

				if( digitWidth == 0 )
					doFill = true;
			}

			buffer[ bi ] = c;
		}
		bi += fieldWidth;
	}	
}
end:	
buffer[ bi ] = 0; // null terminate
return bi;
}



float
atof(char* str)
{
if((str == 0) || (strlen(str) == 0)) return 0;

int hl = 0;
int fr = 0;
float retVal = 0;

char* pntPtr = strrchr(str,'.');
for(char ct = strlen(str) - 1; ct >=0; --ct)
{
	if(str[ct] != '0') break;

	str[ct] = 0;
}


if(pntPtr != 0)
{
	fr = atoi(pntPtr + 1);
	*pntPtr = 0;
}

hl = atoi(str);

char foo[32];
if(pntPtr != 0)
{
	*pntPtr = '.';

	float pwr = 1.0;
	for(unsigned char ct = 0; ct < strlen(pntPtr + 1); ++ct)
	{
		pwr = float32_mul(pwr, 10.0);
	}

	retVal = float32_div(float32_from_int32((int32) fr), pwr);
}


if(hl < 0)
{
	retVal = float32_sub(float32_from_int32((int32) hl), retVal);
}
else
{
	retVal = float32_add(float32_from_int32((int32) hl), retVal);
}

return retVal;
}

Edited by Neil

Share this post


Link to post
Share on other sites

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(nv, float32_mul(float32_from_int32((int32)cnt), float32_from_int32((int32)(cnt + 1))));

 

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

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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(x, x);
for(unsigned char cnt = 2; cnt < 50; cnt+=2)

 

Cheers

 

Reynard

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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(float x)
{
float_rounding_mode = float_round_down;

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

x = float32_mul(x, x);
for(unsigned char cnt = 2; cnt < 50; cnt+=2)
{
	ac = float32_add(ac, float32_div(xv, nv));
	nv = float32_mul(nv, float32_from_int32((int32)((unsigned int)(cnt * (cnt + 1)))));
	xv = float32_mul(xv, x);
	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

Share this post


Link to post
Share on other sites

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!!

Share this post


Link to post
Share on other sites

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(float x)
{
float p[4] = {
	-0.499999993585,
	 0.041666636258,
	-0.0013888361399,
	 0.00002476016134
};
float frac, t, y, t2 = 1.0;
unsigned char quad, i;

if (float32_lt(x,  0)) {
	x = x ^ 0x80000000l;			 // absolute value of input
}

quad = (unsigned char)float32_to_int32(float32_div(x , PI_DIV_BY_TWO));
frac = float32_sub(float32_div(x , PI_DIV_BY_TWO), float32_from_int32((int32)quad));
quad = quad % 4;					// quadrant (0 to 3)

if (quad == 0 || quad == 2)
	t = float32_mul(frac, PI_DIV_BY_TWO);
else if (quad == 1)
	t = float32_mul(float32_sub(1.0, frac), PI_DIV_BY_TWO);
else // should be 3
	t = float32_mul(float32_sub(frac, 1.0), PI_DIV_BY_TWO);
y = 0.999999999781;
t = float32_mul(t, t);

for (i = 0; i <= 3; i++)
{
	t2 = float32_mul(t2, t);
	y = float32_add(y, float32_mul(p[i], t2));
}

if (quad == 2 || quad == 1)
y = y ^ 0x80000000l;  // correct sign

return y;
}

float
sin(float x)
{
return cos(float32_sub(x, PI_DIV_BY_TWO));
}

 

Cheers

 

Reynard

Share this post


Link to post
Share on other sites

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.

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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(2)

float exp(float x)
{
float y, r, twopwrn;
int n;
unsigned char sign;

float_rounding_mode = float_round_down;

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

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

if (float32_lt(x, 0.0))		// 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 (2^n).
*((long*)&twopwrn) = (long)((n + 127) << 23);

// y = (y/ln(2)) - n
y = float32_sub(float32_div(y, LN2), float32_from_int32(n));

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

// Calc. res = 2^n(1 + y)
y = float32_mul(twopwrn, float32_add(1.0, y));

if (sign)		// was x negative ?
	y = float32_div(1.0, y);	// yes, invert. y = 1/y

return(y);
}

 

Cheers

 

Reynard

Share this post


Link to post
Share on other sites

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.

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites

Your content will need to be approved by a moderator

Guest
You are commenting as a guest. If you have an account, please sign in.
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoticons maximum are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...

×