SRC: Определитель матрицы
От: Vampire Россия  
Дата: 25.09.02 20:03
Оценка: -2
Определитель матрицы (determinant) применяеться при решении системы линейных уровнений.
Полезен будет для студентов.

Матрица не ограничена в размере.

#include "stdlib.h"
#include "math.h"

void form_minor(double (*matr)[], double (**minor)[],
                int n, int m, int elem_i, int elem_j)
{
    int i, j;
    int k=0, l=0;
    for(j=0; j<m; j++)
    {
        if (j != elem_j)
        {
            for(i=0; i<n; i++)
            {
                if (i != elem_i)
                {
                    (**minor)[(n-1)*k+l] = (*matr)[n*j+i];
                    l++;
                }
            }
            k++;
            l=0;
        }
    }
    return;
}

double determinant(double (*matr) [], int n, int m)
{
    double res=0;
    double (*minor)[];
    int i, j;
    int znak;

    minor = NULL;
    minor = (double (*)[])calloc(sizeof(double), (n-1)*(m-1));
    j = 0;
    if ((n<3) &&(m<3))
    {
        res = (*matr)[0] * (*matr)[3] - (*matr)[2] * (*matr)[1];
        return res;
    }
    for (i=0; i<n; i++)
    {
        form_minor(matr, &minor, n, m, i, j);
        znak = (int)pow(-1, i+j);
        res += znak * (*matr)[i + j * n] *
            determinant(minor, n-1, m-1);
    }
    return res;
}


Тестовая часть

int _tmain(int argc, _TCHAR* argv[])
{
//
    int i, j;
    int n, m;
    double det;
    double (*matr)[];
    double in;

    matr = NULL;
    n=3; m=3;

    matr = (double(*)[])calloc(sizeof(double), n*m);
    for (j=0; j<m; j++)
        for (i=0; i<n; i++)
        {
            printf(" Input %d, %d : ", i, j);
            scanf("%f", &in);
            (*matr)[j*n+i] = in;
        }
        printf("\n %f ", determinant(matr, n, m));
        det=9*3-2-5*(9-8)+7*(3-9*4);
        printf("\n Control = %f", det);

        free(matr);

    return 0;
}
Если долго мучиться что нибудь получится
Re: SRC: Определитель матрицы
От: Аноним  
Дата: 26.09.02 05:31
Оценка:
Здравствуйте Vampire, Вы писали:

Крайне потрясающе изумительно неоптимальный алгоритм. Его сложность O(n!). А это жутко.

Кстати, почему есть m и n отдельно? Ты можешь вычислить определитель неквадратной матрицы?
Re: SRC: Определитель матрицы
От: Bell Россия  
Дата: 26.09.02 06:44
Оценка:
Здравствуйте Vampire, Вы писали:

Мдааа....
Мож стоит хотя бы вспомнить метод Гаусса? Приведение матрицы к верхне(нижне)треугольному виду с последующим перемножением элементов главной диагонали гораздо эффективнее.
Любите книгу — источник знаний (с) М.Горький
Re[2]: SRC: Определитель матрицы
От: Кодт Россия  
Дата: 26.09.02 07:15
Оценка:
Здравствуйте Аноним, Вы писали:

А>Здравствуйте Vampire, Вы писали:


А>Крайне потрясающе изумительно неоптимальный алгоритм. Его сложность O(n!). А это жутко.


А>Кстати, почему есть m и n отдельно? Ты можешь вычислить определитель неквадратной матрицы?


У Vampire подпись изумительная: Если долго мучиться, что-нибудь получится
Перекуём баги на фичи!
Re[2]: SRC: Определитель матрицы
От: Vampire Россия  
Дата: 01.10.02 15:12
Оценка:
Здравствуйте Bell, Вы писали:

B>Здравствуйте Vampire, Вы писали:


B>Мдааа....

B>Мож стоит хотя бы вспомнить метод Гаусса? Приведение матрицы к верхне(нижне)треугольному виду с последующим перемножением элементов главной диагонали гораздо эффективнее.

