What is the fastest (known) algorithm to find the n-th Catalan number mod m?

Question

The problem is to find the n-th Catalan number mod m, where m is NOT prime, m = (10^14 + 7). Here are the list of methods that I have tried: (max N = 10,000)

Dynamic programming for table look-up, too slow
Use Catalan formula ncr(2*n, n)/(n + 1), again it wasn't fast enough due to the ncr function, can't speed up using exponentiation squaring because m is not prime.
Hardcode a table of pre-generated Catalans, but it failed due to the file size limit.
Recurrence relation C(i,k) = C(i-1,k-1) + C(i-1,k), this is way too slow

So I wonder is there any other faster algorithm to find the n-th Catalan number that I'm not aware of?

Using Dynamic Programming

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;
    }
}

Using original formula

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;
}

Using recurrence relation

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];
}

Exponentiation squaring has nothing to do with m being prime. It can be done on all modulus. — Mysticial, Aug 08 '12 at 21:23
Q: Homework? What does "dynamic programming for table look-up" mean (and how does it differ from 3))? Have you looked it up in your copies of Sedgewick or Knuth, or "Numerical Algorithms"? Also check here: http://math.stackexchange.com/questions/161243/ballot-numbers-sum-up-to-catalan-numbers — paulsm4, Aug 08 '12 at 21:23
@Mysticial: Sorry for the confusion, I meant to say the **Modular Multiplicative Inverse** using Euler theorem require `m` to be a prime which can be calculated as `fast_pow(n, m - 2, m)` — roxrook, Aug 08 '12 at 21:26
Ah. Modular multiplicative inverse doesn't require `m` to be prime. It only requires that `n` and `m` be relatively prime. (`GCD(n,m) == 1`, where `n` is the number you're trying to invert) — Mysticial, Aug 08 '12 at 21:28
@paulsm4: 1) is pre-computed, while 3) is hardcode, i.e. `ull catalan[MAX] = { xxx, yyy, zzz .... };`. Thanks for the reference. — roxrook, Aug 08 '12 at 21:33
Can you leverage the fact that (x*y) mod n = ((x mod n) * (y mod n)) mod n? Essentially, write your own factorial_mod_n and corresponding nCr_mod_n functions that (ab)use that fact to use the explicit formula. — Ben Thul, Aug 08 '12 at 21:36
Do you think O(N lg N) {with O(N) space usage} would be fast enough? — Ben Voigt, Aug 08 '12 at 21:41
@BenVoigt: There are about 10000 test cases and max `N = 10000` while the time limit is 3 seconds; however I'm not sure how they measure this time limit. Could you share your idea? any `log(N)` algorithm sounds appealing to me. — roxrook, Aug 08 '12 at 21:52
@Mysticial: Factorials of N (mod M) can be simplified if M and N are not relatively prime. — Wug, Aug 08 '12 at 21:54

score 14 · Answer 1 · edited May 23 '17 at 11:45

14

From an answer I wrote just as this question about computation of nCr got closed, which ended up in the comments instead:

I'm not sure this is the absolute fastest, but it should be efficient enough. The key is that modulo multiplication decomposes, but division doesn't, so one first must reduce the fraction, as follows:

Since n <= 10000, it's very possible to build an array of size 2*n.
Use the Sieve of Eratosthenes to find and store pairs of factors for all composite numbers up to 20000). This step only needs to be done once regardless of how many Catalan numbers are to be calculated.
Make another table of size 2*n which represents the exponent of each factor.
Now, iterate the product in the Catalan formula .
Break each factor into prime factors using the sieve table, incrementing the exponent table for each term in the numerator and decrementing it for each term in the denominator.
No entry will ever end up negative.
Now, use modulo arithmetic to multiply together the non-cancelled factors.
No division operations are required at any point. Nor any fractions.
Demonstration of my method applied to multi-nCr: http://ideone.com/Weeg6

To use this for Catalan numbers, you can use this instead of the loops inside calc_combinations:

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

The code then looks like this: http://ideone.com/ZZApk

Solution

#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);
    }
}

Test Driver

#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";
}

Completed demo: http://ideone.com/FDWfB

Here's another related question, which I answered with code and demonstration: Number of combinations (N choose R) in C++

edited May 23 '17 at 11:45

Community

1
1

answered Aug 08 '12 at 22:43

Ben Voigt

277,958
43
419
720

