第二章数值计算

麦子葡萄糖

已于 2024-01-26 21:15:16 修改

阅读量838

点赞数 24

文章标签：算法数据结构 c语言

于 2024-01-06 22:00:51 首次发布

本文链接：https://blog.csdn.net/2302_79995889/article/details/135432567

版权

2-1 随机数

例题8 均匀分布随机数

⭐️背景：

伪随机数是用确定性的算法计算出来自[0,1]均匀分布的随机数序列，并不真正的随机，但具有类似于随机数的统计特征，如均匀性、独立性等。在计算伪随机数时，若使用的初值（种子）不变，那么伪随机数的数序也不变。伪随机数可以用计算机大量生成，在模拟研究中为了提高模拟效率，一般采用伪随机数代替真正的随机数。有些周期函数的系数可以调整，之后它们的周期非常大，基本上与非周期的函数效果一样。

研究伪随机数的计算有很多，这边先介绍线性同余发生器（Linear congruential generator），简称LCG，是一种能产生具有不连续计算的伪随机序列的分段线性方程的算法，它代表了最古老和最知名的伪随机序列生成器算法之一，其理论相对容易理解，并且易于实现和快速。

使用初值 x_0 ， $x_n=(Ax_{n-1}+C)mod M$ ，A、C、M为合适的整数，可产生0~M范围的值。

若 M=2^n ， Amod8=5 ，C为奇数，那么在0~M-1的范围内值分别会出现一次，周期为M。像这种每个值的出现频度一致的时候叫做均匀分布随机数。

以下的算法里令A=109，C=1021，M=32768， x_0 = 13；

⭐️运行：

#include <stdio.h>
#define A 109
#define C 1021
#define M 32768
/*生成初值*/
unsigned rndnum = 13;
/*生成0~M的随机数*/
unsigned irnd()
{
    rndnum=(rndnum*A+C)%M;
    return rndnum;
}
/*生成100个随机数*/
int main()
{
    int i;
    for(i=0;i<100;i++){
        printf("%8d",irnd());
    }
    printf("\n");
    return 0;
}

运行结果如下：

2438 4619 12972 5945 26434 31511 27848 21797 17598 18659 3236 26065 24058 1903 11840 13629 12022 699 11676 28521 29618 18119 9912 85 10286 8083 30100 5121 2154 6431 13872 5741 4198 32619 17548 13209 31778 24183 15528 22405 18334 579 31364 11825 11994 30415 6688 9117 11734 2075 30588 25545 146 16935 11928 23221 8974 28915 7028 13409 20810 8319 23056 23757 1862 7371 18028 32761 258 29143 31880 2533 14974 27555 22628 9873 28602 5679 30208 16893 7350 15739 12636 2089 32114 28039 9848 25877 3566 29267 12628 1217 2602 22495 28144 21293 28198 27179 14412 31833

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习8-1 均匀性的检定

⭐️背景：

为了检测例题8中产生的随机数的均匀性，我们使用 $\chi^2$ 检定手法来计算。

产生1~M范围内的随机数时，这里让i值产生次数为 f_i 、期望值为 F_i ， $\chi ^2=\sum_{i=1}^{M}\frac{(f_i-F_i)^2}{F_i}$ ， $\chi^2$ 的值越小越能证明均匀性。

令M=10，那么产生1000个随机数时，期望值为1000/10=100。

⭐️运行：

#include <stdio.h>
#define N 1000         //随机数产生次数
#define M 10           //整数随机数的范围
#define F (N/M)        //期望值
#define SCALE (40.0/F) //直方图的高

unsigned rndnum = 13;
/*产生0~32767范围内的随机数*/
unsigned irnd()
{
    rndnum=(rndnum*109+1021)%32768;
    return rndnum;
}
/*产生0~1之间的随机数*/
double rnd()
{
    return irnd()/32767.1;
}
int main()
{
    int i , j ,rank ,hist[M+1];
    double e = 0.0;
    for(i=1;i<=M;i++)
        hist[i]=0;
    /*产生1~M之间的随机数*/
    for(i=0;i<N;i++){
       rank=(int)(rnd()*M+1);
       hist[rank]++;
    }
    for(i=1;i<=M;i++){
        printf("%3d : %3d",i,hist[i]);
        for(j=0;j<=hist[i]*SCALE;j++)
            printf("*");
        printf("\n");
        e=e+(double)(hist[i]-F)*(hist[i]-F)/F;
    }
    printf("x2=%f",e);
    return 0;
}

