Форум на Исходниках.RU https://forum.sources.ru Форум на Исходниках.RU Форум на Исходниках.RU https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905787 Sun, 16 Jun 2024 16:39:44 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905787 H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905645 Wed, 12 Jun 2024 13:49:29 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905645 https://ru.wikipedia.org/wiki/GNU_Scientific_Library

Там по документации есть функция gsl_linalg_HH_solve()

Документация:

https://www.gnu.org/software/gsl/doc/latex/gsl-ref.pdf]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905581 Tue, 11 Jun 2024 06:26:16 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905581 2024_06_11___002.GIF (, : 412)
MathCAD как и ожидалось считает Result с большой точностью после запятой
Хотелось бы сделать такую же функцию lsolve(A,B)
моя выдает 9999, ... можно конечно округлить, но где гарантия что для другого случая не затрется нужная цифра.]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905570 Mon, 10 Jun 2024 16:19:21 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905570 Qraizer, я в той же теме не так же плотно! Надо просто замутить стресс-тесты на числа "до запятой" порядка 5000-7000, и после запятой. Провести деление и умножение е на три. И оценить потом то, что преобразуется.]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905555 Mon, 10 Jun 2024 13:42:46 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905555 P.S. Я gmp тоже в глаза не видел до этой темы. Сужу исключительно по мануалам.

Добавлено
Ну так, положу тут просто.
    #include <iostream>
    #include <iomanip>
    #include <boost/multiprecision/cpp_bin_float.hpp>
     
    namespace bmp = boost::multiprecision;
    using longlongFloat = bmp::number<bmp::cpp_bin_float<400>>;
     
    int main()
    {
      longlongFloat a, b, result;
     
      a = longlongFloat("234777440777777777777777777777777777777.7324856732495879347856398564739865473986739658739867398657349856734589673463371");
      b = 17.0;
     
      result = a / b;
      result*= 17.0;
     
      std::cout << "Are you fun now, Majestio? " << std::setprecision(120) << std::scientific
                << result << std::endl;
    }
]]> Qraizer C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905554 Mon, 10 Jun 2024 13:29:10 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905554 Цитата Qraizer @
Ну так ты заюзал целочисленный класс. Юзай вещественный.

Заюзал, но ... от перемены мест слагаемых - легче не становится! :lol:

    #include <gmp.h>
     
    int main() {
        mpf_t a, b, result;
        mpf_init(a);
        mpf_init(b);
        mpf_init(result);
     
        mpf_set_d(a, 234777440777777777777777777777777777777.7324856732495879347856398564739865473986739658739867398657349856734589673463371);
        mpf_set_d(b, 17.0);
     
        mpf_div(result, a, b);
        mpf_mul_ui(result, result, 17.0);
     
        gmp_printf("Qraizer, you should be sad: %.1Ff\n", result);
     
        mpf_clear(a);
        mpf_clear(b);
        mpf_clear(result);
     
        return 0;
    }

    Qraizer, you should be sad: 234777440777777784841000000000000000000.0

Не, ну может я где и накосячил - заюзал впервые либу вчера! Но как-то не впечатляет пока :rolleyes:]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905546 Mon, 10 Jun 2024 12:11:55 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905546 Цитата Majestio @
И так. Число 10 делим на 3, делим, делим, делим, а потом умножаем на то же 3. Барабаны, фейерверки ... и получаем число 9. Аплодисменты и занавес!
Ну так ты заюзал целочисленный класс. Юзай вещественный.]]> Qraizer C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905536 Mon, 10 Jun 2024 10:29:48 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905536 Цитата H g @
Вопрос как увеличить точность решения ?

В таких случаях я всегда пытаюсь преобразовать
формулу вычислений таким образом, чтобы:
1. деление выполнялось один раз и было последним действием.
даже если это приводит к увеличению количества вычислений
2. произвести последнее деление с "правильным" округлением.
---
Если есть такая возможность, перейти к целым значениям.
Например, если в результате ожидаются дробные миллиметры,
считайте в тысячных долях миллиметров. Главное, чтобы
разрядной сетки целых хватило для всех преобразований.
---
Да, переход к double правильное решение.
Но надо понимать, что использование float/double это
принципиально приближённые вычисления. И если мы начали это
использовать, значит и результат будет приближённый.
Ничего ужасного в этом нет.]]> ЫукпШ C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905535 Mon, 10 Jun 2024 10:27:14 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905535 Цитата
Попробую еще метод Крамера

