Changeset 2639 in Sophya

Timestamp:

Nov 22, 2004, 5:57:40 PM (21 years ago)

Author:

cmv

Message:

factoriel et approx cmv 22/11/04

Location:

trunk/SophyaLib/NTools

Files:

• nbmath.c (modified) (1 diff)
• nbmath.h (modified) (1 diff)

trunk/SophyaLib/NTools/nbmath.c

@@ -2734 +2734 @@
 return(nsuse);
 }
+
+/*=========================================== cmv 22/11/04 ===*/
+/*
+++
+double log_factoriel(unsigned int n)
+
+        Compute the neperian log of n!
+--
+*/
+double log_factoriel(unsigned int n)
+{
+  unsigned int i;
+  double sum=0.;
+  if(n<2) return sum;
+  for(i=n;i>1;i--) sum += log((double)i);
+  return sum;
+}
+
+/*
+++
+double log_stirling(unsigned int n)
+
+        Compute the neperian log of the Stirling approx of n!
+--
+*/
+double log_stirling(unsigned int n)
+/*
+Approx: n! ~ sqrt(2Pi n) n^n exp(-n) = sqrt(2Pi) n^(n+1/2) exp(-n)
+        log(n!) ~ (n+1/2) log(n) -n + 1/2 log(2Pi)
+(ordre 0: log(n!) ~ n (log(n)-1)
+*/
+{
+  if(n<1) return 0.;
+  return (n+0.5)*log(n) - n + Hln2pi;
+}
+
+/*
+++
+double log_gosper(unsigned int n);
+
+        Compute the neperian log of the Gosper approx of n!
+--
+*/
+double log_gosper(unsigned int n)
+/*
+Approx: n! ~ sqrt( (2n+1/3) Pi ) n^n exp(-n)
+        log(n!) ~ 1/2 log( (2n+1/3) Pi ) + n log(n) -n
+*/
+{
+  if(n<1) return 0.;
+  return 0.5*log((2.*n+1./3.)*M_PI) + n*log(n) -n;
+}
+
+/*
+++
+double give_log_factoriel(unsigned int n);
+
+        Compute the neperian log of the approx of n!
+--
+*/
+static short give_log_factoriel_OK = 0;
+#define give_log_factoriel_LIM 21
+static double give_log_factoriel_TAB[give_log_factoriel_LIM];
+double give_log_factoriel(unsigned int n)
+{
+  int i;
+  if(give_log_factoriel_OK==0) {
+    give_log_factoriel_OK=1;
+    for(i=0;i<give_log_factoriel_LIM;i++)
+      give_log_factoriel_TAB[i]=log_factoriel(i);
+  }
+
+  if(n<give_log_factoriel_LIM) return give_log_factoriel_TAB[n];
+  return log_gosper(n);
+}
+#undef give_log_factoriel_LIM

trunk/SophyaLib/NTools/nbmath.h

@@ -90 +90 @@
              ,int npass,float perclean,int lp);
 
+double log_factoriel(unsigned int n);
+double log_stirling(unsigned int n);
+double log_gosper(unsigned int n);
+double give_log_factoriel(unsigned int n);
+
 #ifdef __cplusplus
 }