1
0
mirror of https://github.com/cc65/cc65.git synced 2024-06-08 15:29:37 +00:00

some cleanup, updated softmath with some google results :), updated readme

This commit is contained in:
mrdudz 2023-09-04 01:05:12 +02:00
parent d771e86cdf
commit ee46184c62
18 changed files with 480 additions and 195 deletions

View File

@ -38,7 +38,13 @@ float __fastcall__ roundf(float x); /* C99 */
float __fastcall__ truncf(float x); /* C99 */
/* double ceil(double x) */
float ceilf(float x); /* C99 */
float __fastcall__ ceilf(float x); /* C99 */
/* double fmod(double x, double y); */
float __fastcall__ fmodf(float x, float y);
/* double floor(double x); */
float __fastcall__ floorf(float x);
/* beware, this is not standard */
#ifndef M_PI
@ -48,8 +54,5 @@ float ceilf(float x); /* C99 */
/* FIXME */
float __fastcall__ _fatan2(float x, float y);
float ffloor(float x);
#endif /* _MATH_H_ */

View File

@ -1,26 +0,0 @@
#include <_float.h>
#include <math.h>
#include "cbmfp.h"
#if 0
#define _fcmplt(_d, _s) (_fcmp((_d), (_s)) == 1)
#define _fcmpgt(_d, _s) (_fcmp((_d), (_s)) == 255)
#define _fcmpeq(_d, _s) (_fcmp((_d), (_s)) == 0)
/* FIXME: this is really too simple */
float ffloor(float x)
{
signed long n;
float d;
n = _ftoi(x); /* FIXME: long */
d = _itof(n); /* FIXME: long */
if (_fcmpeq(d, x) || _fcmplt(_itof(0), x)) {
return d;
}
return _fsub(d, 1);
}
#endif

View File

@ -1,136 +0,0 @@
#include <_float.h>
#include <math.h>
#include "cbmfp.h"
#if 0
// convert float into a string. this is surprisingly complex, so we just use
// the kernal function, and then fix up the result
char *ftoa(char *buf, float n)
{
char i, ii, epos = 0, ex;
char tempbuf[0x20];
_ftostr(tempbuf, n);
// find position of the 'e'
i=0;while(tempbuf[i]) {
if (tempbuf[i] == 69) { /* 'e' */
epos = i;
break;
}
i++;
}
if (epos == 0) {
i = ii = 0;
// no exponent, we can return the number as is
if (tempbuf[i] == '-') {
buf[ii] = tempbuf[i];
i++;ii++;
}
// else {
// buf[ii] = ' ';
// ii++;
// }
// skip space
if (tempbuf[i] == ' ') {
i++;
}
if (tempbuf[i] == '.') {
i++;
buf[ii] = '0'; ii++;
buf[ii] = '.'; ii++;
}
while (tempbuf[i]) {
buf[ii] = tempbuf[i];
i++;ii++;
}
buf[ii] = 0;
} else {
// we have an exponent, get rid of it
i = ii = 0;
if (tempbuf[i] == '-') {
buf[ii] = tempbuf[i];
i++;ii++;
}
// else {
// buf[ii] = ' ';
// ii++;
// }
// skip space
if (tempbuf[i] == ' ') {
i++;
}
ex = ((tempbuf[epos+2] - '0') * 10) + (tempbuf[epos+3] - '0');
if (tempbuf[epos+1] == '+') {
// positive exponent, move decimal point to right, add zeros to the right
// first copy until we see either the decimal point, or the 'e'
while (tempbuf[i] && tempbuf[i] != '.' && tempbuf[i] != 69) {
buf[ii] = tempbuf[i];
i++;ii++;
}
// 'e' found, add as many zeros as in the exponent
if (tempbuf[i] == 69) {
while(ex) {
buf[ii] = '0';
ii++;
ex--;
}
}
if (tempbuf[i] == '.') {
i++; // skip '.'
// copy digits until 'e'
while(ex && tempbuf[i] != 69) {
buf[ii] = tempbuf[i];
ii++;i++;
ex--;
}
// add zeros
while(ex) {
buf[ii] = '0';
ii++;
ex--;
}
}
} else {
// negative exponent, move decimal point to left, add zeros to left
buf[ii] = '0'; ii++;
buf[ii] = '.'; ii++;
ex--;
// add zeros
while(ex) {
buf[ii] = '0';
ii++;
ex--;
}
while (tempbuf[i] && tempbuf[i] != 69) {
if (tempbuf[i] == '.') {
i++;
}
buf[ii] = tempbuf[i];
i++;ii++;
}
}
buf[ii] = 0;
}
return buf;
}
#endif