运行结果如下：

1 : 101*****************************************
2 : 98****************************************
3 : 96***************************************
4 : 107*******************************************
5 : 100*****************************************
6 : 94**************************************
7 : 101*****************************************
8 : 104******************************************
9 : 106*******************************************
10 : 93**************************************
x2=2.080000%

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习8-2 正态分布随机数（Box-Muller）

⭐️背景：

除了之前介绍过的均匀分布随机数，正态分布随机数也一样重要。在这里，使用Box-Muller 算法产生服从正态分布的随机数。基本思想是先得到服从均匀分布的随机数、再将服从均匀分布的随机数转变为服从正态分布。

我们已知m为平均值， $\sigma$ 为标准偏差，正态分布为 $N(m,\sigma ^2)=\frac{1}{\sigma \sqrt{2\pi}}e-\frac{(x-m)^2}{2\sigma ^2}$ 。

在Box-Muller 算法里，我们使用2个均匀分布的随机数r1和r2，可得到：

$x=\sigma \sqrt{-2logr_1}cos2{\pi}r_2+m$

$y=\sigma \sqrt{-2logr_1}sin2{\pi}r_2+m$

⭐️运行：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
void brnd(double sig,double m,double* x,double* y)
{
    double r1 ,r2;
    r1=rand()/(RAND_MAX+0.1);r2=rand()/(RAND_MAX+0.1);//0~1的随机数
    *x=sig*sqrt(-2*log(r1))*cos(2*3.14159*r2)+m;
    *y=sig*sqrt(-2*log(r1))*sin(2*3.14159*r2)+m;
}
int main()
{
    int i,j,hist[100];//设置为100可以更好地覆盖函数生成的随机数的范围
    double x,y;
    for(i=0;i<20;i++)
        hist[i]=0;
    for(i=0;i<1000;i++){
        brnd(2.5,10.0,&x,&y);
        hist[(int)x]++;
        hist[(int)y]++;
    }
    for(i=0;i<=20;i++){
        printf("%3d : I",i);
        for(j=0;j<=hist[i]/10;j++){
            printf("*");
        }
        printf("\n");
    }
    return 0;
}

运行结果如下：

0 : I*
1 : I*
2 : I*
3 : I*
4 : I***
5 : I******
6 : I*************
7 : I**********************
8 : I****************************
9 : I******************************
10 : I********************************
11 : I****************************
12 : I********************
13 : I************
14 : I*******
15 : I***
16 : I**
17 : I*
18 : I*
19 : I*
20 : I**

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-2 数值积分

例题9 定积分（梯形法则）

⭐️背景：

梯形法则是采用梯形来估计曲线下方面积，这等同将被积函数近似为直线函数，被积的部分近似为梯形，要求得较准确的数值，可以将要求积的区间分成多个小区间。

比如，当我们要求a~b范围内f(x)的定积分，我们把[a,b]等间距分成n份，那么定积分近似于n个梯形的面积。在这里，令h=(b-a)/n。

$\int_{a}^{b}f(x)dx=\frac{h}{2}(f(a)+f(a+h))+\frac{h}{2}(f(a+h)+f(a+2h))+\frac{h}{2}(f(a+2h)+f(a+3h))+...+\frac{h}{2}(f(a+(n-1)h)+f(b)) =h(\frac{1}{2}(f(a)+f(b))+f(a+h)+f(a+2h)+...+f(a+(n-1)h))$

⭐️运行：

#include <stdio.h>
#include <math.h>
#define f(x) (sqrt(4-(x)*(x)))
int main()
{
    int k;
    double a , b , n , h , x ,s ,sum;
    printf("请输入积分区间 A , B ? ");
    scanf("%lf %lf",&a,&b);
    n = 50;
    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,%f]区间内,sqrt(4-x*x)的定积分为%f\n",a,b,sum);
    return 0;
}

