找到第n个Catalan数对m取模的最快算法是什么?(已知)

17
问题是找到模 m = (10^14 + 7) 的第 n个卡特兰数,其中m不是质数。这里是我尝试的方法列表:(最大N = 10,000
  1. 使用动态规划进行表查找,速度太慢
  2. 使用卡特兰公式ncr(2 * n,n)/(n + 1),由于ncr函数而不够快,不能使用指数平方加速,因为m不是质数。
  3. 硬编码预生成的Catalans表,但由于文件大小限制而失败。
  4. 递归关系C(i,k) = C(i-1,k-1)+C(i-1,k),这太慢了

那么,我想知道是否有其他更快的算法来找到我不知道的第n个卡特兰数?

使用动态规划

void generate_catalan_numbers() {
    catalan[1] = 1;
    for (int i = 2; i <= MAX_NUMBERS; i++) {
        for (int j = 1; j <= i - 1; j++) {
            catalan[i] = (catalan[i] + ((catalan[j]) * catalan[i - j]) % MODULO) % MODULO;
        }
        catalan[i] = catalan[i] % MODULO;
    }
}

使用原始公式

ull n_choose_r(ull n, ull r) {
    if (n < r)
        return 0;

    if (r > n/2) {
        r = n - r;
    }

    ull result = 1;
    ull common_divisor;
    for (int i = 1; i <= r; ++i) {
        common_divisor = gcd(result, i);
        result /= common_divisor;
        result *= (n - i + 1) / (i / common_divisor);
    }

    return result;
}

使用递归关系

ull n_choose_r_relation(ull n, ull r) {
    for (int i = 0; i <= n + 1; ++i) {
        for (int k = 0; k <= r && k <= i; ++k) {
            if (k == 0 || k == i) {
                ncr[i][k] = 1;
            }
            else {
                ncr[i][k] = (ncr[i - 1][k - 1] + ncr[i - 1][k]) % MODULO;
            }
        }
    }

    return ncr[n][r];
}

@Mysticial:抱歉让你感到困惑,我想说的是使用欧拉定理计算模反元素需要m为质数,可以通过fast_pow(n, m - 2, m)来计算。 - roxrook
啊,模乘逆并不需要m是质数。它只需要nm互质(GCD(n,m) == 1,其中n是你要求逆的数字)。 - Mysticial
@paulsm4:1)是预计算的,而3)是硬编码,即ull catalan[MAX] = { xxx,yyy,zzz .... };。感谢提供参考。 - roxrook
你能利用 (x*y) mod n = ((x mod n) * (y mod n)) mod n 这个事实吗?本质上,编写自己的 factorial_mod_n 和相应的 nCr_mod_n 函数,滥用这个事实来使用显式公式。 - Ben Thul
你认为 O(N lg N) {使用 O(N) 空间} 足够快吗? - Ben Voigt
显示剩余4条评论
2个回答

14

我写的答案是在这个有关计算nCr的问题关闭后写的,最终只出现在评论中:

我不确定这是绝对最快的方法,但应该足够高效。关键在于模乘法可以分解,但除法不能,因此首先必须将分数化简为以下形式:

  • 由于 n <= 10000,因此很可能构建大小为 2*n 的数组。

  • 使用埃拉托色尼筛法查找并存储所有小于 20000 的合数的因子对。这一步只需要执行一次,无论要计算多少个卡特兰数。

  • 创建另一个大小为 2*n 的表,表示每个因子的指数。

  • 现在,迭代卡特兰公式中的乘积enter image description here

  • 使用筛选表将每个因子分解为质因数,并为分子中的每个项增加指数表,在分母中的每个项减少指数表。

  • 没有任何条目会变成负数。

  • 现在,使用模运算将未被抵消的因子相乘。

  • 在任何时候都不需要进行除法运算。也没有分数。

  • 我的方法应用于 multi-nCr 的演示:http://ideone.com/Weeg6

如果您想用这个方法计算卡特兰数,可以使用这个方法替代 calc_combinations 中的循环:

    for( unsigned k = 2; k <= N; ++k ) {
         factor<+1>(k+N);
         factor<-1>(k);
    }

代码看起来像这样:http://ideone.com/ZZApk

解决方案