View File

@ -13,12 +13,12 @@ You can however try the current state of development. You should be able to
build small (and slow...) programs that use floats on any supported target.
- Build the compiler/toolchain/libs from this fptest branch.
- Now you can build the samples and/or tests.
- Now you can build and run the samples and/or tests.
```
samples/floattest.c
samples/mandelfloat.c
samples/mathtest.c (requires full math.h)
samples/tgisincos.c (requires sin/cos from math.h)
samples/mathtest.c
samples/tgisincos.c
```
full math.h is available for C64 when linking agains fp754kernal.o (see below)
@ -31,6 +31,9 @@ full math.h is available for C64 when linking agains fp754kernal.o (see below)
simulator for running test programs, and test changes in the compiler using
our test bench, and that against a library that is known to somewhat work
correctly :)
- The default library also contains a collection of math functions (which i
dubbed "softmath"), which may fill the gap until more targets have specific
wrappers and/or a generic math library (in assembly) was written.
- The default library can be overridden by linking an override file, similar to
how you can use the soft80 implementation for conio. Right now such override
files are provided for the C64 (c64-fp754kernal.o) and VIC20
@ -42,20 +45,17 @@ full math.h is available for C64 when linking agains fp754kernal.o (see below)
Please see below for more info on what is needed to do this, should you be
interested in doing this. Please contact me before you are putting work into
this, so we can discuss some things and prevent anyone wasting time :)
- It might be possible to produce a similar kernal- or OS- wrapper override file
as the C64 one for other targets (or port the C64 one to other CBM targets).
- If you create a new one, keep in mind that the compiler *right now* will
currently work with IEEE754 floats, which your library calls must also work
with (which will involve converting forth and back to whatever other format
at runtime), and there is no easy way tp change that.
- Similar to the softfloat library, it would be nice to have a complete
reference implementation for the math stuff in softmath as well.
- Also the math stuff should be implemented in assembly at some point :)
### Roadmap
- Test/Fix using the Softfloat lib some more, fix as much tests as possible
- Find some more generic math functions that we can use in softmath
- When all obvious tests have been created and work OK, we can merge
- for the failing tests, create "negative" cases
@ -78,8 +78,6 @@ each other (not necessarily by me :)):
- also we must *implement* the native format on the host in fp.c
- if the native format uses a different number of bytes than 4 for its native
format, we must add support for this in the compiler at various places
- add a cmdline option to the compiler to switch the float binary type (754,
cbm, woz, ...)
- last not least a wrapper library that uses the native format must be created
- it is not unlikely that we will need extra tests for the native format
@ -235,7 +233,7 @@ default for all targets.
#### softmath
Contains a collection of math functions, hopefully some day enough to completely
Contains a collection of math functions, (hopefully) enough to completely
implement math.h in C. This is currently used by default for all targets.
#### cbmkernal
@ -294,18 +292,20 @@ These are optional, required for standard libm.
```
func description softmath cbmfp 754
float powf(float f, float a) - * -
float powf(float f, float a) * * -
float sinf(float s) * * -
float cosf(float s) * * -
float logf(float x) - * -
float expf(float x) - * -
float sqrtf(float x) - * -
float tanf(float x) - * -
float atanf(float x) - * -
float logf(float x) * * -
float expf(float x) * * -
float sqrtf(float x) * * -
float tanf(float x) * * -
float atanf(float x) * * -
float fabsf(float x) * * -
float roundf(float x) * * -
float truncf(float x) * * -
float ceilf(float x) * - -
float floorf(float x) * - -
float fmodf(float x, float y) * - -
```
### extra functions

View File

@ -2,11 +2,3 @@
This is a collection of various math related functions found in math.h. All
functions are single precision (only), and are rather selected for their
compactness than speed.
STILL MISSING:
float logf(float x)
float expf(float x)
float sqrtf(float x)
float tanf(float x)
float atanf(float x)

View File

