Kryptografie: Barrett-Reduction



  • Hallo!

    ich beschäftige mich zur Zeit mit Kryptografie (Keygenerierung für DSA z.B.) und brauche hierbei natürlich Algorithmen zur Berechnung vom Modulo für große Zahlen. Hierbei bin ich auf die Barrett-Reduction gestoßen und fand diese recht gut zu implementieren. Leider funktioniert mein Code nicht immer (z.B. rechnet der Code 601 mod 5 = 4). Ich habe leider garkeine Ahnung mehr, woran es liegen könnte. Der Code:

    // returns r := x mod m
    template <std::uint16_t N, typename AtomicType>
    big_int<N, AtomicType> barrett_reduce(const big_int<N, AtomicType>& x, const big_int<N, AtomicType>& m) {
    	if(2 * m.hibit_() >= N * sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits)
    		throw math_error("invalid sizes for barrett reduction");
    	int k = m.hibit_() + 1;
    	const big_int<N, AtomicType>& b2k = big_int<N, AtomicType>(1) << (2 * k); // b2k = b^(2k)
    	const big_int<N, AtomicType>& mu = long_division_(b2k, m).first; // µ = b2k/m
    	big_int<N, AtomicType> q1 = x >> ((k - 1) * sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits), // q1 = x / b^(k-1)
    						   q2 = q1 * mu,
    						   q3 = q2 >> ((k + 1) * sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits), // q3 = q2 / b^(k+1)
    						   bkp1 = big_int<N, AtomicType>(1) << ((k + 1) * sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits), // bkp1 = b^(k+1)
    						   r1 = x - ((x >> (k + 1)) << (k + 1)), // r1 = x mod (k + 1)
    						   mq3 = q3 * m,
    						   r2 = mq3 - ((mq3 >> (k + 1)) << (k + 1)), // r2 = mq3 mod (k + 1)
    						   r = r1 - r2;
    	if(r < big_int<N, AtomicType>(0)) r += bkp1;
    	while(r >= m) r -= m;
    	return r;
    }
    

    Ich habe die big_int-Klasse selber geschrieben, aber jeder einzelne Rechenschritt funktioniert, also glaube ich nicht, dass es an der Klasse liegt. Sie stellt Zahlen im Zweierkomplement zur Basis zwei dar, also ist im Speicher wie ein std::uint32_t-Array.
    Hat jemand von euch eine Idee?

    Edit: Ich habe mich an dieses Paper gehalten: http://130.203.133.150/viewdoc/download;jsessionid=FCBCD05AB106B1251256F7AA8D63431C?doi=10.1.1.156.5934&rep=rep1&type=pdf





  • Hab's hingekriegt!

    Ich hab einfach schlecht gelesen. Die Basis für den Algorithmus ist natürlich nicht 2 sondern die Wortbreite (in meinem Code also sizeof(AtomType) * std::numeric_limits<unsigned char>::max() = 32 meistens). Den Code habe ich dementsprechend angepasst:

    // returns r := x mod m
    	template <std::uint16_t N, typename AtomicType>
    	big_int<N, AtomicType> barrett_reduce(const big_int<N, AtomicType>& x, const big_int<N, AtomicType>& m) {
    		const int atomic_bitsize = sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits;
    		if((2 * m.hibit_() / atomic_bitsize) >= N)
    			throw math_error("invalid sizes for barrett reduction");
    		int k = (m.hibit_() / atomic_bitsize) + 1;
    		const big_int<N, AtomicType>& b2k = big_int<N, AtomicType>(1) << (2 * k * atomic_bitsize); // b2k = b^(2k)
    		const big_int<N, AtomicType>& mu = long_division_(b2k, m).first; // µ = b2k/m
    		big_int<N, AtomicType> q1 = x >> ((k - 1) * atomic_bitsize), // q1 = x / b^(k-1)
    							   q2 = q1 * mu,
    							   q3 = q2 >> ((k + 1) * atomic_bitsize), // q3 = q2 / b^(k+1)
    							   bkp1 = big_int<N, AtomicType>(1) << ((k + 1) * atomic_bitsize), // bkp1 = b^(k+1)
    							   r1 = x - ((x >> ((k + 1) * atomic_bitsize)) << ((k + 1) * atomic_bitsize)), // r1 = x mod (k + 1)
    							   mq3 = q3 * m,
    							   r2 = mq3 - ((mq3 >> ((k + 1) * atomic_bitsize)) << ((k + 1) * atomic_bitsize)), // r2 = mq3 mod (k + 1)
    							   r = r1 - r2;
    		if(r < big_int<N, AtomicType>(0)) r += bkp1;
    		while(r >= m) r -= m;
    		return r;
    	}
    

    Dadurch, dass mein big_int im Zweierkomplement implementiert ist und intern Arrays benutzt, die genau die Wortbreite haben, ist das zu implementieren fast ein Kinderspiel. Für alle Nachfolger, die sich dafür vielleicht interessieren: der Code ist nicht perfekt, es fehlt (noch) eine Option, das µ vorrechnen zu lassen. Bei Zeiten kann ich das vielleicht hier ergänzen.

    @ cookie: Danke, dass du dir Gedanken gemacht hast und tut mir leid, falls ich hier gespoilert habe. 😉



  • Hm... zu früh gefreut. Wenn ich den Algorithmus auf "große Zahlen" modulo "kleine Zahlen loslasse", geht es irgendwie schief, undzwar die Berechnung von r1/r2:

    r1 beinhaltet doch einfach nur die niederwertigsten (k + 1) Stellen von x in der Basis atomic_bitsize (hier 32). Genauso beinhaltet doch r2 die niederwertigsten Stellen von q3 * m. Hierbei berechne ich q3 * m einfach mit einem Karatsuba-Algorithmus. Doch die Differenz von r1 und r2 ist am Ende oft einfach riesig, sodass das Programm dann in der Endschleife Zeile 19 ewig hängenbleibt.

    Merkwürdigerweise passiert das nicht, wenn der "Divisor" m eine Zweierpotenz ist. Woran kann das nur liegen? ...

    EDIT: Es muss ja gelten q3 = x/m. Wenn ich z.B. durch 3 modulo rechnen will, stimmt das allerdings bei mir nicht. Jetzt setzte ich q3 einfach mal per long division auf x / m, dann haut alles hin. Also muss entweder µ oder q1 falsch sein.



  • Ohne mich weiter damit beschäftigt zu haben: der erstbeste Wikipedia Artikel artikuliert gleich in einer der ersten Zeilen die notwendigen Vorbedingungen für den Algorithmus. Dort lautet es: "c = a mod n" unter der Bedingung, dass "a < n^2" ist. Nun ist das ja gerade für besonders große a und besonders kleine n meinem Verständnis nach nicht wirklich gegeben.

    Edit:
    Ich hätte eine Anmerkung zu Deinem Quellcode: Die Funktion ist richtig vollgepflastert mit wiederkehrenden langen template-Typen und sizeofs, die Variablen sind alle einbuchstabig und gehen total darin unter. Als ich die ersten paar Male auf den Code geschaut habe, dachte ich, das wäre kryptischer compile-time-Metaprogramming-Code, aber jetzt sehe ich, dass es eine stinknormale Funktion ist.
    Eventuell würden da einige const size_t's und funktionslokale typedefs helfen, sowie die Neubenennung der Variablen!



  • Die englische Wikipedia ist bisher die einzige Quelle, wo ich das gesehen habe. Ich halte es für falsch. Es muss lediglch gelten, dass bei der Berechnung von x mod m das x mindestens 2 mal so viele Stellen (in der entsprechenden Basis) hat wie die größte, nicht verschwindene Stelle von m. Das wird locker erfüllt.

    Bzgl. des Codes: der zweite hat weniger Redundanzen. Die Benennung der Variablen folgt der Konvention des Algorithmus aus Handbook of Applied Cryptography.



  • Ich habe jetzt schon mehrere Deubgläufe gemacht, aber ich sehe einfach keinen Widerspruch zu dem Pseudocode-Algorithmus.

    Folgender Fall führt z.B. in die Zwickmühle:

    Suche x mod m
    x = (24482, 99, 24482) (Radix-Darstellung zur Basis b = 2^32, also ein 32-Bitblock bzw. jede Zahl ist ein std::uint32_t-Integer)
    m = (0, 0, 3)
    -> k = 1
    -> b^(2k) = (1, 0, 0)
    -> µ = b^(2k)/m = (0, 1431655765, 1431655765)
    -> q1 = (29832, 99, 24482) (Bitshift von x nach rechts mit (k - 1) * 32 = 0
    -> q2 = (4294965512, 2863311497, 2863303370)
    -> q3 = (0, 0, 4294965512)
    -> b^(k+1) = (1, 0, 0)
    -> r1 = (0, 99, 24482)
    -> mq3 = (0, 2, 4294961944)
    -> r2 = (0, 2, 4294961944)
    -> r = r1 - r2 = (0, 96, 29834) ([b]!!!!!!!!!!!!![/b])
    

    Und da diese letzte Zahl zu groß ist, verrennt sich der Algorithmus dann in der while-Schleife so lange, bis er von (0, 96, 29834) (in Dezimal 412.316.890.250) auf den 3er-Interval kommt. Was läuft hier nur schief? 😕

    Edit: Oh ja, und der Code:

    // returns r := x mod m
    template <std::uint16_t N, typename AtomicType>
    big_int<N, AtomicType> barrett_reduce(const big_int<N, AtomicType>& x, const big_int<N, AtomicType>& m, big_int<N, AtomicType> mu) {
    	const int atomic_bitsize = sizeof(AtomicType) * std::numeric_limits<unsigned char>::digits;
    	if((2 * m.hibit_() / atomic_bitsize) >= N)
    		throw math_error("invalid sizes for barrett reduction");
    	const int k = (m.hibit_() / atomic_bitsize) + 1;
    	if(mu < big_int<N, AtomicType>(0)) {
    		const big_int<N, AtomicType> b2k = big_int<N, AtomicType>(1) << (2 * k * atomic_bitsize); // b2k = b^(2k)
    		mu = long_division_(b2k, m).first; // µ = b2k/m
    	}
    	//big_int<N, AtomicType> Q = x / m, R = x % m;
    	int shift = (k + 1) * atomic_bitsize;
    	big_int<N, AtomicType> q1 = x >> ((k - 1) * atomic_bitsize), // q1 = x / b^(k-1)
    						   q2 = q1 * mu,
    						   q3 = q2 >> shift, // q3 = q2 / b^(k+1) =kongruent= floor(x / m),
    						   bkp1 = big_int<N, AtomicType>(1) << shift, // bkp1 = b^(k+1)
    						   r1 = x - ((x >> shift) << shift), // r1 = x mod (k + 1)
    						   mq3 = q3 * m, // mq3 = floor(x / m) * m (remember: x mod m = x - floor(x / m) * m
    						   r2 = mq3 - ((mq3 >> shift) << shift), // r2 = mq3 mod (k + 1)
    						   r = r1 - r2;
    	if(r < big_int<N, AtomicType>(0)) r += bkp1;
    	while(r >= m) r -= m;
    	return r;
    }
    


  • Wie es aussieht, ist die Bedingung nicht, dass der Dividend mindestens, sondern höchstens doppelt so lang ist wie der Divisor. Der Algorithmus funktioniert nämlich für große Dividenden (und beliebige Zweierpotenzen, warum, muss ich erst mal mathematisch ergründen).
    Wie auch immer, er hat scheinbar immer funktioniert.



  • Das ist doch das, was a < n^2 ausdrückt. Logarithmieren zu Deiner Basis ergibt:
    log_b(a) < log_b(n^2) <=> log_b(a) < 2*log_b(n) <=> "Anzahl Ziffern für a" < 2*"Anzahl Ziffern für n"


Log in to reply