Thanks a lot, by the way the link to `multi-nCr` doesn't have any `ncr` function? Are you pointing to the wrong link? or I misunderstood somewhere? – roxrook Aug 08 '12 at 22:55
@Chan: `calc_combinations` is a multi-bin nCr function, splitting N elements among K bins of size r1, r2, r3, ... rk. You'll have only two bins, making it quite a bit easier, but you'll have to add the modulo multiplication step at the end (I just printed the result out in text as a product of prime factors). – Ben Voigt Aug 08 '12 at 22:57
I'm quite confused with the code, could you add a quick minimal example of `ncr` says `ncr(20, 5)`? – roxrook Aug 08 '12 at 23:00
@Chan: Sure. I just changed the line in `main` that sets `bin_sizes`.... http://ideone.com/5L9LU – Ben Voigt Aug 08 '12 at 23:01
Thanks a lot. I'm eager to try it out, I'll let you know. By the way, your bin elements are `15, 5`, is it equivalent to `ncr(15, 5)`? – roxrook Aug 08 '12 at 23:02
@Chan: No. If there's one bin of size 5 and one bin of size 15, there's 20 total items being distributed. Anyway, I added a Catalan demonstration to the answer also. – Ben Voigt Aug 08 '12 at 23:05
That would be great! Thanks again. – roxrook Aug 08 '12 at 23:06
@Chan: Ok, all done. Please run some tests and let me know if it works. – Ben Voigt Aug 08 '12 at 23:11
I got vector subscript out of range exception, something is wrong. Could you double check? – roxrook Aug 09 '12 at 00:02
The one above, I'm using Visual Studio C++ 2010. – roxrook Aug 09 '12 at 00:07
@Chan: Ok, I did find a problem. It needs to be `powers.resize(2*N+1);` instead of `N+1`. – Ben Voigt Aug 09 '12 at 00:09
Thanks once again, by the way it still got TLE though. – roxrook Aug 09 '12 at 00:16
How does the runtime compare to your earlier attempts? – Ben Voigt Aug 09 '12 at 00:19
I just realized it hung up when I called `calc_catalan(n)` the second time. I ran a max-test from 1->10000, after printing the first one, it was freeze. If I understand your method correctly, I only need to call `fill_sieve` once right? – roxrook Aug 09 '12 at 00:24
@Chan: Yes, only call `fill_sieve` once. I'll look and see if anything obviously isn't being reset. – Ben Voigt Aug 09 '12 at 00:29
Try `powers.resize(0)` in between runs. – Ben Voigt Aug 09 '12 at 00:31
The `powers` vector got messed up, there are some very big numbers in that vector after the first call. Anyway, I will implement `exponentiation squaring` for that part as well. – roxrook Aug 09 '12 at 00:31
The DP solution took only 2.72 seconds while your solution is 32.xx seconds. The allocation/de-allocation time for vector would cause the issue. – roxrook Aug 09 '12 at 00:37
Hmmm, on ideone, it takes only 2.75 seconds for 10001 runs: http://ideone.com/S8TKE. You could alternatively make `powers` the maximum size always and just zero it between calls. But that makes almost no difference on ideone: http://ideone.com/RhvhF – Ben Voigt Aug 09 '12 at 00:38
Switching to an array made an improvement though: http://ideone.com/OK1SD and getting rid of all vectors: http://ideone.com/bJSUs Slight additional improvement from reversing the multiply loop: http://ideone.com/hWQuT – Ben Voigt Aug 09 '12 at 00:41
Thanks for the effort, let's see if it works, but let me try to add exponentiation squaring instead of the using the `while (powers[i])` loop. – roxrook Aug 09 '12 at 00:46
@Chan: Yeah, that sounds like a big win. But in my testing it is worse. http://ideone.com/iBD6M – Ben Voigt Aug 09 '12 at 00:55
Can get some speed back by only doing exponentiation squaring when `i<50` (since high factors are not often repeated): http://ideone.com/DOvrp, but still slower than the naive repeated multiplication, at least on ideone. – Ben Voigt Aug 09 '12 at 01:00

score 4 · Accepted Answer · 2012-08-10T00:25:28.627

Easy, peasy. Compute the prime factors of the binomial coefficient. A simple task using a sieve. I won't get into the rest of it, but a powermod computation is trivial, and you don't even need a divide.

For N = 10000, I get 42224403014400 in short order.

But, if you want the full number itself, the 10000'th Catalan number itself is ...

    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

Thanks, after reviewing my Number Theory notes, I can see it much clearer now. — roxrook, Aug 09 '12 at 16:44

What is the fastest (known) algorithm to find the n-th Catalan number mod m?

2 Answers2

Solution

Test Driver

Linked

Related