KONKURS - Maksymalne przyspieszenie kodu.

Witam!

Konkurs polega na maksymalnym przyspieszeniu podanego poniżej kodu. Był on napisany przeze mnie dawno temu, dlatego proszę o wyrozumiałość. W każdym razie, działa.

REGULAMIN:
0) Organizatorem konkursu jest użytkownik merlinnot.

1. Regulamin może być w każdej chwili zmieniony/poszerzony przez Organizatora.
2. Wielowątkowość jest niedozwolona.
3. Wstawki ASM są dozwolone.
4. Programy są testowane przez organizatora na czystej instalacji systemu Ubuntu Linux 12.10, z najnowszą wersją kompilatora GCC.
5. Kody mają być przygotowane w formie pliku KONKURS.cpp, pierwszy post pod tym tematem jest jednocześnie kodem do przerobienia i wzorem postu z rozwiązaniem.
6. Można wysyłać wiele rozwiązań (w granicach rozsądku).

Powodzenia!

KONKURS.cpp

#include <cstdio>
#include <math.h>
#include <vector>
#include <algorithm>
#define lli long long int
std::vector<unsigned int>MI;
std::vector<unsigned int>LI;
bool tmp;
int ROZ(unsigned long, bool &);
int WYNIK();
int WN(int, int);
inline void RETR();
int main ()
{
	int a, b, ltestow;
	scanf("%d", &ltestow);
	while(ltestow--)
	{
		scanf("%d%d", &a, &b);
		MI.clear();
		LI.clear();
		WN(a, b);
		std::sort(MI.begin(), MI.end());
		std::sort(LI.begin(), LI.end());
		RETR();
		printf("%d\n", WYNIK());
	}
}

int ROZ(unsigned long n, bool & ML)
{
	if(n>1) 
	{
		unsigned int i=2; 
		while (i<=(unsigned long)sqrt((double)n))
		{
			while(!(n%i))
			{
				n/=i; 
				if(ML)
					MI.push_back(i);
				else
					LI.push_back(i);
			}
			if(n==1) break;
				i++;
		}
		if(n>1) 
			if(ML)
				MI.push_back(n);
			else
				LI.push_back(n);
	}
	return 0;
}
int WN(int n, int k)
{
	int li=1, mi=1;
	for(int i=1; i<=k; i++)
	{
		tmp=false;
    		ROZ(n+1-i, tmp);
		tmp=true;
			ROZ(i, tmp);
	}
	return li/mi%10000;
}
inline void RETR()
{
	for(int i=MI.size(); i>0; i--)
	{
		LI.erase(lower_bound(LI.begin(), LI.end(), MI[0]));
        	MI.erase(MI.begin());
	}
}
int WYNIK()
{
	int retr=1;
	for(int i=0; i<LI.size(); i++)
	{
		retr*=LI[i];
		retr%=10000;
	}
	return retr;
} 

Dział Edukacja??? To już flame było lepsze.

Z ciekawości - co robi ten kod - jakiś konkretny algorytm? Zamierzam go przepisać na coś bardziej odpowiadającego mojej wrażliwej estetyce :>.

Wydaje mi sie czy long long int dziala tylko w Visualu?

$ echo "123434 54234" | time ./wzor
8000
36.75user 0.00system 0:36.76elapsed 99%CPU (0avgtext+0avgdata 10080maxresident)k
0inputs+0outputs (0major+1125minor)pagefaults 0swaps
$ echo "123434 54234" | time ./moj
8000
0.44user 0.01system 0:00.46elapsed 100%CPU (0avgtext+0avgdata 18848maxresident)k
0inputs+0outputs (0major+1418minor)pagefaults 0swaps

chyba niezła optymalizacja ;) gotowe programy mamy wrzucać na forum czy gdzieś wysyłać?

IdeOne zameldowało: C++0x (gcc-4.5.1)

wynik: sukces	     czas: 0s    pamięć: 6852 kB     zwrócona wartość: 0

