2016-09-30 14 views
0

Мне нужно вычислить с помощью биномиальных коэффициентов выражения (x + y) ** n, причем n очень велико (порядка 500-1000). Первый алгоритм для вычисления биномиальных коэффициентов, который пришел мне на ум, был multiplicative formula. Так что я закодировал в моей программе, какОшибка округления биномиальных коэффициентов

long double binomial(int k, int m) 
{ 
    int i,j; 
    long double num=1, den=1; 
    j=m<(k-m)?m:(k-m); 
    for(i=1;i<=j;i++) 
    { 
     num*=(k+1-i); 
     den*=i; 
    } 
    return num/den; 
} 

Этот код очень быстро на одной основной нити, по сравнению, например, с recursive formula, хотя последний один меньше подвержен ошибкам округления, так как включает в себя только суммы, а не подразделения. Так что я хотел протестировать эти альгоны для больших ценностей и попытался оценить 500 выбрать 250 (заказ 10^160). Я обнаружил, что «относительная ошибка» меньше 10^(- 19), поэтому в основном они одного и того же числа, хотя они отличаются примерно от 10^141.

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

Я искал Google таблицы для биномиальных коэффициентов, поэтому я мог копировать их, но лучший один я нашел остановки при п = 100 ...

+0

Вам будет лучше с большой библиотекой чисел и целым умножением/делением. –

+0

