MAPM & ExpressionParser: Support long long int and cleanup Stirling's

Change-Id: I4a2bd3f36eb891a5f3167e93057ded7427d6cf26
Reviewed-on: https://review.haiku-os.org/c/haiku/+/8996
Reviewed-by: John Scipione <[email protected]>
Tested-by: Commit checker robot <[email protected]>

headers/libs/mapm/m_apm.h

@@ -225,6 +225,7 @@ extern	char	*m_apm_lib_short_version(char *);
 extern	void	m_apm_set_string(M_APM, char *);
 extern	void	m_apm_set_double(M_APM, double);
 extern	void	m_apm_set_long(M_APM, long);
+extern	void	m_apm_set_long_long(M_APM, long long);
 
 extern	void	m_apm_to_string(char *, int, M_APM);
 extern  void	m_apm_to_fixpt_string(char *, int, M_APM);
@@ -421,6 +422,8 @@ public:
 		{create();m_apm_set_long(val(),l);}
 	MAPM(long l) /* Constructor from long int */
 		{create();m_apm_set_long(val(),l);}
+	MAPM(long long ll) /* Constructor from long long int */
+		{create();m_apm_set_long_long(val(),ll);}
 	/* Destructor */
 	~MAPM() {destroy();}
 	

src/kits/shared/ExpressionParser.cpp

@@ -748,14 +748,12 @@ ExpressionParser::_ParseFactorial(MAPM value)
 	if (fTokenizer->NextToken().type == TOKEN_FACTORIAL) {
 		fTokenizer->RewindToken();
 		_EatToken(TOKEN_FACTORIAL);
-		if (value < 1000)
+		if (value.abs() < 1000) {
 			return value.factorial();
-		else {
+		} else {
 			// Use Stirling's approximation (9 term expansion)
 			// http://en.wikipedia.org/wiki/Stirling%27s_approximation
 			// http://www.wolframalpha.com/input/?i=stirling%27s+series
-			// all constants must fit in a signed long for MAPM
-			// (LONG_MAX = 2147483647)
 			return value.pow(value) / value.exp()
 				* (MAPM(2) * MAPM(MM_PI) * value).sqrt()
 				* (MAPM(1) + (MAPM(1) / (MAPM(12) * value))
@@ -763,24 +761,10 @@ ExpressionParser::_ParseFactorial(MAPM value)
 					- (MAPM(139) / (MAPM(51840) * value.pow(3)))
 					- (MAPM(571) / (MAPM(2488320) * value.pow(4)))
 					+ (MAPM(163879) / (MAPM(209018880) * value.pow(5)))
-						// 2147483647 * 35 + 84869155 = 75246796800
+					+ (MAPM(5246819) / (MAPM(75246796800LL) * value.pow(6)))
-					+ (MAPM(5246819) / ((MAPM(2147483647) * MAPM(35)
+					- (MAPM(534703531) / (MAPM(902961561600LL) * value.pow(7)))
-						+ MAPM(84869155)) * value.pow(6)))
+					- (MAPM(4483131259LL) / (MAPM(86686309913600LL) * value.pow(8)))
-						// 2147483647 * 420 + 1018429860 = 902961561600
+					+ (MAPM(432261921612371LL) / (MAPM(514904800886784000LL) * value.pow(9))));
-					- (MAPM(534703531) / ((MAPM(2147483647) * MAPM(420)
-						+ MAPM(1018429860)) * value.pow(7)))
-						// 2147483647 * 2 + 188163965 = 4483131259
-						// 2147483647 * 40366 + 985018798 = 86686309913600
-					- ((MAPM(2147483647) * MAPM(2) + MAPM(188163965))
-						/ ((MAPM(2147483647) * MAPM(40366) + MAPM(985018798))
-							* value.pow(8)))
-						// 2147483647 * 201287 + 1380758682 = 432261921612371
-						// 2147483647 * 239771232 + 1145740896
-						// = 514904800886784000
-					+ ((MAPM(2147483647) * MAPM(201287) + MAPM(1380758682))
-						/ ((MAPM(2147483647) * MAPM(239771232)
-								+ MAPM(1145740896))
-							* value.pow(9))));
 		}
 	}
 

src/libs/mapm/m_apm_lc.h

@@ -22,6 +22,9 @@
  *      This is the local header file needed to build the library
  *
  *      $Log: m_apm_lc.h,v $
+ *      Revision 1.46  2025/02/16 12:22:24  jscipione
+ *      add new prototype for M_long_long_2_ascii
+ *
  *      Revision 1.45  2007/12/04 01:26:02  mike
  *      add support for Digital Mars compiler
  *
@@ -388,6 +391,7 @@ extern	void	M_apm_normalize(M_APM);
 extern	void	M_apm_scale(M_APM, int);
 extern	void	M_apm_pad(M_APM, int);
 extern  void    M_long_2_ascii(char *, long);