wejście: brak
wyjście:
(123434,54234) = 8000

inny wynik (1234567,617283) = 0 czas czas: 0.02s
gdy zmieniłem tak, że liczy dla b<=a<=100 000 000
policzenie (123434,54234) zajęło 0,02s, natomiast (100000000,50000000) = 7040 trwało 2.31s
(a Wolfram się poddał :-) )

"Ten wątek został przeniesiony..." może ktoś wyjaśni ten dowcip? bo jakoś nie chwytam
forum dla zainteresowanych programowaniem, zwykle amatorsko, ale mamy i kilku zawodowców i ...
topic dotyczący algorytmiki, programowania i innych takich jest "off topic"
kretynieje społeczeństwo czy tylko ja?

Twoim zdaniem wątek o konkursie pasował do działu Edukacja?

Opis działu Edukacja: Tutaj możesz podyskutować o książkach, nauce programowania, kursach
Opis działu Off-Topic: Miejsce na dyskusje niepasujące do pozostałych kategorii forum
Tak naprawdę jeszcze lepiej byłoby w dziale C++, ale tam mogłoby się zagubić...

Zapraszam do przysyłania rozwiązań na
01101101 01100101 01110010 01101100 01101001 01101110 01101110 01101111 01110100 01000000 01100111 01101101 01100001 01101001 01101100 00101110 01100011 01101111 01101101, pamiętajcie, że możecie wysłać kilka zgłoszeń :)

Dla bardziej wnikliwych:
MD5 maila: c401e8ae9798deb87574dfa90878963b;
crc32: 1357118514;
sha1: 564487b98f9cdac2d78fa02c55eba97166d9a19f;
;]

Może kup jakiegoś VPSa za te ciężko zarobione 10 złotych, postaw mały frontend do wszystkiego i szerzej rozreklamuj konkurs? Za opłatą mogę pomóc, kilka razy zajmowałem się organizacją takich rzeczy (konkurs na najbardziej poryty kod i eksperyment fizzbuzz).

merlinnot napisał(a):

KONKURS.cpp

#include <cstdio>
//CIACH

Kod z załącznika różni się od kodu w poście. To w końcu która wersja jest właściwa?

I jaki jest zakres danych wejściowych?

Kod w załączniku wykonuje obliczenia dla jednej pary liczb, natomiast zamieszczony w poście może liczyć dla wielu.
Same obliczenia wykonywane są identycznie.

Należy przyjąć, że: 0 <= b <= a, chyba też, że a*10000 < 2^64, ale liczyłoby się to koszmarnie długo.

Częścią zagadki jest odgadnięcie co właściwie liczymy, a to może być trudne dla naszych młodszych koleżanek i kolegów. Proponuję by nie kończyć konkursu w niedzielę, a tylko powiedzmy zamknąć pierwszy etap i ujawnić co to za obliczenia, bo w końcu problem jest wszystkim doskonale znany.

autor "konkursu" nieźle zakręcił problem. Jeśli już się zrozumie co to robi, to optymalne rozwiązanie można znaleźć tu na forum, było to wałkowane kilka razy, a jakieś 6mc temu to ktoś robił poważnie testy jaka metoda jest najszybsza. Spokojnie można to wyliczyć w czasie o(a-b).

Drodzy forumowicze, na maila dostałem pierwsze zgłoszenie, wieczorem testuję i wrzucam wyniki! Widać, że dla niektórych zagadka nie stanowi najmniejszego problemu :)

bez podania zakresu danych wejściowych, to nie może być jeden konkurs tylko wiele.

Testów jest nie więcej niż 10^6, 0 <= b <= a dla pojedynczego testu. I mieści się w incie.

Wyniki na dzień dzisiejszy:

Kod konkursowy:
8m 11.751s

Kod użytkownika dawidgarus:
0m 0.250s

Kod użytkownika Xitami:
0m 0.988s

Zostały dwie godziny, lecz jeżeli ktoś jeszcze miałby ochotę się za to zabrać, zawsze mogę przedłużyć trochę konkurs :)