Не очень. Если при Гауссе погрешность возникала за счет деления, то в этом методе процесс лезет в область больших чисел. Да и определитель 100 x 100 тоже не подарок. Нужно что то "свое" придумать. Как-то использовать сильную разреженность...]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905530 Mon, 10 Jun 2024 09:12:18 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905530 Цитата MBo @
Итеративные методы пробовали? (например, Гаусса-Зейделя)

Нет, пока с прямыми "воюю" ))
В моем представлении, если ожидается результат какой-то иррациональный, тогда приближаться к нему нужно с помощью итераций.
У меня же наоборот, результат чаще всего целое число, или уж точно из области рациональных чисел.
Я к итерационным способам отношусь с недоверием. Где гарантия того, что на данных, которые введет пользователь метод не улетит черти зная куда, сложнее проверить.

Добавлено
Смотрю определители, нужные для метода Крамера хорошо считаются когда много нулей ...
Так может и не плохо получится с моим "велосипедом" )]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905514 Mon, 10 Jun 2024 09:00:20 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905514 H g
Итеративные методы пробовали? (например, Гаусса-Зейделя)]]> MBo C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905513 Mon, 10 Jun 2024 08:19:00 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905513
Так как у меня решение СЛАУ требуется исключительно в прикладных целях, наверное, нужно уточнить пределы задачи. Число уравнений (матрица) максимум до 100 x 100, если будет медленно работать - не критично, такой случай маловероятен. Обычно порядка 20 x 20.
Кроме того матрица сильно разреженная. Вот примерно такая, (решение в MathCAD)
2024_06_10___001.GIF (, : 551) ]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905478 Mon, 10 Jun 2024 05:56:00 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905478 Qraizer Если эта тема больше в "Алгоритмы" подходит, можете перенести туда. Правда активность на форуме не столь велика, поэтому перенос в "Алгоритмы" разве что для порядка. По теме : Направления, в которых можно посмотреть, по советам Majestio, amk, Qraizer я получил. Всем спасибо.]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905474 Mon, 10 Jun 2024 00:18:47 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905474
Скрытый текст
Цитата Qraizer @
Ребят, ещё немного и в Холивары скину.

Да пусть тут останется. Тема-то не про алгоритм как таковой, а об обеспечение точности.
А холивара не будет. amk, залетел, дан ценные указания использовать GMP и не "велосипедить"
Но произошел обломчег - его либа после двух действий дает 10 == 9 ...
that_s_fine.jpg (, : 505)
Вообще не вижу смысла холиварить :lol:
]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905472 Sun, 09 Jun 2024 17:28:45 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905472 H g, тут тема, как я смотрю, сильно переплетена и с инструментами (сиречь Прочим), и с Алгоритмами. Конкретно по Общим вопросам могу лишь предложить double вместо float, а также соптимизировать расчёты, избавивишись от потери точности, где возможно. (Если что, архитектура x86 и некоторых других умеет отлеживать потерю точности при вычисления и в нужным местах кидать эксепшны или хотя бы отражать факт в статусах. Но могут понадобиться инструменты, например, отладчик.) Типа вычитания близких по величине чисел, сложения малых и больших итп. По другим аспектам... в общем, реши, плз, как лучше: отзеркалить, перенести и куда.

Добавлено
P.S. Могу порекомендовать boost::multiprecision, когда-то был нужен коллеге, я мне было просто интересно повозиться. За производительность не скажу, но 1000-битные и даже -значные мантиссы легко держала. Интерфейс у неё, правда, не столь удобен, как хотелось бы, так что конкуренты в этом плане у неё есть, думаю.]]> Qraizer C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905465 Sun, 09 Jun 2024 16:14:45 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905465 Цитата amk @
Где-нибудь при обработке 20 строки числа перестанут помещаться в отведённые для них ячейки.