Возможный дубликат [Количество комбинаций (N выберите R) в C++] (http://stackoverflow.com/questions/9330915/number-of-combinations-n-choose-r-in-c) –

+0

Я не думаю, что это дубликат, так как я знаю, что формула правильная, но расчеты округлены и поэтому появляются ошибки. –

ответ

1

Если вы только вычисления индивидуальных биномиальных коэффициентов С (п, к) с n довольно большой, но не больше, чем около 1750 , то вам лучше всего с достойной библиотеки C, чтобы использовать функцию tgammal стандартной библиотеки:

tgammal(n+1)/(tgammal(n-k+1) * tgammal(k+1)) 

Испытано с реализацией ГНУ libm, что последовательно полученные результаты в течение нескольких U LP точного значения и, как правило, лучше, чем решения, основанные на умножении и делении.

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

Если n настолько велик, что tgammal(n+1) превышает диапазон длинные двойной (более 1754), но не настолько велик, что числитель переполняется, то мультипликативное решение лучшее, что вы можете получить без bignum библиотеки. Однако вы также можете использовать

expl(lgammal(n+1) - lgammal(n-k+1) - lgammal(k+1)) 

, который менее точен, но проще в кодировке. (Кроме того, если для вас полезен логарифм коэффициента, приведенная выше формула будет работать на довольно большом диапазоне n и k. Не использовать expl повысит точность.)

Если вам нужен диапазон биномиальных коэффициентов с тем же значением n, то вам лучше всего является итеративным дополнение:

void binoms(unsigned n, long double* res) { 
    // res must have (n+3)/2 elements 
    res[0] = 1; 
    for (unsigned i = 2, half = 0; i <= n; ++i) { 
    res[half + 1] = res[half] * 2; 
    for (int k = half; k > 0; --k) 
     res[k] += res[k-1]; 
    if (i % 2 == 0) 
     ++half; 
    } 
} 

выше производит только коэффициенты с к от 0 до п/2. Он имеет немного большую ошибку округления, чем мультипликативный алгоритм (по крайней мере, когда k приближается к n/2), но это намного быстрее, если вам нужны все коэффициенты и у него есть больший диапазон допустимых входных данных.

+0

Первое - действительно хорошее решение для того, что мне нужно сделать. Спасибо –

+0

Я просто попытался вычислить 500 выбрать 250. Я столкнулся с результатами моей биномиальной функции и предлагаемого метода tgammal с помощью некоторого онлайн-инструмента, и я подтверждаю, что использование tgammal более точно. –

1

Чтобы получить точные результаты целочисленные для малых k и m, лучшим решением может быть (небольшое изменение вашего кода):

unsigned long binomial(int k, int m) 
    { 
    int i,j; unsigned long num=1; 
    j=m<(k-m)?m:(k-m); 
    for(i=1;i<=j;i++) 
    { 
    num*=(k+1-i); 
    num/=i; 
    } 
    return num; 
    } 

Каждый раз, когда вы получаете комбинаторное число после деления num/=i, вы не получите усечения. Чтобы получить приблизительные результаты для более крупных k и m, ваше решение может быть хорошим. Но будьте осторожны, что умножение long double происходит намного медленнее, чем умножение и деление целых чисел (unsigned long или size_t). Если вы хотите, чтобы большие цифры были точными, вероятно, большое целое число class должно быть закодировано или включено в библиотеку. Вы также можете использовать Google, если есть быстрый факториальный алгоритм для n! чрезвычайно большого целого n. Это может помочь и с комбинаторикой. Формула Стирлинга является хорошим приближением для ln(n!), когда n большой. Все зависит от того, насколько вы точны.

+0

Лучше ли делиться внутри цикла? Я думал, что это было менее точно, чем просто сделать одно разделение, но, возможно, для больших чисел лучше разделить шаг за шагом. –

+0

'unsigned long' - это 64-разрядный целочисленный тип, поэтому ошибки не будет. Это всего лишь вопрос избежать переполнения, насколько мы можем. –

+0

Если вы хотите использовать 'long double', то лучше разделите вне цикла. Это происходит главным образом потому, что деление происходит медленнее, чем умножение. Точность должна быть одинаковой. –

0

Возможно, это не то, что ищет ОП, но можно аналитически аппроксимировать nCr для больших n с функцией двоичной энтропии. Он упоминается в

+0

Бинарная функция энтропии получена из формулы Стирлинга без коррекции асимптотического расширения O (1/n). Если то, что ищет ОП, является просто функцией вероятности биномиального распределения, лучшим решением было бы это сделать именно для малых п и использовать гауссовское приближение для «n> 100, p = 0,05-0,95' и приближение Пуассона для 'n> 100, p <0.05 or p> 0.95', я думаю. –

+0

Да, одно решение может делиться и приближаться, когда n выбирает k «большой». Проблема снова в том, насколько я усекаюсь в таких случаях, и, поскольку время выполнения моей реализации формулы умножения на данный момент не является проблемой, если оно стоит –

1

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

  1. Реализовать формулу с большими целыми числами (долго долго, например)
  2. операций Попытка деления как можно скорее (как это было предложено Zhuoran)
  3. Добавьте код, чтобы проверить правильность каждого деления и умножения
  4. Разрешить неправильные деления или умножения, например
    • попробовать разделение в петле, предложенной Zhuoran, но если это не прибегать к начальному алгоритму (накапливает произведение делителя в берлоге)
    • магазина неразрешенного мультипликатора, делители в дополнительных длинных целых числах и попытаться решить их в следующих итерационных циклах
  5. Если вы действительно используете большие числа, то ваш результат может не соответствовать длине целого числа. то в этом случае вы можете переключиться на длинный двойной или использовать свое личное хранилище LongInteger.

Это скелет кода, чтобы дать вам идею:

long long binomial_l(int k, int m) 
{ 
    int i,j; 
    long long num=1, den=1; 
    j=m<(k-m)?m:(k-m); 
    for(i=1;i<=j;i++) 
    { 
     int multiplier=(k+1-i); 
     int divisor=i; 
     long long candidate_num=num*multiplier; 
     //check multiplication 
     if((candidate_num/multiplier)!=num) 
     { 
      //resolve exception... 
     } 
     else 
     { 
      num=candidate_num; 
     } 

     candidate_num=num/divisor; 
     //check division 
     if((candidate_num*divisor)==num) 
     { 
      num=candidate_num; 
     } 
     else 
     { 
      //resolve exception 
      den*=divisor; 
      //this multiplication should also be checked... 
     } 
    } 
    long long candidate_result= num/den; 
    if((candidate_result*den)==num) 
    { 
     return candidate_result; 
    } 
    // you should not get here if all exceptions are resolved 
    return 0; 
} 
+0

Я не понимаю, как разрешить неправильные разделы. Я имею в виду, если я сделаю один и тот же расчет за каждый цикл, я бы ожидал такого же неправильного результата? Не могли бы вы немного разобраться? Спасибо –

+0

Ну, идея состоит в том, чтобы реализовать исключение как отложенный алгоритм. Например, в моем снимке кода я пытаюсь сделать деление как можно скорее («candid_num = num/divisor» в моем коде). Если деление правильное, я беру кандидата, иначе я «храню» делитель в логарифме, чтобы отложить деление, когда это возможно, в этом случае в «long long candid_result = num/den». – SorMun

 Смежные вопросы

  • Нет связанных вопросов^_^