@ -0,0 +1,72 @@
#include <math.h>
#if 1
/* atan(x)= x(c1 + c2*x**2 + c3*x**4)/(c4 + c5*x**2 + c6*x**4 + x**6)
Accurate to about 13.7 decimal digits over the range [0, pi/12]. */
float _atan(float x)
{
const float c1= 48.70107004404898384f;
const float c2= 49.5326263772254345f;
const float c3= 9.40604244231624f;
const float c4= 48.70107004404996166f;
const float c5= 65.7663163908956299f;
const float c6= 21.587934067020262f;
float x2;
x2 = x * x;
return (x * (c1 + x2 * (c2 + x2 * c3)) / (c4 + x2 * (c5 + x2 * (c6 + x2))));
}
#define TAN_SIXTHPI 0.009138776996f
#define TAN_TWELFTHPI 0.004569293096f
float atanf(float x)
{
float y;
int complement= 0; // true if arg was >1
int region= 0; // true depending on region arg is in
int sign= 0; // true if arg was < 0
if (x < 0.0f ) {
// x = -x;
// x = x * -1.0f;
x = -1.0f * x;
sign = 1; // arctan(-x)=-arctan(x)
}
if (x > 1.0f) {
x = 1.0f / x; // keep arg between 0 and 1
complement = 1;
}
if (x > TAN_TWELFTHPI) {
/* FIXME: reduce arg to under tan(pi/12) */
#if 1
float n, m;
n = (TAN_SIXTHPI * x) + 1.0f;
m = (x - TAN_SIXTHPI);
// x = (x - TAN_SIXTHPI) / (1.0f + (TAN_SIXTHPI * x));
x = m / n;
#else
x = fmodf(x, TAN_SIXTHPI);
#endif
region = 1;
}
y = _atan(x);
if (region) { y += (M_PI / 6.0f); } /* correct for region we're in */
if (complement) { y= (M_PI / 2.0f) - y; } /* correct for 1/x */
if (sign) { y =- y; } /* correct for negative arg */
return y;
}
#endif
#if 0
float atanf(float x)
{
x = x;
return 0.0f;
}
#endif

View File

