초심(初心)으로

/*
 * ------------------------------------
 *     사다리꼴 공식에 의한 정적분    *
 * ------------------------------------
 */

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

#define f(x) (sqrt(4-(x)*(x)))   /* 피적분함수 */

int main(void)
{
    int k;
    double a,b,n,h,x,s,sum;

    printf("적분구간 A,B ? ");
    scanf("%lf %lf",&a,&b);

    n=50;                /* a ~ b사이의 구간수 */
    h=(b-a)/n;        /* 구간폭 */
    x=a; s=0;
    for (k=1;k<=n-1;k++)
 {
        x=x+h;
        s=s+f(x);
 }
    sum=h*((f(a)+f(b))/2+s);
    printf("   /%f\n",b);
    printf("   |  sqrt(4-x*x) =%f\n",sum);
    printf("   /%f\n",a);

    return 0;
}

// 함수 f(x)의 정적분을 해석적으로 수식화해서 구하는 것이 아니라 구간을 아주 작게 분할

해서 근사값을 구하는 방법을 '수치적분'이라 한다.

a 와 b 구간을 n개의 아주 작은 구간으로 분할해서 각 구간 면적의 합계를 구한다.

/*
 * --------------------------------
 *      테일러 전개 (cos (x))     *
 * --------------------------------
 */

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

double mycos(double);

int main(void)
{
    double x,rd=3.14159/180;
    printf("    x      mycos(x)        cos(x)\n");
    for (x=0;x<=180;x=x+10)
        printf("%5.1f%14.6f%14.6f\n",x,mycos(x*rd),cos(x*rd));

 return 0;
}
double mycos(double x)
{
    double EPS=1e-08;
    double s=1.0,e=1.0,d=1.0;
    int k;

    x=fmod(x,2*3.14159265358979);      /* x값을 0 ~ 2π에 맞춘다. */
    for (k=1;k<=200;k=k+2) {
        d=s;
        e=-e*x*x/(k*(k+1));
        s=s+e;
        if (fabs(s-d)<EPS*fabs(d))     /* 중단오차 */
            return s;
    }
    return 9999.0;                     /* 수렴하지 않을 때 */
}

// 테일러 전개는 중심에 가까울수록 좋은 근사를 나타내므로 x값은 0~2pi의 범위에 속하도록

도를 라이안(radian)으로 변환해서 계산한다.

/*
 * ------------------
 *      이분법      *
 * ------------------
 */

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

#define f(x) ((x)*(x)*(x)-(x)+1)
#define EPS 1e-8                   /* 중단오차 */
#define LIMIT 50                    /* 중단회수 */

int main(void)
{
    double low,high,x;
    int k=1;

    low=-2.0; high=2.0;
    for (k=1;k<=LIMIT;k++) {
        x=(low+high)/2;
        if (f(x)>0)
            high=x;
        else
            low=x;                // 수렴조건으로 Abs(f(x))<EPS를 사용해도 된다.
        if (f(x)==0 || fabs(high-low)<fabs(low)*EPS){
            printf("x=%f\n",x);
            break;
        }
    }
    if (k>LIMIT)
        printf("수렴하지 않는다.\n");

 return 0;
}

// 1차 방정식(그래프 내에서 직선)외의 방정식을 비선형 방정식이라 한다. 이와 같은 방정식의

근을 구하는 방법에는 이분법이 있다. 이분법에서는 데이터 범위를 반으로 나누어 근이 어는 쪽에

있는지, 조사를 반복해서 조사범위를 근을 향해서 점차 좁혀간다.

/*
 * -----------------------
 *      Lagrange보간     *
 * -----------------------
 */

#include <stdio.h>

double lagrange(double *,double *,int,double);

int main(void)
{
    static double x[]={0.0,1.0,3.0,6.0,7.0},
                  y[]={0.8,3.1,4.5,3.9,2.8};

    double t;

    printf("      x      y\n");
    for (t=0.0;t<=7.0;t=t+.5)
        printf("%7.2f%7.2f\n",t,lagrange(x,y,5,t));

 return 0;
}
double lagrange(double x[],double y[],int n,double t)
{
    int i,j;
    double s,p;

    s=0.0;
    for (i=0;i<n;i++){
        p=y[i];
        for (j=0;j<n;j++){
            if (i!=j)
                p=p*(t-x[j])/(x[i]-x[j]);
        }
        s=s+p;
    }
    return s;
}

//  Lagrange보간   :  몇 쌍의 x, y 데이터가 주어졌을 때, 이 점들을 지나는 보간 다항식을

  Lagrange보간으로 구하고 주어진 점 이외의 점 데이터를 구한다. n-1차의 다항식이 된다.

/*
 * ---------------------------
 *      무한 자리수 계산     *
 * ---------------------------
 */

#include <stdio.h>

#define CIPHER 20                  /* 자리수 */
#define N ((CIPHER-1)/4+1)     /* 배열크기 */

void ladd(short *,short *,short *);
void lsub(short *,short *,short *);
void print(short *);

int main(void)
{
    static short a[N+2]={1999,4444,7777,2222,9999},
                 b[N+2]={ 111,6666,3333,8888,1111},
                 c[N+2];

    ladd(a,b,c); print(c);
    lsub(a,b,c); print(c);

 return 0;
}
void ladd(short a[],short b[],short c[])    /* 긴 자리수 덧셈 */
{
    short i,cy=0;
    for (i=N-1;i>=0;i--){
        c[i]=a[i]+b[i]+cy;
        if (c[i]<10000)
            cy=0;
        else {
            c[i]=c[i]-10000;
            cy=1;
        }
    }
}
void lsub(short a[],short b[],short c[])  /* 긴 자리수 뺄셈 */
{
    short i,brrw=0;
    for (i=N-1;i>=0;i--){
        c[i]=a[i]-b[i]-brrw;
        if (c[i]>=0)
            brrw=0;
        else {
            c[i]=c[i]+10000;
            brrw=1;
        }
    }
}
void print(short c[])        /* 긴 자리수 출력 */
{
    short i;
    for (i=0;i<N;i++)
        printf("%04d ",c[i]);
    printf("\n");
}