możesz podać czasy dla (2 000 000 000, 1 000 000 000)?

Oczywiście :)

Xitami:

merlinnot@merlinnot:~/Pobrane$ time ./KONKURS < hinput.in > input2.out

real	1m19.630s
user	1m19.553s
sys	0m0.056s

dawidgarus:

merlinnot@merlinnot:~/Pobrane$ time ./KONKURS < hinput.in > input2.out

real	0m33.641s
user	0m33.630s
sys	0m0.000s
 

tak sobie myślałem, że kogoś to obchodzi
pewnie nikogo poza nami trzema, chętnie zobaczyłbym. Mój kod to http://ideone.com/CWLwL

Chętnie zamieściłbym kod Dawida, ale nie zrobię tego bez jego zgody.

Tak paczając w ten jego kod, to *le mindfuck totale.

Chętnie spotkam się z konstruktywną krytyką mojego
zamieść
może ktoś powiesi jakieś mądre psy

Mam nadzieję, że nie będzie miał mi tego za złe, w końcu już wygrał, a inni mogą się z tego wiele nauczyć.

#include <stdio.h>
#include <stdlib.h>

/*int invMod(unsigned a, unsigned b) {
        const unsigned m = a;
        register unsigned c, quot;
        register int p=1, q=0, r=0, s=1;
        while(b) {
                quot = a/b;
                c = p - quot*r;
                p = r;
                r = c;
                c = q - quot*s;
                q = s;
                s = c;
                c = a % b;
                a = b;
                b = c;
        }
        return (q+m) % m;
        //return q>0?q:(q+m);
}*/

struct data {
        unsigned n, k;
        unsigned d2, d5;
        unsigned res1, res2;
};

struct datacontainer {
        unsigned k;
        data *d;
};

int dccmp(const void *a, const void *b) {
        return (long)((datacontainer*)a)->k - (long)((datacontainer*)b)->k;
}