运行结果如下：

请输入积分区间 A , B ? 0 2
[0.000000,2.000000]区间内,sqrt(4-x*x)的定积分为3.138269

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习9 定积分（辛普森积分法）

⭐️背景：

辛普森法则是一种数值积分方法，以五次曲线逼近的方式取代矩形或梯形积分公式，以求得定积分的数值近似解。

例如，通过 (x_0,f(x_0)) ， (x_1,f(x_1)) ， (x_2,f(x_2)) 3个点的二次曲线的方程式g(x)为：

$g(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2)$

这个g(x)在 x_0 ~ x_2 范围内的积分值可以用以下式子表示。

$\int_{x_0}^{x_2}g(x)dx=\frac{h}{3}(f(x_0)+4f(x_1)+f(x_2))$

令 h=x_1-x_0=x_2-x_1 ，运用到a~b区间内，

$\int_{a}^{b}=\frac{h}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+...+2f(x_{2n-2})+4f(x{2n-1})+f(x_{2n}))=\frac{h}{3}(f(x_0)+f(x_{2n})+4(f(x_1)+f(x_3)+...f(x_{2n-3})+f(x_{2n-1}))+2(f(x_2)+f(x_4)+...+f(x_{2n-2}))$

令 $h=\frac{b-a}{2}$ ， f_0 为奇数项的和， f_e 为偶数项的和

$\int_{a}^{b}f(x)dx=\frac{h}{3}(f(a)+f(b)+4(f_0+f(b-h))+2f_e)$

⭐️运行：

#include <stdio.h>
#include <math.h>
#define f(x) (sqrt(4-(x)*(x)))
int main()
{
    int k;
    double a,b,n,h,fo,fe,sum;//fo为奇数项的和，fe为偶数项的和
    printf("请输入积分区间 A B ?");
    scanf("%lf %lf",&a,&b);
    n=50;
    h=(b-a)/(2*n);
    fo=0,fe=0;
    for(k=1;k<=2*n-3;k=k+2){
        fo=fo+f(a+h*k);
        fe=fe+f(a+h*(k+1));
    }
    sum=h/3*(f(a)+f(b)+4*(fo+f(b-h))+2*fe);
    printf("[%f,%f]区间内,sqrt(4-x*x)的定积分为%f\n",a,b,sum);
    return 0;
}

运行结果如下：

请输入积分区间 A B ?0 2
[0.000000,2.000000]区间内,sqrt(4-x*x)的定积分为3.141133

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-3 泰勒展开

例题10 泰勒展开(exp(x))

⭐️背景：

$e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+...+\frac{x^{k-1}}{(k-1)!}+\frac{x^k}{k!}+...$

<——————d——————>

<————————s————————>

在迭代算法或级数求和中，EPS 用于判断结果是否已经足够接近真实值。当连续两次迭代的结果之差小于 EPS 时，可以认为算法已经收敛，因此可以停止进一步计算。这种方法特别适用于无法准确预测迭代次数就能达到期望精度的情况。EPS 被设定为 1e-08，表示任何小于 0.00000001 的数值在计算中都可以被视为无足轻重的。这在判断指数级数是否已经足够接近其真实值时非常有用，特别是考虑到级数可能需要很多项才能收敛，尤其是对于较大的输入值。

⭐️运行：

#include <stdio.h>
#include <math.h>
double myexp(double x)
{
    double EPS=1e-08;
    double s=1.0,e=1.0,d;
    int k;
    for(k=1;k<=200;k++){
        d=s;
        e=e*x/k;
        s=s+e;
        if(fabs(s-d)<EPS*fabs(d))
            return s;
    }
    return 0.0;
}
int main()
{
    double x;
    printf("   x         myexp(x)       exp(x)\n");
    for(x=0;x<=40;x=x+10)
        printf("%5.1f %14.6g %14.6g\n",x,myexp(x),exp(x));
    return 0;
}

运行结果如下：
x myexp(x) exp(x)
0.0 1 1
10.0 22026.5 22026.5
20.0 4.85165e+08 4.85165e+08
30.0 1.06865e+13 1.06865e+13
40.0 2.35385e+17 2.35385e+17

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习10-1 负值版 $e^{-x}$

⭐️背景：

$e^{-x}=\frac{1}{e^x}$

⭐️运行：

#include <stdio.h>
#include <math.h>
double myexp(double x)
{
    /*定义了一个很小的值作为误差容限，用于判断序列是否收敛*/
    double EPS=1e-08;
    double s=1.0,e=1.0,d,a;
    int k;
    a=fabs(x);
    /*如果在 200 次迭代之内收敛，则返回计算的指数值*/
    for(k=1;k<=200;k++){
        d=s;
        e=e*a/k;
        s=s+e;
        if(fabs(s-d)<EPS*fabs(d)){
            if(x>0)
                return s;
            else
                return 1.0/s;
        }
    }
    return 0.0;
}
int main()
{
    double x;
    printf("   x         myexp(x)       exp(x)\n");
    for(x=-40;x<=40;x=x+10)
        printf("%5.1f %14.6g %14.6g\n",x,myexp(x),exp(x));
    return 0;
}

运行结果如下：
x myexp(x) exp(x)
-40.0 4.24835e-18 4.24835e-18
-30.0 9.35762e-14 9.35762e-14
-20.0 2.06115e-09 2.06115e-09
-10.0 4.53999e-05 4.53999e-05
0.0 1 1
10.0 22026.5 22026.5
20.0 4.85165e+08 4.85165e+08
30.0 1.06865e+13 1.06865e+13
40.0 2.35385e+17 2.35385e+17

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习10-2 cosx

⭐️背景：

$cosx=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^4}{6!}+...$

