Структура данных для хранения коэффициентов многомерного полинома?

Можно смотреть на это как на структуру данных, которая должна реализовать угол D-мерной кубической сетки отрезанный плоскостью проходящей через ближайшие к нижнему левому углы вершины

Понял только в общих чертах что нужно, выглядит так что подойдёт обычный многоменный массив на векторе, только смещение нужно считать чуть хитрее чтобы не хранить отсутствующие элементы.

как std::map<uint32_t, double>, есть ощущение что доступ к этой штуке сильно тормозит

Конечно тормозит, не нужно map использовать если тебе не нужна упорядоченность или нет security implications. Если лень думать над вектором бери unordered_map и всё, о коллизиях не надо думать.

Сейчас это реализовано тупо, как std::map<uint32_t, double>, есть ощущение что доступ к этой штуке сильно тормозит.

Не совсем очевидно что Вы называете «доступом». Вставлять и удалять надо редко, а искать часто? И насколько «sparse» набор ключей? И сколько их вообще? И много ли таких структурок будет плодиться (может overhead по памяти вообще не критичен)?

P.S. Что-то мне подсказывает что от поисков совсем можно избавиться и на одном проходе обсчитывать всё …

А можно юзкейс или пример кода как это используется?

Если рассмотреть двумерный случай, кажется что можно представить многочлен как массив коэффициентов и массив степеней. Если чуть-чуть перестать экономить память и хранить степени тоже в double интуитивно кажется, что можно паковать такое в обычный vector с любым уровнем вложенности, предваряя его количеством коэффициентов в текущем уровне.

Т.е. для двумерного 2(3y + 1)x^2 + 5x + 0 это будет что-то типа:

2.0 - степень "внешнего" многочлена (от X)
1.0 - степень "внутреннего" многочлена (от Y)
3.0
1.0
2.0 - коэффициент первого члена X

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

подойдёт обычный многоменный массив на векторе, только смещение нужно считать чуть хитрее чтобы не хранить отсутствующие элементы.

Вот в этом и вопрос - как считать смещения;-) Собственно я уже догадался, но получается довольно сложный рекурсиный алгоритм расчета смещения - мб обычный метод половинного деления будет быстрее.

о коллизиях не надо думать.

Можно не думать, если не лень потом переделывать.

Не совсем очевидно что Вы называете «доступом».

В основном вставка и поиск, удаление неактуально. Но общее число элементов можно посчитать заранее, так что только поиск.

насколько «sparse» набор ключей?

Тут все плохо. Заняты первые биты в каждом байте ключа.

И много ли таких структурок будет плодиться (может overhead по памяти вообще не критичен)?

К сожалению критичен. Дело даже не в их количестве, а в том что при больших D и N оверхед (при хранении в обычном D мерном массиве) может быть настолько чудовищным что не хватит памяти узла кластера. Мне доступны узлы с памятью 256…768Гб.

от поисков совсем можно избавиться и на одном проходе обсчитывать всё

Увы нет. Эта структура накапливает данные от одного прохода по более сложной структуре, вот та вообще никуда не лезет и поэтому воcпроизводится на лету за один проход;-)

А можно юзкейс или пример кода как это используется?

Я последний месяц увлекся компьютерной алгеброй для полуаналитического расчета одного интеграла. В лоб он берется плохо. Под интегралом экспонента от полинома, экспонента раскладывается в ряд, дальше интеграл берется аналитически и остается полином от параметров. Ну вот это для хранения суммы вот таких штук (a+b+…)^N, только перед каждой буковкой еще коэффициент - эти коэффициенты и надо хранить

Я наверное плохо сформулировал. Сами степени целочисленные и хранить их не нужно - степень задается позицией в массиве. В одномерном случае полином это просто вектор его коэффициентов. В двумерном это вот

a 0 0 0
b c 0 0
d e f 0
g h j k 

В многомерном все плохо, ну то есть я могу себе это представить, но объяснить получается не очень.

То что ты предлагаешь это некоей факторизованное представление как я понял. Как ни странно да, я вчера перед сном о чем то похожем подумал - там на каждом уровне число коэффциентов можно точно посчитать (заранее). Получается некое хитрое дерево упакованное в вектор, но теперь непонятно - стоит ли заморачиваться, не будет ли обычное половинное деление проще и быстрее?

Видимо пока не попробую не узнаю…

Если в двумерном примере выше все (или почти все) a..k не нулевые (приемлемо хранить их все) - то они однозначно запихиваются в вектор с довольно несложной (x,y) -> pos. Несложной даже для произвольного (но известного заранее) числа измерений.