+extern  void    M_long_long_2_ascii(char *, long long);
 extern	void	M_check_PI_places(int);
 extern  void    M_calculate_PI_AGM(M_APM, int);
 extern  void    M_set_to_zero(M_APM);

src/libs/mapm/mapm_exp.c

@@ -365,7 +365,7 @@ void	M_raw_exp(M_APM rr, int places, M_APM xx)
 {
 M_APM   tmp0, digit, term;
 int	tolerance,  local_precision, prev_exp;
-long    m1;
+long long    m1;
 
 tmp0  = M_get_stack_var();
 term  = M_get_stack_var();
@@ -378,11 +378,11 @@ prev_exp        = 0;
 m_apm_add(rr, MM_One, xx);
 m_apm_copy(term, xx);
 
-m1 = 2L;
+m1 = 2LL;
 
 while (TRUE)
   {
-   m_apm_set_long(digit, m1);
+   m_apm_set_long_long(digit, m1);
    m_apm_multiply(tmp0, term, xx);
    m_apm_divide(term, local_precision, tmp0, digit);
    m_apm_add(tmp0, rr, term);

src/libs/mapm/mapm_gcd.c

@@ -62,7 +62,7 @@ void	m_apm_gcd(M_APM r, M_APM u, M_APM v)
 {
 M_APM   tmpM, tmpN, tmpT, tmpU, tmpV;
 int	kk, kr, mm;
-long    pow_2;
+long long    pow_2;
 
 /* 'is_integer' will return 0 || 1 */
 
@@ -174,11 +174,11 @@ else
      {
       mm = kk / 28;
       kr = kk % 28;
-      pow_2 = 1L << kr;
+      pow_2 = 1LL << kr;
 
       if (mm == 0)
         {
-	 m_apm_set_long(tmpN, pow_2);
+	 m_apm_set_long_long(tmpN, pow_2);
          m_apm_multiply(r, tmpU, tmpN);
 	}
       else
@@ -201,7 +201,7 @@ else
 	   }
 	 else
 	   {
-	    m_apm_set_long(tmpM, pow_2);
+	    m_apm_set_long_long(tmpM, pow_2);
             m_apm_multiply(tmpT, tmpN, tmpM);
             m_apm_multiply(r, tmpU, tmpT);
 	   }

src/libs/mapm/mapm_lg4.c

@@ -54,7 +54,7 @@ void	M_log_near_1(M_APM rr, int places, M_APM xx)
 {
 M_APM   tmp0, tmp1, tmp2, tmpS, term;
 int	tolerance, dplaces, local_precision;
-long    m1;
+long long    m1;
 
 tmp0 = M_get_stack_var();
 tmp1 = M_get_stack_var();
@@ -72,7 +72,7 @@ m_apm_copy(term, tmpS);
 m_apm_multiply(tmp0, tmpS, tmpS);
 m_apm_round(tmp2, (dplaces + 6), tmp0);
 
-m1 = 3L;
+m1 = 3LL;
 
 while (TRUE)
   {
@@ -86,7 +86,7 @@ while (TRUE)
    if (local_precision < 20)
      local_precision = 20;
 
-   m_apm_set_long(tmp1, m1);
+   m_apm_set_long_long(tmp1, m1);
    m_apm_round(term, local_precision, tmp0);
    m_apm_divide(tmp0, local_precision, term, tmp1);
    m_apm_add(tmp1, tmpS, tmp0);

src/libs/mapm/mapm_log.c

@@ -146,7 +146,8 @@ M_restore_stack(2);                    /* restore the 2 locals we used here */
 void	m_apm_log(M_APM r, int places, M_APM a)
 {
 M_APM   tmp0, tmp1, tmp2;
-int	mexp, dplaces;
+long long     mexp;
+int     dplaces;
 
 if (a->m_apm_sign <= 0)
   {
@@ -217,7 +218,7 @@ else
    
    M_log_basic_iteration(tmp0, dplaces, tmp2);
    
-   m_apm_set_long(tmp1, (long)mexp);
+   m_apm_set_long_long(tmp1, mexp);
    m_apm_multiply(tmp2, tmp1, MM_lc_log10);
    m_apm_add(tmp1, tmp2, tmp0);
    

src/libs/mapm/mapm_set.c

@@ -100,6 +100,11 @@ if (M_lbuf != 0)
 }
 /****************************************************************************/
 void	m_apm_set_long(M_APM atmp, long mm)
+{
+  m_apm_set_long_long(atmp, (long long)mm);
+}
+/****************************************************************************/
+void	m_apm_set_long_long(M_APM atmp, long long mm)
 {
 int     len, ii, nbytes;
 char	*p, *buf, ch, buf2[64];
@@ -112,7 +117,7 @@ if (mm == 0)
    return;
   }
 
-M_long_2_ascii(buf2, mm);     /* convert long -> ascii in base 10 */
+M_long_long_2_ascii(buf2, mm); /* convert long long int -> ascii in base 10 */
 buf = buf2;
 
 if (mm < 0)

src/libs/mapm/mapmasn0.c

@@ -109,7 +109,7 @@ void	M_arctan_near_0(M_APM rr, int places, M_APM aa)
 {
 M_APM   tmp0, tmp2, tmpR, tmpS, digit, term;
 int	tolerance, dplaces, local_precision;
-long    m1;
+long long    m1;
 
 tmp0  = M_get_stack_var();
 tmp2  = M_get_stack_var();
@@ -126,7 +126,7 @@ m_apm_copy(tmpS, aa);
 m_apm_multiply(tmp0, aa, aa);
 m_apm_round(tmp2, (dplaces + 8), tmp0);
 
-m1 = 1L;
+m1 = 1LL;
 
 while (TRUE)
   {
@@ -148,7 +148,7 @@ while (TRUE)
      local_precision = 20;
 
    m1 += 2;
-   m_apm_set_long(digit, m1);
+   m_apm_set_long_long(digit, m1);
    m_apm_round(term, local_precision, tmp0);
    m_apm_divide(tmp0, local_precision, term, digit);
    m_apm_subtract(tmpR, tmpS, tmp0);
@@ -171,7 +171,7 @@ while (TRUE)
      local_precision = 20;
 
    m1 += 2;
-   m_apm_set_long(digit, m1);
+   m_apm_set_long_long(digit, m1);
    m_apm_round(term, local_precision, tmp0);
    m_apm_divide(tmp0, local_precision, term, digit);
    m_apm_add(tmpS, tmpR, tmp0);

src/libs/mapm/mapmrsin.c

@@ -64,7 +64,7 @@ void	M_raw_sin(M_APM rr, int places, M_APM xx)
 {
 M_APM	sum, term, tmp2, tmp7, tmp8;
 int     tolerance, flag, local_precision, dplaces;
-long	m1, m2;
+long long	m1, m2;
 
 sum  = M_get_stack_var();
 term = M_get_stack_var();
@@ -80,7 +80,7 @@ m_apm_round(tmp2, (places + 6), tmp8);
 dplaces   = (places + 8) - xx->m_apm_exponent;
 tolerance = xx->m_apm_exponent - (places + 4);
 
-m1   = 2L;
+m1   = 2LL;
 flag = 0;
 
 while (TRUE)
@@ -128,7 +128,7 @@ void	M_raw_cos(M_APM rr, int places, M_APM xx)
 {
 M_APM	sum, term, tmp7, tmp8, tmp9;
 int     tolerance, flag, local_precision, prev_exp;
-long	m1, m2;
+long long	m1, m2;
 
 sum  = M_get_stack_var();
 term = M_get_stack_var();
@@ -146,13 +146,13 @@ local_precision = places + 8;
 tolerance       = -(places + 4);
 prev_exp        = 0;
 
-m1   = 1L;
+m1   = 1LL;
 flag = 0;
 
 while (TRUE)
   {
    m2 = m1 * (m1 + 1);
-   m_apm_set_long(tmp7, m2);
+   m_apm_set_long_long(tmp7, m2);
 
    m_apm_multiply(tmp8, term, tmp9);
    m_apm_divide(term, local_precision, tmp8, tmp7);

src/libs/mapm/mapmutl2.c

@@ -244,17 +244,41 @@ else
 void    M_long_2_ascii(char *output, long input)
 {
 long    t, m;
+char    *p;
+
+m = input;
+p = output;
+t = 2147000000L;          /* something < 2^31 */
+
+if ((m > t) || (m < -t))  /* handle the bigger numbers with 'sprintf'. */
+  {			  /* let them worry about wrap-around problems */
+   sprintf(p, "%ld", m);  /* at 'LONG_MIN', etc. */
+  }
+else
+  {
+   M_long_long_2_ascii(output, ((long long)input));
+  }
+}
+/****************************************************************************/
+/*
+ *
+ *	convert a signed long long to ASCII in base 10
+ *
+ */
+void    M_long_long_2_ascii(char *output, long long input)
+{
+long long t, m;
 int     i, j;
 char    *p, tbuf[64];
 
 m = input;
 p = output;
 i = 0;
-t = 2147000000L;          /* something < 2^31 */
+t = 9223372036854000000LL;  /* something < 2^63 */
 
 if ((m > t) || (m < -t))  /* handle the bigger numbers with 'sprintf'. */
   {			  /* let them worry about wrap-around problems */
-   sprintf(p, "%ld", m);  /* at 'LONG_MIN', etc.                       */
+   sprintf(p, "%lld", m);  /* at 'LLONG_MIN', etc. */
   }
 else
   {