Треугольный вид мне никогда не нравился
А что до оптимизации, тут конечно все правы Я об этом забочусь в последнюю очередь.
Алгоритм приведен в качестве сырца для студентов.
Если долго мучиться что нибудь получится
Re[2]: SRC: Определитель матрицы
От: Vampire Россия  
Дата: 01.10.02 15:15
Оценка: :)
Здравствуйте Аноним, Вы писали:

А>Здравствуйте Vampire, Вы писали:


А>Крайне потрясающе изумительно неоптимальный алгоритм. Его сложность O(n!). А это жутко.


А что в нем сложного
Решенее в лоб. Прямая реализация математического алгоритма.

А>Кстати, почему есть m и n отдельно? Ты можешь вычислить определитель неквадратной матрицы?


Не пробовал
Если долго мучиться что нибудь получится
Re[3]: SRC: Определитель матрицы
От: Кодт Россия  
Дата: 01.10.02 15:39
Оценка:
Здравствуйте Vampire, Вы писали:

V>А что до оптимизации, тут конечно все правы Я об этом забочусь в последнюю очередь.

V>Алгоритм приведен в качестве сырца для студентов.

Студенту-математику не нужно, он и так все прочтет в виде "лемма 123 теоремы 4 об определителях" в каком-нибудь талмуде.
А студенту-программисту подсовывать алгоритм со сложностью O(n!) — это подлянка. Потому что он ведь запомнит...

Мой любимый пример (автор — Дейкстра или Хоар, — не помню) — это функция вычисления числа Фибоначчи.
unsigned F(unsigned n)
{
  return (n < 2) ? 1 : F(n-1) + F(n-2);
}

При всей его изящности он требует O(n) стека и O(F(n)) времени.

Хотя итерационное решение занимает O(1) стека и O(n) времени.
А еще существует формула (увы, не помню) с квадратными корнями и возведением в n степень (что занимает O(log n)) времени.
Перекуём баги на фичи!
Re: SRC: Определитель матрицы
От: Аноним  
Дата: 23.03.04 18:25
Оценка:
Здравствуйте, Vampire, Вы писали:

V>Определитель матрицы (determinant) применяеться при решении системы линейных уровнений.

V>Полезен будет для студентов.

V>Матрица не ограничена в размере.


V>
V>#include "stdlib.h"
V>#include "math.h"

V>void form_minor(double (*matr)[], double (**minor)[],
V>                int n, int m, int elem_i, int elem_j)
V>{
V>    int i, j;
V>    int k=0, l=0;
V>    for(j=0; j<m; j++)
V>    {
V>        if (j != elem_j)
V>        {
V>            for(i=0; i<n; i++)
V>            {
V>                if (i != elem_i)
V>                {
V>                    (**minor)[(n-1)*k+l] = (*matr)[n*j+i];
V>                    l++;
V>                }
V>            }
V>            k++;
V>            l=0;
V>        }
V>    }
V>    return;
V>}

V>double determinant(double (*matr) [], int n, int m)
V>{
V>    double res=0;
V>    double (*minor)[];
V>    int i, j;
V>    int znak;

V>    minor = NULL;
V>    minor = (double (*)[])calloc(sizeof(double), (n-1)*(m-1));
V>    j = 0;
V>    if ((n<3) &&(m<3))
V>    {
V>        res = (*matr)[0] * (*matr)[3] - (*matr)[2] * (*matr)[1];
V>        return res;
V>    }
V>    for (i=0; i<n; i++)
V>    {
V>        form_minor(matr, &minor, n, m, i, j);
V>        znak = (int)pow(-1, i+j);
V>        res += znak * (*matr)[i + j * n] *
V>            determinant(minor, n-1, m-1);
V>    }
V>    return res;
V>}
V>


V>Тестовая часть