В двумерном нулей чуть меньше половины, дальше с ростом D все становится гораздо хуже. При N=10 вот такая картина маслом:

from math import *
N = 10
for D in range(1,10): print D, sum(factorial(l+D-1)/factorial(D-1)/factorial(l) for l in range(N+1)), N**D
... 
1 11 10
2 66 100
3 286 1000
4 1001 10000
5 3003 100000
6 8008 1000000
7 19448 10000000
8 43758 100000000
9 92378 1000000000

скажем до 8ми

Может завести структуру с битовой маской и динамическим массивом в котором будут последовательно без нулей лежать все коэффициенты? Соответственно чтобы вынуть число проверяешь по маске отлично ли оно от нуля, если да то считаешь единицы перед ним (по той же маске), так получаешь индекс в массиве.

Ничего более плотного в голову не приходит.

Я не очень понимаю как быстро посчитать единицы перед числом. Вчера я придумал рекуррентную формулу, но она довольно сложная получается.

Ничего не понятно. Может тебе нужен обычный двухмерный массив размера DxN?

Можно смотреть на это как на структуру данных, которая должна реализовать угол

Структура данных, которая реализует угол... :(

Или как на результат возведения в степень

Так твоя структура должна «реализовать угол» или «как на результат возведения в степень»?

Выражайся яснее.

https://ru.wikipedia.org/wiki/%D0%A2%D1%80%D0%B5%D1%83%D0%B3%D0%BE%D0%BB%D1%8C%D0%BD%D0%B8%D0%BA_%D0%9F%D0%B0%D1%81%D0%BA%D0%B0%D0%BB%D1%8F

Речь идет сумме треугольников Паскаля со степенями 0…N обобщенной на случай произвольной размерности D.

яснее я выразится не смогу.

Алгоритм можно придумать, но лично я бы просто сгенерил метод с развесистым ифом от 2х чисел. Там меньше 64х кейсов.

интересный подход. Хотя рекурсия мне как то ближе, но такой иф будет скорее всего эффективнее?

Если рассматривать двумерный случай (треугольник паскаля это же двумерный случай?), то я предложил бы линейный массив аля (std::vector<double>). Нужный индекс в массиве вычислять по формуле.

Думаю с многомерным случаем можно сделать тоже самое.

Думаю с многомерным случаем можно сделать тоже самое.

Вот именно об этом я в исходном посте и писал, я тоже так думаю. Вопрос - КАК именно сделать то же самое;-)

интересный подход. Хотя рекурсия мне как то ближе

Ты ученый, я погромист.

Так у тебя будет проверка 2х условий самое большее 16 раз. Мне кажется это проще чем рекурсия.

Если D максимум 8, а N максимум 256, 8*256 - это не так уж много.

Думаю можно замутить обращение по индексу линейного массива с предварительно подготовленной таблицей смещений.

КТо на что учился, да. Рекурсию руками можно написать в пару строчек, а ифы эти придется скриптом делать.

Но пожалуй тут прав pathfinder - можно фигакнуть табличку и развернуть рекурсию в цикл. О табличке я как то вчера невнятно думал, сходу казалось сложным…

Единственное что, там степени по маске выдирать из ключа один фиг придется. А половинное деление я уже набросал, там вообще ничего делать не надо, только сравнивать ключи.

Ща еще подумаю и погляжу что половинное деление дает.

Попробую привести пример для треугольника паскаля(когда D=2 ?).

Пускай у нас есть координаты x и y - номер строки и номер элемента относительно строки. В каждой последующей строке на один элемент больше.

Пускай у нас массив индексов первого элемента для каждой строки M. Тогда i=M[x]+y.

Теперь это надо перенести на многомерный случай.

Ну да, как то так.

Небольшой пример. Приведем индексы для некоего треугольника Паскаля.

0
1 2
3 4 5
6 7 8 9

Пускай у нас есть таблица для каждого первого элемента строки.

M= [0,1,3,6]

Тогда можно посчитать индекс по формуле i=M[x]+y

Теперь вопрос, а как быть с трехмерным случаем?

Примерно так же. Вектор разбивается на блоки по первой степени, внутри каждого блока своя индексация на основе оставшихся степеней. Размер таблицы в любом случае будет меньше чем размер вектора коэффициентов.

Хотя если там чуть подумать - скорее всего все эти таблицы по степеням очень похожи.

Если я ничего не путаю, то точно так же (мне трудно в голове представлять многомерный вариант треугольника Паскаля).

Пускай у нас трехмерный случай D=3. Координаты x,y,z. Если зафиксировать одну из координат равную константе, то оставшееся множество будет двухмерным случаем D=2.

