Skip to content

Commit 38e8866

Browse files
author
camilo
committed
added tgamma and lgamma functions to ffmath. Small refactor and checks. bump to 1.1.3
1 parent d4c4396 commit 38e8866

File tree

13 files changed

+578
-224
lines changed

13 files changed

+578
-224
lines changed

‎check/main.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
using namespace qlibs;
44

55
int main() {
6-
6+
/*
77
mat<2,2> x(
88
1.0f, 2.0f,
99
3.0f, 4.0f
@@ -21,6 +21,6 @@ int main() {
2121
x*=x;
2222
result.display();
2323
x.display();
24-
24+
*/
2525
return 0;
2626
}

‎check/qlibs_cpp_test.cpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -304,6 +304,8 @@ void test_ffmath(void)
304304
cout << ffmath::wrapToPi( -4.5 ) << endl;
305305
cout << ffmath::getNan() << endl;
306306
cout << ffmath::getInf() << endl;
307+
cout << ffmath::tgamma( 12.5 ) << endl;
308+
cout << ffmath::lgamma( 12.5 ) << endl;
307309
}
308310

309311
void test_mat( void )
@@ -425,11 +427,13 @@ int main()
425427
test_fp16();
426428
test_ltisys();
427429
test_tdl();
430+
431+
test_mat();
432+
test_interp1();
428433
test_fis();
429434
test_fis2();
430435
test_fis3();
431436
test_ffmath();
432-
test_mat();
433-
test_interp1();
437+
434438
return 0;
435439
}

‎library.json‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
"maintainer": true
1717
}
1818
],
19-
"version": "1.1.2",
19+
"version": "1.1.3",
2020
"license": "MIT",
2121
"frameworks": "arduino",
2222
"platforms": "*"

‎library.properties‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
name=qlibs
2-
version=1.1.2
2+
version=1.1.3
33
license=MIT
44
author=J. Camilo Gomez C. <[email protected]>
55
maintainer=J. Camilo Gomez C. <[email protected]>

‎src/ffmath.cpp‎

Lines changed: 321 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ static inline void cast_reinterpret( T1 &f, const T2 &u )
1111

1212
static float getAbnormal( const int i );
1313
static float compute_cbrt( float x , bool r );
14+
static float lgamma_positive( float x );
1415

