Sharp PC-G850VS

[robve]

(02-22-2021 10:29 PM)rprosperi Wrote:  This is quite surprising, but thanks for insights on the additional C function support for the built-in BASIC functions. I recall he told me there was a certain entire class of functions that were not provided, which I took to mean something like 'all the trig Fns', 'all the Stats Fns', etc. but perhaps he merely meant that basic functions such as you've provided in your 2nd post, which can be easily assembled from available core functions, were not provided? Also, ROM space limits was a factor (it always is!) so best to spend ROM space on the much more difficult transcendental functions.

Very good points!

I should add that most students learning C typically do not need any of the advanced math functions. They spend a lot of time on other things, like getting to grips with the concept of C pointers to build data structures (or depending what PL you like best).

It might be useful to reproduce most of the math.h functions for the PC-G850(V)(S), even if it is just an exercise. This is a professional calculator forum after all

The C math functions are carefully crafted to prevent catastrophic cancellation, roundoff, and overflow issues. These issues are not for beginners. Entire books have been written on this subject.

Here is my attempt at reproducing the math.h functions and constants for the PC-G850(V)(S). EDIT: added new functions atan2, fmod, erf, erfc, tgamma, and lgamma with everything tested and verified:

/* math.h for the PC-G850(V)(S) - Dr. Robert van Engelen */
/* Written in compact form to save space */
/* WARNING: <,<=,>,>= on floats have a bug on some G850's, instead of x<y use x-y<0 */
/* Note: erf and erfc precision is at least 1.2e-7 (erfcc) */
/* Note: lgamma and tgamma precision is at least 2e-10 (Lanczos) */
/* License: BSD-3 open source */

#define M_E 2.718281828
#define M_LOG2E 1.366512921
#define M_LOG10E .4342944819
#define M_LN2 .693147181
#define M_LN10 2.302585093
#define M_PI 3.141592654
#define M_PI_2 1.570796327
#define M_PI_4 .7853981634
#define M_1_PI .3183098862
#define M_2_PI .6366197724
#define M_2_SQRTPI 1.128379167
#define M_SQRT2 1.414213562
#define M_SQRT1_2 .7071067812

double fabs(double x){return x<0?-x:x;}
double copysign(double x,double y){return y<0?-fabs(x):fabs(x);}
double log2(double x){return log(x)/M_LN2;}
double exp2(double x){return exp(M_LN2*x);}
double cbrt(double x){return pow(x,1./3);}
double floor(double x){return fabs(x)-1e9<0?x<0?(long)(x-.9999999999):(long)x:x;}
double ceil(double x){return -floor(-x);}
double trunc(double x){return fabs(x)-1e9<0?(long)x:x;}
double round(double x){return floor(x+.5);}
double modf(double x,double*y){return x-(*y=fabs(x)-1e9<0?x-(long)x:0);}
double fdim(double x,double y){return x-y>0?x-y:0;}
double fmax(double x,double y){return x-y>0?x:y;}
double fmin(double x,double y){return x-y<0?x:y;}
double hypot(double x,double y){double a=fabs(x),b=fabs(y);return a?b?a-b<0?b*sqrt(1+x/y*x/y):a*sqrt(1+y/x*y/x):a:b;}
double atan2(double y,double x){return x>0?atan(y/x):y?asin(y<0?-1:1)-atan(x/y):acos(-(x<0));}
double fmod(double x,double y){return x-y*trunc(x/y);}

double erfc(double x){double p=1/(1+fabs(x)/2),q=-.18628806+p*(.27886807+p*(-1.13520398+p*(1.4885187+p*(-.82215223+p*.17087277))));p*=exp(-x*x-1.26551223+p*(1.00002368+p*(.37409196+p*(.09678418+p*q))));return x<0?2-p: p;}
double erf(double x){return 1-erfc(x);}
double lgamma(double x){return log((1+76.18009173/(x+1)-86.50532033/(x+2)+24.01409824/(x+3)-1.231739572/(x+4)+1.208650974e-3/(x+5)-5.395239385e-6/(x+6))*2.506628275/x)+(x+.5)*log(x+5.5)-x-5.5;}
double tgamma(double x){return exp(lgamma(x));}

Note that hypot avoids catastrophic overflow by rearranging terms. Macros are used only if the argument(s) is/are guaranteed to be evaluated just once in the macro body. This avoids the common problem with calls like fmax(x++,y++) that fail as a macro. Otherwise, I tried to be as efficient as possible. Unlike the BASIC INT() function posted earlier, the floor(), ceil(), trunc(), round(), and modf() all work on unlimited doubles by using the fact that a long integer can hold up to 9 digits and G850(V)(S) doubles have 10 digit mantissas. This means that values >=1e9 have no fractional part. I used this observation to trick the code to check for |x|<1e9 which means the value has a fraction and can use a long cast in that case to obtain the integer part, otherwise just return the value unchanged.

If you spot an issue or an opportunity for further optimization then let us know!

- Rob

EDIT: replaced #defines with functions and added warning.