From 4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gwenol=C3=A9=20Beauchesne?= Date: Wed, 4 Jun 2003 18:44:09 +0000 Subject: Import dietlibc 0.22 + other fixes for AMD64 --- mdk-stage1/dietlibc/libm/acosh.c | 6 ++ mdk-stage1/dietlibc/libm/asinh.c | 6 ++ mdk-stage1/dietlibc/libm/atanh.c | 8 ++ mdk-stage1/dietlibc/libm/bessel.c | 171 ++++++++++++++++++++++++++++++++++++++ mdk-stage1/dietlibc/libm/cosh.c | 9 ++ mdk-stage1/dietlibc/libm/erf.c | 95 +++++++++++++++++++++ mdk-stage1/dietlibc/libm/gamma.c | 98 ++++++++++++++++++++++ mdk-stage1/dietlibc/libm/ipow.c | 29 +++++++ mdk-stage1/dietlibc/libm/poly.c | 41 +++++++++ mdk-stage1/dietlibc/libm/pow.c | 42 ++++++++++ mdk-stage1/dietlibc/libm/rint.c | 5 ++ mdk-stage1/dietlibc/libm/sinh.c | 9 ++ mdk-stage1/dietlibc/libm/tanh.c | 7 ++ 13 files changed, 526 insertions(+) create mode 100644 mdk-stage1/dietlibc/libm/acosh.c create mode 100644 mdk-stage1/dietlibc/libm/asinh.c create mode 100644 mdk-stage1/dietlibc/libm/atanh.c create mode 100644 mdk-stage1/dietlibc/libm/bessel.c create mode 100644 mdk-stage1/dietlibc/libm/cosh.c create mode 100644 mdk-stage1/dietlibc/libm/erf.c create mode 100644 mdk-stage1/dietlibc/libm/gamma.c create mode 100644 mdk-stage1/dietlibc/libm/ipow.c create mode 100644 mdk-stage1/dietlibc/libm/poly.c create mode 100644 mdk-stage1/dietlibc/libm/pow.c create mode 100644 mdk-stage1/dietlibc/libm/rint.c create mode 100644 mdk-stage1/dietlibc/libm/sinh.c create mode 100644 mdk-stage1/dietlibc/libm/tanh.c (limited to 'mdk-stage1/dietlibc/libm') 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 + +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 + +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 + +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 + +#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 + +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 +#include + +#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 +#include + +#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 +/* + * 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 +#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 +#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 + +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 + +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 + +double tanh ( double x ) +{ + long double y = exp (x + x); + return (y - 1.) / (y + 1.); +} -- cgit v1.2.1