⭐️运行：

#include <stdio.h>
#include <math.h>
double mycos(double x)
{
    /*定义了一个很小的值作为误差容限，用于判断序列是否收敛*/
    double EPS=1e-08;
    double s=1.0,e=1.0,d;
    int k;
    x=fmod(x,2*3.14159265358979);//角度归一化:使用fmod函数将x规范化到[0,2pi]范围内
    /*如果在 200 次迭代之内收敛，则返回计算的指数值*/
    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;
}
int main()
{
    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;
}

运行结果如下：

x mycos(x) cos(x)
0.0 1.000000 1.000000
10.0 0.984808 0.984808
20.0 0.939693 0.939693
30.0 0.866026 0.866026
40.0 0.766045 0.766045
50.0 0.642788 0.642788
60.0 0.500001 0.500001
70.0 0.342021 0.342021
80.0 0.173649 0.173649
90.0 0.000001 0.000001
100.0 -0.173647 -0.173647
110.0 -0.342019 -0.342019
120.0 -0.499998 -0.499998
130.0 -0.642786 -0.642786
140.0 -0.766043 -0.766043
150.0 -0.866024 -0.866024
160.0 -0.939692 -0.939692
170.0 -0.984807 -0.984807
180.0 -1.000000 -1.000000

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-4 非线性方程的解法

例题11 2分法

⭐️背景：

1次方程以外的方程叫做非线性方程，其中一种解非线性方程的方法叫做2分法。做法如下：

①解的左右设置两点a和b分别作为low和high的初始值；

②求low和high的中点x的值，x=(low+high)/2；

③若f(x)>0，说明根在x左侧，所以令high=x；

若f(x)<0，说明根在x右侧，所以令low=x；

④当f(x)为0或者|high-low|/|low|<EPS的时候x为解，否则重复②以后的步骤。

本题我们尝试用利用2分法求方程式 $x^3-x+1=0$ 的解。

⭐️运行：

#include <stdio.h>
#include <math.h>
#define f(x) (x)*(x)*(x)-(x)+1
#define EPS 1e-8
#define LIMIT 50
int main()
{
    double low,high,x;
    int k;
    low=-2.0,high=2.0;
    for(k=1;k<=LIMIT;k++){
        x=(low+high)/2;
        if(f(x)>0)
            high=x;
        if(f(x)<0)
            low=x;
        if(f(x)==0||fabs(high-low)<fabs(low)*EPS){
            printf(" x = %f\n",x);
            break;
        }
    }
    if(k>LIMIT)
        printf("在此范围内无解\n");
    return 0;
}