int main(void) {
        unsigned n, k, l;
        unsigned num;
        unsigned invMod16[] = {1, 11, 13, 7, 9, 3, 5, 15};
        unsigned invMod625[] = {1,
   1,  313,  417,  469,  0,  521,  268,  547,  139,  0,  341,  573,  577,  134,  0,  586,  478,  382,  329,  0,
 506,  483,  462,  599,  0,  601,  463,   67,  194,  0,  121,  293,  322,  239,  0,  191,  473,  477,  609,  0,
  61,  253,  407,  554,  0,  231,  133,  612,  574,  0,  576,  613,  342,  544,  0,  346,  318,   97,  339,  0,
  41,  373,  377,  459,  0,  161,   28,  432,  154,  0,  581,  408,  137,  549,  0,  551,  138,  617,  269,  0,
 571,  343,  497,  439,  0,  516,  273,  277,  309,  0,  261,  428,  457,  379,  0,  306,   58,  287,  524,  0,
 526,  288,  267,  619,  0,  171,  368,  272,  539,  0,  366,  173,  177,  159,  0,  361,  203,  482,  604,  0,
  31,  333,  437,  499,  0,  501,  438,  542,  344,  0,  396,  393,   47,   14,  0,  216,   73,   77,    9,  0,
 461,  603,  507,  204,  0,  381,  608,  587,  474,  0,  476,  588,  192,   69,  0,  621,  418,  447,  114,  0,
  66,  598,  602,  484,  0,  561,  378,  532,  429,  0,  106,  258,  112,  449,  0,  451,  113,  467,  419,  0,
 221,  443,  222,  214,  0,  541,  498,  502,  334,  0,   36,  153,  557,   29,  0,  456,  533,  262,  424,  0,
 426,  263,  117,  144,  0,  446,  468,  622,  314,  0,  391,  398,  402,  184,  0,  136,  553,  582,  254,  0,
 181,  183,  412,  399,  0,  401,  413,  392,  494,  0,   46,  493,  397,  414,  0,  241,  298,  302,   34,  0,
 236,  328,  607,  479,  0,  531,  458,  562,  374,  0,  376,  563,   42,  219,  0,  271,  518,  172,  514,  0,
  91,  198,  202,  509,  0,  336,  103,    7,   79,  0,  256,  108,   87,  349,  0,  351,   88,  317,  569,  0,
 496,  543,  572,  614,  0,  566,   98,  102,  359,  0,  436,  503,   32,  304,  0,  606,  383,  237,  324,  0,
 326,  238,  592,  294,  0,   96,  568,  347,   89,  0,  416,  623,    2,  209,  0,  536,  278,   57,  529,  0,
 331,   33,  387,  299,  0,  301,  388,  242,   19,  0,  321,  593,  122,  189,  0,  266,  523,  527,   59,  0,
  11,   53,   82,  129,  0,   56,  308,  537,  274,  0,  276,  538,  517,  369,  0,  546,  618,  522,  289,  0,
 116,  423,  427,  534,  0,  111,  453,  107,  354,  0,  406,  583,   62,  249,  0,  251,   63,  167,   94,  0,
 146,   18,  297,  389,  0,  591,  323,  327,  384,  0,  211,  228,  132,  579,  0,  131,  233,  212,  224,  0,
 226,  213,  442,  444,  0,  371,   43,   72,  489,  0,  441,  223,  227,  234,  0,  311,    3,  157,  179,  0,
 481,  508,  362,  199,  0,  201,  363,   92,  169,  0,  596,   68,  472,  589,  0,  291,  123,  127,   84,  0,
 411,  403,  182,  404,  0,  206,  158,  512,  174,  0,  176,  513,  367,  519,  0,  196,   93,  247,   64,  0,
 141,   23,   27,  559,  0,  511,  178,  207,    4,  0,  556,  433,   37,  149,  0,  151,   38,   17,  244,  0,
 421,  118,   22,  164,  0,  616,  548,  552,  409,  0,  611,  578,  232,  229,  0,  281,   83,  187,  124,  0,
 126,  188,  292,  594,  0,   21,  143,  422,  264,  0,  466,  448,  452,  259,  0,   86,  353,  257,  454,  0,
   6,  358,  337,   99,  0,  101,  338,  567,  319,  0,  246,  168,  197,  364,  0,  316,  348,  352,  109,  0,
 186,  128,  282,   54,  0,  356,    8,  487,   74,  0,   76,  488,  217,   44,  0,  471,  193,  597,  464,  0,
 166,  248,  252,  584,  0,  286,  528,  307,  279,  0,   81,  283,   12,   49,  0,   51,   13,  492,  394,  0,
  71,  218,  372,  564,  0,   16,  148,  152,  434,  0,  386,  303,  332,  504,  0,  431,  558,  162,   24,  0,
  26,  163,  142,  119,  0,  296,  243,  147,   39,  0,  491,   48,   52,  284,  0,  486,   78,  357,  104,  0,
 156,  208,  312,  624};
 // tablica z odwrotnosciami modulo 625, zeby generowac dynamicznie odkomentuj kod ponizej i funkcje invMod
        register unsigned i, j;
        register unsigned d2, d5;
        register unsigned result1, result2;
        unsigned maxk;
        /*int invMod625[625];
        for (i=1; i<625; ++i) {
                if (i % 5) {
                        invMod625[i] = invMod(625, i);
                //printf("% 4d, ", invMod625[i]);
                //} else {
                        //printf(" 0, ");
                }
                //if (i % 20 == 0) printf("\n");
        }*/

        scanf("%d", &num);
        data *d = (data*)malloc(sizeof(data)*num);
        datacontainer *dc = (datacontainer*)malloc(sizeof(datacontainer)*num);
	if (!d || !dc) {
		fprintf(stderr, "Nie mozna zaalokowac pamieci\n");
		exit(0);
	}
        maxk = 0;
        for (l=0; l<num; ++l) {
                scanf("%d%d", &n, &k);
                if (k*2 > n) k = n-k;
                if (k > maxk) maxk = k;
                d[l] = (data){n, k};
                dc[l] = (datacontainer){k, d+l};
        }
        qsort(dc, num, sizeof(datacontainer), dccmp);
        d2 = d5 = 0;
        result1 = result2 = 1;
	l = 0;
	while(dc[l].k == 0) {
		dc[l].d->d2 = 0;
		dc[l].d->d5 = 0;
		dc[l].d->res1 = 1;
		dc[l].d->res2 = 1;
		++l;
	}
        for (i=1; i<=maxk; ++i) {
                for(j=i; ~j&1; j=j>>1, ++d2);
                result1 = (result1 * invMod16[(j&15)>>1]) &15;
                for(j=i; j%5==0; j/=5, ++d5);
                result2 = (result2 * invMod625[j%625]) % 625;
                while (dc[l].k == i) {
                        dc[l].d->d2 = d2;
                        dc[l].d->d5 = d5;
                        dc[l].d->res1 = result1;
                        dc[l].d->res2 = result2;
                        ++l;
                }
        }
        for(l=0; l<num; ++l) {
                n = d[l].n;
                k = d[l].k;
                d2 = d[l].d2;
                d5 = d[l].d5;
                result1 = d[l].res1;
                result2 = d[l].res2;
                for (i=0; i<k; ++i) {
                        for(j=n-i; d2 && ~j&1; j=j>>1, --d2);
                        result1 = (result1 * (j&15)) &15;
                        for(j=n-i; d5 && j%5==0; j/=5, --d5);
                        result2 = (result2 * (j % 625)) % 625;
                }
                d2 = ((result2 - result1) % 625) * 586;
                result2 = (result1 + (d2<<4)) % 10000;
                printf("%d\n", result2);
        }
        return 0;
}
 

