From 1071a6c6ae98b2f1017cc0bc323860ac0d145d8c Mon Sep 17 00:00:00 2001 From: johanknol Date: Fri, 11 Jan 2002 09:17:33 +0000 Subject: [PATCH] added math libs of Jesus Calvino-Fraga git-svn-id: https://sdcc.svn.sourceforge.net/svnroot/sdcc/trunk/sdcc@1787 4a8a32a2-be11-0410-ad9d-d568d2c75423 --- device/include/math.h | 102 +++++++++++++++++++++++++++++----------- device/lib/Makefile.in | 6 ++- device/lib/acosf.c | 31 ++++++++++++ device/lib/asincosf.c | 80 +++++++++++++++++++++++++++++++ device/lib/asinf.c | 32 +++++++++++++ device/lib/atan2f.c | 45 ++++++++++++++++++ device/lib/atanf.c | 70 +++++++++++++++++++++++++++ device/lib/ceilf.c | 31 ++++++++++++ device/lib/cosf.c | 29 ++++++++++++ device/lib/coshf.c | 28 +++++++++++ device/lib/cotf.c | 41 ++++++++++++++++ device/lib/expf.c | 86 +++++++++++++++++++++++++++++++++ device/lib/fabsf.c | 31 ++++++++++++ device/lib/floorf.c | 31 ++++++++++++ device/lib/frexpf.c | 37 +++++++++++++++ device/lib/ldexpf.c | 36 ++++++++++++++ device/lib/libfloat.lib | 26 ++++++++++ device/lib/log10f.c | 27 +++++++++++ device/lib/logf.c | 65 +++++++++++++++++++++++++ device/lib/modff.c | 27 +++++++++++ device/lib/powf.c | 31 ++++++++++++ device/lib/sincosf.c | 85 +++++++++++++++++++++++++++++++++ device/lib/sincoshf.c | 89 +++++++++++++++++++++++++++++++++++ device/lib/sinf.c | 29 ++++++++++++ device/lib/sinhf.c | 28 +++++++++++ device/lib/sqrtf.c | 51 ++++++++++++++++++++ device/lib/tancotf.c | 85 +++++++++++++++++++++++++++++++++ device/lib/tanf.c | 29 ++++++++++++ device/lib/tanhf.c | 60 +++++++++++++++++++++++ 29 files changed, 1319 insertions(+), 29 deletions(-) create mode 100644 device/lib/acosf.c create mode 100644 device/lib/asincosf.c create mode 100644 device/lib/asinf.c create mode 100644 device/lib/atan2f.c create mode 100644 device/lib/atanf.c create mode 100644 device/lib/ceilf.c create mode 100644 device/lib/cosf.c create mode 100644 device/lib/coshf.c create mode 100644 device/lib/cotf.c create mode 100644 device/lib/expf.c create mode 100644 device/lib/fabsf.c create mode 100644 device/lib/floorf.c create mode 100644 device/lib/frexpf.c create mode 100644 device/lib/ldexpf.c create mode 100644 device/lib/log10f.c create mode 100644 device/lib/logf.c create mode 100644 device/lib/modff.c create mode 100644 device/lib/powf.c create mode 100644 device/lib/sincosf.c create mode 100644 device/lib/sincoshf.c create mode 100644 device/lib/sinf.c create mode 100644 device/lib/sinhf.c create mode 100644 device/lib/sqrtf.c create mode 100644 device/lib/tancotf.c create mode 100644 device/lib/tanf.c create mode 100644 device/lib/tanhf.c diff --git a/device/include/math.h b/device/include/math.h index dc237d1c..326ab00e 100644 --- a/device/include/math.h +++ b/device/include/math.h @@ -1,33 +1,79 @@ -/*------------------------------------------------------------------------- - math.h - ANSI functions forward declarations - - Written By - Sandeep Dutta . sandeep.dutta@usa.net (1998) - - This program is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2, or (at your option) any - later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program; if not, write to the Free Software - Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - - In other words, you are welcome to use, share and improve this program. - You are forbidden to forbid anyone else to use, share and improve - what you give them. Help stamp out software-hoarding! --------------------------------------------------------------------------*/ - -#ifndef __SDC51_MATH_H -#define __SDC51_MATH_H 1 -#error Floating point not yet completely implemented -#endif +/* math.h: Floating point math function declarations + Copyright (C) 2001 Jesus Calvino-Fraga, jesusc@ieee.org + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +/* Version 1.0 - Initial release */ + +#ifndef _INC_MATH +#define _INC_MATH + +#define PI 3.1415926536 +#define TWO_PI 6.2831853071 +#define HALF_PI 1.5707963268 +#define QUART_PI 0.7853981634 +#define iPI 0.3183098862 +#define iTWO_PI 0.1591549431 +#define TWO_O_PI 0.6366197724 + +// EPS=B**(-t/2), where B is the radix of the floating-point representation +// and there are t base-B digits in the significand. Therefore, for floats +// EPS=2**(-12). Also define EPS2=EPS*EPS. +#define EPS 244.14062E-6 +#define EPS2 59.6046E-9 +#define XMAX 3.402823466E+38 + +union float_long +{ + float f; + long l; +}; + +/********************************************** + * Prototypes for float ANSI C math functions * + **********************************************/ + +/* Trigonometric functions */ +float sinf(const float x); +float cosf(const float x); +float tanf(const float x); +float cotf(const float x); +float asinf(const float x); +float acosf(const float x); +float atanf(const float x); +float atan2f(const float x, const float y); + +/* Hyperbolic functions */ +float sinhf(const float x); +float coshf(const float x); +float tanhf(const float x); + +/* Exponential, logarithmic and power functions */ +float expf(const float x); +float logf(const float x); +float log10f(const float x); +float powf(const float x, const float y); +float sqrtf(const float a); + +/* Nearest integer, absolute value, and remainder functions */ +float fabsf(const float x); +float frexpf(const float x, int *pw2); +float ldexpf(const float x, const int pw2); +float ceilf(float x); +float floorf(float x); +float modff(float x, float * y); + +#endif /* _INC_MATH */ diff --git a/device/lib/Makefile.in b/device/lib/Makefile.in index b9b62da8..6927200d 100644 --- a/device/lib/Makefile.in +++ b/device/lib/Makefile.in @@ -57,7 +57,11 @@ SOURCES = _atoi.c _atol.c _autobaud.c _bp.c _schar2fs.c \ _strstr.c _strtok.c _uchar2fs.c _uint2fs.c \ _ulong2fs.c malloc.c serial.c ser_ir.c printfl.c \ printf_large.c vprintf.c puts.c gets.c \ - assert.c _strcat.c time.c printf_fast.c bpx.c + assert.c _strcat.c time.c printf_fast.c bpx.c \ + fabsf.c frexpf.c ldexpf.c expf.c powf.c sincosf.c sinf.c \ + cosf.c logf.c log10f.c sqrtf.c tancotf.c tanf.c cotf.c \ + asincosf.c asinf.c acosf.c atanf.c atan2f.c sincoshf.c \ + sinhf.c coshf.c tanhf.c floorf.c ceilf.c modff.c OBJECTS = $(patsubst %.c,$(PORTDIR)/%.rel,$(SOURCES)) diff --git a/device/lib/acosf.c b/device/lib/acosf.c new file mode 100644 index 00000000..110ffa40 --- /dev/null +++ b/device/lib/acosf.c @@ -0,0 +1,31 @@ +/* acosf.c: Computes arc cosine of a 32-bit float + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float asincosf(const float x, const int isacos); + +float acosf(const float x) reentrant +{ + if(x== 1.0) return 0.0; + else if(x==-1.0) return PI; + else if(x== 0.0) return HALF_PI; + return asincosf(x,1); +} diff --git a/device/lib/asincosf.c b/device/lib/asincosf.c new file mode 100644 index 00000000..b8005ecb --- /dev/null +++ b/device/lib/asincosf.c @@ -0,0 +1,80 @@ +/* asincosf.c: Computes asin or acos of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P1 0.933935835E+0 +#define P2 -0.504400557E+0 +#define Q0 0.560363004E+1 +#define Q1 -0.554846723E+1 +#define Q2 0.100000000E+1 + +#define P(g) (P2*g+P1) +#define Q(g) ((Q2*g+Q1)*g+Q0) + +float asincosf(const float x, const int isacos) +{ + float y, g, r; + int i; + float a[2]={ 0.0, QUART_PI }; + float b[2]={ HALF_PI, QUART_PI }; + + y=fabsf(x); + i=isacos; + if (y < EPS) r=y; + else + { + if (y > 0.5) + { + i=1-i; + if (y > 1.0) + { + errno=EDOM; + return 0.0; + } + g=(0.5-y)+0.5; + g=ldexpf(g,-1); + y=sqrtf(g); + y=-(y+y); + } + else + { + g=y*y; + } + r=y+y*((P(g)*g)/Q(g)); + } + if (isacos) + { + if (x < 0.0) + r=(b[i]+r)+b[i]; + else + r=(a[i]-r)+a[i]; + } + else + { + r=(a[i]+r)+a[i]; + if (x<0.0) r=-r; + } + return r; +} diff --git a/device/lib/asinf.c b/device/lib/asinf.c new file mode 100644 index 00000000..29291364 --- /dev/null +++ b/device/lib/asinf.c @@ -0,0 +1,32 @@ +/* asinf.c: Computes asin(x) + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float asincosf(const float x, const int isacos); + +float asinf(const float x) reentrant +{ + if(x== 1.0) return HALF_PI; + else if(x==-1.0) return -HALF_PI; + else if(x== 0.0) return 0.0; + else return asincosf(x,0); +} + diff --git a/device/lib/atan2f.c b/device/lib/atan2f.c new file mode 100644 index 00000000..09d6ff96 --- /dev/null +++ b/device/lib/atan2f.c @@ -0,0 +1,45 @@ +/* atan2f.c: Computes atan2(x) where x is a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float atan2f(const float x, const float y) +{ + float r; + + if ((x==0.0) && (y==0.0)) + { + errno=EDOM; + return 0.0; + } + + if(fabsf(y)>=fabsf(x)) + { + r=atanf(x/y); + if(y<0.0) r+=(x>=0?PI:-PI); + } + else + { + r=-atanf(y/x); + r+=(x<0.0?-HALF_PI:HALF_PI); + } + return r; +} diff --git a/device/lib/atanf.c b/device/lib/atanf.c new file mode 100644 index 00000000..a558948a --- /dev/null +++ b/device/lib/atanf.c @@ -0,0 +1,70 @@ +/* atanf.c: Computes arctan of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P0 -0.4708325141E+0 +#define P1 -0.5090958253E-1 +#define Q0 0.1412500740E+1 +#define Q1 0.1000000000E+1 + +#define P(g,f) ((P1*g+P0)*g*f) +#define Q(g) (Q1*g+Q0) + +#define K1 0.2679491924 /* 2-sqrt(3) */ +#define K2 0.7320508076 /* sqrt(3)-1 */ +#define K3 1.7320508076 /* sqrt(3) */ + +float atanf(const float x) reentrant +{ + float f, r, g; + int n=0; + static float a[]={ 0.0, 0.5235987756, 1.5707963268, 1.0471975512 }; + + f=fabsf(x); + if(f>1.0) + { + f=1.0/f; + n=2; + } + if(f>K1) + { + f=((K2*f-1.0)+f)/(K3+f); + // What it is actually wanted is this more accurate formula, + // but SDCC optimizes it and then it does not work: + // f=(((K2*f-0.5)-0.5)+f)/(K3+f); + n++; + } + if(fabsf(f)1) r=-r; + r+=a[n]; + if(x<0.0) r=-r; + return r; +} + diff --git a/device/lib/ceilf.c b/device/lib/ceilf.c new file mode 100644 index 00000000..f3599f8a --- /dev/null +++ b/device/lib/ceilf.c @@ -0,0 +1,31 @@ +/* ceilf.c: Returns the integer larger or equal than x + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float ceilf(float x) reentrant +{ + long r; + r=x; + if (r<0) + return r; + else + return (r+((r + +float sincosf(const float x, const int iscos); + +float cosf(const float x) reentrant +{ + if (x==0.0) return 1.0; + return sincosf(x, 1); +} diff --git a/device/lib/coshf.c b/device/lib/coshf.c new file mode 100644 index 00000000..3a12eda8 --- /dev/null +++ b/device/lib/coshf.c @@ -0,0 +1,28 @@ +/* coshf.c: Computes cosh(x) where x is a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float sincoshf(const float x, const int iscosh); + +float coshf(const float x) reentrant +{ + return sincoshf(x, 1); +} diff --git a/device/lib/cotf.c b/device/lib/cotf.c new file mode 100644 index 00000000..9e8329a5 --- /dev/null +++ b/device/lib/cotf.c @@ -0,0 +1,41 @@ +/* cotf.c: Computes cot(x) where x is a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float tancotf(const float x, const int iscot); + +float cotf(const float x) reentrant +{ + float y; + + y=fabsf(x); + if (y<1.0E-30) //This one requires more thinking... + { + errno = ERANGE; + if (x<0.0) + return -XMAX; + else + return XMAX; + } + return tancotf(x, 1); +} + diff --git a/device/lib/expf.c b/device/lib/expf.c new file mode 100644 index 00000000..083de0c2 --- /dev/null +++ b/device/lib/expf.c @@ -0,0 +1,86 @@ +/* expf.c: Computes e**x of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P0 0.2499999995E+0 +#define P1 0.4160288626E-2 +#define Q0 0.5000000000E+0 +#define Q1 0.4998717877E-1 + +#define P(z) ((P1*z)+P0) +#define Q(z) ((Q1*z)+Q0) + +#define C1 0.693359375 +#define C2 -2.1219444005469058277e-4 + +#define BIGX 88.72283911 /* ln(XMAX) */ +#define EXPEPS 1.0E-7 /* exp(1.0E-7)=0.0000001 */ +#define K1 1.4426950409 /* 1/ln(2) */ + +float expf(const float x) +{ + int n; + float xn, g, r, z, y; + int sign; + + if(x>=0.0) + { y=x; sign=0; } + else + { y=-x; sign=1; } + + if(yBIGX) + { + if(sign) + { + errno=ERANGE; + return XMAX; + } + else + { + return 0.0; + } + } + + z=y*K1; + n=z; + + if(n<0) --n; + if(z-n>=0.5) ++n; + xn=n; + g=((y-xn*C1))-xn*C2; + z=g*g; + r=P(z)*g; + r=0.5+(r/(Q(z)-r)); + + n++; + z=ldexpf(r, n); + if(sign) + return 1.0/z; + else + return z; +} + diff --git a/device/lib/fabsf.c b/device/lib/fabsf.c new file mode 100644 index 00000000..827a99a1 --- /dev/null +++ b/device/lib/fabsf.c @@ -0,0 +1,31 @@ +/* fabsf.c: Returns the absolute value of a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float fabsf(const float x) reentrant +{ + union float_long fl; + + fl.f = x; + fl.l &= 0x7fffffff; + return fl.f; +} diff --git a/device/lib/floorf.c b/device/lib/floorf.c new file mode 100644 index 00000000..b0f474f3 --- /dev/null +++ b/device/lib/floorf.c @@ -0,0 +1,31 @@ +/* floorf.c: Returns the integer smaller or equal than x + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float floorf (float x) reentrant +{ + long r; + r=x; + if (r<=0) + return (r+((r>x)?-1:0)); + else + return r; +} diff --git a/device/lib/frexpf.c b/device/lib/frexpf.c new file mode 100644 index 00000000..df2e3815 --- /dev/null +++ b/device/lib/frexpf.c @@ -0,0 +1,37 @@ +/* frexpf.c: Returns the exponent and mantisa of a 32 bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float frexpf(const float x, int *pw2) +{ + union float_long fl; + long int i; + + fl.f=x; + /* Find the exponent (power of 2) */ + i = ( fl.l >> 23) & 0x000000ff; + i -= 0x7e; + *pw2 = i; + fl.l &= 0x807fffff; /* strip all exponent bits */ + fl.l |= 0x3f000000; /* mantissa between 0.5 and 1 */ + return(fl.f); +} diff --git a/device/lib/ldexpf.c b/device/lib/ldexpf.c new file mode 100644 index 00000000..549fff63 --- /dev/null +++ b/device/lib/ldexpf.c @@ -0,0 +1,36 @@ +/* ldexpf.c: Build a float from a mantisa and exponent. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float ldexpf(const float x, const int pw2) +{ + union float_long fl; + long e; + + fl.f = x; + + e=(fl.l >> 23) & 0x000000ff; + e+=pw2; + fl.l= ((e & 0xff) << 23) | (fl.l & 0x807fffff); + + return(fl.f); +} diff --git a/device/lib/libfloat.lib b/device/lib/libfloat.lib index ad35b908..20f89112 100644 --- a/device/lib/libfloat.lib +++ b/device/lib/libfloat.lib @@ -18,3 +18,29 @@ _fs2ulong _fs2schar _fs2sint _fs2slong +fabsf +frexpf +ldexpf +expf +powf +sincosf +sinf +cosf +logf +log10f +sqrtf +tancotf +tanf +cotf +asincosf +asinf +acosf +atanf +atan2f +sincoshf +sinhf +coshf +tanhf +floorf +ceilf +modff diff --git a/device/lib/log10f.c b/device/lib/log10f.c new file mode 100644 index 00000000..13f77cda --- /dev/null +++ b/device/lib/log10f.c @@ -0,0 +1,27 @@ +/* log10f.c: Computes the base 10 log of a 32 bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float log10f(const float x) reentrant +{ + return logf(x)*0.4342944819; +} diff --git a/device/lib/logf.c b/device/lib/logf.c new file mode 100644 index 00000000..2d3599ce --- /dev/null +++ b/device/lib/logf.c @@ -0,0 +1,65 @@ +/* logf.c: Computes the natural log of a 32 bit float as outlined in [1]. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +/*Constans for 24 bits or less (8 decimal digits)*/ +#define A0 -0.5527074855E+0 +#define B0 -0.6632718214E+1 +#define A(w) (A0) +#define B(w) (w+B0) + +#define C0 0.70710678118654752440 +#define C1 0.693359375 /*355.0/512.0*/ +#define C2 -2.121944400546905827679E-4 + +float logf(const float x) reentrant +{ + float Rz, f, z, w, znum, zden, xn; + int n; + + if (x<=0.0) + { + errno=EDOM; + return 0.0; + } + f=frexpf(x, &n); + znum=f-0.5; + if (f>C0) + { + znum-=0.5; + zden=(f*0.5)+0.5; + } + else + { + n--; + zden=znum*0.5+0.5; + } + z=znum/zden; + w=z*z; + + Rz=z+z*(w*A(w)/B(w)); + xn=n; + return ((xn*C2+Rz)+xn*C1); +} diff --git a/device/lib/modff.c b/device/lib/modff.c new file mode 100644 index 00000000..ace1565b --- /dev/null +++ b/device/lib/modff.c @@ -0,0 +1,27 @@ +/* modff.c: Returns both the integer and fraction of a float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float modff(float x, float * y) +{ + *y=((int)x); + return (x-*y); +} diff --git a/device/lib/powf.c b/device/lib/powf.c new file mode 100644 index 00000000..ef78f586 --- /dev/null +++ b/device/lib/powf.c @@ -0,0 +1,31 @@ +/* powf.c: Computes x**y where x and y are 32-bit floats. + WARNING: less that 6 digits accuracy. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float powf(const float x, const float y) +{ + if(y == 0.0) return 1.0; + if(y==1.0) return x; + if(x <= 0.0) return 0.0; + return expf(logf(x) * y); +} diff --git a/device/lib/sincosf.c b/device/lib/sincosf.c new file mode 100644 index 00000000..e4ed08c0 --- /dev/null +++ b/device/lib/sincosf.c @@ -0,0 +1,85 @@ +/* sincosf.c: Computes sin or cos of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define r1 -0.1666665668E+0 +#define r2 0.8333025139E-2 +#define r3 -0.1980741872E-3 +#define r4 0.2601903036E-5 + +/* PI=C1+C2 */ +#define C1 3.140625 +#define C2 9.676535897E-4 + +/*A reasonable value for YMAX is the int part of PI*B**(t/2)=3.1416*2**(12)*/ +#define YMAX 12867.0 + +float sincosf(const float x, const int iscos) +{ + float y, f, r, g, XN; + int N, sign; + + if(iscos) + { + y=fabsf(x)+HALF_PI; + sign=0; + } + else + { + if(x<0.0) + { y=-x; sign=1; } + else + { y=x; sign=0; } + } + + if(y>YMAX) + { + errno=ERANGE; + return 0.0; + } + + /*Round y/PI to the nearest integer*/ + N=((y*iPI)+0.5); /*y is positive*/ + + /*If N is odd change sign*/ + if(N&1) sign=~sign; + + XN=N; + /*Cosine required? (is done here to keep accuracy)*/ + if(iscos) XN-=0.5; + + y=fabsf(x); + r=(int)y; + g=y-r; + f=((r-XN*C1)+g)-XN*C2; + + g=f*f; + if(g>EPS2) //Used to be if(fabsf(f)>EPS) + { + r=(((r4*g+r3)*g+r2)*g+r1)*g; + f+=f*r; + } + return (sign?-f:f); +} diff --git a/device/lib/sincoshf.c b/device/lib/sincoshf.c new file mode 100644 index 00000000..4339ae98 --- /dev/null +++ b/device/lib/sincoshf.c @@ -0,0 +1,89 @@ +/* sincoshf.c: Computes sinh or cosh of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P0 -0.713793159E+1 +#define P1 -0.190333999E+0 +#define Q0 -0.428277109E+2 +#define Q1 0.100000000E+1 + +#define P(z) (P1*z+P0) +#define Q(z) (Q1*z+Q0) + +#define K1 0.69316101074218750000E+0 /* ln(v) */ +#define K2 0.24999308500451499336E+0 /* v**(-2) */ +#define K3 0.13830277879601902638E-4 /* v/2-1 */ + +//WMAX is defined as ln(XMAX)-ln(v)+0.69 +#define WMAX 44.93535952E+0 +//WBAR 0.35*(b+1) +#define WBAR 1.05 +#define YBAR 9.0 /*Works for me*/ + +float sincoshf(const float x, const int iscosh) +{ + float y, w, z; + int sign; + + if (x<0.0) { y=-x; sign=1; } + else { y=x; sign=0; } + + if ((y>1.0) || iscosh) + { + if(y>YBAR) + { + w=y-K1; + if (w>WMAX) + { + errno=ERANGE; + z=XMAX; + } + else + { + z=expf(w); + z+=K3*z; + } + } + else + { + z=expf(y); + w=1.0/z; + if(!iscosh) w=-w; + z=(z+w)*0.5; + } + if(sign) z=-z; + } + else + { + if (y + +float sincosf(const float x, const int iscos); + +float sinf(const float x) reentrant +{ + if (x==0.0) return 0.0; + return sincosf(x, 0); +} diff --git a/device/lib/sinhf.c b/device/lib/sinhf.c new file mode 100644 index 00000000..3741d401 --- /dev/null +++ b/device/lib/sinhf.c @@ -0,0 +1,28 @@ +/* sinhf.c: Computes sinh(x) where x is a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float sincoshf(const float x, const int iscosh); + +float sinhf(const float x) reentrant +{ + return sincoshf(x, 0); +} diff --git a/device/lib/sqrtf.c b/device/lib/sqrtf.c new file mode 100644 index 00000000..477d535a --- /dev/null +++ b/device/lib/sqrtf.c @@ -0,0 +1,51 @@ +/* sqrtf.c: Computes square root of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +float sqrtf(const float x) reentrant +{ + float f, y; + int n; + + if (x==0.0) return x; + else if (x==1.0) return 1.0; + else if (x<0.0) + { + errno=EDOM; + return 0.0; + } + f=frexpf(x, &n); + y=0.41731+0.59016*f; /*Educated guess*/ + /*For a 24 bit mantisa (float), two iterations are sufficient*/ + y+=f/y; + y=ldexpf(y, -2) + f/y; /*Faster version of 0.25 * y + f/y*/ + + if (n&1) + { + y*=0.7071067812; + ++n; + } + return ldexpf(y, n/2); +} diff --git a/device/lib/tancotf.c b/device/lib/tancotf.c new file mode 100644 index 00000000..72696127 --- /dev/null +++ b/device/lib/tancotf.c @@ -0,0 +1,85 @@ +/* tancotf.c: Computes tan or cot of a 32-bit float as outlined in [1] + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P0 0.100000000E+1 +#define P1 -0.958017723E-1 +#define Q0 0.100000000E+1 +#define Q1 -0.429135777E+0 +#define Q2 0.971685835E-2 + +#define C1 1.5703125 +#define C2 4.83826794897E-4 + +#define P(f,g) (P1*g*f+f) +#define Q(g) ((Q2*g+Q1)*g+Q0) + +//A reasonable choice for YMAX is the integer part of B**(t/2)*PI/2: +#define YMAX 6433.0 + +float tancotf(const float x, const int iscotan) +{ + float f, g, xn, xnum, xden; + int n; + + if (fabsf(x) > YMAX) + { + errno = ERANGE; + return 0.0; + } + + /*Round x*2*PI to the nearest integer*/ + n=(x*TWO_O_PI+(x>0.0?0.5:-0.5)); /*works for +-x*/ + xn=n; + + xnum=(int)x; + xden=x-xnum; + f=((xnum-xn*C1)+xden)-xn*C2; + + if (fabsf(f) < EPS) + { + xnum = f; + xden = 1.0; + } + else + { + g = f*f; + xnum = P(f,g); + xden = Q(g); + } + + if(n&1) + //xn is odd + { + if(iscotan) return (-xnum/xden); + else return (-xden/xnum); + } + else + { + if(iscotan) return (xden/xnum); + else return (xnum/xden); + } +} + diff --git a/device/lib/tanf.c b/device/lib/tanf.c new file mode 100644 index 00000000..48c18bff --- /dev/null +++ b/device/lib/tanf.c @@ -0,0 +1,29 @@ +/* tanf.c: Computes tan(x) where x is a 32-bit float. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* Version 1.0 - Initial release */ + +#include + +float tancotf(const float x, const int iscot); + +float tanf(const float x) reentrant +{ + return tancotf(x, 0); +} + diff --git a/device/lib/tanhf.c b/device/lib/tanhf.c new file mode 100644 index 00000000..5737c001 --- /dev/null +++ b/device/lib/tanhf.c @@ -0,0 +1,60 @@ +/* tanhf.c: Computes tanh(x) where x is a 32-bit float as outlined in [1]. + + Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +/* [1] William James Cody and W. M. Waite. _Software manual for the + elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */ + +/* Version 1.0 - Initial release */ + +#include +#include + +#define P0 -0.8237728127E+0 +#define P1 -0.3831010665E-2 +#define Q0 0.2471319654E+1 +#define Q1 0.1000000000E+1 + +/* ln(3)/2 */ +#define K1 0.5493061443E+0 +/* SBIG=[ln(2)+(t+1)*ln(B)]/2 */ +#define SBIG 9.01091 + +#define P(g) ((P1*g+P0)*g) +#define Q(g) (Q1*g+Q0) + +float tanhf(const float x) reentrant +{ + float f, g, r; + + f=fabsf(x); + if(f>SBIG) r=1.0; + else if(f>K1) + { + r=0.5-1.0/(expf(f+f)+1.0); + r+=r; + } + else if(f