V>
V>int _tmain(int argc, _TCHAR* argv[])
V>{
V>//
V>    int i, j;
V>    int n, m;
V>    double det;
V>    double (*matr)[];
V>    double in;

V>    matr = NULL;
V>    n=3; m=3;

V>    matr = (double(*)[])calloc(sizeof(double), n*m);
V>    for (j=0; j<m; j++)
V>        for (i=0; i<n; i++)
V>        {
V>            printf(" Input %d, %d : ", i, j);
V>            scanf("%f", &in);
V>            (*matr)[j*n+i] = in;
V>        }
V>        printf("\n %f ", determinant(matr, n, m));
V>        det=9*3-2-5*(9-8)+7*(3-9*4);
V>        printf("\n Control = %f", det);

V>        free(matr);

V>    return 0;
V>}
V>

Не могли бы вы для тех же студентов написать ту же хню на Java, а то после ее (Java) недельного изучения (до этого к програмированию не имел отношения) слабо.
Заранее спасибо.
Был бы рад получить на burnnick_u@rambler.ru
Re[3]: SRC: Определитель матрицы
От: Andy77 Ниоткуда  
Дата: 23.03.04 19:51
Оценка:
Здравствуйте, Vampire, Вы писали:

А>>Крайне потрясающе изумительно неоптимальный алгоритм. Его сложность O(n!). А это жутко.


V>А что в нем сложного

V>Решенее в лоб. Прямая реализация математического алгоритма.

Re[4]: SRC: Определитель матрицы
От: McSeem2 США http://www.antigrain.com
Дата: 23.03.04 20:21
Оценка:
К>Мой любимый пример (автор — Дейкстра или Хоар, — не помню) — это функция вычисления числа Фибоначчи.
К>
К>unsigned F(unsigned n)
К>{
К>  return (n < 2) ? 1 : F(n-1) + F(n-2);
К>}
К>

К>При всей его изящности он требует O(n) стека и O(F(n)) времени.

Фисла Чибоначи — так их ведь мало совсем и для int32 и для int64 и даже для int128. Зачем вообще вычислять — можно просто взять из таблицы. O(1) времени и O(1) стека для i-того числа или bsearch (lower_bound) для числа, скажем, не превышающего заданное. O(log N).
McSeem
Я жертва цепи несчастных случайностей. Как и все мы.
Re[2]: SRC: Определитель матрицы
От: Кодт Россия  
Дата: 24.03.04 08:47
Оценка:
Здравствуйте, Аноним, Вы писали:

А>Не могли бы вы для тех же студентов написать ту же хню на Java, а то после ее (Java) недельного изучения (до этого к програмированию не имел отношения) слабо.


Вот "ту же ню" точно нет смысла портировать, тем более для студентов.
Крайне неэффективный алгоритм: требует O(N) стека и O(N!) времени.

В то же время метод Гаусса — O(N^2) памяти (дубликат матрицы) и O(N^2) времени.
Перекуём баги на фичи!
Re[4]: SRC: Определитель матрицы
От: Аноним  
Дата: 24.03.04 08:49
Оценка:
Здравствуйте, Кодт, Вы писали:

Ну да, а ведь есть еще алгоритмы и с более быстро растущей сложностью, их бы студентам-программистам
И стоя за спиной ухмыляться, "Ну что, не оптимизируется, сынок?"
Re[3]: SRC: Определитель матрицы
От: Кодт Россия  
Дата: 24.03.04 09:01
Оценка:
К>В то же время метод Гаусса — O(N^2) памяти (дубликат матрицы) и O(N^2) времени.

И, кстати говоря, метод Гаусса (модифицированный) обеспечивает наибольшую точность.
Перекуём баги на фичи!
Re[5]: SRC: Определитель матрицы
От: Kir. Россия  
Дата: 24.03.04 09:04
Оценка:
Здравствуйте, Аноним, Вы писали:

А>И стоя за спиной ухмыляться, "Ну что, не оптимизируется, сынок?"

Просто писать нужно на ASM`е
Писание же твое принято бысть и уразумлено внятельно. (С) Иван IV
Re[6]: SRC: Определитель матрицы
От: Кодт Россия  
Дата: 24.03.04 09:53
Оценка:
Здравствуйте, Kir., Вы писали:

K>Здравствуйте, Аноним, Вы писали:


А>>И стоя за спиной ухмыляться, "Ну что, не оптимизируется, сынок?"

K>Просто писать нужно на ASM`е