mvt8 napisał(a):

1. ludzie naucza sie pisac fajny kod

Tja, wręcz zjawiskowy:

        unsigned invMod625[] = {1,
   1,  313,  417,  469,  0,  521,  268,  547,  139,  0,  341,  573,  577,  134,  0,  586,  478,  382,  329,  0,
 506,  483,  462,  599,  0,  601,  463,   67,  194,  0,  121,  293,  322,  239,  0,  191,  473,  477,  609,  0,
  61,  253,  407,  554,  0,  231,  133,  612,  574,  0,  576,  613,  342,  544,  0,  346,  318,   97,  339,  0,
  41,  373,  377,  459,  0,  161,   28,  432,  154,  0,  581,  408,  137,  549,  0,  551,  138,  617,  269,  0,
 571,  343,  497,  439,  0,  516,  273,  277,  309,  0,  261,  428,  457,  379,  0,  306,   58,  287,  524,  0,
 526,  288,  267,  619,  0,  171,  368,  272,  539,  0,  366,  173,  177,  159,  0,  361,  203,  482,  604,  0,
  31,  333,  437,  499,  0,  501,  438,  542,  344,  0,  396,  393,   47,   14,  0,  216,   73,   77,    9,  0,
 461,  603,  507,  204,  0,  381,  608,  587,  474,  0,  476,  588,  192,   69,  0,  621,  418,  447,  114,  0,
  66,  598,  602,  484,  0,  561,  378,  532,  429,  0,  106,  258,  112,  449,  0,  451,  113,  467,  419,  0,
 221,  443,  222,  214,  0,  541,  498,  502,  334,  0,   36,  153,  557,   29,  0,  456,  533,  262,  424,  0,
 426,  263,  117,  144,  0,  446,  468,  622,  314,  0,  391,  398,  402,  184,  0,  136,  553,  582,  254,  0,
 181,  183,  412,  399,  0,  401,  413,  392,  494,  0,   46,  493,  397,  414,  0,  241,  298,  302,   34,  0,
 236,  328,  607,  479,  0,  531,  458,  562,  374,  0,  376,  563,   42,  219,  0,  271,  518,  172,  514,  0,
  91,  198,  202,  509,  0,  336,  103,    7,   79,  0,  256,  108,   87,  349,  0,  351,   88,  317,  569,  0,
 496,  543,  572,  614,  0,  566,   98,  102,  359,  0,  436,  503,   32,  304,  0,  606,  383,  237,  324,  0,
 326,  238,  592,  294,  0,   96,  568,  347,   89,  0,  416,  623,    2,  209,  0,  536,  278,   57,  529,  0,
 331,   33,  387,  299,  0,  301,  388,  242,   19,  0,  321,  593,  122,  189,  0,  266,  523,  527,   59,  0,
  11,   53,   82,  129,  0,   56,  308,  537,  274,  0,  276,  538,  517,  369,  0,  546,  618,  522,  289,  0,
 116,  423,  427,  534,  0,  111,  453,  107,  354,  0,  406,  583,   62,  249,  0,  251,   63,  167,   94,  0,
 146,   18,  297,  389,  0,  591,  323,  327,  384,  0,  211,  228,  132,  579,  0,  131,  233,  212,  224,  0,
 226,  213,  442,  444,  0,  371,   43,   72,  489,  0,  441,  223,  227,  234,  0,  311,    3,  157,  179,  0,
 481,  508,  362,  199,  0,  201,  363,   92,  169,  0,  596,   68,  472,  589,  0,  291,  123,  127,   84,  0,
 411,  403,  182,  404,  0,  206,  158,  512,  174,  0,  176,  513,  367,  519,  0,  196,   93,  247,   64,  0,
 141,   23,   27,  559,  0,  511,  178,  207,    4,  0,  556,  433,   37,  149,  0,  151,   38,   17,  244,  0,
 421,  118,   22,  164,  0,  616,  548,  552,  409,  0,  611,  578,  232,  229,  0,  281,   83,  187,  124,  0,
 126,  188,  292,  594,  0,   21,  143,  422,  264,  0,  466,  448,  452,  259,  0,   86,  353,  257,  454,  0,
   6,  358,  337,   99,  0,  101,  338,  567,  319,  0,  246,  168,  197,  364,  0,  316,  348,  352,  109,  0,
 186,  128,  282,   54,  0,  356,    8,  487,   74,  0,   76,  488,  217,   44,  0,  471,  193,  597,  464,  0,
 166,  248,  252,  584,  0,  286,  528,  307,  279,  0,   81,  283,   12,   49,  0,   51,   13,  492,  394,  0,
  71,  218,  372,  564,  0,   16,  148,  152,  434,  0,  386,  303,  332,  504,  0,  431,  558,  162,   24,  0,
  26,  163,  142,  119,  0,  296,  243,  147,   39,  0,  491,   48,   52,  284,  0,  486,   78,  357,  104,  0,
 156,  208,  312,  624};
 