运行结果如下：

x = -1.324718

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习11 牛顿法

⭐️背景：

牛顿法的概念如下：

①令根附近的适当的值 x_0 为初始值；

②连接y=f(x)在x= x_0 上的接线，与x轴的交点设为 x_1 ，按照同样手法求出 $x_2,x_3,...,x_{n-1},x_n$ ；

③求出满足 $\frac{|x_n-x_{n-1}|}{|x_{n-1}|}<EPS$ 的 x_n 的值。

利用 $f'(x_{n-1})\approx \frac{f(x_{n-1})}{x_{n-1}-x_n}$ 的关系式，可以求出 $x_n=x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}$ 。其实，牛顿法会比2分法收敛的速度更快。

⭐️运行：

#include <stdio.h>
#include <math.h>
#define f(x) ((x)*(x)*(x)-(x)+1)
#define g(x) (3*(x)*(x)-1)
#define EPS 1e-8
#define LIMIT 50

int main()
{
    double x=-2.0,dx;
    int k;
    for(k=1;k<=LIMIT;k++){
        dx=x;
        x=x-f(x)/g(x);
        if(fabs(x-dx)<fabs(dx)*EPS){
            printf(" x = %f\n",x);
            break;
        }
    }
    if(k>LIMIT)
        printf("不收敛");
    return 0;
}

运行结果如下：

x = -1.324718

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-5 插值

例题12 拉格朗日插值法

⭐️背景：

当给若干个点 $(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})$ ，可以通过拉格朗日插值求出通过这些点的插值多项式。通过这些点的函数f(x)可以表示为：

$f(x)=\sum_{i=0}^{n-1}(\prod_{j=0}^{n-1}\frac{x-x_j}{x_i-x_j})y_i$ （不包含i=j的项）

⭐️运行：

#include <stdio.h>
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;
}
int main()
{
    double x[]={0.0,1.0,3.0,6.0,7.0};
    double y[]={0.8,3.1,4.5,3.9,2.8};
    double t;
    printf("    x       y\n");
    for(t=0;t<=7.0;t=t+.5)
        printf("%7.2f%7.2f\n",t,lagrange(x,y,5,t));
    return 0;
}

运行结果如下：

x y
0.00 0.80
0.50 2.15
1.00 3.10
1.50 3.74
2.00 4.14
2.50 4.38
3.00 4.50
3.50 4.54
4.00 4.53
4.50 4.48
5.00 4.37
5.50 4.19
6.00 3.90
6.50 3.46
7.00 2.80

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习12 牛顿插值法

⭐️背景：

$f(x)=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_{n-1}(x-x_0)(x-x_1)..(x-x_{n-2})$

这是牛顿插值多项式，也可以使用Honer方法（练习1-1）变形。

$f(x)=(...((a_{n-1}(x-x_{n-2})+a_{n-2})(x-x_{n-3})+a_{n-3})(x-x_{n-4})...+a_1)(x-x_0)+a_0$

⭐️运行：

#include <stdio.h>
double newton(double x[],double y[],int n,double t)
{
    static int flag=1;
    static double a[100];
    double w[100],s;
    int i,j;
    if(flag==1){
        for(i=0;i<n;i++){
            w[i]=y[i];
            for(j=i-1;j>=0;j--)
                w[j]=(w[j+1]-w[j])/(x[i]-x[j]);
            a[i]=w[0];
        }
        flag=-1;
    }
    s=a[n-1];
    for(i=n-2;i>=0;i--)
        s=s*(t-x[i])+a[i];
    return s;
}
int main()
{
    double x[]={0.0,1.0,3.0,6.0,7.0};
    double y[]={0.8,3.1,4.5,3.9,2.8};
    double t;
    printf("    x       y\n");
    for(t=0;t<=7.0;t=t+.5)
        printf("%7.2f%7.2f\n",t,newton(x,y,5,t));
    return 0;
}

运行结果如下：

x y
0.00 0.80
0.50 2.15
1.00 3.10
1.50 3.74
2.00 4.14
2.50 4.38
3.00 4.50
3.50 4.54
4.00 4.53
4.50 4.48
5.00 4.37
5.50 4.19
6.00 3.90
6.50 3.46
7.00 2.80

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-6 多位数运算

