Обращение полностью заполненной матрицы 4000х4000 (16 млн. элементов) на современных процессорах занимает в районе секунды или меньше. Конечно неизвестны все нюансы, но выглядит как пример чудовищно неэффективной реализации.
В этом случае, для меня проще следующий способ: т.к. преобразование поворота всегда поворачивает относительно точки (0, 0), чтобы выполнить поворот относительно точки (x, y), нам нужно выполнить преобразование, чтобы эта точка переместилась в точку (0, 0), повернуть, после чего выполнить преобразование, чтобы точка вернулась на свое место.
40 секунд для матрицы 7600х7600 что то долго, тем более на GPU
Проверил свою реализацию 15 летней давности для обращения float64 матриц через LUP разложение с распараллеливанием для cpu. 19 секунд на i7-6800К в 12 потоков. Но 12 потоков это с hyper threading. Без него 23 секунды в 6 потоков/38 секунд в 3 потока, думаю масштабироваться должно близко к линейному, на современных процессорах по идее в несколько секунд должно укладываться. Реализацией могу поделиться, если интересно.
Область неквадратная - в ней количество столбцов больше чем количество строк, поэтому у нее нет одной диагонали, будет как минимум два варианта диагонали, каждый их которых заранее определяет, где будет проход. Можете сами проверить на примере области 4x3. Если вы поставите "монстров" на верхней диагонали (в координаты (x, y): (2, 2,), (3,1)), то правила расстановки не позволяют вам перекрыть проход справа, а если на нижней диагонали ((2,3), (3,2)), то не сможете перекрыть проход слева. Для зеркальных диагоналей то же самое.
Знаем в какой именно. Допустим у нас конфигурация, что закрашена клетка слева вверху от той, куда мы уперлись. Тогда у нас есть две области: слева вверху от конфигурации 2х2 из этих двух клеток и справа внизу. Для этих областей известен проход как к линии целиком над ними, так и выход с линии целиком под ними. Т.е. вопрос только в проходе через эти области, при этом для каждой области ситуация аналогична исходной. Очевидно, что проход точно есть в той области, где число столбов больше, чем число строк. Это условие будет всегда выполняться как минимум для одной из этих двух областей. Для конфигурации, когда закрашена клетка справа сверху, ситуация просто зеркальна.
За логарифмическое количество попыток путь находится. Спускаемся по среднему столбцу, пока не уткнемся. Исходя из локальной конфигурации 3х2, на определение которой нужно не больше двух попыток, можно либо понять, что мы нашли проход, либо узнать в левой или в правой части (либо в обеих) точно есть проход, тогда задача сводится к исходной, только с шириной в два раза меньше. Таким образом путь находится бинарным поиском.
А вот если взять F = 34, то R = 636291452, а 795364315 = 636291452 * 5/4. Суть метода для 27 я до конца не понял, конкретно что дает умножение на этот коэффициент и как потом от него избавляются.
Не знаю, сам вывел. Написал тестовую программку для проверки, формула работает. Рассуждал примерно так:
Перевод дробного числа в представление с фиксированной точкой производится путем его умножения на 2^F (обозначим M = 2^F) и отбрасывания дробной части. В нашем случае мы умножаем на обратное значение для нашего делителя (D) равное 1/D. Пусть его представление с фиксированной точкой будет Rfixed = [1/D * M] ([] - взятие целой части). Реальное дробное значение этого представления будет R = Rfixed / M Таким образом 1/D = R + e, где (e) ошибка округления, причем 0 <= e < 1/M Соответественно R = 1/D - e, а e = 1/D - R Когда мы умножаем любое число X на R то получаем: X * R = X/D - e * X Таким образом вычисленное значение будет всегда меньше реального и при округлении у нас может получаться значение меньше на 1 чем должно быть. Поэтому прибавим к Rfixed единицу: Rcfixed = Rfixed + 1 Его дробное значение будет соответственно Rc = Rcfixed / M = R + 1/M = X / D + X * (1/M - e) обозначим ec = 1/M - e, при этом 0 < ec <= 1/M Соответственно X * Rс = X/D + eс * X Нам нужно обеспечить выполнение условия: X * Rс < [X/D] + 1 иначе при умножении X на Rc результат будет больше на 1, чем должно быть. Отсюда: X/D + ec * X < [X/D] + 1 Заметим, что максимальным значением для X/D - [X/D] будет (D-1)/D Таким образом: (D-1)/D + eс * X < 1, откуда: eс * X < 1/D Это условие будет выполняться для любого X, если оно выполняется для максимального значения X = Xmax, т.е. eс * Xmax < 1/D, откуда: eс < 1/(D * Xmax), обозначив Z = D * Xmax eс < 1/Z Здесь можно поступить просто, используя факт, что ec <= 1/M, обеспечить выполнение условия 1/M < 1/Z, откуда: M > Z, или 2^F > Z Для примера из статьи, где D = 10, а Xmax = 255, Z = D * Xmax = 2550, получаем F=12: 2^12 = 4096 > 2550 Но, в этом случае будет нужно больше 2*N бит для промежуточного результата, это в общем случае проблема, поэтому нужно проверять более строгое ограничение: eс < 1/Z 1/M - e < 1/Z M * (1/Z + e) > 1 M > Z/(1 + e * Z) M > Z/(1 + (1/D - R) * Z) (это неравенство можно преобразовать в (M + (M - Rfixed * D) * Xmax) > Z и в коде лучше использовать его, т.к. тут все вычисления целочисленные) то видим что, для F=11 M = 2048 R = 204 Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 204/2048)*2550) = 1277.4951076320849 это условие выполняется, а для F = 10 M = 1024 R = 102 Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 102/1024)*2550) = 1277.4951076320849 уже не выполняется, т.ч. минимальное значение F=11
Как я уже писал в предыдущем комментарии, есть числа, при оптимизации деления на которые указанным способом, нужно больше чем 2 * N бит (для них формула дает F = N + Nd). Для таких чисел компилятор использует другую оптимизацию. Вот пример для деления 32-битных чисел на 25 и 27:
#include <cstdint>
uint32_t div25(uint32_t x)
{
return x / 25;
}
uint32_t div27(uint32_t x)
{
return x / 27;
}
clang генерит:
div25(unsigned int):
mov eax, edi
imul rax, rax, 1374389535
shr rax, 35
ret
div27(unsigned int):
mov eax, edi
imul rax, rax, 795364315
shr rax, 32
sub edi, eax
shr edi
add eax, edi
shr eax, 4
ret
Ну вообще я не прав. Для D = 11 приведенный алгоритм не работает. Точное условие, при котором деление будет точным для всех X, это выполнение неравенства:
2^F > Z/(1 + (1/D - (R - 1)/2^F) * Z), где Z = D * (2^N - 1) (вывод приводить здесь не буду)
для N = 8, D =11, это неравенство выполняется для F=12, а для F=11 уже не выполняется. Но тут получается интересная ситуация, в этом случае нам потребуется больше чем 2*N бит для представления R * X, а это может уже приводить к тому, что такую оптимизацию делением для некоторых чисел выполнить просто невозможно, потому что, если например для 32 битных чисел потребуется больше 64 бит для промежуточного результата, то на 64-битной машине за одно умножение результат уже не получишь.
Какие то сложные вычисления. На сколько я могу судить обратное значение делителя и сдвиг вычисляются достаточно просто. Для деления N-битного числа X на число D (при условии, что D не является степенью двойки), нам нужно F бит дробной части fixed point представления обратного значения D:
F = N + Nd - 1, где Nd - число бит необходимое для представления числа D
fixed point представление для обратного значения делителя D
R = [2^F / D] + 1, где [] - взятие целой части
соответственно для примера из статьи, N = 8, D =10
Ну даже если данные реально с плавающей точкой, но диапазон их значений известен, то возможно проще их квантовать и считать статистику по квантованным целочисленным значениям. Конечно квантование будет вносить свою ошибку, но она должна быть постоянной и ее величиной можно управлять, меняя шаг квантования. Зато от накопительной ошибки можно полностью избавиться.
А откуда берутся значения с плавающей точкой? Если речь о данных из "реального мира", то обычно они приходят с АЦП, т.е. по факту целочисленны. В этом случае можно все накопление/обновление сумм выполнять в целочисленной арифметике и никаких проблем с точностью не будет.
Книга действительно читается с таким же интересом как художественная, сложно оторваться. Я не физик, но на мой взгляд изложено все понятно и логично, в отличие от всякого научпопа по той же теории инфляции, где при чтении чувствуешь себя дураком.
Хотелось бы узнать, как обстоят дела с данной теорией в научном сообществе в текущий момент? Какая критика присутствует? Какие основные аргументы оппонентов?
Потому что на основе одного графика никаких выводов сделать нельзя. А выводы о наличии эффекта делаются на основе совместного анализа двух графиков, конкретно их разницы: если разница больше нуля - переоценка, меньше нуля - недооценка.
Вообще обычно нужны ближайшие точки двух прямых. Пересечение - это частный случай, когда эти точки совпадают.
Решается достаточно просто. Даны две прямые:
a(t) = a0 + a * t
b(s) = b0 + b * s
Вектор между ближайшими точками прямых должен быть перпендикулярен обеим прямым. Т.е. имеем систему двух линейных уравнений:
dot(a, a(t) - b(s)) = 0
dot(b, a(t) - b(s)) = 0
Если уравнения линейно зависимы - прямые параллельны. Если нет, решаем, находим t и s. Если расстояние между соответствующим им точками равно нулю - прямые пересекаются, если не равно - это кратчайшее расстояние между прямыми.
Обращение полностью заполненной матрицы 4000х4000 (16 млн. элементов) на современных процессорах занимает в районе секунды или меньше. Конечно неизвестны все нюансы, но выглядит как пример чудовищно неэффективной реализации.
Ну точку вращения всегда нужно учитывать.
В этом случае, для меня проще следующий способ: т.к. преобразование поворота всегда поворачивает относительно точки (0, 0), чтобы выполнить поворот относительно точки (x, y), нам нужно выполнить преобразование, чтобы эта точка переместилась в точку (0, 0), повернуть, после чего выполнить преобразование, чтобы точка вернулась на свое место.
Т.е. итоговое преобразование рассчитывается так:
QTrasnform trot = ... // Преобразование поворота
QTransform tcenter = QTransform::fromTranslate(-x, -y);
QTransform t = tcenter * trot * tcenter.inverted();
Не проще использовать QPainter::setTransform() вместо QPixmap::transformed(), чтобы не воевать с подобными проблемами?
Не знаю только как лучше это сделать :) Код старый, лежит в виде нескольких файликов. Могу скинуть архив, но не знаю куда.
40 секунд для матрицы 7600х7600 что то долго, тем более на GPU
Проверил свою реализацию 15 летней давности для обращения float64 матриц через LUP разложение с распараллеливанием для cpu. 19 секунд на i7-6800К в 12 потоков. Но 12 потоков это с hyper threading. Без него 23 секунды в 6 потоков/38 секунд в 3 потока, думаю масштабироваться должно близко к линейному, на современных процессорах по идее в несколько секунд должно укладываться. Реализацией могу поделиться, если интересно.
Так это же, судя по всему, цитата из документации с сайта microsoft. А там вся документация это и есть машинный перевод с английского.
Область неквадратная - в ней количество столбцов больше чем количество строк, поэтому у нее нет одной диагонали, будет как минимум два варианта диагонали, каждый их которых заранее определяет, где будет проход. Можете сами проверить на примере области 4x3. Если вы поставите "монстров" на верхней диагонали (в координаты (x, y): (2, 2,), (3,1)), то правила расстановки не позволяют вам перекрыть проход справа, а если на нижней диагонали ((2,3), (3,2)), то не сможете перекрыть проход слева. Для зеркальных диагоналей то же самое.
Знаем в какой именно. Допустим у нас конфигурация, что закрашена клетка слева вверху от той, куда мы уперлись. Тогда у нас есть две области: слева вверху от конфигурации 2х2 из этих двух клеток и справа внизу. Для этих областей известен проход как к линии целиком над ними, так и выход с линии целиком под ними. Т.е. вопрос только в проходе через эти области, при этом для каждой области ситуация аналогична исходной. Очевидно, что проход точно есть в той области, где число столбов больше, чем число строк. Это условие будет всегда выполняться как минимум для одной из этих двух областей. Для конфигурации, когда закрашена клетка справа сверху, ситуация просто зеркальна.
За логарифмическое количество попыток путь находится. Спускаемся по среднему столбцу, пока не уткнемся. Исходя из локальной конфигурации 3х2, на определение которой нужно не больше двух попыток, можно либо понять, что мы нашли проход, либо узнать в левой или в правой части (либо в обеих) точно есть проход, тогда задача сводится к исходной, только с шириной в два раза меньше. Таким образом путь находится бинарным поиском.
Для 27 не влезает. Для 27
F = 37
R = 5090331611
А вот если взять F = 34, то R = 636291452, а 795364315 = 636291452 * 5/4. Суть метода для 27 я до конца не понял, конкретно что дает умножение на этот коэффициент и как потом от него избавляются.
Не знаю, сам вывел. Написал тестовую программку для проверки, формула работает. Рассуждал примерно так:
Перевод дробного числа в представление с фиксированной точкой производится путем его умножения на 2^F (обозначим M = 2^F) и отбрасывания дробной части. В нашем случае мы умножаем на обратное значение для нашего делителя (D) равное 1/D. Пусть его представление с фиксированной точкой будет Rfixed = [1/D * M] ([] - взятие целой части).
Реальное дробное значение этого представления будет R = Rfixed / M
Таким образом 1/D = R + e, где (e) ошибка округления, причем 0 <= e < 1/M
Соответественно R = 1/D - e, а e = 1/D - R
Когда мы умножаем любое число X на R то получаем:
X * R = X/D - e * X
Таким образом вычисленное значение будет всегда меньше реального и при округлении у нас может получаться значение меньше на 1 чем должно быть.
Поэтому прибавим к Rfixed единицу:
Rcfixed = Rfixed + 1
Его дробное значение будет соответственно
Rc = Rcfixed / M = R + 1/M = X / D + X * (1/M - e)
обозначим
ec = 1/M - e, при этом 0 < ec <= 1/M
Соответственно
X * Rс = X/D + eс * X
Нам нужно обеспечить выполнение условия:
X * Rс < [X/D] + 1
иначе при умножении X на Rc результат будет больше на 1, чем должно быть. Отсюда:
X/D + ec * X < [X/D] + 1
Заметим, что максимальным значением для X/D - [X/D] будет (D-1)/D
Таким образом:
(D-1)/D + eс * X < 1, откуда:
eс * X < 1/D
Это условие будет выполняться для любого X, если оно выполняется для максимального значения X = Xmax, т.е.
eс * Xmax < 1/D, откуда:
eс < 1/(D * Xmax), обозначив Z = D * Xmax
eс < 1/Z
Здесь можно поступить просто, используя факт, что ec <= 1/M, обеспечить выполнение условия
1/M < 1/Z, откуда:
M > Z, или
2^F > Z
Для примера из статьи, где D = 10, а Xmax = 255, Z = D * Xmax = 2550, получаем F=12:
2^12 = 4096 > 2550
Но, в этом случае будет нужно больше 2*N бит для промежуточного результата, это в общем случае проблема, поэтому нужно проверять более строгое ограничение:
eс < 1/Z
1/M - e < 1/Z
M * (1/Z + e) > 1
M > Z/(1 + e * Z)
M > Z/(1 + (1/D - R) * Z) (это неравенство можно преобразовать в (M + (M - Rfixed * D) * Xmax) > Z и в коде лучше использовать его, т.к. тут все вычисления целочисленные)
то видим что, для F=11
M = 2048
R = 204
Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 204/2048)*2550) = 1277.4951076320849
это условие выполняется, а для F = 10
M = 1024
R = 102
Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 102/1024)*2550) = 1277.4951076320849
уже не выполняется, т.ч. минимальное значение F=11
Как я уже писал в предыдущем комментарии, есть числа, при оптимизации деления на которые указанным способом, нужно больше чем 2 * N бит (для них формула дает F = N + Nd). Для таких чисел компилятор использует другую оптимизацию. Вот пример для деления 32-битных чисел на 25 и 27:
clang генерит:
Ну вообще я не прав. Для D = 11 приведенный алгоритм не работает. Точное условие, при котором деление будет точным для всех X, это выполнение неравенства:
2^F > Z/(1 + (1/D - (R - 1)/2^F) * Z), где Z = D * (2^N - 1) (вывод приводить здесь не буду)
для N = 8, D =11, это неравенство выполняется для F=12, а для F=11 уже не выполняется. Но тут получается интересная ситуация, в этом случае нам потребуется больше чем 2*N бит для представления R * X, а это может уже приводить к тому, что такую оптимизацию делением для некоторых чисел выполнить просто невозможно, потому что, если например для 32 битных чисел потребуется больше 64 бит для промежуточного результата, то на 64-битной машине за одно умножение результат уже не получишь.
Какие то сложные вычисления. На сколько я могу судить обратное значение делителя и сдвиг вычисляются достаточно просто. Для деления N-битного числа X на число D (при условии, что D не является степенью двойки), нам нужно F бит дробной части fixed point представления обратного значения D:
F = N + Nd - 1, где Nd - число бит необходимое для представления числа D
fixed point представление для обратного значения делителя D
R = [2^F / D] + 1, где [] - взятие целой части
соответственно для примера из статьи, N = 8, D =10
Nd = 4
F = 8 + 4 - 1 = 11
R = [2048/10] + 1 = [204.8] + 1 = 205
Через SVD это просто делается:
https://www.ltu.se/cms_fs/1.51590!/svd-fitting.pdf
Часть 3: Fitting Planes and Lines by Orthogonal Distance Regression
Ну даже если данные реально с плавающей точкой, но диапазон их значений известен, то возможно проще их квантовать и считать статистику по квантованным целочисленным значениям. Конечно квантование будет вносить свою ошибку, но она должна быть постоянной и ее величиной можно управлять, меняя шаг квантования. Зато от накопительной ошибки можно полностью избавиться.
А откуда берутся значения с плавающей точкой? Если речь о данных из "реального мира", то обычно они приходят с АЦП, т.е. по факту целочисленны. В этом случае можно все накопление/обновление сумм выполнять в целочисленной арифметике и никаких проблем с точностью не будет.
Книга действительно читается с таким же интересом как художественная, сложно оторваться. Я не физик, но на мой взгляд изложено все понятно и логично, в отличие от всякого научпопа по той же теории инфляции, где при чтении чувствуешь себя дураком.
Хотелось бы узнать, как обстоят дела с данной теорией в научном сообществе в текущий момент? Какая критика присутствует? Какие основные аргументы оппонентов?
Потому что на основе одного графика никаких выводов сделать нельзя. А выводы о наличии эффекта делаются на основе совместного анализа двух графиков, конкретно их разницы: если разница больше нуля - переоценка, меньше нуля - недооценка.
Вообще обычно нужны ближайшие точки двух прямых. Пересечение - это частный случай, когда эти точки совпадают.
Решается достаточно просто. Даны две прямые:
a(t) = a0 + a * t
b(s) = b0 + b * s
Вектор между ближайшими точками прямых должен быть перпендикулярен обеим прямым. Т.е. имеем систему двух линейных уравнений:
dot(a, a(t) - b(s)) = 0
dot(b, a(t) - b(s)) = 0
Если уравнения линейно зависимы - прямые параллельны. Если нет, решаем, находим t и s. Если расстояние между соответствующим им точками равно нулю - прямые пересекаются, если не равно - это кратчайшее расстояние между прямыми.
Я так понимаю, там в вычислениях участвуют числа гораздо большего порядка.