Szybkie pewno to jest, ale na pewno nie fajne.

to ja może jutro

to ja może wytłumaczę działanie mojego kodu dla zainteresowany. chcemy policzyć n po k mod 10000.
10000 = 2^{4 * 5}4 = 16 * 625
liczę osobno n po k mod 16 i n po k mod 625 i korzystam z chińskiego twierdzenia o resztach aby wyliczyć z tego n po k mod 10000.
dlaczego tak dzielę? gdy liczę n po k mod p^m (pierwsza do potęgi m) to mogę zapisać to jako:
[n(n-1)...(n-k+1)] / [12...k] mod p^{m. dla mianownika dla wszystkich liczb niepodzielnych przez p mogę policzyć odwrotność modulo p}m:
np.:
4 po 3 mod 2 = (4 * 3 * 2) / (1 * 2 * 3) mod 2 = [4 * 3 * 2 * (1^{-1 mod 2) * (3}-1 mod 2)]/2 mod 2
pozostaje tylko pozbyć się tych 2 z mianownika. liczymy ile ich jest i w liczniku znajdujemy liczby podzielne przez 2 i je dzielimy (np. 4)
232(1^{-1 mod 2)*(3}-1 mod 2) mod 2 a to już prosto wyliczyć.

w tablicy struktur data trzymam dane wejściowe. nie muszę za każdym razem liczyć ile wynosi iloczyn odwrotności modulo dla liczby od 1 do k. liczę raz dla od 1 do kmax (największy k z danych wejściowych) i pośrednie wyniki zapisuje w strukturze data odpowiadającej danemu k w result1 dla mod 16 i result2 dla mod 625. w d2 i d5 trzymam też ile razy liczba 2 i 5 mieściła się w liczniku, żeby potem liczby z mianownika przez nie podzielić. datacontainer służy po to, aby dane posortować i żeby po kolei wstawiać wyliczony iloczyn, ale zachować kolejność gdy potem będę liczył iloczyn w mianowniku i drukował wyniki.