Напишите на ассемблере функцю Аккермана
Перекуём баги на фичи!
Re: SRC: Определитель матрицы
От: screw_cms Россия ICQ: 168185721
Дата: 24.03.04 10:21
Оценка:
Вот более оптимальная реализация (тут есть и определитель и ещё куча всего) — pure c, not c++

mathlib.h

#ifndef __MATHLIB_H_ALREADY_INCLUDED__
#define __MATHLIB_H_ALREADY_INCLUDED__

typedef double CELL_TYPE;

#ifndef bool
#define bool int
#endif
#ifndef true
#define true 1
#endif
#ifndef false
#define false 0
#endif

#define USE_FAST_KRAMER

// определим функцию AL_ADDITION как макрос
// алгебраическое дополнение к эл-ту (i,j) = -1^(i+j)
#define AL_ADDITION(x,y) ((x+y)%2==0?1:-1)

// Умножаем матрицы: mm(m x p) = lm(m x n) * rm(n x p)
void MatrixMult( CELL_TYPE** lm, CELL_TYPE** rm, CELL_TYPE** mm, int m, int n, int p);

// транспонируем кв. матрицу размера (size x size)
void MatrixTrans( CELL_TYPE** m, int size);

// Приведение матрицы в верхнетреугольному виду
void Gauss( CELL_TYPE** m, int size );

// считаем определитель кв. матрицы размера (size x size)
double MatrixDeterminant( CELL_TYPE** m, int size);

// быстрый определитель (методом Гаусса)
double GaussFastDeterm( CELL_TYPE** m, int size );

// m- матрица размера size x size, 
// v- вектор размера со свободными членами
// x- вектор размера для результата
bool Kramer( CELL_TYPE** m, int size, CELL_TYPE* b, CELL_TYPE* x);

// макроподстановки выделения памяти под вектор
#define VectorAlloc(size) (CELL_TYPE*)malloc( size*sizeof(CELL_TYPE) )
#define VectorFree(pVector) free(pVector)

// функции выделения памяти под матрицу
CELL_TYPE** MatrixAlloc( int size );
void MatrixFree( CELL_TYPE** ppMatrix, int size );


#endif


mathlib.cpp


#include "mathlib.h"
// Выделение памяти из кучи
#include <malloc.h>

// internal routines
CELL_TYPE** MatrixAlloc( int size )
{
    int i;
    CELL_TYPE** ppMatrix;
    
    // создаём временную матрицу размеров (size-1)x(size-1)
    // т.к. мы не собираемся ограничивать макс. размер, то
    // придётся выделять память из кучи
    ppMatrix = (CELL_TYPE**)malloc( size*sizeof(CELL_TYPE*) );
    for (i=0;i<size;i++)
        ppMatrix[i] = VectorAlloc( size );

    return ppMatrix;
}

void MatrixFree( CELL_TYPE** ppMatrix, int size )
{
    int i;
    
    // освобождаём память
    for (i=0;i<size;i++)
        VectorFree( ppMatrix[i] );
    free( ppMatrix );
}

// Умножаем матрицы: mm(m x p) = lm(m x n) * rm(n x p)
void MatrixMult( CELL_TYPE** lm, CELL_TYPE** rm, CELL_TYPE** mm, int m, int n, int p)
{
    int i,j,k;    
    for (i=0; i<m; i++)
        for (j=0; j<p; j++)
        {
            mm[i][j] = 0;
            for (k=0; k<n; k++)
                mm[i][j] += lm[i][k]*rm[k][j];
        }
}

// транспонируем кв. матрицу размера (size x size)
void MatrixTrans( CELL_TYPE** m, int size)
{
    int i,j;
    CELL_TYPE cTmp;

    for(i=1; i<size;i++)
        for(j=0;j<i;j++)
        {
            cTmp = m[i][j];
            m[i][j] = m[j][i];
            m[j][i] = cTmp;
        }
}