#include <utility>
#include <vector>

std::vector< std::pair<int, int> > factor_table;
void fill_sieve( int n )
{
        factor_table.resize(n+1);
        for( int i = 1; i <= n; ++i )
                factor_table[i] = std::pair<int, int>(i, 1);
        for( int j = 2, j2 = 4; j2 <= n; (j2 += j), (j2 += ++j) ) {
                if (factor_table[j].second == 1) {
                        int i = j;
                        int ij = j2;
                        while (ij <= n) {
                                factor_table[ij] = std::pair<int, int>(j, i);
                                ++i;
                                ij += j;
                        }
                }
        }
}

std::vector<unsigned> powers;

template<int dir>
void factor( int num )
{
        while (num != 1) {
                powers[factor_table[num].first] += dir;
                num = factor_table[num].second;
        }
}

void calc_catalan(unsigned N)
{
    powers.resize(0);
    powers.resize(2*N+1);
    for( unsigned k = 2; k <= N; ++k ) {
         factor<+1>(k+N);
         factor<-1>(k);
    }
}

测试驱动程序

#include <iostream>
#include <cmath>
int main(void)
{
    fill_sieve(20000);
    unsigned N = 9913;
    unsigned long long M = 1000000000007LL;
    calc_catalan(N);
    unsigned long long result = 1;
    for( unsigned i = 0; i < powers.size(); ++i ) {
        while (powers[i]--) {
            result *= i;
            result %= M;
        }
    }
    std::cout << "Catalan(" << N << ") modulo " << M << " = " << result << "\n\n";
}

完成的演示:http://ideone.com/FDWfB

这里有另一个相关问题,我用代码和演示进行了回答:C++中的组合数(N选择R)数量


顺便说一下,multi-nCr 的链接没有任何 ncr 函数?你是指错了链接吗?还是我哪里理解有误?非常感谢。 - roxrook
DP解决方案只需2.72秒,而您的解决方案需要32.xx秒。向量的分配/释放时间会导致问题。 - roxrook
嗯,在 ideone 上,10001 次运行只需要 2.75 秒:http://ideone.com/S8TKE。你也可以将 powers 的大小设置为最大,并在调用之间将其归零。但是在 ideone 上几乎没有什么区别:http://ideone.com/RhvhF。 - Ben Voigt
转换为数组后有所改善:http://ideone.com/OK1SD 并且摆脱所有向量:http://ideone.com/bJSUs 通过反转乘法循环略微提高了性能:http://ideone.com/hWQuT - Ben Voigt
可以通过仅在 i<50 时执行指数平方运算来提高一些速度(因为高因子不经常重复):http://ideone.com/DOvrp,但至少在 ideone 上仍然比朴素的重复乘法慢。 - Ben Voigt
显示剩余18条评论

4

简单易行。计算二项式系数的质因数是使用筛法的简单任务。我不会深入探讨其余部分,但幂模运算是微不足道的,甚至不需要除法。

对于N = 10000,我很快得到了42224403014400。

但是,如果您想要完整的数字本身,则第10000个Catalan数本身是...

    22453781249338521563359358425736057870110358621936588777329371383585