zrobię przykład. liczymy 10 po 6 mod 16
(10 * 9 * 8 * 7 * 6 * 5) / (1 * 2 * 3 * 4 * 5 * 6) mod 16
zajmijmy się mianownikiem:
result1 = 1
d2 = 0
1 ^-1 mod 16 = 1 | result1 = (result1 * 1) % 16
2 ^-1 mod 16 - nie da się polczyć | d2++
3 ^-1 mod 16 = 11 | result1 = (result1 * 11) % 16
4 ^-1 mod 16 - nie da się policzyć | d2+=2
5 ^-1 mod 16 = 13 | result1 = (result1 * 13) % 16
6 ^{-1 mod 16 - nie da się policzyć, ale 6 = 2*3, a 3}-1 mod 16 = 11 | d2++, result = (result1 * 11) % 16
czyli d2 = 4, result1 = 5
teraz zajmujemy się licznikiem
10) d2>0 i 10 jest podzielne przez 2. d2-- i zostajemy z 5 | result1 = (result1 * 5) % 16
9) 9 nie jest podzielne przez 2 | result1 = (result1 * 9) mod 16
8) d2>0 i 8 jest podzielne przez 2^3. d2-=3 i zostajemy z 1 | result1 = (result1 * 1) % 16
// d2 == 0 czyli juz nie musimy dzielic:
result1 = (result1 * 7) % 16
result1 = (result1 * 6) % 16
result1 = (result1 * 5) % 16
result1 = 2

10 po 6 mod 16 = 2
10 po 6 mod 625 = 210 (nie chce mi sie liczyć)
naszym rozwiązaniem jest liczba x spełniająca:
x mod 16 = 2
x mod 625 = 210
pierwsze równanie zapisuje jako:
x = 2+16i
podstawiam do drugiego:
2+16i mod 625 = 210
16i mod 625 = 208 | mnożę stronami przez 16^-1 mod 625 = 586
i mod 625 = 208586 mod 625
i = 208586 mod 625
i = 13
wstawiam do pierwszego i mam:
x = 2 + 16*13 = 210
10 po 6 mod 10000 = 210

pierwszy pomysł to liczenie takie jak dla małych liczb (jak tu newton()), ale nie wszystko da się odwrócić (podzielić) mod 10000 to to odłożyłem, a wystarczyło pomyśleć o krok dalej. Jeszcze do tego wrócę (z Olem).

a to moja zabawka

#include <iostream>
#include <bitset>
using namespace std;
 
#define N  2147483647
bitset <N/2+1> t;
// wolę samodzielne rozwiązania,
// trzeba by sprawdzić o ile gorszy byłby robiony na piechotę
// w sumie to kilka linijek
#define M 10000
 
int newton( int n, int k ) { // n>0, k>1, n<32?, binomial(n,k) % M
        int i, r = n;
        for( i=2; i<=k; i++ )
                r = r * --n / i;
        return r % M; }
 
int bpow( int r, long long int x, int n ) { //(r*x^n) % M
        x %= M;   
        while( 1 ) {
                if( n & 1 )
                        r = (r*x) % M; 
                if( n >>= 1 ) 
                        x = (x*x) % M;
                else 
                        return r; }}
 
