二元域的高斯消元法C语言

发布网友发布时间：2022-08-14 08:16

共2个回答

热心网友时间：2024-10-27 12:42

// 自己写的，二元域没啥挑战性，还是弄正常的
// 有逆矩阵时返回，无时返回NULL

/*
    matrix.h
*/
#ifndef HEAD_matrix
#define HEAD_matrix

typedef struct matrix *Matrix;

Matrix createMatrix(long double *data, unsigned long length);
Matrix createIdentityMatrix(unsigned long length);
Matrix plicateMatrix(const Matrix trx);
void freeMatrix(Matrix trx);
Matrix caculateInversedMatrix(const Matrix trx);
const char *allocateDiscriptionOfMatrix(Matrix trx);

#endif

/*
    matrix.c
*/
#include <stdlib.h>
#include <assert.h>
#include "matrix.h"

#define __nullable
#define __notnull

#define lineOfMatrix(trx, x) (trx->matrixData+trx->length*(x))
#define elementOfMatrix(trx, x, y) (lineOfMatrix(trx, x)[y])
#define BytesOfMatrixSize (sizeof(long double)*length*length)

struct matrix
{
unsigned long length;
long double matrixData[];
};

static void manipulateMatrixMutiplyK(Matrix trx, unsigned long lineIndex, long double __notnull k);
static void manipulateMatrixLineXMinusKTimesOfY(Matrix trx, unsigned long xLineIndex, long double __notnull k, unsigned long yLineIndex);
static unsigned long matrixLineNoneZeroCount(const Matrix trx, unsigned long columnIndex, unsigned long searchFromLineIndex);
static void manipulateMatrixSwapLineXY(Matrix trx, unsigned long xLineIndex, unsigned long yLineIndex);

const wchar_t *allocateDiscriptionOfMatrix(Matrix trx)
{
FILE *tempFilePointer = tmpfile();
static const wchar_t *format = L"Matrix %lux%lu\n";
unsigned long length = trx->length;
unsigned long stringLength = fwprintf(tempFilePointer, format, length, length);
for(unsigned long lineIndex = 0; lineIndex<length; lineIndex++)
{
//wcscat(result, L"[");
for(unsigned long elementIndexOfLine = 0; elementIndexOfLine<length; elementIndexOfLine++)
if(elementIndexOfLine+1==length)
stringLength += fwprintf(tempFilePointer, L"%5.4Lg", elementOfMatrix(trx, lineIndex, elementIndexOfLine));
else
stringLength += fwprintf(tempFilePointer, L"%5.4Lg, ", elementOfMatrix(trx, lineIndex, elementIndexOfLine));
//if(lineIndex+1 == length)
//wcscat(result, L"]");
//else
//wcscat(result, L"]\n");
}
stringLength += length*3-1;
char *result;
if((result = malloc(stringLength+1))!=NULL)
{
unsigned long indexOfResultOutput = 0;
indexOfResultOutput += swprintf(result, stringLength+1, format, length, length);
for(unsigned long lineIndex = 0; lineIndex<length; lineIndex++)
{
memcpy(result+indexOfResultOutput, L"[", sizeof(wchar_t)*2);
indexOfResultOutput += 2;
for(unsigned long elementIndexOfLine = 0; elementIndexOfLine<length; elementIndexOfLine++)
if(elementIndexOfLine+1==length)
indexOfResultOutput += swprintf(result+indexOfResultOutput, stringLength+1-indexOfResultOutput, L"%5.4Lg", elementOfMatrix(trx, lineIndex, elementIndexOfLine));
else
indexOfResultOutput += swprintf(result+indexOfResultOutput, stringLength+1-indexOfResultOutput, L"%5.4Lg, ", elementOfMatrix(trx, lineIndex, elementIndexOfLine));
if(lineIndex+1 == length)
{
memcpy(result+indexOfResultOutput, L"]", sizeof(wchar_t)*2);
indexOfResultOutput += 2;
}
else
{
memcpy(result+indexOfResultOutput, L"]\n", sizeof(wchar_t)*3);
indexOfResultOutput += 3;
}
}
}
fclose(tempFilePointer);
return result;
}

static unsigned long

