'프로그래밍' 카테고리의 글 목록 (23 Page)

감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 Java 언어 소스

프로그래밍/Java 2012. 12. 13. 00:12

Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp, sin, pow 함수들의 구현에 따라 오차가 더 있을 수 있다. C 언어에는 long double 타입이 있어서 좀 더 정밀한 계산을 할 수 있지만. Java 언어에는 그렇지 못하여 그냥 double 타입을 사용하였다. 비교를 위해 아파치의 commons의 Math 라이브러라에서 지원하는 Gamma.logGamma 함수를 사용한 결과도 함께 출력하도록 하였다. (commons 라이브러리는 http://commons.apache.org/math/download_math.cgi 에서 다운로드할 수 있다.

commons 2.2와 common 3.0은 차이가 있으며 서로 호환되지 않는다.

우선 import 구문이 서로 다르다.

import org.apache.commons.math.special.*; // commons 2.2의 경우

import org.apache.commons.math3.special.*; // commons 3.0의 경우

common 2.2의 경우에는 org.apache.commons.math.special.Gamma 클래스에 gamma 라는함수가 있어서 이를 직접 호출하면 되지만, commons 3.0의 경우에는 org.apache.commons.math3.special.Gamma 클래스에 gamma 라는 함수가 없으므로, logGamma 함수와 exp 함수의 합성함수로 계산하였다.

exp(Gamma.logGamma(double)) 함수와 (자체 구현된) gamma2(double) 함수에 의한 계산 결과를 비교하기 바란다.

/**
* Filename: TestLanczos_02.java
*
* Compile: javac -d . -cp .:commons-math3-3.0.jar TestLanczos_02.java
* Execute: java -classpath .:commons-math3-3.0.jar TestLanczos_02
*
* Note: You can download commons-math3-3.0-bin.tar.gz or commons-math3-3.0-bin.zip
* at http://commons.apache.org/math/download_math.cgi
*
* ----------------------------------------------------------------------------
* See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function
* See: http://en.wikipedia.org/wiki/Gamma_function
* See: http://http://jany.tistory.com/53
* See: http://support.realsoftware.com/listarchives/realbasic-nug/2003-08/msg03123.html
* ----------------------------------------------------------------------------
* Date: 2012. 12. 12.
* Copyright (C) Pil Ho, Kim (pkim (AT) scripts.pe.kr)
*/

import java.io.*;
import org.apache.commons.math3.special.Gamma;

public class TestLanczos_02 {

 private static double frac(double x) {
 double result;
 if (x >= 0.0)
 result = Math.abs(x) - Math.floor(Math.abs(x));
 else
 result = -(Math.abs(x) - Math.floor(Math.abs(x)));
 return result;
 }

 public static double gamma2(double y) {
 double pi = 3.1415926535897932385;
 double sq2pi = 2.50662827463100050241577; // sqrt(2Pi)
 double g = 607./128; // best resu'ts when 4<=g<=5
 double t, w, gam, x = y - 1.0;
 int i, cg = 14;
 double[] c = new double[16];

 // Lanczos approximation for the complex plane
 // calculated using vpa digits(256)
 // the best set of coeffs was selected from
 // a solution space of g=0 to 32 with 1 to 32 terms
 // these coeffs really give superb performance
 // of 15 sig. digits for 0<=real(z)<=171
 // coeffs should sum to about g*g/2+23/24

 // http://www.numericana.com/answer/info/godfrey.htm

 c[ 1] = 0.99999999999999709182; // thiBasic arrays start at 1 ?
 c[ 2] = 57.156235665862923517;
 c[ 3] = -59.597960355475491248;
 c[ 4] = 14.136097974741747174;
 c[ 5] = -0.49191381609762019978;
 c[ 6] = 0.33994649984811888699/10000;
 c[ 7] = 0.46523628927048575665/10000;
 c[ 8] = -0.98374475304879564677/10000;
 c[ 9] = 0.15808870322491248884/1000;
 c[10] = -0.21026444172410488319/1000;
 c[11] = 0.21743961811521264320/1000;
 c[12] = -0.16431810653676389022/1000;
 c[13] = 0.84418223983852743293/10000;
 c[14] = -0.26190838401581408670/10000;
 c[15] = 0.36899182659531622704/100000;

 if ( x < 0.0 ) {
 x = -x;
 if (frac(x) == 0.0) {
 gam = Math.pow(10.0, 4932);
 }
 else {
 t = c[1];
 for (i = 1; i <= cg; i++) {
 t = t + c[i+1]/(x+i);
 }
 w = x + g + 0.5;
 gam = Math.pow(w, x+0.5)*Math.exp(-w)*sq2pi*t;
 gam = pi*x/(gam*Math.sin(pi*x));
 }
 }
 else {
 t = c[1];
 for (i = 1; i <= cg; i++) {
 t = t + c[i+1]/(x+i);
 }
 w = x + g + 0.5;
 gam = Math.pow(w, x+0.5)*Math.exp(-w)*sq2pi*t;
 }
 return gam;
 }

 public static void main(String[] args) {
 double y, x = 1.0;

 while (x > 0.0) {
 BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
 String line = null;

 try {
 System.out.printf("Input x(0 to exit): ");
 line = br.readLine();
 x = Double.parseDouble((line));
 }
 catch (IOException ioe) {
 System.err.println("IO error trying to read your input number!");
 System.err.println("Try to input again.");
 }

 if (x == 0.0) {
 break;
 }
 y = gamma2(x);
 System.out.printf(" Gamma(%f) = (%f)! = %30f\n", x, x - 1, y);

 System.out.printf(" Gamma.logGamma(%f) = log((%f)!) = %f\n", x, x - 1, Gamma.logGamma(x));
 System.out.printf(" Math.exp(Gamm.logGamma(%f)) = (%f)! = %30f\n", x, x - 1, Math.exp(Gamma.logGamma(x)));
 }
 }
}

/*
100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Execute: java -classpath .:commons-math3-3.0.jar TestLanczos_02
Result:
Input x(0 to exit): 101
Gamma(101.000000) = (100.000000)! = 93326215443944150000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
Gamma.logGamma(101.000000) = log((100.000000)!) = 363.739376
Math.exp(Gamm.logGamma(101.000000)) = (100.000000)! = 93326215443942250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
Input x(0 to exit): 0
*/

'프로그래밍 > Java' 카테고리의 다른 글

Java 언어로 평방근, 입방근, n제곱근 구하는 함수를 구현하고 테스트하기 (1)	2013.01.11
Java 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.01
Java 와 근사공식을 이용한 PI/sqrt(12) 의 근사값 계산 (0)	2012.01.21
Java 언어로 함수 프로그래밍 미리 해보기 (0)	2010.07.02
스트링 리스트에서 스트링 찾기(find) with Java (0)	2009.04.22

Posted by Scripter

,

감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 C 언어 소스

프로그래밍/C 2012. 12. 12. 18:23

Lanczos 알고리즘은 Stirling 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp, sin, pow 함수들의 구현에 따라 오차가 더 있을 수 있다. 소스를 보면 (조금이라도 유효수자의 개수가 더 많은 계산을 하기 위해) double 타입이 아니라 long double 타입으로 처리하고 있음을 알 수 있다.

long double frac(long double x)

long double gamma2(long double y)

Unix/Linux의 C 컴파일러 gcc 를 사용하면 math 라이브러리에 tgamma 라는 함수가 이미 있다. 이 함수와 (자체 구현된) gamma2 함수에 의한 계산 결과를 비교하기 바란다.

// Filename: testLanczos-02.c
//
// Compile: gcc -o testLanczos-02 -lm testLanczos-02.c
// Execute: ./testLanczos-02
// Output:
// Gamma(10) = (9)! = 362880
// tgamma(10) = (9)! = 362880
//
// See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function
// See: http://en.wikipedia.org/wiki/Gamma_function
// See: http://http://jany.tistory.com/53
// See: http://support.realsoftware.com/listarchives/realbasic-nug/2003-08/msg03123.html

#include <stdio.h>
#include <math.h>

long double frac(long double x)
{
long double result;

if (x >= 0.0)
 result = fabs(x) - floor(x);
else
 result = -(fabs(x) - floor(x));
return result;
}

long double gamma2(long double y)
{
 long double pi = 3.1415926535897932385;
 long double sq2pi = 2.50662827463100050241577; // sqrt(2Pi)
 long double g = 607./128; // best resu'ts when 4<=g<=5
 long double t, w, gam, x = y - 1.0;
 int i, cg = 14;
 long double c[16];

 // Lanczos approximation for the complex plane
 // calculated using vpa digits(256)
 // the best set of coeffs was selected from
 // a solution space of g=0 to 32 with 1 to 32 terms
 // these coeffs really give superb performance
 // of 15 sig. digits for 0<=real(z)<=171
 // coeffs should sum to about g*g/2+23/24

 // http://www.numericana.com/answer/info/godfrey.htm

 c[ 1] = 0.99999999999999709182; // thiBasic arrays start at 1 ?
 c[ 2] = 57.156235665862923517;
 c[ 3] = -59.597960355475491248;
 c[ 4] = 14.136097974741747174;
 c[ 5] = -0.49191381609762019978;
 c[ 6] = 0.33994649984811888699/10000;
 c[ 7] = 0.46523628927048575665/10000;
 c[ 8] = -0.98374475304879564677/10000;
 c[ 9] = 0.15808870322491248884/1000;
 c[10] = -0.21026444172410488319/1000;
 c[11] = 0.21743961811521264320/1000;
 c[12] = -0.16431810653676389022/1000;
 c[13] = 0.84418223983852743293/10000;
 c[14] = -0.26190838401581408670/10000;
 c[15] = 0.36899182659531622704/100000;

 // printf(" y = %LG\n", y);
 // printf(" x = %LG\n", x);

 if ( x < 0.0 ) {
 x = -x;
 if (frac(x) == 0.0) {
 gam = pow(10.0, 4932);
 }
 else {
 t = c[1];
 for (i = 1; i <= cg; i++) {
 t = t + c[i+1]/(x+i);
 }
 w = x + g + 0.5;
 gam = pow(w, x+0.5)*exp(-w)*sq2pi*t;
 gam = pi*x/(gam*sin(pi*x));
 }
 }
 else {
 t = c[1];
 for (i = 1; i <= cg; i++) {
 t = t + c[i+1]/(x+i);
 }
 w = x + g + 0.5;
 // printf(" w = %LG\n", w);
 gam = pow(w, x+0.5)*exp(-w)*sq2pi*t;
 }
 return gam;
}

int main()
{
 long double y, x = 1.0;
 while (x > 0) {
 printf("Input x(0 to exit): ");
 scanf("%LF", &x);
 // printf("x = %LG\n", x);
 if (x == 0.0) {
 break;
 }
 y = gamma2(x);
 printf(" Gamma(%LG) = (%LG)! = %30LF\n", x, x - 1, y);
 printf("tgamma(%LG) = (%LG)! = %30lf\n", x, x - 1, tgamma(x));
 }
 // printf("All done. Press any key to finish\n");
 return 0;
}

/*
100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Execute: ./testLanczos-02
Result:
Input x(0 to exit): 101
Gamma(101) = (100)! = 93326215443944119652178484825626229006662095017918287040779978859652504703806235528941436474245838774007726779438736471616107541007146050712395913562333118464.000000
tgamma(101) = (100)! = 93326215443947553183793338240612302366006877769275129529769934567734921397913148008797092443225880267680360266584467988186500321788890593731926090238149001216.000000
Input x(0 to exit): 0
*/

'프로그래밍 > C' 카테고리의 다른 글

C 언어와 GSL 라이브러리로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.10
C 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.01
리눅스 환경에 MPFR 설치와 테스트 (0)	2012.04.30
long long 타입의 정수 계산을 100경 까지 테스트하는 C 소스 (0)	2012.04.28
ISOC99 를 지원하는 C 컴파일러로 복소수 계산하는 C 예제 (0)	2012.04.28

Posted by Scripter

,

감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 Python 언어 소스

프로그래밍/Python 2012. 12. 12. 17:59

Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp 함수의 구현에 따라 오차가 더 있을 수 있다.

#!/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
"""

'프로그래밍 > Python' 카테고리의 다른 글

Python 언어로 scipy, numpy, mpmath 모듈을 이용하여 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.09
Python 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.01
Python 2.6, 2.6 및 3.2 로 복소수의 포맷(format) 출력 (0)	2012.12.07
pyglet 에 한글 사용하는 예제 (0)	2011.08.26
Python 3.2 를 이용한 간단한 수학식 계산하기 (0)	2011.08.18

Posted by Scripter

,

Python 2.6, 2.6 및 3.2 로 복소수의 포맷(format) 출력

프로그래밍/Python 2012. 12. 7. 11:04

* 리눅스의 Python 2.6 으로 실행한 결과

Python 2.6.5 (r265:79063, Oct 1 2012, 22:04:36)
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> x = complex(1, 2./3)
>>> x
(1+0.66666666666666663j)
>>> "{0.real:.5}+{0.imag:.5}j".format(x)
'1.0+0.66667j'
>>> "{0.real:.5f}{0.imag:+.5f}j".format(x)
'1.00000+0.66667j'
>>> print "{0.real:.5f}{0.imag:+.5f}j".format(x)
1.00000+0.66667j
>>> "{0.real:.5f}{0.imag:+.5f}j".format(complex(1, -2./3))
'1.00000-0.66667j'
>>> print "{0.real:.5f}{0.imag:+.5f}j".format(complex(1, -2./3))
1.00000-0.66667j
>>> format(1+1j, '')
'(1+1j)'
>>> print format(1+1j, '')
(1+1j)

* 윈도우의 Python 2.7 로 실행한 결과 (Python 2.6의 경우와 동일함)

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> x = complex(1, 2./3)
>>> x
(1+0.6666666666666666j)
>>> "{0.real:.5}+{0.imag:.5}j".format(x)
'1.0+0.66667j'
>>> "{0.real:.5f}{0.imag:+.5f}j".format(x)
'1.00000+0.66667j'
>>> print "{0.real:.5f}{0.imag:+.5f}j".format(x)
1.00000+0.66667j
>>> "{0.real:.5f}{0.imag:+.5f}j".format(complex(1, -2./3))
'1.00000-0.66667j'
>>> print "{0.real:.5f}{0.imag:+.5f}j".format(complex(1, -2./3))
1.00000-0.66667j
>>> format(1+1j, '')
'(1+1j)'
>>> print format(1+1j, '')
(1+1j)

* 윈도우의 Python 3.2 로 실행한 결과 (결과는 print() 함수의 사용법 차이 외에는 Python 2.7의 경우와 동일함)

Python 3.2 (r32:88445, Feb 20 2011, 21:29:02) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> x = complex(1, 2./3)
>>> x
(1+0.6666666666666666j)
>>> "{0.real:.5}+{0.imag:.5}j".format(x)
'1.0+0.66667j'
>>> "{0.real:.5f}{0.imag:+.5f}j".format(x)
'1.00000+0.66667j'
>>> print( "{0.real:.5f}{0.imag:+.5f}j".format(x) )
1.00000+0.66667j
>>> print( "{0.real:.5f}{0.imag:+.5f}j".format(complex(1, -2./3)) )
1.00000-0.66667j
>>> format(1+1j, '')
'(1+1j)'
>>> print( format(1+1j, '') )
(1+1j)

'프로그래밍 > Python' 카테고리의 다른 글

Python 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.01
감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 Python 언어 소스 (0)	2012.12.12
pyglet 에 한글 사용하는 예제 (0)	2011.08.26
Python 3.2 를 이용한 간단한 수학식 계산하기 (0)	2011.08.18
Python 2.7 을 이용한 간단한 수학식 계산하기 (0)	2011.08.18

Posted by Scripter

,

높은 버전(1.9.2 이상)의 Ruby 언어에서 (UTF-8 인코딩) 한글 문자열 거꾸로 하기

프로그래밍/Ruby 2012. 12. 6. 09:57

* 아래는 높은 버전(1.9.2 이상)의 Ruby를 사용할 때 적용된다.

UTF--8 인코딩이면 영문이든 한글이든 모두 문자열을 거꾸로 하기가 잘된다.

String.reverse 또는 Strig.reverse() 하면 거꾸로 된 문자열을 얻는다.

예제 1. UTF-8 한글 문자열을 잘 처리하는 예제

# -*- encoding: utf-8 -*-

a = "Hello, world!"
b = "안녕하세요? 아햏햏"
puts "%s --> %s" % [a, a.reverse()]
puts "%s --> %s" % [b, b.reverse()]
####################
# Expected:
#   Hello, world! --> !dlrow ,olleH
#   안녕하세요? 아햏햏 --> 햏햏아 ?요세하녕안
# Result:
#   Hello, world! --> !dlrow ,olleH
#   안녕하세요? 아햏햏 --> 햏햏아 ?요세하녕안
####################

버전 1.8.6 이하의 Ruby 언어에서는 유니코드를 거의 지원하지 않아서 iconv 를 이용하여 문자열 인코딩 변환을 하였지만. 버전 1.9.2 이상의 Ruby 언어에서는 유니코드를 지원하므로 reverse 뿐만 아니라 split를 써도 된다.

다음은 reverse 대신 split를 적용하여 한글은 거꾸로 하는 에제이다.

예제 2. 높은 버전의 Ruby 언어에서 UTF-8 한글 문자열을 분리하여 거꾸로 하는 예제

# -*- encoding: utf-8 -*-

a = "Hello, world!"
b = "안녕하세요? 아햏햏"

# 루비 언어에서 문자열 분리하여 거꾸로 하기
for i in b.split(//u).reverse do
print "%s \n" % i
end
for i in b.upcase.gsub(/([ㄱ-ㅎㅏ-ㅣ가-힣])/u,'\1.').split('.').reverse do
print "%s \n" % i
end

####################
# Expected:
# 햏 
# 햏 
# 아 
# 
# ? 
# 요 
# 세 
# 하 
# 녕 
# 안 
# 햏 
# 햏 
# ? 아 
# 요 
# 세 
# 하 
# 녕 
# 안 
# Result:
# 햏 
# 햏 
# 아 
# 
# ? 
# 요 
# 세 
# 하 
# 녕 
# 안 
# 햏 
# 햏 
# ? 아 
# 요 
# 세 
# 하 
# 녕 
# 안 
####################

'프로그래밍 > Ruby' 카테고리의 다른 글

Ruby 언어로 평방근, 입방근, n제곱근 구하는 함수를 구현하고 테스트하기 (0)	2013.01.11
Ruby 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0)	2013.01.01
Mac OS X Lion 에서 Ruby 1.9.2 설치하여 Tk GUI 사용하기 (0)	2011.08.25
숫자 맞추기 게임 with Ruby (0)	2009.11.05
스트링 리스트에서 스트링 찾기(find) with Ruby (0)	2009.04.22

Posted by Scripter

,

리눅스용 계산기 bc를 이용한 간단한 계산

프로그래밍/BC 2012. 12. 5. 12:18

bc는 유닉스/리눅스가 설치된 한경이면 어디든 있는 계산기 프로그램이다.

bc 사용 메뉴얼

$ man bc

NAME
       bc - An arbitrary precision calculator language

SYNTAX
       bc [ -hlwsqv ] [long-options] [ file ... ]

DESCRIPTION
       bc is a language that supports arbitrary precision numbers with inter‐
       active execution of statements. There are some similarities in the
       syntax to the C programming language.   A standard math library is
       available by command line option. If requested, the math library is
       defined before processing any files. bc starts by processing code from
       all the files listed on the command line in the order listed.   After
       all files have been processed, bc reads from the standard input. All
       code is executed as it is read. (If a file contains a command to halt
       the processor, bc will never read from the standard input.)

       This version of bc contains several extensions beyond traditional bc
       implementations and the POSIX draft standard. Command line options can
       cause these extensions to print a warning or to be rejected. This doc‐
       ument describes the language accepted by this processor.   Extensions

bc - 무한 정확도 계산기 언어
설명
   bc 는 대화형으로 문장을 실행하는 무한 정확도의 숫자를 지원하는언어이며 약간 C 언어와 비슷한 문법을 가지고 있다. 명령행 옵션을 주면 표준 수
   학 라이브러리를 사용할 수 있다. 옵션을 주면 화일들을 처리하기에 앞서 수학 라이브러리가 정의된다. bc 는 우선 명령에서 주어진 화일 순서대 로
   처 리한다. 화일을 모두 처리한 후 bc 는 표준 입력을 읽는다. 모든 코드는 읽는 즉시 실행된다. ( 만약 화일 내의 코드에 처리를 중지하라는명령이
   있다면 bc는 표준입력에서 읽지 않을 것이다. )

bc 현재 버전에서는 전통적인 bc 기능과 POSIX 표준 이외의확장 기능을 포함하고 있다. 명령행 옵션을 주면 확장 기능에 대한 경고 메세지를보여주 고
처리를 무시하게 할 수 있다. 이 문서에서는 GNU 버전의 처리기에서사용하는 언어를 설명한다. 확장 기능도 같이 설명한다.

옵션 -q (welcome 메세지 출력 방지)

옵션 -l (수학 함수 라이브러리 사용) 아래에서는 e(3/0/4) 즉 e^(3/4) 계산에 사용됨

옵션 -ql (-q와 -l을 합친 옵션)

quit; (실행 후 즉시 bc 종료하기)

아래는 이항분포 B(6, 1-e^(3/4)) 에서 확률 P(X>=4) 를 구하는 bc 스크립트 소스이다.

* 소스 파일명: test-01.bc

#!/usr/bin/bc -ql
# Execute: bc -l test-01.bc

define f(x) {
 if (x <= 1)
 return (1);
 return (f(x - 1) * x);
}
define p() {
 return 1.0 - e(-3.0/4.0);
}
define b(x) {
 return (f(6)/f(x)/f(6-x))*(p()^x)*((1-p())^(6-x));
}
1 - (b(0) + b(1) + b(2) + b(3));
quit;

* 실행:
$ bc -ql test-01.bc

또는
$ cat test-01.bc | bc -ql

* 실행 결과:

.39688469988268586698

* 셀 스크립트로 bc를 실행하는 간단한 예제들

$ echo "1+2+3+4+5" | bc -l
15
$ echo "1+1/2+1/3+1/4+1/5" | bc -l
2.28333333333333333333
$ echo "1+1/2+1/3+1/4+1/5" | bc
1
$ echo "scale=5; 1+1/2+1/3+1/4+1/5" | bc
2.28333

$ cat for.txt
for (i=0; i<10; i++) {
print i; print "\n"
}
$ cat for.txt | bc -l
0
1
2
3
4
5
6
7
8
9

* 직각을 낀 두 변의 길이가 0.00000001 인 직각이등변삼각형의 빗변의 길이 구하기

$ echo "sqrt(0.1^8^2 + 0.1^8^2)" | bc -l
0
$ echo "sqrt(1^2 + 1^2)*0.1^8" | bc -l
.00000001414213562373

* 원주율 PI 값 계산하기

$ time echo "scale=10000; 4*a(1)" | bc -l
3.141592653589793238462643383279502884197169399375105820974944592307\
81640628620899862803482534211706798214808651328230664709384460955058\
22317253594081284811174502841027019385211055596446229489549303819644\
28810975665933446128475648233786783165271201909145648566923460348610\
45432664821339360726024914127372458700660631558817488152092096282925\
40917153643678925903600113305305488204665213841469519415116094330572\
70365759591953092186117381932611793105118548074462379962749567351885\
75272489122793818301194912983367336244065664308602139494639522473719\
07021798609437027705392171762931767523846748184676694051320005681271\
45263560827785771342757789609173637178721468440901224953430146549585\
37105079227968925892354201995611212902196086403441815981362977477130\
99605187072113499999983729780499510597317328160963185950244594553469\
08302642522308253344685035261931188171010003137838752886587533208381\
42061717766914730359825349042875546873115956286388235378759375195778\
18577805321712268066130019278766111959092164201989380952572010654858\
63278865936153381827968230301952035301852968995773622599413891249721\
77528347913151557485724245415069595082953311686172785588907509838175\
46374649393192550604009277016711390098488240128583616035637076601047\
10181942955596198946767837449448255379774726847104047534646208046684\
25906949129331367702898915210475216205696602405803815019351125338243\
00355876402474964732639141992726042699227967823547816360093417216412\
19924586315030286182974555706749838505494588586926995690927210797509\
30295532116534498720275596023648066549911988183479775356636980742654\
25278625518184175746728909777727938000816470600161452491921732172147\
72350141441973568548161361157352552133475741849468438523323907394143\
33454776241686251898356948556209921922218427255025425688767179049460\
16534668049886272327917860857843838279679766814541009538837863609506\
80064225125205117392984896084128488626945604241965285022210661186306\
74427862203919494504712371378696095636437191728746776465757396241389\
08658326459958133904780275900994657640789512694683983525957098258226\
20522489407726719478268482601476990902640136394437455305068203496252\
45174939965143142980919065925093722169646151570985838741059788595977\
29754989301617539284681382686838689427741559918559252459539594310499\
72524680845987273644695848653836736222626099124608051243884390451244\
13654976278079771569143599770012961608944169486855584840635342207222\
58284886481584560285060168427394522674676788952521385225499546667278\
23986456596116354886230577456498035593634568174324112515076069479451\
09659609402522887971089314566913686722874894056010150330861792868092\
08747609178249385890097149096759852613655497818931297848216829989487\
22658804857564014270477555132379641451523746234364542858444795265867\
82105114135473573952311342716610213596953623144295248493718711014576\
54035902799344037420073105785390621983874478084784896833214457138687\
51943506430218453191048481005370614680674919278191197939952061419663\
42875444064374512371819217999839101591956181467514269123974894090718\
64942319615679452080951465502252316038819301420937621378559566389377\
87083039069792077346722182562599661501421503068038447734549202605414\
66592520149744285073251866600213243408819071048633173464965145390579\
62685610055081066587969981635747363840525714591028970641401109712062\
80439039759515677157700420337869936007230558763176359421873125147120\
53292819182618612586732157919841484882916447060957527069572209175671\
16722910981690915280173506712748583222871835209353965725121083579151\
36988209144421006751033467110314126711136990865851639831501970165151\
16851714376576183515565088490998985998238734552833163550764791853589\
32261854896321329330898570642046752590709154814165498594616371802709\
81994309924488957571282890592323326097299712084433573265489382391193\
25974636673058360414281388303203824903758985243744170291327656180937\
73444030707469211201913020330380197621101100449293215160842444859637\
66983895228684783123552658213144957685726243344189303968642624341077\
32269780280731891544110104468232527162010526522721116603966655730925\
47110557853763466820653109896526918620564769312570586356620185581007\
29360659876486117910453348850346113657686753249441668039626579787718\
55608455296541266540853061434443185867697514566140680070023787765913\
44017127494704205622305389945613140711270004078547332699390814546646\
45880797270826683063432858785698305235808933065757406795457163775254\
20211495576158140025012622859413021647155097925923099079654737612551\
76567513575178296664547791745011299614890304639947132962107340437518\
95735961458901938971311179042978285647503203198691514028708085990480\
10941214722131794764777262241425485454033215718530614228813758504306\
33217518297986622371721591607716692547487389866549494501146540628433\
66393790039769265672146385306736096571209180763832716641627488880078\
69256029022847210403172118608204190004229661711963779213375751149595\
01566049631862947265473642523081770367515906735023507283540567040386\
74351362222477158915049530984448933309634087807693259939780541934144\
73774418426312986080998886874132604721569516239658645730216315981931\
95167353812974167729478672422924654366800980676928238280689964004824\
35403701416314965897940924323789690706977942236250822168895738379862\
30015937764716512289357860158816175578297352334460428151262720373431\
46531977774160319906655418763979293344195215413418994854447345673831\
62499341913181480927777103863877343177207545654532207770921201905166\
09628049092636019759882816133231666365286193266863360627356763035447\
76280350450777235547105859548702790814356240145171806246436267945612\
75318134078330336254232783944975382437205835311477119926063813346776\
87969597030983391307710987040859133746414428227726346594704745878477\
87201927715280731767907707157213444730605700733492436931138350493163\
12840425121925651798069411352801314701304781643788518529092854520116\
58393419656213491434159562586586557055269049652098580338507224264829\
39728584783163057777560688876446248246857926039535277348030480290058\
76075825104747091643961362676044925627420420832085661190625454337213\
15359584506877246029016187667952406163425225771954291629919306455377\
99140373404328752628889639958794757291746426357455254079091451357111\
36941091193932519107602082520261879853188770584297259167781314969900\
90192116971737278476847268608490033770242429165130050051683233643503\
89517029893922334517220138128069650117844087451960121228599371623130\
17114448464090389064495444006198690754851602632750529834918740786680\
88183385102283345085048608250393021332197155184306354550076682829493\
04137765527939751754613953984683393638304746119966538581538420568533\
86218672523340283087112328278921250771262946322956398989893582116745\
62701021835646220134967151881909730381198004973407239610368540664319\
39509790190699639552453005450580685501956730229219139339185680344903\
98205955100226353536192041994745538593810234395544959778377902374216\
17271117236434354394782218185286240851400666044332588856986705431547\
06965747458550332323342107301545940516553790686627333799585115625784\
32298827372319898757141595781119635833005940873068121602876496286744\
60477464915995054973742562690104903778198683593814657412680492564879\
85561453723478673303904688383436346553794986419270563872931748723320\
83760112302991136793862708943879936201629515413371424892830722012690\
14754668476535761647737946752004907571555278196536213239264061601363\
58155907422020203187277605277219005561484255518792530343513984425322\
34157623361064250639049750086562710953591946589751413103482276930624\
74353632569160781547818115284366795706110861533150445212747392454494\
54236828860613408414863776700961207151249140430272538607648236341433\
46235189757664521641376796903149501910857598442391986291642193994907\
23623464684411739403265918404437805133389452574239950829659122850855\
58215725031071257012668302402929525220118726767562204154205161841634\
84756516999811614101002996078386909291603028840026910414079288621507\
84245167090870006992821206604183718065355672525325675328612910424877\
61825829765157959847035622262934860034158722980534989650226291748788\
20273420922224533985626476691490556284250391275771028402799806636582\
54889264880254566101729670266407655904290994568150652653053718294127\
03369313785178609040708667114965583434347693385781711386455873678123\
01458768712660348913909562009939361031029161615288138437909904231747\
33639480457593149314052976347574811935670911013775172100803155902485\
30906692037671922033229094334676851422144773793937517034436619910403\
37511173547191855046449026365512816228824462575916333039107225383742\
18214088350865739177150968288747826569959957449066175834413752239709\
68340800535598491754173818839994469748676265516582765848358845314277\
56879002909517028352971634456212964043523117600665101241200659755851\
27617858382920419748442360800719304576189323492292796501987518721272\
67507981255470958904556357921221033346697499235630254947802490114195\
21238281530911407907386025152274299581807247162591668545133312394804\
94707911915326734302824418604142636395480004480026704962482017928964\
76697583183271314251702969234889627668440323260927524960357996469256\
50493681836090032380929345958897069536534940603402166544375589004563\
28822505452556405644824651518754711962184439658253375438856909411303\
15095261793780029741207665147939425902989695946995565761218656196733\
78623625612521632086286922210327488921865436480229678070576561514463\
20469279068212073883778142335628236089632080682224680122482611771858\
96381409183903673672220888321513755600372798394004152970028783076670\
94447456013455641725437090697939612257142989467154357846878861444581\
23145935719849225284716050492212424701412147805734551050080190869960\
33027634787081081754501193071412233908663938339529425786905076431006\
38351983438934159613185434754649556978103829309716465143840700707360\
41123735998434522516105070270562352660127648483084076118301305279320\
54274628654036036745328651057065874882256981579367897669742205750596\
83440869735020141020672358502007245225632651341055924019027421624843\
91403599895353945909440704691209140938700126456001623742880210927645\
79310657922955249887275846101264836999892256959688159205600101655256\
375676

real    2m52.269s
user    2m52.200s
sys    0m0.000s

[참고 사이트]

1) http://ydux.tistory.com/27

2) UNIX 기본 명령어 - 강원대학교 컴퓨터

3) 수치계산기로서의 Linux

4) 영문 위키피디어에서 설명하는 bc 프로그래밍

5) Online version of GNU bc

6) Bc for Windows

Posted by Scripter

,

Sage로 행렬의 LU분해 구하기

프로그래밍/Sage 2012. 12. 4. 17:59

sage: A = matrix(CDF, [[1, 2, 3], [4, 5, 6]])
sage: C = transpose(A)*A
sage: P, L, U = C.LU()
sage: P
[0.0 0.0 1.0]
[1.0 0.0 0.0]
[0.0 1.0 0.0]
sage: L
[           1.0            0.0            0.0]
[ 0.62962962963            1.0            0.0]
[0.814814814815            0.5            1.0]
sage: U
[           27.0            36.0            45.0]
[            0.0 -0.666666666667 -1.33333333333]
[            0.0             0.0            -0.0]
sage: L*U
[27.0 36.0 45.0]
[17.0 22.0 27.0]
[22.0 29.0 36.0]
sage: P*C
[27.0 36.0 45.0]
[17.0 22.0 27.0]
[22.0 29.0 36.0]
sage: transpose(P)*L*U
[17.0 22.0 27.0]
[22.0 29.0 36.0]
[27.0 36.0 45.0]
sage: C
[17.0 22.0 27.0]
[22.0 29.0 36.0]
[27.0 36.0 45.0]
sage: transpose(P)*L*U == C
False
sage: transpose(P)*L*U - C
[-3.5527136788e-15               0.0               0.0]
[              0.0               0.0               0.0]
[              0.0               0.0               0.0]
sage: O = matrix(CDF, 3, 3, 0)
sage: transpose(P)*L*U - C == O
False
sage: transpose(P)*L*U - C - O
[-3.5527136788e-15               0.0               0.0]
[              0.0               0.0               0.0]
[              0.0               0.0               0.0]

Posted by Scripter

,

SciLab 으로 행렬의 LU분해 구하기

프로그래밍/SciLab 2012. 12. 4. 17:29

-->A = [1, 2, 3; 4, 5, 6]; C = A'*A;

-->[L, U, E] = lu(C)
E =

    0.    0.    1.
    1.    0.    0.
    0.    1.    0.
U =

    27.    36.          45.
    0.   - 0.6666667 - 1.3333333
    0.     0.           0.
L =

    1.           0.     0.
    0.6296296    1.     0.
    0.8148148    0.5    1.

-->L*U
ans =

    27.    36.    45.
    17.    22.    27.
    22.    29.    36.

-->E*C
ans =

    27.    36.    45.
    17.    22.    27.
    22.    29.    36.

-->E'*L*U
ans =

    17.    22.    27.
    22.    29.    36.
    27.    36.    45.

-->C
C =

    17.    22.    27.
    22.    29.    36.
    27.    36.    45.

-->E'*L*U - C
ans =

   1.0D-14 *

- 0.3552714    0.    0.
    0.           0.    0.
    0.           0.    0.

-->E'*L*U - C == zeros(2, 3)
ans =

F

-->E'*L*U - C == zeros(3, 3)
ans =

F T T
T T T
T T T

-->E'*L*U == C
ans =

F T T
T T T
T T T

Posted by Scripter

,

Jython의 PyServlet에서 한글 파라미터 처리하기

프로그래밍/Jytthon 2012. 11. 29. 21:39

* web.xml 파일에 추가될 내용

<servlet-name>PyServlet</servlet-name>

<servlet-class>org.python.util.PyServlet</servlet-class>

<init-param>

<param-name>python/home</param-name>

<param-value>/opt/usr/local/jython2.5.3</param-value>

</init-param>

</servlet>

<servlet-mapping>

<servlet-name>PyServlet</servlet-name>

<url-pattern>*.py</url-pattern>

</servlet-mapping>

* 수정 전 PyServlet: enter.py

import javax.servlet.http as http class enter(http.HttpServlet): def doPost(self, request, response): session = request.getSession(1) counter = session.getAttribute("counter") try: counter = int(counter) except: # counter is None counter = 0 counter += 1 session.setAttribute("counter", counter) response.contentType = "text/html; charset=utf-8" out = response.outputStream print >> out, """The string you entered is: %s. You have played this game %s times.\n<hr>""" % ( request.getParameter("thing"), counter) self.doGet(request, response) def doGet(self, request, response): response.contentType = "text/html; charset=utf-8" out = response.outputStream print >> out, """Enter a string: <form method="post"> <input type="text" name="thing"> <input type="submit"> </form>"""

* 수정 후 PyServlet: enter.py

# -*- encoding: utf-8 -*- import javax.servlet.http as http class enter(http.HttpServlet): def doPost(self, request, response): session = request.getSession(1) counter = session.getAttribute("counter") try: counter = int(counter) except: # counter is None counter = 0 counter += 1 session.setAttribute("counter", counter) request.setCharacterEncoding("utf-8") response.contentType = "text/html; charset=utf-8" out = response.outputStream print >> out, """The string you entered is: %s. You have played this game %s times.\n<hr>""" % ( request.getParameter("thing").encode("utf-8"), counter) self.doGet(request, response) def doGet(self, request, response): response.contentType = "text/html; charset=utf-8" out = response.outputStream print >> out, """Enter a string: <form method="post"> <input type="text" name="thing"> <input type="submit"> </form>"""

* Mac OS X Lion에서 Tomcat 7 & Jython 2.5.3을 사용하여 enter.py를 실행한 장면