1516
static inline float absolute( const float x )
1617
{
@@ -31,6 +32,13 @@ static float getAbnormal( const int i )
3132
return f_ab[ i ];
3233
}
3334
/*============================================================================*/
35+
bool ffmath::isEqual( const float a,
36+
const float b,
37+
const float tol ) noexcept
38+
{
39+
return ( ffmath::absf( a - b ) <= ffmath::absf( tol ) );
40+
}
41+
/*============================================================================*/
3442
float ffmath::getInf( void )
3543
{
3644
return getAbnormal( 0 );
@@ -568,3 +576,316 @@ float ffmath::nextAfter( float x, float y )
568576
return retVal;
569577
}
570578
/*============================================================================*/
579+
float ffmath::tgamma( float x )
580+
{
581+
float result;
582+
583+
const auto fClass = ffmath::classify( x );
584+
if ( classification::FFP_NAN == fClass ) {
585+
result = getNan();
586+
}
587+
else if ( classification::FFP_ZERO == fClass ) {
588+
result = getInf(); /* a huge value */
589+
}
590+
else if ( classification::FFP_INFINITE == fClass ) {
591+
if ( x > 0.0f ) {
592+
result = getInf(); /* a huge value */
593+
}
594+
else {
595+
result = getNan();
596+
}
597+
}
598+
else {
599+
bool parity = false;
600+
float fact = 1.0F;
601+
float y = x;
602+
float y1;
603+
604+
if ( y <= 0.0F ) {
605+
float isItAnInt;
606+
607+
y = -x;
608+
y1 = ffmath::trunc( y );
609+
isItAnInt = y - y1;
610+
if ( !ffmath::isEqual( 0.0F, isItAnInt ) ) {
611+
const float tmp = 2.0F*ffmath::trunc( y1*0.5F );
612+
613+
if ( !ffmath::isEqual( y1, tmp ) ) {
614+
parity = true;
615+
}
616+
fact = -FFP_PI/ffmath::sin( FFP_PI*isItAnInt );
617+
y += 1.0F;
618+
}
619+
else {
620+
result = ffmath::getNan();
621+
goto EXIT_TGAMMA;
622+
}
623+
}
624+
if ( y < 1.19209290E-07F ) { /* y < eps */
625+
if ( y >= 1.175494351e-38F ) { /* y >= MinimumX */
626+
result = 1.0F/y;
627+
}
628+
else {
629+
result = getInf();
630+
}
631+
}
632+
else if ( y < 12.0F ) {
633+
float num = 0.0F, den = 1.0F, z;
634+
int n = 0;
635+
636+
y1 = y;
637+
if ( y < 1.0F ) {
638+
z = y;
639+
y += 1.0F;
640+
}
641+
else {
642+
n = static_cast<int>( y ) - 1;
643+
/*cstat -CERT-FLP36-C */
644+
y -= static_cast<float>( n );
645+
/*cstat +CERT-FLP36-C */
646+
z = y - 1.0F;
647+
}
648+
649+
num = z*( num + -1.71618513886549492533811e+0F );
650+
den = ( den*z ) -3.08402300119738975254353e+1F;
651+
num = z*( num + 2.47656508055759199108314e+1F );
652+
den = ( den*z ) + 3.15350626979604161529144e+2F;
653+
num = z*( num - 3.79804256470945635097577e+2F );
654+
den = ( den*z ) - 1.01515636749021914166146e+3F;
655+
num = z*( num + 6.29331155312818442661052e+2F );
656+
den = ( den*z ) - 3.10777167157231109440444e+3F;
657+
num = z*( num + 8.66966202790413211295064e+2F );
658+
den = ( den*z ) + 2.25381184209801510330112e+4F;
659+
num = z*( num - 3.14512729688483675254357e+4F );
660+
den = ( den*z ) + 4.75584627752788110767815e+3F;
661+
num = z*( num - 3.61444134186911729807069e+4F );
662+
den = ( den*z ) - 1.34659959864969306392456e+5F;
663+
num = z*( num + 6.64561438202405440627855e+4F );
664+
den = ( den*z ) - 1.15132259675553483497211e+5F;
665+
666+
result = ( num/den ) + 1.0F;
667+
if ( y1 < y ) {
668+
result /= y1;
669+
}
670+
else if ( y1 > y ) {
671+
for ( int i = 0; i < n ; ++i ) {
672+
result *= y;
673+
y += 1.0F;
674+
}
675+
}
676+
}
677+
else {
678+
if ( x <= 171.624F ) { /* x <= xBig */
679+
const float yy = y*y;
680+
float sum = 5.7083835261e-03F;
681+
682+
sum = ( sum/yy ) - 1.910444077728e-03F;
683+
sum = ( sum/yy ) + 8.4171387781295e-04F;
684+
sum = ( sum/yy ) - 5.952379913043012e-04F;
685+
sum = ( sum/yy ) + 7.93650793500350248e-04F;
686+
sum = ( sum/yy ) - 2.777777777777681622553e-03F;
687+
sum = ( sum/yy ) + 8.333333333333333331554247e-02F;
688+
689+
sum = ( sum /y ) - y + FFP_LN_SQRT_2PI;
690+
sum += ( y - 0.5F )*ffmath::log( y );
691+
result = ffmath::exp( sum );
692+
}
693+
else {
694+
result = ffmath::getInf();
695+
}
696+
}
697+
if ( parity ) {
698+
result = -result;
699+
}
700+
if ( !ffmath::isEqual( fact, 1.0F ) ) {
701+
result = fact/result;
702+
}
703+
}
704+
705+
EXIT_TGAMMA:
706+
return result;
707+
}
708+
/*============================================================================*/
709+
static float lgamma_positive( float x )
710+
{
711+
constexpr float d1 = -5.772156649015328605195174e-1F;
712+
constexpr float d2 = 4.227843350984671393993777e-1F;
713+
constexpr float d4 = 1.791759469228055000094023e+0F;
714+
constexpr float p1[ 8 ] = { 4.945235359296727046734888e+0F,
715+
2.018112620856775083915565e+2F,
716+
2.290838373831346393026739e+3F,
717+
1.131967205903380828685045e+4F,
718+
2.855724635671635335736389e+4F,
719+
3.848496228443793359990269e+4F,
720+
2.637748787624195437963534e+4F,
721+
7.225813979700288197698961e+3F };
722+
constexpr float q1[ 8 ] = { 6.748212550303777196073036e+1F,
723+
1.113332393857199323513008e+3F,
724+
7.738757056935398733233834e+3F,
725+
2.763987074403340708898585e+4F,
726+
5.499310206226157329794414e+4F,
727+
6.161122180066002127833352e+4F,
728+
3.635127591501940507276287e+4F,
729+
8.785536302431013170870835e+3F };
730+
constexpr float p2[ 8 ] = { 4.974607845568932035012064e+0F,
731+
5.424138599891070494101986e+2F,
732+
1.550693864978364947665077e+4F,
733+
1.847932904445632425417223e+5F,
734+
1.088204769468828767498470e+6F,
735+
3.338152967987029735917223e+6F,
736+
5.106661678927352456275255e+6F,
737+
3.074109054850539556250927e+6F };
738+
constexpr float q2[ 8 ] = { 1.830328399370592604055942e+2F,
739+
7.765049321445005871323047e+3F,
740+
1.331903827966074194402448e+5F,
741+
1.136705821321969608938755e+6F,
742+
5.267964117437946917577538e+6F,
743+
1.346701454311101692290052e+7F,
744+
1.782736530353274213975932e+7F,
745+
9.533095591844353613395747e+6F };
746+
constexpr float p4[ 8 ] = { 1.474502166059939948905062e+04F,
747+
2.426813369486704502836312e+06F,
748+
1.214755574045093227939592e+08F,
749+
2.663432449630976949898078e+09F,
750+
2.940378956634553899906876e+10F,
751+
1.702665737765398868392998e+11F,
752+
4.926125793377430887588120e+11F,
753+
5.606251856223951465078242e+11F };
754+
constexpr float q4[ 8 ] = { 2.690530175870899333379843e+03F,
755+
6.393885654300092398984238e+05F,
756+
4.135599930241388052042842e+07F,
757+
1.120872109616147941376570e+09F,
758+
1.488613728678813811542398e+10F,
759+
1.016803586272438228077304e+11F,
760+
3.417476345507377132798597e+11F,
761+
4.463158187419713286462081e+11F };
762+
constexpr float pnt68 = 0.6796875F;
763+
float result, y;
764+
765+
if ( x > 171.624F ) {
766+
result = ffmath::getInf(); /* a huge value */
767+
}
768+
else {
769+
float corrector, num, den;
770+
771+
y = x;
772+
if ( y <= 1.19209290E-07F ) { /* y < eps */
773+
result = -ffmath::log( y );
774+
}
775+
else if ( y <= 1.5F ) {
776+
float xMinus;
777+
778+
if ( y < pnt68 ) {
779+
corrector = -ffmath::log( y );
780+
xMinus = y;
781+
}
782+
else {
783+
corrector = 0.0F;
784+
xMinus = ( y - 0.5F ) - 0.5F;
785+
}
786+
if ( ( y <= 0.5F ) || ( y >= pnt68 ) ) {
787+
den = 1.0F;
788+
num = 0.0F;
789+
for ( size_t i = 0U ; i < 8U ; ++i ) {
790+
num = ( num*xMinus ) + p1[ i ];
791+
den = ( den*xMinus ) + q1[ i ];
792+
}
793+
result = corrector + ( xMinus*( d1 + xMinus*( num/den ) ) );
794+
}
795+
else {
796+
xMinus = ( y - 0.5F ) - 0.5F;
797+
den = 1.0F;
798+
num = 0.0F;
799+
for ( size_t i = 0U; i < 8U ; ++i ) {
800+
num = num*xMinus + p2[ i ];
801+
den = den*xMinus + q2[ i ];
802+
}
803+
result = corrector + ( xMinus*( d2 + xMinus*( num/den ) ) );
804+
}
805+
}
806+
else if ( y <= 4.0F ) {
807+
const float xMinus = y - 2.0F;
808+
den = 1.0F;
809+
num = 0.0F;
810+
for ( size_t i = 0U; i < 8U ; ++i ) {
811+
num = num*xMinus + p2[ i ];
812+
den = den*xMinus + q2[ i ];
813+
}
814+
result = xMinus*( d2 + xMinus*( num/den ) );
815+
}
816+
else if ( y <= 12.0F ) {
817+
const float xMinus = y - 4.0F;
818+
den = -1.0F;
819+
num = 0.0F;
820+
for ( size_t i = 0U; i < 8U ; ++i ) {
821+
num = num*xMinus + p4[ i ];
822+
den = den*xMinus + q4[ i ];
823+
}
824+
result = d4 + xMinus*( num/den );
825+
}
826+
else {
827+
result = 0.0F;
828+
if ( y <= 4294967296.87842273712158203125F ) { /* y < xBig^(1/4)*/
829+
const float yy = y*y;
830+
result = 5.7083835261e-03F;
831+
result = ( result/yy ) - 1.910444077728e-03F;
832+
result = ( result/yy ) + 8.4171387781295e-04F;
833+
result = ( result/yy ) - 5.952379913043012e-04F;
834+
result = ( result/yy ) + 7.93650793500350248e-04F;
835+
result = ( result/yy ) - 2.777777777777681622553e-03F;
836+
result = ( result/yy ) + 8.333333333333333331554247e-02F;
837+
}
838+
result /= y;
839+
corrector = ffmath::log( y );
840+
result += ffmath::FFP_LN_SQRT_2PI - ( 0.5F*corrector );
841+
result += y*( corrector - 1.0F );
842+
}
843+
}
844+
845+
return result;
846+
}
847+
/*============================================================================*/
848+
float ffmath::lgamma( float x )
849+
{
850+
float result;
851+
852+
const auto fClass = ffmath::classify( x );
853+
if ( classification::FFP_NAN == fClass ) {
854+
result = getNan();
855+
}
856+
else if ( classification::FFP_ZERO == fClass ) {
857+
result = getInf();
858+
}
859+
else if ( classification::FFP_INFINITE == fClass ) {
860+
result = getInf();
861+
}
862+
else {
863+
float y;
864+
865+
if ( x < 0.0F ) {
866+
if ( x <= -4503599627370496.0F ) { /* x < 2^52 */
867+
result = getInf();
868+
}
869+
else {
870+
float a, y1, isItAnInt;
871+
872+
y = -x;
873+
y1 = ffmath::trunc( y );
874+
isItAnInt = y - y1;
875+
if ( ffmath::isEqual( 0.0F, isItAnInt ) ) {
876+
result = getInf();
877+
}
878+
else {
879+
a = sin( FFP_PI*isItAnInt );
880+
result = ffmath::log( FFP_PI/ffmath::absf( a*x ) ) - lgamma_positive( -x );
881+
}
882+
}
883+
}
884+
else {
885+
result = lgamma_positive( x );
886+
}
887+
}
888+
889+
return result;
890+
}
891+
/*============================================================================*/

0 commit comments

Comments
 (0)