Matrix createMatrix(long double *data, unsigned long length)
{
if(length==0||data==NULL) return NULL;
Matrix result;
if((result=malloc(sizeof(struct matrix)+BytesOfMatrixSize)==NULL) return NULL;
memcpy(result->matrixData, data, BytesOfMatrixSize);
return result;
}

Matrix createIdentityMatrix(unsigned long length)
{
if(length==0) return NULL;
Matrix result;
if((result=malloc(sizeof(struct matrix)+BytesOfMatrixSize)==NULL) return NULL;
for(unsigned long index = 0; index<length; index++)
elementOfMatrix(result, index, index) = 1.0L;
return result;
}

Matrix plicateMatrix(const Matrix trx)
{
if(trx==NULL) return NULL;
Matrix result;
unsigned long length = trx->length;
if((result=malloc(sizeof(struct matrix)+BytesOfMatrixSize)==NULL) return NULL;
memcpy(result, trx, sizeof(struct matrix)+BytesOfMatrixSize);
return result;
}

void freeMatrix(Matrix trx)
{
free(trx);
}

/*FUNCTION: caculateInversedMatrix
* SUCCEED: inversedMatrix
*  FAILED: NULL
*    INFO: if can't find it inversedMatrix return NULL
*/
Matrix caculateInversedMatrix(const Matrix trx)
{
if(trx==NULL) return NULL;
const unsigned long length = trx->length;
Matrix result = createIdentityMatrix(length);
Matrix modifiedMatrix = plicateMatrix(trx);

for(unsigned long lineIndex = 0; lineIndex < length; lineIndex++)
{
unsigned long searchForCurrentLineCount;
if((searchForCurrentLineCount = searchForCurrentLineCount(modifiedMatrix, lineIndex, lineIndex))==0)
{
freeMatrix(modifiedMatrix);
freeMatrix(result);
return NULL;
}
if(searchForCurrentLineCount-1!=lineIndex)
{
manipulateMatrixSwapLineXY(modifiedMatrix, lineIndex, searchForCurrentLineCount-1);
manipulateMatrixSwapLineXY(result, lineIndex, searchForCurrentLineCount-1);
}
for(unsigned long restLineIndex = lineIndex+1; restLineIndex<length; restLineIndex++)
if(elementOfMatrix(modifiedMatrix, restLineIndex, lineIndex)!=0.0L)
{
long double k = elementOfMatrix(modifiedMatrix, restLineIndex, lineIndex)/elementOfMatrix(modifiedMatrix, lineIndex, lineIndex);
manipulateMatrixLineXMinusKTimesOfY(modifiedMatrix, restLineIndex, k, lineIndex);
manipulateMatrixLineXMinusKTimesOfY(result, restLineIndex, k, lineIndex);
}
}
for(unsigned long lineCount = length; lineCount>0; lineCount--)
{
#ifdef lineIndex
#error lineIndex plicated definition as macro
#else
#define lineIndex (lineCount-1)

for(unsigned long restLineIndex = 0; restLineIndex<lineIndex; restLineIndex++)
if(elementOfMatrix(modifiedMatrix, restLineIndex, lineIndex)!=0.0L)
{
long double k = elementOfMatrix(modifiedMatrix, restLineIndex, lineIndex)/elementOfMatrix(modifiedMatrix, lineIndex, lineIndex);
manipulateMatrixLineXMinusKTimesOfY(modifiedMatrix, restLineIndex, k, lineIndex);
manipulateMatrixLineXMinusKTimesOfY(result, restLineIndex, k, lineIndex);
}
#undef lineIndex
#endif
}
for(unsigned long lineIndex = 0; lineIndex<length; lineIndex++)
if(elementOfMatrix(modifiedMatrix, lineIndex, lineIndex)!=1.0L)
{
long double k = 1/elementOfMatrix(modifiedMatrix, lineIndex, lineIndex);
//manipulateMatrixMutiplyK(modifiedMatrix, lineIndex, k);
manipulateMatrixMutiplyK(result, lineIndex, k);
}
freeMatrix(modifiedMatrix);
return result;
}

/*FUNCTION: matrixLineNoneZeroCount
* SUCCEED: Count of line returned (not index)
*  FAILED: Zero
*    INFO:
*/
static unsigned long matrixLineNoneZeroCount(const Matrix trx, unsigned long columnIndex, unsigned long searchFromLineIndex)
{
assert(trx!=NULL);
assert(columnIndex < trx->length);
assert(searchFromLineIndex < trx->length);
for(unsigned long lineIndex = searchFromLineIndex; lineIndex < trx->length; lineIndex++)
if(elementOfMatrix(trx, lineIndex, columnIndex)!=0.0L) return lineIndex+1;
return 0;
}

static void manipulateMatrixMutiplyK(Matrix trx, unsigned long lineIndex, long double __notnull k)
{
assert(trx!=NULL);
assert(lineIndex < trx->length);
assert(k != 0.0L);
long double *thisLine = lineOfMatrix(trx, lineIndex);
for(unsigned long index = 0; index<trx->length; index++)
thisLine[index] = thisLine[index]*k;
}

static void manipulateMatrixLineXMinusKTimesOfY(Matrix trx, unsigned long xLineIndex, long double __notnull k, unsigned long yLineIndex)
{
assert(trx!=NULL);
assert(k!=0.0L);
assert(xLineIndex < trx->length);
assert(yLineIndex < trx->length);
long double *xLine = lineOfMatrix(trx, xLineIndex),
*yLine = lineOfMatrix(trx, yLineIndex);
for(unsigned long index = 0; index<trx->length; index++)
xLine[index] -= k*yLine[index];
}

static void manipulateMatrixSwapLineXY(Matrix trx, unsigned long xLineIndex, unsigned long yLineIndex)
{
assert(trx!=NULL);
assert(xLineIndex < trx->length);
assert(yLineIndex < trx->length);
long double *xLine = lineOfMatrix(trx, xLineIndex),
*yLine = lineOfMatrix(trx, yLineIndex);
for(unsigned long index = 0; index<trx->length; index++)
{
xLine[index] ^= yLine[index];
yLine[index] ^= xLine[index];
xLine[index] ^= yLine[index];
}
}

热心网友时间：2024-10-27 12:42

//
//Created By Kevin Feng
//

#include<stdio.h>
#define MAX 10

double A[MAX][MAX]; //系数矩阵
double b[MAX];      //右端项
double X[MAX];      //迭代向量
int NUM;            //A的阶数
int size;           //最大迭代次数



int main(void)
{
    int i,j,k;     //计数器
    float Aik;     //正消过程用到的变量名
    float S;       //回代过程用到的变量名
    //以下代码输入系数矩阵A,右端项b
    printf("请输入系数矩阵A的阶数：");
    scanf("%d",&NUM);;
    size=NUM;
    for(i=1;i<=size;i++)
    {
        printf("请输入A的第%d行元素,各元素间以空格间隔:\n",i);
        for(j=1;j<=size;j++)
            scanf("%lf",&A[i-1][j-1]);
    }

    printf("输入右端项b,各元素间以空格间隔:\n");
    for(i=1;i<=size;i++)
    {
        scanf("%lf",&b[i-1]);
    }


    //在屏幕中输出用户输入的系数矩阵A和矩阵B
    printf("\n亲爱的，您输入的维度是%d!\n您输入的矩阵A[][]:\n\n",NUM);    //在屏幕中输出用户输入的矩阵A
    for(i=0;i<size;i++)
    {
        for(j=0;j<size;j++)
            printf("%f\t",A[i][j]);
        printf("\n\n");
    }
    printf("\n您输入的矩阵b[]：\n\n");
    for(i=0;i<size;i++)
        printf("\t%f\t\n\n",b[i]);
    printf("\n\v");//打印矩阵b[]


    //以下代码是高斯消去法的主要步骤
    for(k=0;k<size-1;k++)
    {
        if(!A[k][k])
            return -1;
        for(i=k+1;i<size;i++)
        {
            Aik=A[i][k]/A[k][k];
            for(j=k;j<size;j++)
            {
                A[i][j]=A[i][j]-Aik*A[k][j];
            }
            b[i]=b[i]-Aik*b[k];
        }
    }//首先通过正消过程，参见书本137页
    printf("A[]\n");
    for(i=0;i<size;i++)
    {
        for(j=0;j<size;j++)
            printf("%f\t",A[i][j]);
        printf("\n\n");
    }
    printf("\nb[]\n");
    for(i=0;i<size;i++)
        printf("\t%f\t\n\n",b[i]);
    printf("\n\n"); //此处返回正消过程的结果
    X[size-1]=b[size-1]/A[size-1][size-1];
    for(k=size-2;k>=0;k--)
    {
        S=b[k];
        for(j=k+1;j<size;j++)
        {
            S=S-A[k][j]*X[j];
        }
        X[k]=S/A[k][k];
    }    //回代

    //此处打印出向量X
    printf("The solution x[]=\n\n");
    for(i=0;i<size;i++)
        printf("\t%f\t\n\n",X[i]);
    return 0;
}

from 网页链接

追问第一，这不是二元域的求解
第二，无解和无限解情况没给出啊