summaryrefslogtreecommitdiffstats
path: root/mdk-stage1/dietlibc/libm
diff options
context:
space:
mode:
Diffstat (limited to 'mdk-stage1/dietlibc/libm')
-rw-r--r--mdk-stage1/dietlibc/libm/acosh.c6
-rw-r--r--mdk-stage1/dietlibc/libm/asinh.c6
-rw-r--r--mdk-stage1/dietlibc/libm/atanh.c8
-rw-r--r--mdk-stage1/dietlibc/libm/bessel.c171
-rw-r--r--mdk-stage1/dietlibc/libm/cosh.c9
-rw-r--r--mdk-stage1/dietlibc/libm/erf.c95
-rw-r--r--mdk-stage1/dietlibc/libm/gamma.c98
-rw-r--r--mdk-stage1/dietlibc/libm/ipow.c29
-rw-r--r--mdk-stage1/dietlibc/libm/poly.c41
-rw-r--r--mdk-stage1/dietlibc/libm/pow.c42
-rw-r--r--mdk-stage1/dietlibc/libm/rint.c5
-rw-r--r--mdk-stage1/dietlibc/libm/sinh.c9
-rw-r--r--mdk-stage1/dietlibc/libm/tanh.c7
13 files changed, 526 insertions, 0 deletions
diff --git a/mdk-stage1/dietlibc/libm/acosh.c b/mdk-stage1/dietlibc/libm/acosh.c
new file mode 100644
index 000000000..a09a4c9f0
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/acosh.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+double acosh ( double x )
+{
+ return log ( x + sqrt (x*x - 1.) );
+}
diff --git a/mdk-stage1/dietlibc/libm/asinh.c b/mdk-stage1/dietlibc/libm/asinh.c
new file mode 100644
index 000000000..49c6b467f
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/asinh.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+double asinh ( double x )
+{
+ return log ( x + sqrt (x*x + 1.) );
+}
diff --git a/mdk-stage1/dietlibc/libm/atanh.c b/mdk-stage1/dietlibc/libm/atanh.c
new file mode 100644
index 000000000..bdb3367be
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/atanh.c
@@ -0,0 +1,8 @@
+#include <math.h>
+
+extern const float __half;
+
+double atanh ( double x )
+{
+ return __half * log ( (1.+x) / (1.-x) );
+}
diff --git a/mdk-stage1/dietlibc/libm/bessel.c b/mdk-stage1/dietlibc/libm/bessel.c
new file mode 100644
index 000000000..ba8a1afcb
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/bessel.c
@@ -0,0 +1,171 @@
+/*--------------------------------------------------------------------------*
+
+Name j0, j1, jn - Bessel functions
+ y0, y1, yn - Weber functions
+
+Usage double j0 (double x);
+ double j1 (double x);
+ double jn (int n, double x);
+ double y0 (double x);
+ double y1 (double x);
+ double yn (int n, double x);
+
+Prototype in math.h
+
+Description j0, j1 and jn calculate the Bessel function.
+ y0, y1 and yn calcualte the Weber function.
+
+Return value return their return values as doubles.
+
+*---------------------------------------------------------------------------*/
+
+#include <math.h>
+
+#define M_C 0.5772156649015328
+#if 0
+#define M_1_PI 0.318309886183790671538
+#define M_2_PI 0.636619772367581343076
+#define M_PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148
+#endif
+
+
+#define EXPL(x) ((((short *)&x)[4] & 0x7FFF) >> 0)
+#define EXPD(x) ((((short *)&x)[3] & 0x7FF0) >> 4)
+#define EXPF(x) ((((short *)&x)[1] & 0x7F80) >> 7)
+
+#define SQUARE(x) (long) (My - (x) * (x) )
+
+
+static long double P ( int My, double* x )
+{
+ long double Sum = 0.;
+ long double Fact = 1.;
+ long double z182 = -0.015625 / (x[0] * x[0]);
+ register int i;
+
+ for ( i = 1; ; i += 2 ) {
+ Fact *= SQUARE(i+i-1) * SQUARE(i+i+1) * z182 / (i*(i+1));
+ if ( EXPL (Fact) < 0x3FFF-53 )
+ break;
+ Sum += Fact;
+ }
+ return 1. + Sum;
+}
+
+static long double Q ( int My, double* x )
+{
+ long double Fact = (My-1) / x[0] * 0.125;
+ long double Sum = Fact;
+ long double z182 = -0.015625 / (x[0]*x[0]);
+ register int i;
+
+ for ( i = 2; ; i += 2 ) {
+ Fact *= SQUARE(i+i-1) * SQUARE(i+i+1) * z182 / (i*(i+1));
+ if ( EXPL (Fact) < 0x3FFF-53 )
+ break;
+ Sum += Fact;
+ }
+ return Sum;
+}
+
+
+static long double ___jn ( int n, double* x )
+{
+ long double Sum;
+ long double Fact;
+ long double y;
+ register int i;
+ double xx;
+ long double Xi;
+ int My;
+
+ if ( n < 0 )
+ return n & 1 ? ___jn (-n, x) : -___jn (-n, x);
+
+ if ((x[0] >= 17.7+0.0144*(n*n))) {
+ Xi = x[0] - M_PI * (n*0.5 + 0.25);
+ My = n*n << 2;
+
+ return sqrt ( M_2_PI/x[0] ) * ( P(My,x) * cos(Xi) - Q(My,x) * sin(Xi) );
+ }
+ xx = x[0] * 0.5;
+ Sum = 0.;
+ Fact = 1.;
+ y = -xx * xx;
+
+ for ( i = 1; i <= n; i++ )
+ Fact *= xx/i;
+ for ( i = 1; ; i++ ) {
+ Sum += Fact;
+ Fact *= y / (i*(n+i));
+ if ( EXPL (Sum) - EXPL(Fact) > 53 || !EXPL(Fact) )
+ break;
+ }
+ return Sum;
+}
+
+
+static long double ___yn ( int n, double* x )
+{
+ long double Sum1;
+ long double Sum2;
+ long double Fact1;
+ long double Fact2;
+ long double F1;
+ long double F2;
+ long double y;
+ register int i;
+ double xx;
+ long double Xi;
+ unsigned int My;
+
+ if ( EXPD (x[0]) == 0 )
+ return -1./0.; /* ignore the gcc warning, this is intentional */
+
+ if ( (x[0] >= (n>=32 ? 25.8 : (n<8 ? 17.4+0.1*n : 16.2+0.3*n))) ) {
+ Xi = x[0] - M_PI * (n*0.5+0.25);
+ My = n*n << 2;
+
+ return sqrt ( M_2_PI / x[0] ) * ( P(My,x) * sin(Xi) + Q(My,x) * cos(Xi) );
+ }
+
+ Sum1 = Sum2 = F1 = F2 = 0;
+ Fact1 = 1. / (xx = x[0] * 0.5 );
+ Fact2 = 1.;
+ y = xx*xx;
+
+ for ( i = 1; i < n; i++ )
+ Fact1 *= (n-i) / xx;
+
+ for ( i = 1; i <= n; i++ ) {
+ Sum1 += Fact1;
+ if ( i == n )
+ break;
+ Fact1 *= y/(i*(n-i));
+ }
+
+ for (i=1; i<=n; i++) {
+ Fact2 *= xx / i;
+ F1 += 1. / i;
+ }
+
+ for ( i = 1; ; i++ ) {
+ Sum2 += Fact2 * (F1+F2);
+ Fact2 *= -y / (i*(n+i));
+ if ( EXPL (Sum2) - EXPL (Fact2) > 53 || !EXPL (Fact2) )
+ break;
+ F1 += 1. / (n+i);
+ F2 += 1. / i;
+ }
+
+ return M_1_PI * (2. * (M_C + log(xx)) * ___jn (n, x) - Sum1 - Sum2);
+}
+
+
+double j0 ( double x ) { return ___jn ( 0,&x ); }
+double j1 ( double x ) { return ___jn ( 1,&x ); }
+double jn ( int n, double x ) { return ___jn ( n,&x ); }
+double y0 ( double x ) { return ___yn ( 0,&x ); }
+double y1 ( double x ) { return ___yn ( 1,&x ); }
+double yn ( int n, double x ) { return ___yn ( n,&x ); }
+
diff --git a/mdk-stage1/dietlibc/libm/cosh.c b/mdk-stage1/dietlibc/libm/cosh.c
new file mode 100644
index 000000000..f64d59106
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/cosh.c
@@ -0,0 +1,9 @@
+#include <math.h>
+
+extern const float __half;
+
+double cosh ( double x )
+{
+ long double y = exp (x);
+ return (y + 1./y) * __half;
+}
diff --git a/mdk-stage1/dietlibc/libm/erf.c b/mdk-stage1/dietlibc/libm/erf.c
new file mode 100644
index 000000000..63f52d81f
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/erf.c
@@ -0,0 +1,95 @@
+#include "dietlibm.h"
+
+/*--------------------------------------------------------------------------*
+ z
+ 1 | -x²/2
+Name erf(z) = --------- | e dx
+ sqrt(2pi) |
+ 0
+
+ oo
+ 1 | -x²/2
+ erfc(z) = -------- | e dx
+ sqrt(2pi) |
+ z
+
+Usage double erf (double x);
+ double erfc(double x);
+
+Prototype in math.h
+
+Description erf(x) is the probability a normal distributed event occures
+ within the range [0,x]. erfc(x) is the probability a normal
+ distributed event occures within the range [x,oo].
+
+Return value return their respective function value.
+
+*---------------------------------------------------------------------------*/
+
+
+/* even function in (0): Coefficients for gamma(0) */
+
+static const double tab1 [9 + 1] = {
+ 0.398942280401432677926, -0.066490380066905446321, 9.97355701003581694794E-3, -1.18732821548045439878E-3, 1.15434687616155288764E-4, -9.44465625950361453450E-6, 6.65969351631665127484E-7, -4.12266741486268888409E-8, 2.27352982437280636972E-9, -1.13011716416192129505E-10
+};
+
+/* non even or odd function in (x), x>0: Coefficients for gamma(x), x>0 */
+
+static const double tab2 [] [31 + 1] = {
+ { -0.158655253931457051468, +0.241970724519143349823, -0.120985362259571674911, 0, +0.0201642270432619458197, -4.03284540865238916394E-3, -2.01642270432619458197E-3, +7.68161030219502697887E-4, +1.20025160971797296538E-4, -8.80184513793180174807E-5, -1.86705805956129127862E-6, +7.37124220917704609315E-6, -4.72826391707080259142E-7, -4.83395817951682973566E-7, +6.57036391970156141055E-8, +2.5544260402922190768E-8, -5.4292285616752144141E-9, -1.08932444506260820153E-9, +3.44399256708718202774E-10, +3.6021429664641554881E-11, -1.81147204852239925966E-11, -7.66935128389784976374E-13, +8.19047721646461768154E-13, -3.78144699611990981391E-15, -3.24856460059989147863E-14, +1.44438130842455313227E-15, +1.14391687912824634892E-15, -9.38053726039148625184E-17, -3.59908648108845288945E-17, +4.36020846676166022246E-18, +1.01298640134330880603E-18, -1.68640470512244526894E-19 },
+ { -0.0227501319481792072104, +0.0539909665131880519553, -0.0539909665131880519553, +0.0269954832565940259776, -4.49924720943233766301E-3, -2.24962360471616883129E-3, +1.34977416282970129877E-3, -1.17837426913704081544E-4, -1.15159303574756261652E-4, +3.70473728554448438507E-5, +2.82690796888936559912E-6, -3.54513195524355369855E-6, +3.76695631261094890352E-7, +1.92024079214184701051E-7, -5.22690859049557191018E-8, -4.91799344974114749666E-9, +3.66377919234006038965E-9, -1.5981997209104676352E-10, -1.73812379171063320997E-10, +2.62403075313043113473E-11, +5.60918720760414713346E-12, -1.72126983666416144614E-12, -8.63428809787622525331E-14, +7.89441765474563834480E-14, -3.13747960081562321348E-15, -2.77519506625391157547E-15, +3.29321944203493138076E-16, +7.44375150395529134369E-17, -1.66428523299294690222E-17, -1.32735612757620496568E-18, +6.24122437514304644794E-19, +1.12471123532438919306E-21 },
+ { -1.3498980316300945272E-3, +4.43184841193800717687E-3, -6.64777261790701076574E-3, +5.90913121591734290293E-3, -3.32388630895350538287E-3, +1.10796210298450179421E-3, -1.10796210298450179595E-4, -8.44161602273906129349E-5, +4.35270826172482847927E-5, -6.30190085030867423515E-6, -1.9785037553294674925E-6, +1.05520200284238266374E-6, -1.13913852579575399458E-7, -4.81174572974454799623E-8, +1.78216871733806513653E-8, -5.85637697215219690327E-10, -9.29791350219350980904E-10, +1.96377023046901260016E-10, +1.58870373467897094393E-11, -1.22699105512396660364E-11, +1.08794270836433192571E-12, +3.99646995170699427940E-13, -1.01594404465456044793E-13, -3.33469605506835759271E-15, +4.46588935876766499879E-15, -4.08076707607833277747E-16, -1.17808602368979218862E-16, +2.76224909899945482352E-17, +1.09206599392049874162E-18, -1.03145418746203977253E-18, +6.79984672177279963209E-20, +2.55831283729070534712E-20 },
+ { -3.16712418331199212695E-5, +1.33830225764885351832E-4, -2.67660451529770703664E-4, +3.34575564412213379613E-4, -2.89965489157251595673E-4, +1.8178605666396926958E-4, -8.25286392216793003064E-5, +2.55180251904870680833E-5, -3.91665839292075186649E-6, -7.40182052221464123606E-7, +6.44220233592652481453E-7, -1.73701553397390201613E-7, +9.09595464817154590424E-9, +9.44943118114780783705E-9, -3.29957075383376125942E-9, +2.94920746951281580686E-10, +1.18744773902482360274E-10, -4.42039585809856402486E-11, +3.61422484008923382324E-12, +1.43638335494248833511E-12, -4.58476794992724591068E-13, +2.23496663226445199624E-14, +1.57839046076890756440E-14, -3.67258220998453293248E-15, -1.69716269032291432153E-17, +1.43497778353923791279E-16, -2.14499365995613073838E-17, -1.93255135682867953692E-18, +1.01377499752128183701E-18, -7.55713215369572830154E-20, -2.25510650946079103289E-20, +5.26633993110171917109E-21 },
+ { -2.86651571879193912033E-7, +1.48671951473429770924E-6, -3.7167987868357442731E-6, +5.9468780589371908374E-6, -6.81413110919886450076E-6, +5.92209940035828587496E-6, -4.02653201907205629582E-6, +2.17108246596119665457E-6, -9.25512396325170449452E-7, +3.03096091545533908077E-7, -6.92802772105295808398E-8, +6.69226396924248971087E-9, +2.46006252876483997508E-9, -1.41806830376639605249E-9, +3.44251040657349801884E-10, -2.6965166176434937652E-11, -1.16546962748761528049E-11, +4.91490145086991326748E-12, -7.55854519365765424197E-13, -4.53988828124843593484E-14, +4.71533558309731405623E-14, -9.17323049919073092370E-15, +4.35542982587998484108E-17, +3.71238868922011013332E-16, -7.90772907386322623053E-17, +1.58463483904927528072E-18, +2.61503941976309571331E-18, -5.40699423853895351239E-19, +6.61825040533797444037E-21, +1.68378440730394776550E-20, -3.01930850797704474581E-21, -3.80658085177617928332E-23 },
+ { -9.8658764503769814198E-10, +6.07588284982328549581E-9, -1.82276485494698564874E-8, +3.54426499573024987263E-8, -5.01260335110421053478E-8, +5.48348427196551516061E-8, -4.81513715848495375522E-8, +3.47446467489597046263E-8, -2.08994095347716137282E-8, +1.0554987922587771203E-8, -4.4752674615729637229E-9, +1.57746505810079893253E-9, -4.49697115294871911476E-10, +9.63210042443717269402E-11, -1.16300711402336909847E-11, -1.31070037808191623761E-12, +1.16993345829435057496E-12, -3.40636420312606285351E-13, +5.23724821541706939045E-14, +3.93541148139975862961E-16, -2.59886413069218394637E-15, +7.24729556829529838503E-16, -8.51485747763574768020E-17, -7.86503719948806184368E-18, +5.35986191777031053618E-18, -9.84873767617830925356E-19, +2.93759678710573738811E-20, +2.85458592629073152182E-20, -7.12725445137377009753E-21, +5.25419393758902871947E-22, +1.24299023131490990316E-22, -4.04419210566489645405E-23 },
+ { -1.27981254388583500631E-12, +9.1347204083645933588E-12, -3.19715214292760767584E-11, +7.30777632669167468738E-11, -1.22557498812224960902E-10, +1.60618833847077433236E-10, -1.71047639646627010648E-10, +1.51926349902927316213E-10, -1.14609023345779936276E-10, +7.43697341394886835864E-11, -4.18713451557949730558E-11, +2.05606050331840905587E-11, -8.82161466664564577599E-12, +3.30031395277698236679E-12, -1.06851205331295409813E-12, +2.94333808755089195146E-13, -6.64411715537625335642E-14, +1.11264855981436243262E-14, -8.52918435682649455145E-16, -2.38837813662069487819E-16, +1.23994634366691956599E-16, -3.05269770279941723219E-17, +4.34539596489459676621E-18, -5.55819387468189608390E-20, -1.56974672263484202926E-19, +4.60835492190702561464E-20, -6.61112150617493330405E-21, +7.28424268476803924831E-23, +2.09156005934313228089E-22, -5.29080328670107625978E-23, +5.61375000671507211726E-24, +3.82199410465700894394E-25 },
+ { -6.22096057427178413283E-16, +5.05227108353689229741E-15, -2.02090843341475691883E-14, +5.30488463771373691202E-14, -1.02729512031916810045E-13, +1.56409892294496290711E-13, -1.94849254788406146283E-13, +2.04064637342166989709E-13, -1.83187931471980616892E-13, +1.42994099344605424348E-13, -9.8111907789286062426E-14, +5.96545975367403288587E-14, -3.23370114040930933005E-14, +1.56932853967230342257E-14, -6.83548101324218922896E-15, +2.67410077774155118457E-15, -9.38313996431647887562E-16, +2.94090734842381109313E-16, -8.16448235152204729921E-17, +1.9758222496699617607E-17, -4.03590262164308783690E-18, +6.43662361965717426956E-19, -5.93446415094778572090E-20, -6.07164564350191039536E-21, +4.38906686886388095825E-21, -1.17175498170220204828E-21, +1.98482140750318604418E-22, -1.70803571702439545981E-23, -1.94600332107885234554E-24, +1.10477141319981582738E-24, -2.31975718243847439962E-25, +2.54148402104633283670E-26 },
+ { -1.12858840595384064928E-19, +1.02797735716689148111E-18, -4.62589810725101166456E-18, +1.37063647622252197466E-17, -3.0068337697131575822E-17, +5.2067053140503053517E-17, -7.40914680178037035E-17, +8.9062000172830588611E-17, -9.22563786210983011008E-17, +8.35975730487397716492E-17, -6.70372487553237232779E-17, +4.80088566412770650047E-17, -3.09280630297969106245E-17, +1.8026496052333452774E-17, -9.54924880090907168481E-18, +4.61362333444861021959E-18, -2.03812361224098073479E-18, +8.24578860830779678155E-19, -3.0572087552697254564E-19, +1.03827313453936543577E-19, -3.22407758977306397999E-20, +9.12052549039695437376E-21, -2.33541947993595580264E-21, +5.35339963891271164659E-22, -1.07674173853083520575E-22, +1.82413373046113374293E-23, -2.33864726317468746329E-24, +1.29928813344150027051E-25, +3.86668349205203745336E-26, -1.63203452712600670685E-26, +3.65165372186699607411E-27, -5.51243539825332137371E-28 },
+ { -7.61985302416052609616E-24, +7.69459862670641937159E-23, -3.84729931335320968601E-22, +1.26960877340655919637E-21, -3.10990027829384449637E-21, +6.02935924057670511377E-21, -9.6342786971886625897E-21, +1.30454744197246721374E-20, -1.52745988785284834672E-20, +1.57034665186695273938E-20, -1.43457243961336621961E-20, +1.17567385540485497556E-20, -8.7104848256363928121E-21, +5.87137214731944288587E-21, -3.61951956727412561213E-21, +2.04954715001535632502E-21, -1.06982832733527370879E-21, +5.1628428354196120786E-22, -2.30885865897937993512E-22, +9.58556229281154921137E-23, -3.69911125531027884646E-23, +1.32784897023484841369E-23, -4.43433027366044567275E-24, +1.37688611947822111040E-24, -3.96971995397574368025E-25, +1.06008163579031271153E-25, -2.61149430849477426613E-26, +5.89698164189548613154E-27, -1.20793190886658723050E-27, +2.20446342551066852143E-28, -3.46061447029252398335E-29, +4.28913922246949096952E-30 }
+};
+
+static const double tab3 [8] = { +1, -1, +3, -15, +105, -945, +10395, -135135.0 };
+
+
+/*
+ Calculated: oo
+ 1 | -x²/2
+ gauss(z) = --------- | e dx
+ sqrt(2pi) |
+ z
+
+ gauss ( 0) = 0.5
+ gauss ( 1) ~ 0.1586
+ gauss ( 2) ~ 0.02275
+ gauss ( 4) ~ 3.17e-5
+ gauss (10) ~ 7.62e-24
+ gauss (oo) = 0
+
+ Note: only for z>0
+*/
+#include <stdio.h>
+#include <math.h>
+
+#define M_1_SQRT2PI 0.398942280401432686
+
+static long double gauss ( double x )
+{
+ unsigned int i = (unsigned int)(x + 0.5);
+ double y = x * x;
+
+ if ( i > 150 ) return 0.;
+ if ( i > 10 ) return M_1_SQRT2PI * exp (-0.5*y) / x * __poly (1./y, 7, tab3);
+ if ( i > 0 ) return -__poly ((x-i), 31, tab2 [i-1]);
+ return 0.5 - x * __poly (y, 9, tab1);
+ }
+
+double erf ( double x )
+{
+ return x < 0. ? -0.5 + gauss(-x) : 0.5 - gauss(x);
+}
+
+double erfc ( double x )
+{
+ return x < 0. ? 1.0 - gauss(-x) : gauss(x);
+}
+
diff --git a/mdk-stage1/dietlibc/libm/gamma.c b/mdk-stage1/dietlibc/libm/gamma.c
new file mode 100644
index 000000000..d5f3e4275
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/gamma.c
@@ -0,0 +1,98 @@
+#include "dietlibm.h"
+
+/*--------------------------------------------------------------------------*
+
+Name gamma, lgamma - gamma function
+
+Usage double gamma (double x);
+ double lgamma(double x);
+ extern int signgam;
+
+Prototype in math.h
+
+Description gamma returns the logarithm of the absolute value of the
+ gamma function. So it is possible â(x) for very large x.
+ The sign is stored in signgam, a extern variable
+ overwritten during every call to gamma(). lgamma() is
+ a synonym for gamma().
+ You can calculate â(x) by the following sequence:
+
+ double gammafunction(double x)
+ { double y=exp(gamma(x));
+
+ return signgam ? -y : +y;
+ }
+
+Return value gamma returns a value in range (-0.1208, +oo). For a input
+ value of zero, it returns +oo and errno is set to:
+
+ ERANGE Result out of range
+
+*---------------------------------------------------------------------------*/
+
+#include <stdlib.h>
+#include <math.h>
+
+#define B0 + 1.0l/ 6/ 1/ 2
+#define B1 - 1.0l/ 30/ 3/ 4
+#define B2 + 1.0l/ 42/ 5/ 6
+#define B3 - 1.0l/ 30/ 7/ 8
+#define B4 + 5.0l/ 66/ 9/10
+#define B5 - 691.0l/2730/11/12
+#define B6 + 7.0l/ 6/13/14
+#define B7 - 3617.0l/ 510/15/16
+#define B8 + 43867.0l/ 798/17/18
+#define B9 - 174611.0l/ 330/19/20
+#define B10 + 854513.0l/ 138/21/22
+#define B11 - 236364091.0l/2730/23/24
+#define B12 + 8553103.0l/ 6/25/26
+
+static const double coeff[] = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, B10 };
+int signgam;
+
+#define EXPL(x) (((short *)&x)[4] & 0x7FFF)
+
+static double logfact ( long double x )
+{
+ long double z = 2. * M_PI * x;
+ register int e = EXPL (x);
+
+ static unsigned char list [] = { 6, 4, 3, 3, 2, 2 };
+
+ return (log(x) - 1) * x + 0.5*log(z) + __poly (1./(x*x), e<0x4003 ? 10 : (e>0x4008 ? 1 : list [e-0x4003] ), coeff) / x;
+}
+
+
+double lgamma ( double x )
+{
+ register int k = floor (x);
+ long double w;
+ long double y;
+ long double z;
+
+ signgam = 0;
+
+ if ( k >= 7 )
+ return logfact (x-1);
+
+ if ( k == x )
+ switch (k) {
+ case 1 :
+ case 2 : return 0.000000000000000000000000000l;
+ case 3 : return 0.693147180559945309432805516l;
+ case 4 : return 1.791759469228055000858148560l;
+ case 5 : return 3.178053830347945619723759592l;
+ case 6 : return 4.787491742782045994244981560l;
+ default: return 1./0.; /* ignore the gcc warning, this is intentional */
+ }
+
+ z = logfact (y = x - k + 7.0 - 1);
+ w = 1;
+ for ( k = 7 - k; k--; )
+ w *= y, y -= 1.;
+
+ signgam = k >= 0 ? 0 : k & 1;
+ return z - log (w);
+}
+
+double gamma ( double val ) __attribute__ ((weak,alias("lgamma")));
diff --git a/mdk-stage1/dietlibc/libm/ipow.c b/mdk-stage1/dietlibc/libm/ipow.c
new file mode 100644
index 000000000..399986ea1
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/ipow.c
@@ -0,0 +1,29 @@
+#define _GNU_SOURCE
+#include <math.h>
+/*
+ * This is not standard, but often you only need such this function
+ * which is much shorter than the generic pow() function.
+ *
+ * double ipow ( double mant, int expo );
+ */
+
+double ipow ( double mant, int expo )
+{
+ double ret = 1.;
+ unsigned int e = expo; /* Some attention is necessary for expo = 2^31 */
+
+ if ( (int)e < 0 ) {
+ e = -e;
+ mant = 1./mant;
+ }
+
+ while (1) {
+ if ( e & 1 )
+ ret *= mant;
+ if ( (e >>= 1) == 0 )
+ break;
+ mant *= mant;
+ }
+
+ return ret;
+}
diff --git a/mdk-stage1/dietlibc/libm/poly.c b/mdk-stage1/dietlibc/libm/poly.c
new file mode 100644
index 000000000..cdcfb8c5c
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/poly.c
@@ -0,0 +1,41 @@
+/*--------------------------------------------------------------------------*
+
+Name __poly - generates a polynomial from arguments
+
+Usage double __poly ( double x, int n, const double* c );
+
+Prototype in math.h
+
+Description __poly generates a polynomial in x, of degree n, with
+ coefficients c[0], c[1], ..., c[n]. For example, if n=4,
+ the generated polynomial is
+
+ c[4]*x^4 + c[3]*x^3 + c[2]*x^2 + c[1]*x + c[0]
+
+ The polynomial is calculated using Horner's method:
+
+ polynom = (..((x*c[n] + c[n-1])*x + c[n-2])..)*x + c[0]
+
+Return value __poly returns the value of the polynomial as evaluated for
+ the given x.
+ A range error occurs if the result exceeds double range.
+
+*---------------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include "dietlibm.h"
+
+double __poly ( double x, size_t n, const double* c)
+{
+ long double ret;
+ size_t i;
+
+ i = n;
+ c += n;
+ ret = 0;
+ do
+ ret = ret * x + *c--;
+ while ( i-- );
+
+ return ret;
+}
diff --git a/mdk-stage1/dietlibc/libm/pow.c b/mdk-stage1/dietlibc/libm/pow.c
new file mode 100644
index 000000000..e0e5a2983
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/pow.c
@@ -0,0 +1,42 @@
+
+#include <math.h>
+#include "dietlibm.h"
+
+double pow ( double mant, double expo )
+{
+ unsigned int e;
+ long double ret;
+
+ /* special cases 0^x */
+ if ( mant == 0. ) {
+ if ( expo > 0. )
+ return 0.;
+ else if ( expo == 0. )
+ return 1.;
+ else
+ return 1./mant;
+ }
+
+ /* special cases x^n with n is integer */
+ if ( expo == (int) (e = (int) expo) ) {
+
+ if ( (int)e < 0 ) {
+ e = -e;
+ mant = 1./mant;
+ }
+
+ ret = 1.;
+
+ while (1) {
+ if ( e & 1 )
+ ret *= mant;
+ if ( (e >>= 1) == 0 )
+ break;
+ mant *= mant;
+ }
+ return ret;
+ }
+
+ /* normal case */
+ return exp2 ( log2 (mant) * expo );
+}
diff --git a/mdk-stage1/dietlibc/libm/rint.c b/mdk-stage1/dietlibc/libm/rint.c
new file mode 100644
index 000000000..b6f0d8f85
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/rint.c
@@ -0,0 +1,5 @@
+#include <math.h>
+
+double rint(double x) {
+ return floor(x+0.5);
+}
diff --git a/mdk-stage1/dietlibc/libm/sinh.c b/mdk-stage1/dietlibc/libm/sinh.c
new file mode 100644
index 000000000..ae4542d25
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/sinh.c
@@ -0,0 +1,9 @@
+#include <math.h>
+
+extern const float __half;
+
+double sinh ( double x )
+{
+ long double y = exp (x);
+ return (y - 1./y) * __half;
+}
diff --git a/mdk-stage1/dietlibc/libm/tanh.c b/mdk-stage1/dietlibc/libm/tanh.c
new file mode 100644
index 000000000..21dc3d0c0
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/tanh.c
@@ -0,0 +1,7 @@
+#include <math.h>
+
+double tanh ( double x )
+{
+ long double y = exp (x + x);
+ return (y - 1.) / (y + 1.);
+}
generated by cgit v1.2.1 (git 2.21.0) at 2024-06-01 01:02:43 +0000