int mpd(int n, int p){ // p^s | n!
        int s=(n/=p);
        while(n)
                s+=(n/=p);
        return s; }
 
int policz(int n, int k){
        int i, j, w, nk;
        if( 2*k>n ) 
                k = n-k;
        nk=n-k;  // 0 < k <= n-k < n, bo (n,k)==(n, n-k)

        // osobno traktujemy 2 jako jedyną parzystą liczbę pierwszą
        w = bpow( 1, 2, mpd( n, 2 ) - mpd( k, 2 ) - mpd( nk, 2 ));

        for( i=1, j=2*i+1; j <= k;   j += 2, i++ )
                if( !t.test(i) )
                        w = bpow( w, j, mpd( n, j ) - mpd( k, j ) - mpd( nk, j ));

        for(             ; j <= nk;  j += 2, i++ )
                if( !t.test(i) )
                        w = bpow( w, j, mpd( n, j )               - mpd( nk, j ));

        for(             ; j <= n;   j += 2, i++ )
                if( !t.test(i) )
                        w = ((long long)w*j) % M;

        return w; }     
 
int b[1000000][2];
 
int main(){
        int lt, i, n, k, m=0; 
        
        cin >> lt;
        for( i=0; i<lt; i++ ) {
                cin >> b[i][0] >> b[i][1];
                if( b[i][0] > m )
                        m = b[i][0]; }

        for( i=1, k=3; k*k<m; k+=2, i++ ) 
                if( !t.test(i) )
                        for( n=(i*i + i)*2; n+n < m; n += 2*i+1 ) 
                                t.set(n);
        // ! t.test(i) oznacza, że i*2 + 1 jest liczbą pierwszą

        for( i=0; i<lt; i++ )
                if(      b[i][1] == 0 ) cout << 1 << '\n';               // a == 0
                else if( b[i][0] == 0 ) cout << 0 << '\n';               // b == 0
                else if( b[i][1] == b[i][0] ) cout << 1 << '\n';         // a == b
                else if( b[i][0] < 30 ) cout << newton( b[i][0], b[i][1] ) << '\n';
                else cout << policz( b[i][0], b[i][1] ) << '\n'; 
        return 0; } 

newton() liczenie dla małych liczb, można chyba dla n<32, nie pamiętałem nie sprawdziłem, zostało dla <30

bpow() szybkie, od prawej do lewej potęgowanie modularne

mpd() MathWorld formułka (5), maksymalna potęga liczby pierwszej, dzieląca bez reszty silnię liczby naturalnej

policz() liczę z definicji: n!/k!/(n-k)!, rozkładając na czynniki pierwsze wszystkie trzy czynniki, sprawdzam w jakiej potędze występują kolejne liczby pierwsze w tych czynnikach, odejmuję te wykładniki (co odpowiada dzieleniu) i wymnażam z policzonym dla pierwszych mniejszych od bieżącej.
W trzeciej pętelce już wiem, że pierwsze większe od n-k występują w pierwszej potędze.
Wczoraj miałem gynijalnego pomysła by ominąć zakres [sqrt(n) ... n-k] ale jakoś nic to nie pomogło, otrzeźwieję to to przemyślę.

main() zbieram wejście w tablicy, przy okazji znajduję największą liczbę którą będę musiał rozłożyć na czynniki pierwsze
i zapuszczam sito Eratostenesa dla znalezionego maksimum.
Odcedzone pierwsze mam jako bitset liczb nieparzystych, potrzebne sporo pamięci, maks_N / 16 bajtów, 16 bo w bajcie 8 bitów, a że nieparzyste to co druga liczba to i 16.
przed policzeniem sprawdzam czy bieżące współczynniki nie stanowią szczególnego przypadku, jakieś zera czy jedynki, dla małych liczb można to policzyć "na wprost"

nawet mysz powie, że jaj dziecko najładniejsze, a widział ktoś młode myszy? :-)