📄 hilbert.cpp
字号:
/*
* Solve set of linear equations involving
* a Hilbert matrix
* i.e. solves Hx=b, where b is the vector [1,1,1....1]
*
* Requires: flash.cpp
*
* Copyright (c) Shamus Software 1988-1997
*/
#include <iostream>
#include "flash.h"
using namespace std;
Miracl precision=20;
static Flash A[20][20];
static Flash b[20];
BOOL gauss(Flash A[][20],Flash b[],int n)
{ /* solve Ax=b using Gaussian elimination *
* solution returned in b */
int i,j,k,m;
BOOL ok;
Flash s;
ok=TRUE;
for (i=0;i<n;i++) A[i][n]=b[i];
for (i=0;i<n;i++)
{ /* Gaussian elimination */
m=i;
for (j=i+1;j<n;j++) if (fabs(A[j][i])>fabs(A[m][i])) m=j;
if (m!=i) for (k=i;k<=n;k++)
{
s=A[i][k];
A[i][k]=A[m][k];
A[m][k]=s;
}
if (A[i][i]==0)
{
ok=FALSE;
break;
}
for (j=i+1;j<n;j++)
{
s=A[j][i]/A[i][i];
for (k=n;k>=i;k--) A[j][k]-=s*A[i][k];
}
}
if (ok) for (j=n-1;j>=0;j--)
{ /* Backward substitution */
s=0;
for (k=j+1;k<n;k++) s+=b[k]*A[j][k];
if (A[j][j]==0)
{
ok=FALSE;
break;
}
b[j]=(A[j][n]-s)/A[j][j];
}
return ok;
}
int main()
{ /* solve set of linear equations */
int i,j,n;
miracl *mip=&precision;
mip->RPOINT=OFF; /* use fractions for output */
do
{
cout << "Order of Hilbert matrix H= ";
cin >> n;
} while (n<2 || n>19);
for (i=0;i<n;i++)
{
b[i]=1;
for (j=0;j<n;j++)
A[i][j]=(Flash)1/(i+j+1);
}
if (gauss(A,b,n))
{
cout << "\nSolution is\n";
for (i=0;i<n;i++) cout << "x[" << i+1 << "] = " << b[i] << "\n";
if (mip->EXACT) cout << "Result is exact!\n";
}
else cout << "H is singular!\n";
return 0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -