HW 3, 1.22

By Nasser Abbasi, PYS 240. Feb 2, 2002.

Part (b)

Obtain an expression for n!! in terms on n!

n!!= n (n-2) (n-4) (n-6) …

n! =n(n-1)(n-2)(n-3)….

(2n)!! = (2n) (2n-2) (2n-4) (2n-6) ….

(2n)!! = (2n) [2(n-1)] [2(n-2)] [2(n-3)] [2(n-4)] ….

(2n)!! = [2*2*2*2…..] [ n (n-1) (n-2) (n-3) ….]

(2n)!! = (2^n) n!

Log( 2n!!) = log 2^n + log( n!)

Log( 2n!!) = n log 2 + log( n! ) -------- (1)

So, for large number n, we use Sterling equation to find (n)!, then we find 2n!! from the above.

Sterling equation for large n is

Log(n!) = ½ log(2 n pi) + n log(n) –n + log(1+ 1/(12 n) + 1/(288 n^2))

So, for large n, and using equation (1) above

Log( 2n!!) = n log 2 + { ½ log(2 n pi) + n log(n) –n + log(1+ 1/(12 n) + 1/(288 n^2)) }

Let m = 2n

Log( m!! ) = (m/2) log 2 + { ½ log(m pi) + m/2 log( m/2) – m/2 + log(1+ 1/(6 m) + 1/(288 (m/2)^2)) }

The above equation is what I’ll use in part C.

part( C )

Write a program that prints out n!! by evaluating it using Stirling approximation when n>30, compute 10000!!, 314159!!, and (6.02 x 10^23)!!

Source code

function [realPart, powerToTen]= nma_doubleFactorialSterling(n)

% FUNCTION [realPart,powerToTen]=nma_doubleFactorialSterling(n)

% compute double factorial of a number using the Sterling

% approximation.

% INPUT

% n: The number to find double factorial for.

% OUTPUT

% output is returned in the form of [realPart,powerOfTen]

% where n!! = realPart x 10^(powerToTen)

% realPart : the value in the above expression.

% powerToTen : the power to raise 10 to as shown in the above expression.

% This program also prints the final n!! number for

% display purposes.

% by Nasser Abbasi

% HW 3, 1.22. part C. PHY 240, Feb 2, 2002

% design notes:

% Use this equation, derived in part (b)

% Log( n!! ) = (n/2) log 2 + { 1/2 log(n pi) + n/2 log( n/2) - n/2 + log(1+ 1/(6 n) + 1/(288 (n/2)^2)) }

if nargin ~=1

error('missing input argument');

end

if n <= 0

error('positive value is expected');

end

sterlingTerm = sterlingLogNFactorial(n);

v = ((n/2) * log(2) ) + sterlingTerm;

% convert to log base.

% recall, v here is log(n!!)

% so use the conversion log10(x) = log10(e) * log(x)

v = log10(exp(1)) * v;

powerToTen = floor(v);

realPart = v - powerToTen;

% since log10(n!!) = A + B

% so, n!! = 10 ^(A+B) = 10^A * 10^B

% where A is the powerToTen, and 10^B will be the realPart.

realPart = 10^realPart;

disp(sprintf('%d!! = %f x 10^%d',n,realPart,powerToTen));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [v]= sterlingLogNFactorial(n)

% finds sterling term expressed by this equation derived in part (b):

% 1/2 log(n pi) + n/2 log( n/2) - n/2 + log(1+ 1/(6 n) + 1/(288 (n/2)^2)) }

v= (1/2) * log(n*pi);

v= v + ( (n/2) * log(n/2) );

v= v - n/2;

v= v + log( 1 + 1/(6*n) + 1/(288*(n/2)^2) );

OUTPUT

» help nma_doubleFactorialSterling

FUNCTION [realPart,powerToTen]=nma_doubleFactorialSterling(n)

compute double factorial of a number using the Sterling

approximation.

INPUT

n: The number to find double factorial for.

OUTPUT

output is returned in the form of [realPart,powerOfTen]

where n!! = realPart x 10^(powerToTen)

realPart : the value in the above expression.

powerToTen : the power to raise 10 to as shown in the above expression.

This program also prints the final n!! number for

display purposes.

Example 1

» clear all

» [A,B]= nma_doubleFactorialSterling(1000)

1000!! = 3.993984 x 10^1284

A =

3.9940

B =

1284

Example 2

» clear all

» [A,B]= nma_doubleFactorialSterling(2001)

2001!! = 1.928935 x 10^2870

A =

1.9289

B =

2870

Example 3

» clear all

» [A,B]= nma_doubleFactorialSterling(10000)

10000!! = 5.972727 x 10^17830

A =

5.9727

B =

17830

Example 4

» clear all

» [A,B]= nma_doubleFactorialSterling(314159)

314159!! = 5.406120 x 10^795273

A =

5.4061

B =

795273

Example 5

» clear all

» [A,B]= nma_doubleFactorialSterling(6.02*10^23)

6.020000e+023!! = 1.000000 x 10^7.026936e+024

A =

B =

7.0269e+024