lib/distribution/math_extension.rb in distribution-0.2.0 vs lib/distribution/math_extension.rb in distribution-0.3.0

@@ -1,6 +1,18 @@
-require 'prime'
+if RUBY_VERSION<"1.9"
+  require 'mathn'
+  def Prime.each(upper,&block) 
+    @primes=Prime.new
+    @primes.each do |prime|
+      break if prime > upper.to_i
+      block.call(prime)
+    end
+  end
+else
+  require 'prime'
+end
+
 require 'bigdecimal'
 require 'bigdecimal/math'
 
 # Useful additions to Math
 module Distribution
@@ -36,34 +48,35 @@
       end
       def swing(n)
         return SmallOddSwing[n] if (n<33)
         sqrtN = Math.sqrt(n).floor
         count=0
+        
         Prime.each(n/3) do |prime|
           next if prime<3
           if (prime<=sqrtN)
             q=n
             _p=1
-            while((q/=prime)>0) do
-              if ((q&1)==1)
+            
+            while((q=(q/prime).truncate)>0) do
+              if((q%2)==1)
                 _p*=prime
               end
             end
-            
             if _p>1
               @prime_list[count]=_p
               count+=1
             end
             
           else
-            if ((n/prime)&1==1)
+            if ((n/prime).truncate%2==1)
               @prime_list[count]=prime
               count+=1
             end
           end
         end
-        prod=get_primorial(n/2+1,n)
+        prod=get_primorial((n/2).truncate+1,n)
         prod * @prime_list[0,count].inject(1) {|ac,v| ac*v}
       end
       def get_primorial(low,up)
         prod=1;
         Prime.each(up) do |prime|
@@ -150,11 +163,42 @@
     # Source:
     # * http://mathworld.wolfram.com/BetaFunction.html
     def beta(x,y)
       (gamma(x)*gamma(y)).quo(gamma(x+y))
     end
+    # I_x(a,b): Regularized incomplete beta function
+    #   
+    # Source:
+    #
+    def regularized_beta_function(x,a,b)
+      return 1 if x==1
+      #incomplete_beta(x,a,b).quo(beta(a,b))
+      m=a
+      n=b+a-1
+      (m..n).inject(0) {|sum,j|
+        sum+(binomial_coefficient(n,j)* x**j * (1-x)**(n-j))
+      }
+
+    end
+    # B_x(a,b) : Incomplete beta function
+    # Should be replaced by
+    # http://lib.stat.cmu.edu/apstat/63
+    def incomplete_beta(x,a,b)
+      raise "Not work"
+      return beta(a,b) if x==1
+      
+      ((x**a * (1-x)**b).quo(a)) * hyper_f(a+b,1,a+1,x)
+    end
+    def permutations(x,n)
+      factorial(x).quo(factorial(x-n))
+    end
     
+    def rising_factorial(x,n)
+      factorial(x+n-1).quo(factorial(x-1))
+    end
+   
+    
     LOG_2PI = Math.log(2 * Math::PI)# log(2PI)
     N = 8
     B0  = 1.0
     B1  = -1.0 / 2.0
     B2  = 1.0 / 6.0
@@ -190,19 +234,57 @@
       if (x < 0.0)
         return Math::PI / (Math.sin(Math.PI * x) * Math.exp(loggamma(1 - x))) #/
       end
       Math.exp(loggamma(x))
     end
+    # Binomial coeffients, or:
+    # ( n )
+    # ( k )
+    # Gives the number of different k size subsets of a set size n
+    # 
+    # Replaces (n,k) for (n, n-k) if k>n-k
+    #
+    #  (n)   n^k'    (n)..(n-k+1)
+    #  ( ) = ---- =  ------------
+    #  (k)    k!          k!
+    #
+    def binomial_coefficient(n,k)
+      return 1 if (k==0 or k==n)
+      k=[k, n-k].min
+      (((n-k+1)..n).inject(1) {|ac,v| ac * v}).quo(factorial(k))
+      # Other way to calcule binomial is this: 
+      # (1..k).inject(1) {|ac, i| (ac*(n-k+i).quo(i))}
+    end
+    # Approximate binomial coefficient, using gamma function.
+    # The fastest method, until we fall on BigDecimal!
+    def binomial_coefficient_gamma(n,k)
+      return 1 if (k==0 or k==n)
+      k=[k, n-k].min
+      
+      val=gamma(n+1) / (gamma(k+1)*gamma(n-k+1))
+      if (val.nan?)
+        lg=lgamma(n+1) - (lgamma(k+1)+lgamma(n-k+1))
+        val=Math.exp(lg)
+        # Crash again! We require BigDecimals
+        if val.infinite?
+          val=BigMath.exp(BigDecimal(lg.to_s),16)
+        end
+      end
+      
+      val
+    end
   end
 end
 
 module Math
   include Distribution::MathExtension
-  module_function :factorial, :beta, :gamma, :gosper, :loggamma, :fast_factorial
+  alias :lgamma :loggamma 
+
+  module_function :factorial, :beta, :gamma, :gosper, :loggamma, :lgamma, :binomial_coefficient, :binomial_coefficient_gamma, :regularized_beta_function, :incomplete_beta, :permutations, :rising_factorial , :fast_factorial
 end
 
 # Necessary on Ruby 1.9
 module CMath # :nodoc:
   include Distribution::MathExtension
-  module_function :factorial, :beta, :gamma, :gosper, :loggamma, :fast_factorial
+  module_function :factorial, :beta, :gosper, :loggamma,  :binomial_coefficient, :binomial_coefficient_gamma, :regularized_beta_function, :incomplete_beta, :permutations, :rising_factorial, :fast_factorial
 end