Neil 0 Report post Posted November 24, 2008 (edited) 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 November 27, 2008 by Neil Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted November 27, 2008 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 Quote Share this post Link to post Share on other sites

Neil 0 Report post Posted November 27, 2008 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 Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted November 27, 2008 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 Quote Share this post Link to post Share on other sites

Dave 0 Report post Posted November 28, 2008 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 Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted November 28, 2008 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 Quote Share this post Link to post Share on other sites

Neil 0 Report post Posted November 28, 2008 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!! Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted November 28, 2008 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 Quote Share this post Link to post Share on other sites

PeterK 0 Report post Posted April 24, 2010 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. Quote Share this post Link to post Share on other sites

Dave 0 Report post Posted April 25, 2010 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 Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted April 25, 2010 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 Quote Share this post Link to post Share on other sites

Reynard 0 Report post Posted April 25, 2010 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 Quote Share this post Link to post Share on other sites

PeterK 0 Report post Posted April 25, 2010 You are right about "Pade approximant", which is even better then "Chebyshev Approximation". I use the "Polynomial Approximation from Chebyshev Coefﬁcients" 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 Coefﬁcients", 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. Quote Share this post Link to post Share on other sites

Neil 0 Report post Posted June 3, 2010 You are right about "Pade approximant", which is even better then "Chebyshev Approximation". I use the "Polynomial Approximation from Chebyshev Coefﬁcients" 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 Coefﬁcients", 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 Quote Share this post Link to post Share on other sites