Script Programming(스크립트 프로그래밍)

#!/usr/bin/env python
# -*- encoding:utf-8 -*-

# Filename: testLanczos-01.py
#
#           An approximation for the gamma function by using the Lanczos algorithm
#
# Execute: python testLanczos-01.py
#         or
# Execute: ./testLanczos-01.py
#
# See: http://en.wikipedia.org/wiki/Lanczos_approximation
# See:http://www-history.mcs.st-and.ac.uk/Biographies/Lanczos.html
# See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function

from cmath import *
 
# Coefficients used by the GNU Scientific Library
g = 7
p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
     771.32342877765313, -176.61502916214059, 12.507343278686905,
     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
 
def gamma(z):
    z = complex(z)
    # Reflection formula
    if z.real < 0.5:
        return pi / (sin(pi*z)*gamma(1-z))
    else:
        z -= 1
        x = p[0]
        for i in range(1, g+2):
            x += p[i]/(z+i)
        t = z + g + 0.5
        return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x

def factorial(n):
    if n < 2:
        return 1
    else:
        k = 1
        if n % 2 == 0:
            for i in xrange(n/2):
                k *= (i + 1) * (n - i)
        else:
            for i in xrange(n/2):
                k *= (i + 1) * (n - 1 - i)
            k *= n
        return k

def facto(n):
    if n < 2:
        return 1
    else:
        k = 1
        for i in xrange(2, n + 1):
                k *= i
        return k

if __name__ == "__main__":
    print "gamma(10) = 9! = %s asymptotically" % gamma(10)
    # print "gamma(101) = 100! = %16s asymptotically" % gamma(101)
    print "gamma(101) = 100! = %16s asymptotically" % "{0.real:.15f}{0.imag:+.5f}j".format(gamma(101))

    for i in range(11):
        print "%d! = %d" % (i, factorial(i))

    i = 100
    print "factorial(%d) = %d! = %d" % (i, i, factorial(i))
    print "facto(%d) = %d! = %d" % (i, i, facto(i))

"""
Output:
gamma(10) = 9! = (362880+0j) asymptotically
gamma(101) = 100! = 9.33262154439379e+157+0.00000j asymptotically
0! = 1
1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
factorial(100) = 100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
facto(100) = 100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
"""

$ gnuplot

    G N U P L O T
    Version 4.2 patchlevel 6
    last modified Sep 2009
    System: Linux 2.6.32-33-generic

    Copyright (C) 1986 - 1993, 1998, 2004, 2007 - 2009
    Thomas Williams, Colin Kelley and many others

    Type `help` to access the on-line reference manual.
    The gnuplot FAQ is available from http://www.gnuplot.info/faq/

    Send bug reports and suggestions to <http://sourceforge.net/projects/gnuplot>


Terminal type set to 'wxt'
gnuplot> plot [-5:5] [-5:5] gamma(x)

gnuplot> set xlabel "x values"
gnuplot> set ylabel "y values"
gnuplot> set title "Gamma function"
gnuplot> replot

gnuplot> set xtics 1
gnuplot> set ytics 1
gnuplot> set grid
gnuplot> replot

gnuplot> set xzeroaxis linetype 3 linewidth 0.8
gnuplot> set yzeroaxis linetype 3 linewidth 0.8
gnuplot> replot

gnuplot> set samples 1600
gnuplot> replot

/*
 *  Filename: recFactGMPcpp01.cc
 *
 *  Compile: g++ -o recFactGMPcpp01 recFactGMPcpp01.cc -lgmpxx -lgmp
 */

#include <iostream>
#include <gmpxx.h>

using namespace std;

mpz_class factorial(int n) {
    static mpz_class v;
    v = 1;
    for (int i = 1; i <= n; i++) {
     v = v*i;
    }
    return v;
}

int main(int argc, char *argv[]) {
    int n1 = 9;
    int n2 = 30;
    for (int i = n1; i <= n2; i++) {
        cout << i << "! = " << factorial(i) << endl;
        if (i == 13)
            cout << "-- below calc are regular ----" << endl;
    }
    return 0;
}

#include <stdio.h>
#include <gmp.h>

void factorial(mpz_t v, int n) {
 int i;
    mpz_t t;
    mpz_init (t);
    mpz_init_set_str (v, "1", 0);
    for (i = 1; i <= n; i++) {
        mpz_init_set_si (t, i);
        mpz_mul(v, v, t);
    }
}

int main(int argc, char *argv[]) {
    int n1 = 9;
    int n2 = 30;
    int i;
    mpz_t v;
    mpz_init (v);
    for (i = n1; i <= n2; i++) {
        mpz_init_set_si (v, 1);
        factorial(v, i);
        gmp_printf("%d! = %Zd\n", i, v);
        if (i == 13)
        printf("-- below calc are regular ----\n");
    }
    return 0;
}

    Compile: ghc fact.hs
    Execute: main

    Date: 2010/07/20
    Author: phkim  pkim (AT) scripts.pe.kr
-}

module Main where

  factorial :: Integer -> Integer
  factorial n =  recfact 1 n

  recfact :: Integer -> Integer -> Integer
  recfact acc k = case k of
        0 -> acc
        k-> recfact (acc*k) (k - 1)

  fib :: Integer -> Integer
  fib n = fibGen 0 1 n

  fibGen :: Integer -> Integer -> Integer -> Integer
  fibGen a b n = case n of
        0 -> a
        n -> fibGen b (a + b) (n - 1)

  main :: IO ()
  main = do
        print ("Factorial(30000) has " ++ show (length (show (factorial 30000))) ++ " digits")
        print ("Fibonacci(200000) has " ++ show (length (show (fib 200000))) ++ " digits")

{-
    Expected result:
    "Factorial(30000) has 121288 digits"
    "Fibonacci(200000) has 41798 digits"
-}