@ -1,6 +1,7 @@
#include <math.h>
/* FIXME: this is really too simple */
float ceilf(float x)
{
int n = (int) x;

View File

@ -61,7 +61,7 @@ float cosf(float x)
}
#if 0
float c(float x) {
float cosf(float x) {
return sinf(x + (M_PI / 2.0));
}
#endif

View File

@ -0,0 +1,73 @@
#include <math.h>
#if 1
/* natural logarithm of 2 */
#define LN2 0.693147180559945309417f
// #include <stdio.h>
// char abuf[100];
// char abuf2[100];
// char abuf3[100];
float expf(float x)
{
float p;
int i;
int k;
// float i;
// float k;
float r;
float tn;
float x0 = fabsf(x); // FIXME: somehow this doesnt work!
if (x == 0) {
return 1;
}
x0 = fabsf(x);
// printf("x:%s x0:%s\n", _ftostr(abuf3, x), _ftostr(abuf, x0));
k = ceilf((x0 / LN2) - 0.5f);
p = (float)(1 << (int)k);
r = x0 - (LN2 * (float)k); // internal error
tn = 1.0f;
// printf("k:%s tn:%s r:%s\n", _ftostr(abuf3, k), _ftostr(abuf, tn), _ftostr(abuf2, r ));
for (i = 14; i > 0; --i) {
//tn = tn * (r / (float)i) + 1.0f;
float tmp;
tmp = (r / (float)i);
tn = tn * tmp;
tn = tn + 1.0f;
// printf("i:%d tn:%s tmp:%s\n", i, _ftostr(abuf, tn), _ftostr(abuf2, tmp ));
}
p *= tn;
if (x < 0) {
return 1.0f / p;
}
return p;
}
#endif
#if 0
static float expf(float n) {
int a = 0, b = n > 0;
float c = 1, d = 1, e = 1;
for (b || (n = -n); e + .00001 < (e += (d *= n) / (c *= ++a));); // "Floating point type is currently unsupported"
// for (b || (n = -n); e + .00001 < (e = e + (d *= n) / (c *= ++a));); // "Floating point type is currently unsupported"
// approximately 15 iterations
return b ? e : 1 / e;
}
#endif
#if 0
float expf(float x)
{
x = x;
return 0.0f;
}
#endif

View File

@ -2,5 +2,5 @@
#include <math.h>
float fabsf(float x) {
return x < 0.0 ? -x : x;
return x < 0.0f ? -x : x;
}

View File

@ -0,0 +1,17 @@
#include <math.h>
/* FIXME: this is really too simple */
float ffloor(float x)
{
signed long n;
float d;
n = (signed long)x;
d = (float)n;
if (x >= 0) {
return d;
}
return d - 1;
}

View File

@ -0,0 +1,9 @@
#include <math.h>
float fmodf(float x, float y)
{
float res;
res = x / y;
return x - (truncf(res) * y);
}

View File

@ -0,0 +1,130 @@
#include <math.h>
/* natural logarithm */
#if 1
/*#define LOGBASE 10*/
#define LOGBASE 2 /* e=2.7182818... */
float logf(float x)
{
int i;
float alpha;
float save;
float ans;
alpha = (x-1)/(x+1);
ans = alpha;
save = ans * alpha * alpha;
for (i = 2 ; i <= LOGBASE ; i++) {
ans += (1.0/(float)(2*i-1)) * save;
save = save * alpha * alpha;
}
return 2.0*ans;
}
#endif
#if 0
static int msb(int a)
{
unsigned int r = 0;
while(a >>= 1) {
r++;
}
return r;
}
float logf(float y)
{
float result;
int log2;
float divisor, x;
log2 = msb((int)y);
divisor = (float)(1 << log2);
x = y / divisor;
result = -1.7417939f + (2.8212026f + (-1.4699568f + (0.44717955f - 0.056570851f * x) * x) * x) * x;
result += ((float)log2) * 0.69314718f; // ln(2) = 0.69314718
return result;
}
#endif
#if 0
float logf(float x)
{
// ASSUMING:
// - non-denormalized numbers i.e. x > 2^126
// - integer is 32 bit. float is IEEE 32 bit.
// INSPIRED BY:
// - https://stackoverflow.com/a/44232045
// - http://mathonweb.com/help_ebook/html/algorithms.htm#ln
// - https://en.wikipedia.org/wiki/Fast_inverse_square_root
// FORMULA:
// x = m * 2^p =>
// ln(x) = ln(m) + ln(2)p,
// first normalize the value to between 1.0 and 2.0
// assuming normalized IEEE float
// sign exp frac
// 0b 0 [00000000] 00000000000000000000000
// value = (-1)^s * M * 2 ^ (exp-127)
//
// exp = 127 for x = 1,
// so 2^(exp-127) is the multiplier
// evil floating point bit level hacking
unsigned long bx = * (unsigned long *) (&x);
// extract exp, since x>0, sign bit must be 0
unsigned long ex = bx >> 23;
signed long t = (signed long)ex-(signed long)127;
unsigned long s = (t < 0) ? (-t) : t;
// reinterpret back to float
// 127 << 23 = 1065353216
// 0b11111111111111111111111 = 8388607
bx = 1065353216 | (bx & 8388607);
x = * (float *) (&bx);
// use remez algorithm to find approximation between [1,2]
// - see this answer https://stackoverflow.com/a/44232045
// - or this usage of C++/boost's remez implementation
// https://computingandrecording.wordpress.com/2017/04/24/
// e.g.
// boost::math::tools::remez_minimax<double> approx(
// [](const double& x) { return log(x); },
// 4, 0, 1, 2, false, 0, 0, 64);
//
// 4th order is:
// { -1.74178, 2.82117, -1.46994, 0.447178, -0.0565717 }
//
// 3rd order is:
// { -1.49278, 2.11263, -0.729104, 0.10969 }
return
/* less accurate */
-1.49278+(2.11263+(-0.729104+0.10969*x)*x)*x
/* OR more accurate */
// -1.7417939+(2.8212026+(-1.4699568+(0.44717955-0.056570851*x)*x)*x)*x
/* compensate for the ln(2)s. ln(2)=0.6931471806 */
+ 0.6931471806*t;
}
#endif
#if 0
float logf(float y)
{
return 0.0f;
}
#endif

View File

@ -1,10 +1,11 @@
#include <math.h>
/* FIXME: this is really too simple */
float roundf(float x)
{
if (x > 0.0f) {
return (float)((signed)(x + 0.5f));
return (float)((signed long)(x + 0.5f));
}
return (float)((signed)(x - 0.5f));
return (float)((signed long)(x - 0.5f));
}

View File

@ -0,0 +1,90 @@
#if 1
static float powerOfTen(int num)
{
int i;
// float rst = 1.0f;
float rst;
rst = 1.0f;
if(num >= 0) {
for(i = 0; i < num ; ++i) {
// rst *= 10.0;
rst = rst * 10.0f;
}
} else {
for(i = 0; i < (0 - num); ++i) {
// rst *= 0.1;
rst = rst * 0.1f;
}
}
return rst;
}
#define MAXDIGITS 8
float sqrtf(float a)
{
// float rst = 0.0f;
float rst;
// float z = a;
float z;
signed int i;
// float j = 1.0f;
float j;
float power;
rst = 0.0f;
z = a;
j = 1.0f;
for(i = MAXDIGITS ; i > 0 ; i--) {
power = powerOfTen(i);
// value must be bigger then 0
if(z - (( 2.0f * rst ) + ( j * power)) * ( j * power) >= 0) {
while( z - (( 2.0f * rst ) + ( j * power)) * ( j * power) >= 0) {
//j++;
j = j + 1.0f;
if(j >= 10.0f) {
break;
}
}
//j--; //correct the extra value by minus one to j
j = j - 1.0f;
//z -= (( 2.0f * rst ) + ( j * power)) * ( j * power); //find value of z
z = z - (( 2.0f * rst ) + ( j * power)) * ( j * power); //find value of z
//rst += j * power; // find sum of a
rst = rst + (j * power); // find sum of a
j = 1.0f;
}
}
for(i = 0 ; i >= 0 - MAXDIGITS ; i--) {
power = powerOfTen(i);
if(z - (( 2.0f * rst ) + ( j * power))*( j * power) >= 0) {
while( z - (( 2.0f * rst ) + ( j * power))*( j * power) >= 0) {
//j++;
j = j + 1.0f;
if(j >= 10.0f) {
break;
}
}
//j--; //correct the extra value by minus one to j
j = j - 1.0f;
//z -= (( 2.0f * rst ) + ( j * power))*( j * power); //find value of z
z = z - (( 2.0f * rst ) + ( j * power)) * ( j * power); //find value of z
//rst += j * power; // find sum of a
rst = rst + (j * power); // find sum of a
j = 1.0f;
}
}
// find the number on each digit
return rst;
}
#endif
#if 0
float sqrtf(float a)
{
a = a;
return 0.0f;
}
#endif

View File

@ -0,0 +1,52 @@
#include <math.h>
#if 1
/* The input argument is in radians. Note that the function
computes tan(pi*x/4), NOT tan(x); it's up to the range
reduction algorithm that calls this to scale things properly.
Algorithm:
tan(x)= x(c1 + c2*x**2 + c3*x**4)/(c4 + c5*x**2 + c6*x**4 + x**6)
Accurate to about 14 decimal digits over the range [0, pi/4]. */
static float _tan(float x)
{
const float c1=-34287.4662577359568109624f;
const float c2= 2566.7175462315050423295f;
const float c3=- 26.5366371951731325438f;
const float c4=-43656.1579281292375769579f;
const float c5= 12244.4839556747426927793f;
const float c6=- 336.611376245464339493f;
float x2;
x2 = x * x;
return (x * (c1 + x2 * (c2 + x2 * c3)) / (c4 + x2 * (c5 + x2 * (c6 + x2))));
}
float tanf(float x){
int octant;
x = fmodf(x, (2.0f * M_PI));
octant=(int)(x * (4.0f / M_PI));
switch (octant){
case 0: return _tan(x * (4.0f / M_PI));
case 1: return 1.0f / _tan(((M_PI / 2.0f) - x) * (4.0f / M_PI));
case 2: return -1.0f / _tan((x - (M_PI / 2.0f)) * (4.0f / M_PI));
case 3: return - _tan((M_PI - x) * (4.0f / M_PI));
case 4: return _tan((x - M_PI) * (4.0f / M_PI));
case 5: return 1.0f / _tan(((3.0f * M_PI / 2.0f) - x) * (4.0f / M_PI));
case 6: return -1.0f / _tan((x - (3.0f * M_PI / 2.0f)) * (4.0f / M_PI));
case 7: return - _tan(((2.0f * M_PI)-x) * (4.0f / M_PI));
}
}
#endif
#if 0
float tanf(float x)
{
x = x;
return 0.0f;
}
#endif

View File

@ -1,7 +1,8 @@
#include <math.h>
/* FIXME: this is really too simple */
float truncf(float x)
{
return (float)((signed)x);
return (float)((signed long)x);
}

View File

@ -1,7 +1,13 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef __CC65__
#include <conio.h> // for cgetc
#else
#define cgetc()
char *_ftostr(char *d, float s);
#endif
#ifdef __SIM6502__
#define cgetc()