443658870053449099810271911432021020990539379958970114932732650095370271
397751300183876130693653440780258549445459994177372998459176454278220288
679699783327649549651476024591222065426709156831181207130089121989402216
517545144106669143509197596949973192167548893412063804651413496597406903
967719298471463870452875276986356795262033484770727452974197655810423629
386184662262278329466750526865120502476640878488187299740404235631962632
335108916990663560351330901464515744357084282208286669901241545533951877
777078174205283779947690623035078595904048715811899275348402286537327410
009576296851062523691528014340846065120667839872568170381150542379156626
173532955062796771718993285598391346886779480658586379448386923993317934
139425945651509102645665277040984870211604644540699508509248821099873225
565699224344151993874742555422872473424262356666363196825449089721410665
537521519676271082500130505509387186351879731113568837096419481746389018
721284533242225719341420124434480886444987373634542567071582458263380247
628252179873943804465262216365735901268165347321451279736504798992232739
106390706179212626442096326217616178171108663008963682821183764312867791
507672494716865305031842633900748973827504534625795968537648004286087039
823233370550650634239448544304798764239028734674653967478032618882557954
859328131980782727940394400855369003385513208814011609977239377877068501
893633819436630205358663340684840462204867552576509569736390978718963517
869423927523718504671005747648411794527978689778762460237949479732242725
154275831263823307362585789708343583184171797113785187466609433767144371
710845773715328364171910363978492352051901370003068055356444233141131383
192077598317531370925033378421138581148001529316546340657631162629562941
211065221871760353772365014435796695284269667873562415761642871681276498
507492541421942131281008978510862112693424595990036710403533420006771490
575482785612280198742983770649313043583275207213939274300662039637048647
395250014477941359641726047221826652916778311801541491816826072282488555
018173563867058868251361080516013361134986419403377613243853586312008767
909635869692823359899687030213634793656744420820912530014968355236934193
747181786083577435923400955703014812335311495073521773651461701750485101
119310472898683618090898735223665962918372501660743711042258315604294195
583076309209507444333462531858856911411408798540404888967120239682480627
570158137868956844950713279360385273144560292399045892610118082102910880
862332337854786916935223744892537176357434650161037841572213751901947447
479406915511862629144757855890852243043614898752155191154178797427659170
858428903659564218086017881546286273599385917718058276038925354040884258
022546721698832195059172836919416429064599278227491956109630837263590884
232587058023101145921693423507849076470763334833613166731358258440439729
023251976962577737416518794914009277934381234511794730677137605309953636
716963188964230436087118746073758080815722286112796870306754227017546055
347853334923811143440952672436342961180384459596879312187164969968096364
679341577416027452001090523659332406246454292701122715894579618818643071
139925009651888661718404932582731927646801878919152052218535889565319288
284306134970608577076704660104569794464663831193002735423564364371354521
236158069405955372080665906666149641642367693009585743888230289135078928
729184475260174446278915850624301208853693618442212023236924456444468934
014289741543223145235333811594418344798647068944904371005158995839127368
111629241573877617157577569590584624720552246920280151741755137476154967
741272080362312952750328628775530857638646138592895858764915987201920286
661490154786097488396300779244279606416541720716707237058679072236693234
932525387774462125138686406910133757255779021404876020200833761157767584
015369673586027681003369474431448843539054790848335705489738731700240579
310855452462903455809888697753847348175077261616431384533713924568807999
599683993362082982833949280082553659996487889394727840889035163412693106
865702752400579571351436509808650503057036278511515529330634352096987240
087618010503197530225589878764240330302768263496958673020211712107611762
945771002810537812467742009399047607169797035466100221770262334445478074
080845928677855301631860443068261061887109865290453732333638130446973519
286828584088203627113605849939106943614542645022903932947597417823646592
053417189520415596451505598330301782369213897762201629272201936584136036
027455748892667375417522206148332891409959866390232031014358337935412166
499617373308661369292739138448626161089231445046384163766705419698533262
040353901193260661841441922949263756492472641127072018961101915467728184
640938751407261817683231072132781927769994322689591991504965204544928105
747119997826784396172488376877215547707335474490892399544875233372674064
229287210750045834971802632275569822679385098328070604595140732389126327
092826465756212595551194678295464565601548041854366455751504169209131794
100099734293551231149329072243438440125013340293416345726479426178738686
238273833019523777019099811511419301476900607138083408535229058593795242
998150989330379630607152057165593682028276808657989133687600036850256257
973833780907105126134335912174477305526445570101413725539992976023375381
201759604514592679113676113078381084050224814280307372001545194100603017
219283437543128615425515965977881708976796492254901456997277712672653778
789696887633779923567912536882486775488103616173080561347127863398147885
811314120272830343521897029277536628882920301387371334992369039412492040
272569854478601604868543152581104741474604522753521632753090182704058850
525546680379379188800223157168606861776429258407513523623704438333489387
460217759660297923471793682082742722961582765796049294605969530190679149
426065241142453853283673009798518752237906836442958353267589634936329512
043142900668824981800672231156890228835045258196841806861681826866706774
199447245550164975361170844597908233890221446745462710788815648943858461
7793175431865532382711812960546611287516640

谢谢,经过我复习了数论笔记后,现在我能够更清晰地看到它了。 - roxrook

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接