// Приведение матрицы в верхнетреугольному виду
void Gauss( CELL_TYPE** m, int size )
{
    int i, j, k;
    CELL_TYPE cWorkElem;

    for (i=0; i<size-1; i++)
    {
        cWorkElem = m[i][i];

        if ( cWorkElem==0 )
        {
            // пытаемся найти строку с ненулевым элементом,
            // среди строк, лежащих ниже
            for (j=i+1; i<size; i++)
            {
                cWorkElem = m[j][i];
                if (cWorkElem!=0)
                {
                    // добавляем найденную строку к рабочей
                    for (k=0; k<size; k++)
                        m[i][k]+=m[j][k];
                    break;
                }
            }
        }

        if (cWorkElem!=0)
        {
            // рабочий элемент корректен - 
            // обрабатываем все нижлежащие строки
            for (j=i+1; j<size; j++)
            {
                for (k=size-1;k>=i;k--)
                    m[j][k] -= m[i][k] * m[j][i] / cWorkElem;
            }
        }
    }
}

// считаем определитель кв. матрицы размера (size x size)
double MatrixDeterminant( CELL_TYPE** m, int size)
{
    int i,j,k,l;
    double fResult = 0;
    CELL_TYPE** mTemp;

    // simple cases
    if (size<1) return 0;
    if (size==1) return **m;
    if (size==2) return m[0][0]*m[1][1] - m[1][0]*m[0][1];
    
    // создаём временную матрицу размеров (size-1)x(size-1)
    mTemp = MatrixAlloc( size-1 );
    
    for (i=0; i<size; i++)
    {
        // собираем матрицу-минор, вычёркивая первую строку и i-ый столбец
        for(j=1; j<size; j++)
        {
            k=0; l=0;
            while(k<size)
            {
                if (k!=i)
                {
                    mTemp[j-1][l++] = m[j][k];
                }
                k++;
            }
        }
        // для каждой матрицы ищем определители меньших рекурсивно
        fResult += m[0][i] * AL_ADDITION(0,i) * MatrixDeterminant( mTemp, size-1 );
    }
    
    // освобождаём память
    MatrixFree( mTemp, size-1 );
    return fResult;
}


// быстрый определитель (методом Гаусса)
double GaussFastDeterm( CELL_TYPE** m, int size )
{
    double fResult = 1;
    int i;
    
    Gauss( m, size );
    for (i=0; i<size; i++)
        fResult *= m[i][i];
    
    return fResult;
}

// m- матрица размера size x size, 
// v- вектор размера со свободными членами
// x- вектор размера для результата
bool Kramer( CELL_TYPE** m, int size, CELL_TYPE* b, CELL_TYPE* x)
{
    int i,j,k;
    double fDeterm, fX;
    CELL_TYPE **mTemp;

#ifndef USE_FAST_KRAMER
    fDeterm = MatrixDeterminant( m, size );
#else
    fDeterm = GaussFastDeterm( m, size );
#endif

    if ( fDeterm==0 )
    {
        // Корни вычислить не удастся
        return false;
    }

    // создаём матрицу
    mTemp = MatrixAlloc( size );

    for (i=0; i<size; i++)
    {
        // вычисляем по-очереди корни

        // Формируем матрицу, замещая i-ый столбец на 
        // столбец вектора свободных членов
        for (j=0;j<size;j++)
            for (k=0; k<size;k++)
            {
                if (k==i)
                    mTemp[j][k] = b[j];
                else
                    mTemp[j][k] = m[j][k];
            }
        // считаем корень как определитель сгенерированной
        // делённый на определитель матрицы системы
#ifndef USE_FAST_KRAMER
        fX = MatrixDeterminant( mTemp, size );
#else
        fX = GaussFastDeterm( mTemp, size );
#endif
        x[i] = fX / fDeterm;
    }

    MatrixFree( mTemp, size );
    return true;
}
When in doubt, use brute force. © Ken Thompson

 
Подождите ...
Wait...
Пока на собственное сообщение не было ответов, его можно удалить.