Ты не дочитал мой ответ и давай отвечать ... Но ты не на мой текст отвечаешь, а на то, что ты там сам себе придумал! Какие нафик ячейки?!!! :lol:]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905458 Sun, 09 Jun 2024 12:13:31 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905458 Цитата H g @
В целом если умножение/деление будет выполнено на последнем шаге, то все типа /3 и x3 в конце сократятся
Скорее всего до этого не дойдёт. Где-нибудь при обработке 20 строки числа перестанут помещаться в отведённые для них ячейки. Можно конечно для хранения числителя/знаменателя вместо long int взять long long int = int64_t, тогда это произойдёт попозже, строке так в 40-й.
Цитата Majestio @
amk, ты бы лучше не мне советовал как советовать, а посоветовал бы свой совет ТС'у
Я посоветовал. Воспользоваться не тетрадным вариантом исключения Гаусса, который он реализовал, а нормальным, компьютерным.]]> amk C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905450 Sun, 09 Jun 2024 03:06:47 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905450 Цитата amk @
Majestio, ты бы лучше посоветовал ему не изобретать велосипед, а прикрутить к программе какую-нибудь хорошо отлаженную библиотеку вычислений c повышенной точностью, gmp например.

amk, ты бы лучше не мне советовал как советовать, а посоветовал бы свой совет ТС'у. А давай ты посоветуешь предложенную тобою либу? Заодно и покажешь ее качество и "невелосипедность" на примерах? А мы проверим.

Добавлено
Хотя ... чего катать вату. Возьмем и накидаем пример деления и последующего умножения на "три", с использованием библиотеки GMP, очень отлаженной и очень надёжной. Ну и посмотрим качество совета.

    #include <stdio.h>
    #include <gmp.h>
     
    int main(void) {
        mpz_t x, result;
        mpz_init_set_str(x, "10", 10);
        mpz_init(result);
     
        // Деление на 3
        mpz_tdiv_q_ui(x, x, 3);
     
        // Умножение на 3
        mpz_mul_ui(result, x, 3);
     
        gmp_printf("Результат: %Zd\n", result);
     
        mpz_clear(x);
        mpz_clear(result);
        return 0;
    }


И так. Число 10 делим на 3, делим, делим, делим, а потом умножаем на то же 3. Барабаны, фейерверки ... и получаем число 9. Аплодисменты и занавес! :lol:

amk, а вывод простой - советовать другим как советовать - это явно не твое! :jokingly:]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905449 Sat, 08 Jun 2024 20:12:57 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905449 главный элемент. Достаточно искать по столбцу. Главный - это максимальный по величине (без учёта знака). Своим сравнением с эпсилон ты допускаешь, что коэффициент на который ты умножаешь ведущую строку при исключении переменной из остальных строк может оказываться много больше 1. А на этот коэффициент умножаются не только элементы матрицы, но и их погрешности.
Собственно алгоритмы без выбора главного элемента в программах давно уже не применяются. Они остались только на уроках математики при ручном решении СЛАУ, когда ученик глазами видит, что унего происходит.

И незачем сравнивать диагональный элемент с единицей - за исключением очень специальных случаев экономия вычислений мизерная, а программа усложняется. Ты пишешь у тебя накапливается погрешность, значит у тебя система уравнений далёка от такого специального случая.

Majestio, ты бы лучше посоветовал ему не изобретать велосипед, а прикрутить к программе какую-нибудь хорошо отлаженную библиотеку вычислений c повышенной точностью, gmp например.]]> amk C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905435 Sat, 08 Jun 2024 06:27:16 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905435 Majestio Спасибо за подробный ответ.
Цитата
не работайте с десятичными дробями, работайте со смешанными дробями
В этом что-то есть, хорошая мысль. В целом если умножение/деление будет выполнено на последнем шаге, то все типа /3 и x3 в конце сократятся. Нужно подумать как реализовать.]]> H g C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905434 Sat, 08 Jun 2024 02:22:58 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905434 Цитата H g @
Вопрос как увеличить точность решения ?