Пускай фиксируем z=7. Тогда в оставшемся двухмерном множестве смещения будут иметь вид:

w+0,
w+1,w+2
w+3,w+4,w+5
w+6,w+7,w+8, w+9

Где w - это индекс элемента [z=7,x=0,y=0]

Т.е. для трехмерного случая мы можем завести новую таблицу W. Где каждый её элемент - смещение при x=0, y=0.

Тогда индекс элемента трехмерной сущности будет рассчитываться по формуле i=W[z]+M[y]+x.

В итоге нам достаточно предварительно рассчитать таблицу смещений для всех случаев размером DxN.

Ну что, это считается ответом на вопрос?

Не не, я понял, что степени хранятся (в текущем варианте) исключительно как ключ. По-факту я предлагаю этот ключ размазать по массиву, как некие навигационные узлы дерева, что бы было понятно как по нему спуститься.

Идею с половинным делением - я не понял если честно.

Картинку и представление многочленов в голове примерно такое же и держал. Ты попробуй двумерный случай записать не в матрицу а в вектор, возможно появятся идеи как ещё упаковать нули в одно чиселко. Но одно дело компактность представелния а второе это эффективность поиска, кажется что нужно строить n-арное дерево какое-то.

Одним из, да. Спасибо!

Идею с половинным делением - я не понял если честно.

Степени одного монома (слагаемого многочлена) можно собрать в целое число, по 8 бит на степень. Нам нужно хранить коэффициенты при мономах.

Делаем вектор пар степени:коэффициент, сортируем его по степеням. Дальше ищемся половинным делением. Тот же std::map, только оптимизированный с т.з. поиска.

нужна структура данных которая обеспечит эффективный доступ к коэффициенту если известен набор степеней k.

Глянь cxsc – C++ library for Extended Scientific Computing ©, может там что-нибудь полезное есть.

Я для решения данной задачи, если конечно оверхед по памяти действительно критичен, взял бы разреженные матрицы. Либо дефолтные какие-нибудь, типа Eigen::SparseMatrix или из boost::numeric::ublas, либо руками навелосипедил.

По описанию немного напоминает разреженные матрицы.

Вот, выше, как оказалось, их уже предложили.

О, привет;-)

Это не совсем оно, но похоже - размерность выше, дырок нет но граница неровная. Я знаю два вменяемых решения - std::map (его для записи gmm напр. юзает) или построчный формат, который собственно выше pathfinder предложил де-факто

Вот в этом и вопрос - как считать смещения;-) Собственно я уже догадался, но получается довольно сложный рекурсиный алгоритм расчета смещения - мб обычный метод половинного деления будет быстрее.

Нет, никаких рекурсивных алгоритмов там не должно быть - это будет медленнее хэша. Если я, опять таки, правильно понял что примерно нужно, то нужна просто сумма арифметической прогрессии. Т.е. для примера, двумерный случай:

Полная матрица NxN адресуется m[y*N + x], а если нужно хранить только её диагональную половину, то m[y*N - (y*y-y)/2 + x]

Теперь это надо обобщить на твой многомерный.

Можно не думать, если не лень потом переделывать.

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

Кстати я тут прикинул, что вместо ифов можно на старте сгенерить например двумерный массив где по одной координате индекс по другой битовая маска. Тогда всё вычисление это дважды взятие элемента из массива.

В двумерном нулей чуть меньше половины, дальше с ростом D все становится гораздо хуже. При N=10 вот такая картина маслом:

Возможно мы общаемся на слишком разных языках: в итоге Вы пришли к тому же что я сразу сказал - vector<double> плюс функция дающая позицию в векторе по координате в N-мерном пространстве.

Двумерный случай вообще тривиален, и никаких табличек даже не надо. Например:

int calcPos(int x, int y)
{
   int diag = x + y;
   int prevSize = diag * (diag + 1) / 2;
   return prevSize + x;
}

Ну это просто чтобы идею продемонстрировать, не претендую ;)

Мне кажется, если подумать, то на N-мерный случай оно тоже обобщается, и самое главное - остаётся O(1).

На unordered_map общий выигрыш по скорости втрое (по сравнению с map). Хотя на каких то тестах ЕМНИП unordered был на порядок быстрее.

Половинное деление вообще профита не дало.

Ну это просто чтобы идею продемонстрировать, не претендую ;)

Мне кажется, если подумать, то на N-мерный случай оно тоже обобщается, и самое главное - остаётся O(1).

Вроде как мне удалось сформулировать свою идею в [псевдо] математических терминах:

Как только Вы выведете формулу для числа точек на срезе N-мерного куба равноудаленных от его нулевой вершины (таких что sum(x{i}) = M) Вы получите отображение координат которые Вас интересуют на множество натуральных чисел, причём для данного max(M) такое отображение будет (запамятовал правильный термин) максимально компактно(?) / плотно(?). Мне представляется цена этого отображения будет O(D) и при это оно будет абсолютно нечувствительно к M. Но D у Вас маленькое, и, спешу заметить, то что Вы сейчас делаете уже O(D) (Вам же это очевидно?). Лучшего варианта пока не вижу.

На unordered_map общий выигрыш по скорости втрое (по сравнению с map). Половинное деление вообще профита не дало.

Ну ёпт, а как ты хотел константу против логарифма с кучей indirection’ов. Какое нахрен деление.

Хотя на каких то тестах ЕМНИП unordered был на порядок быстрее.

У мапы против вектора минус в дополнительном indirection и в отсутствии вообще какой бы то ни было локальности. Если осилишь формулу для адресации в векторе - это скорее всего будет самое оптимальное, только

• Нужно правильно выбрать главное измерение (т.е. чтобы в памяти оно лежало в том же порядке в котором ты его обходишь, если такой кейс вообще есть).
• Размерность передавай как шаблонный аргумент, чтобы компилятор мог оптимизировать умножения сдвигами.
• Если измерений так много что формула получится с 100500 умножений, узким местом может стать это (т.е. посчитать хэш и перейти по указателю станет быстрее чем посчитать смещение в векторе). Тогда есть гибридный и простой вариант - оставить хэш, но использовать кастомный аллокатор который гарантированно кладёт твои даблы подряд в памяти, а потом заполнить этот хэш в том же порядке в котором будет доступ. Итого бесплатно получаем локальность.

И есть ещё идея - если тебе нужна локальность сразу по всем осям, можно хранить в узлах (вектора или мапы) не один коэффициент, а сразу гиперкуб. Т.е. для 2d матрицы 8x8 это значило бы хранить 4x4 матриц 2x2. Так что если тебе нужно для элемента ходить в ближайших соседей, они сразу все в L1 будут.

Да, я примерно так же рассуждал. Формула есть, хотя она почему то не общеизвестна (N+D-1)!/N!(D-1)! она и правда O(D) - N! можно сократить

Дальше просто, находим срез гиперкуба отвечающий степени и начинаем понижать размерность.

Но ya-betmen говорил что есть вариант с битовой маской - вот это я пока че то не могу сообразить;-(

И есть ещё идея - если тебе нужна локальность сразу по всем осям, можно хранить в узлах (вектора или мапы) не один коэффициент, а сразу гиперкуб.

Да, это называется тайлинг, мы так часто делаем. Нет, в данном случае это не поможет - обход довольно прихотливый… ну или придется его как то сильно вкуривать, там есть закономерности но они сложные.

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

обход довольно прихотливый… ну или придется его как то сильно вкуривать, там есть закономерности но они сложные.

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

Хотя об чём это я - я забыл что тут просто даблы, значит всяко быстрее будет мапа с открытой адресацией. Это не stl, но легко делается на том же векторе.

вариант с битовой маской

Быстро накидал что я имел в виду для степеней <=8. Но метод можно расширить и на большие степени. По сути получение коэффициента это сдвиг побитовое сложение и получение элемента по индексу.

https://pastebin.com/rX65i89R

Спасибо, понял. Да, это хорошая структура данных что бы не хранить нули (одна из), но она не решает проблем индексации в D-мерном случае. Хотя бы потому, что если использовать наивную индексацию (как в D-мерном кубе) то даже размер маски будет какой то чудовищный.

Я перечитал посты (исходно я решал задачу просто хранения коэффициентов полинома из оп). Если я правильно понимаю то у тебя есть Д-мерное пространство, в точках с целочисленными координатами живут даблы (не более 255 штук, расположены Д-мерным треугольником возле начала отсчёта, всё остальное нули), это и есть твои коэффициенты. Ты хочешь имея на входе Д координат получить число которое по этим координатам лежит. И таких структур тебе нужно много?

Или я чего-то недопонял?

Как ты получил 66 для двумерного массива при N=10? 66 было бы при N=11.

1									
2	11								
3	12	20							
4	13	21	28						
5	14	22	29	35					
6	15	23	30	36	41				
7	16	24	31	37	42	46			
8	17	25	32	38	43	47	50		
9	18	26	33	39	44	48	51	53	
10	19	27	34	40	45	49	52	54	55