// 긴 자리수의 덧셈, 뺄셈 : n자리의 수끼리 덧셈, 뺄셈을 수행한다.

a[]b[]의 하위자리수부터 각각 더하되 계산결과가 10000미만이면 그대로 c[]에 저장하고 10000

이상(네 자리를 넘을 경우)이면 결과 값에서 10000을 뺀 값을 저장하고, 네자리를 계산할 때

자리올림수 1을 더한다. 뺄셈도 같은 방법으로 수행하되 결과가 음수인 경우에는 10000을 더해서

c[]에 저장하고 네 자리를 계산할때 1을 뺀다.

/*
 * -----------------------------------
 *      팩토리얼(factorial) 계산     *
 * -----------------------------------
 */

#include <stdio.h>

void ladd(short *,short *,short *);
void lsub(short *,short *,short *);
void lmul(short *,short,short *);
void printresult(short *);

#define L 64               /* 구하고자 하는 자리수      */
#define L2 ((L+3)/4)       /* 배열크기    */

int main(void)
{
    static short s[L2];
    short k;
    for (k=0;k<L2;k++)
        s[k]=0;

    s[L2-1]=1;
    for (k=1;k<=49;k++){
        lmul(s,k,s);
        printf("%2d!=",k);
        printresult(s);
    }

 return 0;
}
void lmul(short a[],short b,short c[])    /* 긴 자리수 곱셈 */
{
    short i;long d,cy=0;
    for (i=L2-1;i>=0;i--){
        d=a[i];
        c[i]=(d*b+cy)%10000;
        cy=(d*b+cy)/10000;
    }
}

void printresult(short c[])       /* 결과 출력 */
{
    short i;
    for (i=0;i<L2;i++)
        printf("%04d",c[i]);
    printf("\n");
}
void ladd(short a[],short b[],short c[])
{
    short i,cy=0;
    for (i=L2-1;i>=0;i--){
        c[i]=a[i]+b[i]+cy;
        if (c[i]<10000)
            cy=0;
        else {
            c[i]=c[i]-10000;
            cy=1;
        }
    }
}
void lsub(short a[],short b[],short c[])
{
    short i,brrw=0;
    for (i=L2-1;i>=0;i--){
        c[i]=a[i]-b[i]-brrw;
        if (c[i]>=0)
            brrw=0;
        else {
            c[i]=c[i]+10000;
            brrw=1;
        }
    }
}

//  1!부터 49!까지의 값을 64자리까지 구한다.

/*
 * -------------------------------
 *      소수점 1000자리의 π     *
 * -------------------------------
 */

#include <stdio.h>

void ladd(short *,short *,short *);
void lsub(short *,short *,short *);
void ldiv(short *,short,short *);
void printresult(short *);

#define L 1000                    /* 구하고자 하는 자리수     */
#define L1 ((L/4)+1)              /* 배열크기   */
#define L2 (L1+1)                 /* 배열크기 + 1 */
#define N (short)(L/1.39794+1)    /* 계산하려는 항의 수   */

int main(void)
{
    static short s[L2+2],w[L2+2],v[L2+2],q[L2+2];
    short k;
    for (k=0;k<=L2;k++)
        s[k]=w[k]=v[k]=q[k]=0;

    w[0]=16*5; v[0]=4*239;            /* 마친의 공식 */
    for (k=1;k<=N;k++){
        ldiv(w,25,w);
        ldiv(v,239,v);ldiv(v,239,v);
        lsub(w,v,q);ldiv(q,2*k-1,q);
        if ((k%2)!=0)                 /* 홀수항인지 짝수항인지 판정 */
            ladd(s,q,s);
        else
            lsub(s,q,s);
    }
    printresult(s);

 return 0;
}
void printresult(short c[])      /* 결과 출력 */
{
    short i;
    printf("%3d. ",c[0]);        /* 최상위 자리 출력 */
    for (i=1;i<L1;i++)
        printf("%04d ",c[i]);
    printf("\n");
}
void ladd(short a[],short b[],short c[])    /* 긴 자리수 덧셈 */
{
    short i,cy=0;
    for (i=L2;i>=0;i--){
        c[i]=a[i]+b[i]+cy;
        if (c[i]<10000)
            cy=0;
        else {
            c[i]=c[i]-10000;
            cy=1;
        }
    }
}
void lsub(short a[],short b[],short c[])    /* 긴 자리수 뺄셈 */
{
    short i,brrw=0;
    for (i=L2;i>=0;i--){
        c[i]=a[i]-b[i]-brrw;
        if (c[i]>=0)
            brrw=0;
        else {
            c[i]=c[i]+10000;
            brrw=1;
        }
    }
}
void ldiv(short a[],short b,short c[])      /* 긴 자리수 나눗셈 */
{
    short i;long d,rem=0;
    for (i=0;i<=L2;i++){
        d=a[i];
        c[i]=(d+rem)/b;
        rem=((d+rem)%b)*10000;
    }
}

//  pi값을 소수점 이하 1000자리까지 정확하게 구현한다.

부호는 n이 홀수면 양, 짝수면 음이 된다.