Интересный вопрос. Я когда-то им задавался, но не было прикладной задачи - и так руки и не дошли.

Тем не менее, поделюсь своими соображениями ...

Сперва немного рассуждений. Классический пример, приведенный выше (деление, а потом умножение на три), показывает, прямо тыкает носом в то, что использование float и double - тупиковый метод. Возникает законный вопрос, ачо делать-та? И тут возникает ответ: не работайте с десятичными дробями, работайте со смешанными дробями! Естественно, это нужно реализовывать самому (ну или поискать в сети, может уже кто удосужился).

Т.е. любое число у нас представляется примерно следующим образом:

    class BigDecimal {
      string WholePart;   // целая часть
      string Numerator;   // числитель
      string Denominator; // знаменатель
    }

Ну и далее реализация кучи публичных методов и перегрузок операций. Сразу возникает законный вопрос "WTF strings?". Рассказываю, может непоследовательно, но тезисно и с душой :lol:
  1. В своей реализации вообще-то BigDecimal мы должны задействовать какую-либо реализацию BigInt. Я встречал простейшие реализации бигинтов на строках. Возможно эти части и следовало бы объявлять не строками, а сразу большими целыми, но не суть...
  2. В процессе арифметических операций у нас вообще ничего не теряется *, Но за это нужно платить быстро растущими размерами всех полей BigDecimal, равно как и значительным снижением скорости вычислений
  3. В результате нас практически перестает беспокоить потеря точности, но начинает серьёзно напрягать проблема потери производительности и потребления ресурсов
  4. Поэтому, время от времени, по реализованной нами, вернее вами, стратегии, BigDecimal должен производить нормализацию своих данных путем поиска НОД числителя и знаменателя, и сокращать их на этот НОД
  5. Ну и второй вариант - если нормализации явно недостаточно и данные все равно слишком большие, должна быть разработана стратегия управления уменьшением точности - сокращением числителя и знаменателя (одного из) делением на какое-то серьёзное целое нацело, а второго из ... с потерями. Но когда это касается, к примеру 30 000-й цифры после запятой, потери точности буду едва заметны. В любом случае для решения мелких и средних задач по вычислениям хватит и решения из п.4.
  6. И конечно последний шаг - уже после всех вычислений мы можем спокойно вызвать метод, который из BigDecimal нам вернет float или double
* с учетом п.5

Как резюме

В такой реализации мы сохраняли все что можно в процессе вычислений, платили мучительными секундами ожидания и повышения температуры центрального процессора, но получили на выходе максимально возможный по точности результат во float или double.

Перед реализацией своего BigDecimal нужно очень внимательно выбирать либу BigInt, т.к. реализация на строках последней - не лучший вариант. Чуйка подсказывает, что на битовых массивах будет получше.

Ну и последнее. Ложечка дегтя в вашу баночку с вашим мёдом. Ясный перец, что речь идет только об элементарных арифметических операциях. В BigDecimal мы не сможем хранить иррациональные числа, не сможем обеспечивать 100% точность при извлечении корней, вычислении логарифмов & etc. Жизнь - боль!]]> Majestio C/C++: Общие вопросы https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905432 Fri, 07 Jun 2024 21:35:28 +0000 Точность вычислений. Метод Гаусса https://forum.sources.ru/index.php?showtopic=446685&view=findpost&p=3905432
a11x1 + a12x2 + a13x3 + ... a1NxN = b1
a21x1 + a22x2 + a23x3 + ... a2NxN = b2
a31x1 + a32x2 + a33x3 + ... a3NxN = b3
...
aN1x1 + aN2x2 + aN3x3 + ... aNNxN = bN

Получаю марицы a[][], b[], x[]. Решаю с помощью метода Гаусса.
Все работает, но есть одна проблема: накапливается погрешность и в итоге вместо точного ответа в 10000, имею 9999,

Код привожу, но он возможно требует пояснений, так как вырван из контекста