例题13 多位数的加减法

⭐️背景：

当处理不能用long型和double型表示的加减法运算时，我们可以从个位数开始每4位数区分，分别用a[]，b[]容纳。加法的结果用c[]容纳。如果结果不满10000，那么直接放入c[]。如果超过10000，那么将结果减去10000放入c[]，下一位结果+1。减法也同样处理，如果结果是负的，那么结果加上10000放入c[]，下一位结果-1。

⭐️运行：

#include <stdio.h>
#define KETA 20 //定义位数
#define N ((KETA-1)/4+1)

//加法
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]-1;
        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");
}
int main()
{
    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;
}

运行结果如下：

21111111111111111110
18877777444333338887

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

练习13 多位数的乘除法

⭐️背景：

计算long类型和short类型的乘法以及除法，这里的short类型指0~32767的整数。

⭐️运行：

#include <stdio.h>
#define KETA 20 //定义位数
#define N ((KETA-1)/4+1)
//乘法
void lmul(short a[],short b,short c[])
{
    short i;long d,cy=0;
    for(i=N-1;i>=0;i--){
        d=a[i];
        c[i]=(d*b+cy)%10000;
        cy=(d*b+cy)/10000;
    }
}
//除法
void ldiv(short a[],short b,short c[])
{
    short i;long d,rem=0;
    for(i=0;i<N;i++){
        d=a[i];
        c[i]=(short)((d+rem)/b);
        rem=(d+rem)%b*10000;
    }
}
//打印
void print(short c[])
{
    short i;
    for(i=0;i<N;i++)
        printf("%04d",c[i]);
    printf("\n");
}
int main()
{
    short a[N+2]={0,3050,2508,8080,1233},
          c[N+2];
    lmul(a,101,c);print(c);
    ldiv(a,200,c);print(c);
    return 0;
}

运行结果如下：

00308075338960924533
00000015251254404006

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

2-7 多位数计算

例题14 e的1000位计算

⭐️背景：

我们已知e可通过泰勒展开求出近似值，那么也可以使用2-6的除法与加法算出每一项以及和。

⭐️运行：

#include <stdio.h>
#define L 1000//要求到的位数
#define L1 ((L/4)+1)//需要分配的大小
#define L2 (L1+1)//多余1个
#define N 451 //449！<10exp(1000)<451!，为了富余一些算到第451项
//加法
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){
            c[i]=c[i]-10000;
            cy=1;
        }
        else
            cy=0;
    }
}
//除法
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]=(short)((d+rem)/b);
        rem=(d+rem)%b*10000;
    }
}
//打印
void printresult(short c[])
{
    short i;
    printf("%3d. ",c[0]);
    for(i=0;i<L1;i++)
        printf("%04d",c[i]);
    printf("\n");
}
int main()
{
    short s[L2+2],w[L2+2];
    short k;
    for(k=0;k<=L2;k++)
        s[k]=w[k]=0;
    s[0]=w[0]=1;
    for(k=1;k<=N;k++){
        ldiv(w,k,w);
        ladd(s,w,s);
    }
    printresult(s);
    return 0;
}

⭐️运行： 2.00027182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354

☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:* ☆.。.:*・°☆.。.:*・°☆.。.:*・°☆.。.:*

麦子葡萄糖

关注

24
点赞
踩
18

收藏

觉得还不错? 一键收藏
0
评论
第二章数值计算

伪随机数是用确定性的算法计算出来自[0,1]均匀分布的随机数序列，并不真正的随机，但具有类似于随机数的统计特征，如均匀性、独立性等。在计算伪随机数时，若使用的初值（种子）不变，那么伪随机数的数序也不变。（Linear congruential generator），简称LCG，是一种能产生具有不连续计算的伪随机序列的分段线性方程的算法，它代表了最古老和最知名的伪随机序列生成器算法之一，其理论相对容易理解，并且易于实现和快速。，C为奇数，那么在0~M-1的范围内值分别会出现一次，周期为M。
复制链接

扫一扫