view release on metacpan or  search on metacpan

src/rsa/rsa_i31_privexp.c  view on Meta::CPAN

	br_i31_zero(phi, p[0]);
	br_i31_mulacc(phi, p, q);
	len = (phi[0] + 31) >> 5;
	memmove(tmp, phi, (1 + len) * sizeof *phi);
	phi = tmp;
	phi[0] = br_i31_bit_length(phi + 1, len);
	len = (phi[0] + 31) >> 5;

	/*
	 * Divide phi by public exponent e. The final remainder r must be
	 * non-zero (otherwise, the key is invalid). The quotient is k,
	 * which we write over phi, since we don't need phi after that.
	 */
	r = 0;
	for (u = len; u >= 1; u --) {
		/*
		 * Upon entry, r < e, and phi[u] < 2^31; hence,
		 * hi:lo < e*2^31. Thus, the produced word k[u]
		 * must be lower than 2^31, and the new remainder r
		 * is lower than e.
		 */
		uint32_t hi, lo;

		hi = r >> 1;
		lo = (r << 31) + phi[u];
		phi[u] = br_divrem(hi, lo, e, &r);
	}
	if (r == 0) {
		return 0;
	}
	k = phi;

	/*
	 * Compute u and v such that u*e - v*r = GCD(e,r). We use
	 * a binary GCD algorithm, with 6 extra integers a, b,
	 * u0, u1, v0 and v1. Initial values are:
	 *   a = e    u0 = 1   v0 = 0
	 *   b = r    u1 = r   v1 = e-1
	 * The following invariants are maintained:
	 *   a = u0*e - v0*r
	 *   b = u1*e - v1*r
	 *   0 < a <= e
	 *   0 < b <= r
	 *   0 <= u0 <= r
	 *   0 <= v0 <= e
	 *   0 <= u1 <= r
	 *   0 <= v1 <= e
	 *
	 * At each iteration, we reduce either a or b by one bit, and
	 * adjust u0, u1, v0 and v1 to maintain the invariants:
	 *  - if a is even, then a <- a/2
	 *  - otherwise, if b is even, then b <- b/2
	 *  - otherwise, if a > b, then a <- (a-b)/2
	 *  - otherwise, if b > a, then b <- (b-a)/2
	 * Algorithm stops when a = b. At that point, the common value
	 * is the GCD of e and r; it must be 1 (otherwise, the private
	 * key or public exponent is not valid). The (u0,v0) or (u1,v1)
	 * pairs are the solution we are looking for.
	 *
	 * Since either a or b is reduced by at least 1 bit at each
	 * iteration, 62 iterations are enough to reach the end
	 * condition.
	 *
	 * To maintain the invariants, we must compute the same operations
	 * on the u* and v* values that we do on a and b:
	 *  - When a is divided by 2, u0 and v0 must be divided by 2.
	 *  - When b is divided by 2, u1 and v1 must be divided by 2.
	 *  - When b is subtracted from a, u1 and v1 are subtracted from
	 *    u0 and v0, respectively.
	 *  - When a is subtracted from b, u0 and v0 are subtracted from
	 *    u1 and v1, respectively.
	 *
	 * However, we want to keep the u* and v* values in their proper
	 * ranges. The following remarks apply:
	 *
	 *  - When a is divided by 2, then a is even. Therefore:
	 *
	 *     * If r is odd, then u0 and v0 must have the same parity;
	 *       if they are both odd, then adding r to u0 and e to v0
	 *       makes them both even, and the division by 2 brings them
	 *       back to the proper range.
	 *
	 *     * If r is even, then u0 must be even; if v0 is odd, then
	 *       adding r to u0 and e to v0 makes them both even, and the
	 *       division by 2 brings them back to the proper range.
	 *
	 *    Thus, all we need to do is to look at the parity of v0,
	 *    and add (r,e) to (u0,v0) when v0 is odd. In order to avoid
	 *    a 32-bit overflow, we can add ((r+1)/2,(e/2)+1) after the
	 *    division (r+1 does not overflow since r < e; and (e/2)+1
	 *    is equal to (e+1)/2 since e is odd).
	 *
	 *  - When we subtract b from a, three cases may occur:
	 *
	 *     * u1 <= u0 and v1 <= v0: just do the subtractions
	 *
	 *     * u1 > u0 and v1 > v0: compute:
	 *         (u0, v0) <- (u0 + r - u1, v0 + e - v1)
	 *
	 *     * u1 <= u0 and v1 > v0: compute:
	 *         (u0, v0) <- (u0 + r - u1, v0 + e - v1)
	 *
	 *    The fourth case (u1 > u0 and v1 <= v0) is not possible
	 *    because it would contradict "b < a" (which is the reason
	 *    why we subtract b from a).
	 *
	 *    The tricky case is the third one: from the equations, it
	 *    seems that u0 may go out of range. However, the invariants
	 *    and ranges of other values imply that, in that case, the
	 *    new u0 does not actually exceed the range.
	 *
	 *    We can thus handle the subtraction by adding (r,e) based
	 *    solely on the comparison between v0 and v1.
	 */
	a = e;
	b = r;
	u0 = 1;
	v0 = 0;
	u1 = r;
	v1 = e - 1;
	hr = (r + 1) >> 1;