Скрытый текст
    void CalculateMatrixSolve()
    {
        int i, j, k, n, p, r;
        float a[100][100];
        float b[100];
        float x[100];
        float temp, eps;
        n = e::calculatematrixubound; // число уравнений
        eps = 0.0000001f;
        for (i=0; i<e::calculatematrixshift; i++) // это по сути обертка - внешние исходные данные, типа расчет выполнить i раз с разными входными данными.
                                                  // переменная i в самом алгоритме не используется
        {
            for (j=0; j<n; j++)
            {
                for (k=0; k<n; k++)
                {
                    a[j][k] = e::calculatematrixa[i][j][k]; // входная матрица A[j][k]
                }
                b[j] = e::calculatematrixb[i][j]; // входная матрица B[j]
                x[j] = 0.0f;
            }
            for (j=0; j<n; j++)
            {
                if (a[j][j] == 1)
                {
                    for (r=j+1; r<n; r++)
                    {
                        if (fabs(a[r][j]) >= eps)
                        {
                            temp = a[r][j];
                            for (p=j; p<n; p++)
                            {
                                a[r][p] -= (a[j][p] * temp);
                            }
                            b[r] -= (b[j] * temp);
                        }
                    }
                }
                else
                {
                    if (fabs(a[j][j]) < eps)
                    {
                        for (r=j+1; r<n; r++)
                        {
                            if (fabs(a[r][j]) >= eps)
                            {
                                for (p=j; p<n; p++)
                                {
                                    temp = a[j][p];
                                    a[j][p] = a[r][p];
                                    a[r][p] = temp;
                                }
                                temp = b[j];
                                b[j] = b[r];
                                b[r] = temp;
                                break;
                            }
                        }
                        if (fabs(a[j][j]) >= eps)
                        {
                            temp = a[j][j];
                            for (p=j; p<n; p++)
                            {
                                a[j][p] /= temp;
                            }
                            b[j] /= temp;
                            for (r=j+1; r<n; r++)
                            {
                                if (a[r][j] != 0)
                                {
                                    temp = a[r][j];
                                    for (p=j; p<n; p++)
                                    {
                                        a[r][p] -= (a[j][p] * temp);
                                    }
                                    b[r] -= (b[j] * temp);
                                }
                            }
                        }
                        else
                        {
                            return;
                        }
                    }
                    else
                    {
                        temp = a[j][j];
                        for (p=j; p<n; p++)
                        {
                            a[j][p] /= temp;
                        }
                        b[j] /= temp;
                        for (r=j+1; r<n; r++)
                        {
                            if (a[r][j] != 0)
                            {
                                temp = a[r][j];
                                for (p=j; p<n; p++)
                                {
                                    a[r][p] -= (a[j][p] * temp);
                                }
                                b[r] -= (b[j] * temp);
                            }
                        }
                    }
                }
            }
            for (j=(n-1); j>=0; j--)
            {
                x[j] = b[j];
                for (p=0; p<j; p++)
                {
                    b[p] -= (a[p][j] * x[j]);
                }
            }
            for (j=0; j<n; j++)
            {
                e::calculatematrixx[i][j] = x[j]; // выходная матрица X[j] результат
            }
        }
    }

Здесь есть некоторая погрешность (eps), потому что сам метод предполагает деление, чтобы получилась "треугольная" матрица. Так вот, если делитель будет например 3, то 1.0/3.0 = 0.3333... значит деление не может быть выполнено точно и дальше накапливается погрешность. Попробовал поменять (eps) но было еще хуже, при чем при увеличении числа знаков после запятой функция решает неправильно. Путем практики подобрал (eps = 0.0000001f;)
// Если 1.0/3.0 полчим 0.333333, а если потом 0.333333 умножить на 3.0 то получим 0.999999 а не 1.0
// В этом мне кажется загвоздка

Вопрос как увеличить точность решения ?
Предполагаемые варианты:
1) заменить float на double, поможет ли это
2) использовать другой метод не Гаусса, тогда какой..

Может что-нибудь посоветуете ?]]> H g C/C++: Общие вопросы