cyassl-lib - This is a port of cyaSSL 2.7.0.

Users » ashleymills » Code » cyassl-lib

This is a port of cyaSSL 2.7.0.

Dependents: CyaSSL_DTLS_Cellular CyaSSL_DTLS_Ethernet

ctaocrypt/src/tfm.c@0:714293de3836, 2013-09-05 (annotated)

Committer:: ashleymills
Date:: Thu Sep 05 10:33:04 2013 +0000
Revision:: 0:714293de3836

Initial commit

Who changed what in which revision?

User	Revision	Line number	New contents of line
ashleymills	0:714293de3836	1	/* tfm.c
ashleymills	0:714293de3836	2	*
ashleymills	0:714293de3836	3	* Copyright (C) 2006-2013 wolfSSL Inc.
ashleymills	0:714293de3836	4	*
ashleymills	0:714293de3836	5	* This file is part of CyaSSL.
ashleymills	0:714293de3836	6	*
ashleymills	0:714293de3836	7	* CyaSSL is free software; you can redistribute it and/or modify
ashleymills	0:714293de3836	8	* it under the terms of the GNU General Public License as published by
ashleymills	0:714293de3836	9	* the Free Software Foundation; either version 2 of the License, or
ashleymills	0:714293de3836	10	* (at your option) any later version.
ashleymills	0:714293de3836	11	*
ashleymills	0:714293de3836	12	* CyaSSL is distributed in the hope that it will be useful,
ashleymills	0:714293de3836	13	* but WITHOUT ANY WARRANTY; without even the implied warranty of
ashleymills	0:714293de3836	14	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
ashleymills	0:714293de3836	15	* GNU General Public License for more details.
ashleymills	0:714293de3836	16	*
ashleymills	0:714293de3836	17	* You should have received a copy of the GNU General Public License
ashleymills	0:714293de3836	18	* along with this program; if not, write to the Free Software
ashleymills	0:714293de3836	19	* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
ashleymills	0:714293de3836	20	*/
ashleymills	0:714293de3836	21
ashleymills	0:714293de3836	22
ashleymills	0:714293de3836	23	/*
ashleymills	0:714293de3836	24	* Based on public domain TomsFastMath 0.10 by Tom St Denis, tomstdenis@iahu.ca,
ashleymills	0:714293de3836	25	* http://math.libtomcrypt.com
ashleymills	0:714293de3836	26	*/
ashleymills	0:714293de3836	27
ashleymills	0:714293de3836	28	/**
ashleymills	0:714293de3836	29	* Edited by Moisés Guimarães (moisesguimaraesm@gmail.com)
ashleymills	0:714293de3836	30	* to fit CyaSSL's needs.
ashleymills	0:714293de3836	31	*/
ashleymills	0:714293de3836	32
ashleymills	0:714293de3836	33	#ifdef HAVE_CONFIG_H
ashleymills	0:714293de3836	34	#include <config.h>
ashleymills	0:714293de3836	35	#endif
ashleymills	0:714293de3836	36
ashleymills	0:714293de3836	37	/* in case user set USE_FAST_MATH there */
ashleymills	0:714293de3836	38	#include <cyassl/ctaocrypt/settings.h>
ashleymills	0:714293de3836	39
ashleymills	0:714293de3836	40	#ifdef USE_FAST_MATH
ashleymills	0:714293de3836	41
ashleymills	0:714293de3836	42	#include <cyassl/ctaocrypt/tfm.h>
ashleymills	0:714293de3836	43	#include <ctaocrypt/src/asm.c> /* will define asm MACROS or C ones */
ashleymills	0:714293de3836	44
ashleymills	0:714293de3836	45
ashleymills	0:714293de3836	46	/* math settings check */
ashleymills	0:714293de3836	47	word32 CheckRunTimeSettings(void)
ashleymills	0:714293de3836	48	{
ashleymills	0:714293de3836	49	return CTC_SETTINGS;
ashleymills	0:714293de3836	50	}
ashleymills	0:714293de3836	51
ashleymills	0:714293de3836	52
ashleymills	0:714293de3836	53	/* math settings size check */
ashleymills	0:714293de3836	54	word32 CheckRunTimeFastMath(void)
ashleymills	0:714293de3836	55	{
ashleymills	0:714293de3836	56	return FP_SIZE;
ashleymills	0:714293de3836	57	}
ashleymills	0:714293de3836	58
ashleymills	0:714293de3836	59
ashleymills	0:714293de3836	60	/* Functions */
ashleymills	0:714293de3836	61
ashleymills	0:714293de3836	62	void fp_add(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	63	{
ashleymills	0:714293de3836	64	int sa, sb;
ashleymills	0:714293de3836	65
ashleymills	0:714293de3836	66	/* get sign of both inputs */
ashleymills	0:714293de3836	67	sa = a->sign;
ashleymills	0:714293de3836	68	sb = b->sign;
ashleymills	0:714293de3836	69
ashleymills	0:714293de3836	70	/* handle two cases, not four */
ashleymills	0:714293de3836	71	if (sa == sb) {
ashleymills	0:714293de3836	72	/* both positive or both negative */
ashleymills	0:714293de3836	73	/* add their magnitudes, copy the sign */
ashleymills	0:714293de3836	74	c->sign = sa;
ashleymills	0:714293de3836	75	s_fp_add (a, b, c);
ashleymills	0:714293de3836	76	} else {
ashleymills	0:714293de3836	77	/* one positive, the other negative */
ashleymills	0:714293de3836	78	/* subtract the one with the greater magnitude from */
ashleymills	0:714293de3836	79	/* the one of the lesser magnitude. The result gets */
ashleymills	0:714293de3836	80	/* the sign of the one with the greater magnitude. */
ashleymills	0:714293de3836	81	if (fp_cmp_mag (a, b) == FP_LT) {
ashleymills	0:714293de3836	82	c->sign = sb;
ashleymills	0:714293de3836	83	s_fp_sub (b, a, c);
ashleymills	0:714293de3836	84	} else {
ashleymills	0:714293de3836	85	c->sign = sa;
ashleymills	0:714293de3836	86	s_fp_sub (a, b, c);
ashleymills	0:714293de3836	87	}
ashleymills	0:714293de3836	88	}
ashleymills	0:714293de3836	89	}
ashleymills	0:714293de3836	90
ashleymills	0:714293de3836	91	/* unsigned addition */
ashleymills	0:714293de3836	92	void s_fp_add(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	93	{
ashleymills	0:714293de3836	94	int x, y, oldused;
ashleymills	0:714293de3836	95	register fp_word t;
ashleymills	0:714293de3836	96
ashleymills	0:714293de3836	97	y = MAX(a->used, b->used);
ashleymills	0:714293de3836	98	oldused = c->used;
ashleymills	0:714293de3836	99	c->used = y;
ashleymills	0:714293de3836	100
ashleymills	0:714293de3836	101	t = 0;
ashleymills	0:714293de3836	102	for (x = 0; x < y; x++) {
ashleymills	0:714293de3836	103	t += ((fp_word)a->dp[x]) + ((fp_word)b->dp[x]);
ashleymills	0:714293de3836	104	c->dp[x] = (fp_digit)t;
ashleymills	0:714293de3836	105	t >>= DIGIT_BIT;
ashleymills	0:714293de3836	106	}
ashleymills	0:714293de3836	107	if (t != 0 && x < FP_SIZE) {
ashleymills	0:714293de3836	108	c->dp[c->used++] = (fp_digit)t;
ashleymills	0:714293de3836	109	++x;
ashleymills	0:714293de3836	110	}
ashleymills	0:714293de3836	111
ashleymills	0:714293de3836	112	c->used = x;
ashleymills	0:714293de3836	113	for (; x < oldused; x++) {
ashleymills	0:714293de3836	114	c->dp[x] = 0;
ashleymills	0:714293de3836	115	}
ashleymills	0:714293de3836	116	fp_clamp(c);
ashleymills	0:714293de3836	117	}
ashleymills	0:714293de3836	118
ashleymills	0:714293de3836	119	/* c = a - b */
ashleymills	0:714293de3836	120	void fp_sub(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	121	{
ashleymills	0:714293de3836	122	int sa, sb;
ashleymills	0:714293de3836	123
ashleymills	0:714293de3836	124	sa = a->sign;
ashleymills	0:714293de3836	125	sb = b->sign;
ashleymills	0:714293de3836	126
ashleymills	0:714293de3836	127	if (sa != sb) {
ashleymills	0:714293de3836	128	/* subtract a negative from a positive, OR */
ashleymills	0:714293de3836	129	/* subtract a positive from a negative. */
ashleymills	0:714293de3836	130	/* In either case, ADD their magnitudes, */
ashleymills	0:714293de3836	131	/* and use the sign of the first number. */
ashleymills	0:714293de3836	132	c->sign = sa;
ashleymills	0:714293de3836	133	s_fp_add (a, b, c);
ashleymills	0:714293de3836	134	} else {
ashleymills	0:714293de3836	135	/* subtract a positive from a positive, OR */
ashleymills	0:714293de3836	136	/* subtract a negative from a negative. */
ashleymills	0:714293de3836	137	/* First, take the difference between their */
ashleymills	0:714293de3836	138	/* magnitudes, then... */
ashleymills	0:714293de3836	139	if (fp_cmp_mag (a, b) != FP_LT) {
ashleymills	0:714293de3836	140	/* Copy the sign from the first */
ashleymills	0:714293de3836	141	c->sign = sa;
ashleymills	0:714293de3836	142	/* The first has a larger or equal magnitude */
ashleymills	0:714293de3836	143	s_fp_sub (a, b, c);
ashleymills	0:714293de3836	144	} else {
ashleymills	0:714293de3836	145	/* The result has the opposite sign from */
ashleymills	0:714293de3836	146	/* the first number. */
ashleymills	0:714293de3836	147	c->sign = (sa == FP_ZPOS) ? FP_NEG : FP_ZPOS;
ashleymills	0:714293de3836	148	/* The second has a larger magnitude */
ashleymills	0:714293de3836	149	s_fp_sub (b, a, c);
ashleymills	0:714293de3836	150	}
ashleymills	0:714293de3836	151	}
ashleymills	0:714293de3836	152	}
ashleymills	0:714293de3836	153
ashleymills	0:714293de3836	154	/* unsigned subtraction \|\|a\|\| >= \|\|b\|\| ALWAYS! */
ashleymills	0:714293de3836	155	void s_fp_sub(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	156	{
ashleymills	0:714293de3836	157	int x, oldbused, oldused;
ashleymills	0:714293de3836	158	fp_word t;
ashleymills	0:714293de3836	159
ashleymills	0:714293de3836	160	oldused = c->used;
ashleymills	0:714293de3836	161	oldbused = b->used;
ashleymills	0:714293de3836	162	c->used = a->used;
ashleymills	0:714293de3836	163	t = 0;
ashleymills	0:714293de3836	164	for (x = 0; x < oldbused; x++) {
ashleymills	0:714293de3836	165	t = ((fp_word)a->dp[x]) - (((fp_word)b->dp[x]) + t);
ashleymills	0:714293de3836	166	c->dp[x] = (fp_digit)t;
ashleymills	0:714293de3836	167	t = (t >> DIGIT_BIT)&1;
ashleymills	0:714293de3836	168	}
ashleymills	0:714293de3836	169	for (; x < a->used; x++) {
ashleymills	0:714293de3836	170	t = ((fp_word)a->dp[x]) - t;
ashleymills	0:714293de3836	171	c->dp[x] = (fp_digit)t;
ashleymills	0:714293de3836	172	t = (t >> DIGIT_BIT)&1;
ashleymills	0:714293de3836	173	}
ashleymills	0:714293de3836	174	for (; x < oldused; x++) {
ashleymills	0:714293de3836	175	c->dp[x] = 0;
ashleymills	0:714293de3836	176	}
ashleymills	0:714293de3836	177	fp_clamp(c);
ashleymills	0:714293de3836	178	}
ashleymills	0:714293de3836	179
ashleymills	0:714293de3836	180	/* c = a * b */
ashleymills	0:714293de3836	181	void fp_mul(fp_int A, fp_int B, fp_int *C)
ashleymills	0:714293de3836	182	{
ashleymills	0:714293de3836	183	int y, yy;
ashleymills	0:714293de3836	184
ashleymills	0:714293de3836	185	y = MAX(A->used, B->used);
ashleymills	0:714293de3836	186	yy = MIN(A->used, B->used);
ashleymills	0:714293de3836	187
ashleymills	0:714293de3836	188	/* call generic if we're out of range */
ashleymills	0:714293de3836	189	if (y + yy > FP_SIZE) {
ashleymills	0:714293de3836	190	fp_mul_comba(A, B, C);
ashleymills	0:714293de3836	191	return ;
ashleymills	0:714293de3836	192	}
ashleymills	0:714293de3836	193
ashleymills	0:714293de3836	194	/* pick a comba (unrolled 4/8/16/32 x or rolled) based on the size
ashleymills	0:714293de3836	195	of the largest input. We also want to avoid doing excess mults if the
ashleymills	0:714293de3836	196	inputs are not close to the next power of two. That is, for example,
ashleymills	0:714293de3836	197	if say y=17 then we would do (32-17)^2 = 225 unneeded multiplications
ashleymills	0:714293de3836	198	*/
ashleymills	0:714293de3836	199
ashleymills	0:714293de3836	200	#ifdef TFM_MUL3
ashleymills	0:714293de3836	201	if (y <= 3) {
ashleymills	0:714293de3836	202	fp_mul_comba3(A,B,C);
ashleymills	0:714293de3836	203	return;
ashleymills	0:714293de3836	204	}
ashleymills	0:714293de3836	205	#endif
ashleymills	0:714293de3836	206	#ifdef TFM_MUL4
ashleymills	0:714293de3836	207	if (y == 4) {
ashleymills	0:714293de3836	208	fp_mul_comba4(A,B,C);
ashleymills	0:714293de3836	209	return;
ashleymills	0:714293de3836	210	}
ashleymills	0:714293de3836	211	#endif
ashleymills	0:714293de3836	212	#ifdef TFM_MUL6
ashleymills	0:714293de3836	213	if (y <= 6) {
ashleymills	0:714293de3836	214	fp_mul_comba6(A,B,C);
ashleymills	0:714293de3836	215	return;
ashleymills	0:714293de3836	216	}
ashleymills	0:714293de3836	217	#endif
ashleymills	0:714293de3836	218	#ifdef TFM_MUL7
ashleymills	0:714293de3836	219	if (y == 7) {
ashleymills	0:714293de3836	220	fp_mul_comba7(A,B,C);
ashleymills	0:714293de3836	221	return;
ashleymills	0:714293de3836	222	}
ashleymills	0:714293de3836	223	#endif
ashleymills	0:714293de3836	224	#ifdef TFM_MUL8
ashleymills	0:714293de3836	225	if (y == 8) {
ashleymills	0:714293de3836	226	fp_mul_comba8(A,B,C);
ashleymills	0:714293de3836	227	return;
ashleymills	0:714293de3836	228	}
ashleymills	0:714293de3836	229	#endif
ashleymills	0:714293de3836	230	#ifdef TFM_MUL9
ashleymills	0:714293de3836	231	if (y == 9) {
ashleymills	0:714293de3836	232	fp_mul_comba9(A,B,C);
ashleymills	0:714293de3836	233	return;
ashleymills	0:714293de3836	234	}
ashleymills	0:714293de3836	235	#endif
ashleymills	0:714293de3836	236	#ifdef TFM_MUL12
ashleymills	0:714293de3836	237	if (y <= 12) {
ashleymills	0:714293de3836	238	fp_mul_comba12(A,B,C);
ashleymills	0:714293de3836	239	return;
ashleymills	0:714293de3836	240	}
ashleymills	0:714293de3836	241	#endif
ashleymills	0:714293de3836	242	#ifdef TFM_MUL17
ashleymills	0:714293de3836	243	if (y <= 17) {
ashleymills	0:714293de3836	244	fp_mul_comba17(A,B,C);
ashleymills	0:714293de3836	245	return;
ashleymills	0:714293de3836	246	}
ashleymills	0:714293de3836	247	#endif
ashleymills	0:714293de3836	248
ashleymills	0:714293de3836	249	#ifdef TFM_SMALL_SET
ashleymills	0:714293de3836	250	if (y <= 16) {
ashleymills	0:714293de3836	251	fp_mul_comba_small(A,B,C);
ashleymills	0:714293de3836	252	return;
ashleymills	0:714293de3836	253	}
ashleymills	0:714293de3836	254	#endif
ashleymills	0:714293de3836	255	#if defined(TFM_MUL20)
ashleymills	0:714293de3836	256	if (y <= 20) {
ashleymills	0:714293de3836	257	fp_mul_comba20(A,B,C);
ashleymills	0:714293de3836	258	return;
ashleymills	0:714293de3836	259	}
ashleymills	0:714293de3836	260	#endif
ashleymills	0:714293de3836	261	#if defined(TFM_MUL24)
ashleymills	0:714293de3836	262	if (yy >= 16 && y <= 24) {
ashleymills	0:714293de3836	263	fp_mul_comba24(A,B,C);
ashleymills	0:714293de3836	264	return;
ashleymills	0:714293de3836	265	}
ashleymills	0:714293de3836	266	#endif
ashleymills	0:714293de3836	267	#if defined(TFM_MUL28)
ashleymills	0:714293de3836	268	if (yy >= 20 && y <= 28) {
ashleymills	0:714293de3836	269	fp_mul_comba28(A,B,C);
ashleymills	0:714293de3836	270	return;
ashleymills	0:714293de3836	271	}
ashleymills	0:714293de3836	272	#endif
ashleymills	0:714293de3836	273	#if defined(TFM_MUL32)
ashleymills	0:714293de3836	274	if (yy >= 24 && y <= 32) {
ashleymills	0:714293de3836	275	fp_mul_comba32(A,B,C);
ashleymills	0:714293de3836	276	return;
ashleymills	0:714293de3836	277	}
ashleymills	0:714293de3836	278	#endif
ashleymills	0:714293de3836	279	#if defined(TFM_MUL48)
ashleymills	0:714293de3836	280	if (yy >= 40 && y <= 48) {
ashleymills	0:714293de3836	281	fp_mul_comba48(A,B,C);
ashleymills	0:714293de3836	282	return;
ashleymills	0:714293de3836	283	}
ashleymills	0:714293de3836	284	#endif
ashleymills	0:714293de3836	285	#if defined(TFM_MUL64)
ashleymills	0:714293de3836	286	if (yy >= 56 && y <= 64) {
ashleymills	0:714293de3836	287	fp_mul_comba64(A,B,C);
ashleymills	0:714293de3836	288	return;
ashleymills	0:714293de3836	289	}
ashleymills	0:714293de3836	290	#endif
ashleymills	0:714293de3836	291	fp_mul_comba(A,B,C);
ashleymills	0:714293de3836	292	}
ashleymills	0:714293de3836	293
ashleymills	0:714293de3836	294	void fp_mul_2(fp_int * a, fp_int * b)
ashleymills	0:714293de3836	295	{
ashleymills	0:714293de3836	296	int x, oldused;
ashleymills	0:714293de3836	297
ashleymills	0:714293de3836	298	oldused = b->used;
ashleymills	0:714293de3836	299	b->used = a->used;
ashleymills	0:714293de3836	300
ashleymills	0:714293de3836	301	{
ashleymills	0:714293de3836	302	register fp_digit r, rr, tmpa, tmpb;
ashleymills	0:714293de3836	303
ashleymills	0:714293de3836	304	/* alias for source */
ashleymills	0:714293de3836	305	tmpa = a->dp;
ashleymills	0:714293de3836	306
ashleymills	0:714293de3836	307	/* alias for dest */
ashleymills	0:714293de3836	308	tmpb = b->dp;
ashleymills	0:714293de3836	309
ashleymills	0:714293de3836	310	/* carry */
ashleymills	0:714293de3836	311	r = 0;
ashleymills	0:714293de3836	312	for (x = 0; x < a->used; x++) {
ashleymills	0:714293de3836	313
ashleymills	0:714293de3836	314	/* get what will be the next carry bit from the
ashleymills	0:714293de3836	315	* MSB of the current digit
ashleymills	0:714293de3836	316	*/
ashleymills	0:714293de3836	317	rr = *tmpa >> ((fp_digit)(DIGIT_BIT - 1));
ashleymills	0:714293de3836	318
ashleymills	0:714293de3836	319	/* now shift up this digit, add in the carry [from the previous] */
ashleymills	0:714293de3836	320	tmpb++ = ((tmpa++ << ((fp_digit)1)) \| r);
ashleymills	0:714293de3836	321
ashleymills	0:714293de3836	322	/* copy the carry that would be from the source
ashleymills	0:714293de3836	323	* digit into the next iteration
ashleymills	0:714293de3836	324	*/
ashleymills	0:714293de3836	325	r = rr;
ashleymills	0:714293de3836	326	}
ashleymills	0:714293de3836	327
ashleymills	0:714293de3836	328	/* new leading digit? */
ashleymills	0:714293de3836	329	if (r != 0 && b->used != (FP_SIZE-1)) {
ashleymills	0:714293de3836	330	/* add a MSB which is always 1 at this point */
ashleymills	0:714293de3836	331	*tmpb = 1;
ashleymills	0:714293de3836	332	++(b->used);
ashleymills	0:714293de3836	333	}
ashleymills	0:714293de3836	334
ashleymills	0:714293de3836	335	/* now zero any excess digits on the destination
ashleymills	0:714293de3836	336	* that we didn't write to
ashleymills	0:714293de3836	337	*/
ashleymills	0:714293de3836	338	tmpb = b->dp + b->used;
ashleymills	0:714293de3836	339	for (x = b->used; x < oldused; x++) {
ashleymills	0:714293de3836	340	*tmpb++ = 0;
ashleymills	0:714293de3836	341	}
ashleymills	0:714293de3836	342	}
ashleymills	0:714293de3836	343	b->sign = a->sign;
ashleymills	0:714293de3836	344	}
ashleymills	0:714293de3836	345
ashleymills	0:714293de3836	346	/* c = a * b */
ashleymills	0:714293de3836	347	void fp_mul_d(fp_int a, fp_digit b, fp_int c)
ashleymills	0:714293de3836	348	{
ashleymills	0:714293de3836	349	fp_word w;
ashleymills	0:714293de3836	350	int x, oldused;
ashleymills	0:714293de3836	351
ashleymills	0:714293de3836	352	oldused = c->used;
ashleymills	0:714293de3836	353	c->used = a->used;
ashleymills	0:714293de3836	354	c->sign = a->sign;
ashleymills	0:714293de3836	355	w = 0;
ashleymills	0:714293de3836	356	for (x = 0; x < a->used; x++) {
ashleymills	0:714293de3836	357	w = ((fp_word)a->dp[x]) * ((fp_word)b) + w;
ashleymills	0:714293de3836	358	c->dp[x] = (fp_digit)w;
ashleymills	0:714293de3836	359	w = w >> DIGIT_BIT;
ashleymills	0:714293de3836	360	}
ashleymills	0:714293de3836	361	if (w != 0 && (a->used != FP_SIZE)) {
ashleymills	0:714293de3836	362	c->dp[c->used++] = (fp_digit) w;
ashleymills	0:714293de3836	363	++x;
ashleymills	0:714293de3836	364	}
ashleymills	0:714293de3836	365	for (; x < oldused; x++) {
ashleymills	0:714293de3836	366	c->dp[x] = 0;
ashleymills	0:714293de3836	367	}
ashleymills	0:714293de3836	368	fp_clamp(c);
ashleymills	0:714293de3836	369	}
ashleymills	0:714293de3836	370
ashleymills	0:714293de3836	371	/* c = a * 2*d /
ashleymills	0:714293de3836	372	void fp_mul_2d(fp_int a, int b, fp_int c)
ashleymills	0:714293de3836	373	{
ashleymills	0:714293de3836	374	fp_digit carry, carrytmp, shift;
ashleymills	0:714293de3836	375	int x;
ashleymills	0:714293de3836	376
ashleymills	0:714293de3836	377	/* copy it */
ashleymills	0:714293de3836	378	fp_copy(a, c);
ashleymills	0:714293de3836	379
ashleymills	0:714293de3836	380	/* handle whole digits */
ashleymills	0:714293de3836	381	if (b >= DIGIT_BIT) {
ashleymills	0:714293de3836	382	fp_lshd(c, b/DIGIT_BIT);
ashleymills	0:714293de3836	383	}
ashleymills	0:714293de3836	384	b %= DIGIT_BIT;
ashleymills	0:714293de3836	385
ashleymills	0:714293de3836	386	/* shift the digits */
ashleymills	0:714293de3836	387	if (b != 0) {
ashleymills	0:714293de3836	388	carry = 0;
ashleymills	0:714293de3836	389	shift = DIGIT_BIT - b;
ashleymills	0:714293de3836	390	for (x = 0; x < c->used; x++) {
ashleymills	0:714293de3836	391	carrytmp = c->dp[x] >> shift;
ashleymills	0:714293de3836	392	c->dp[x] = (c->dp[x] << b) + carry;
ashleymills	0:714293de3836	393	carry = carrytmp;
ashleymills	0:714293de3836	394	}
ashleymills	0:714293de3836	395	/* store last carry if room */
ashleymills	0:714293de3836	396	if (carry && x < FP_SIZE) {
ashleymills	0:714293de3836	397	c->dp[c->used++] = carry;
ashleymills	0:714293de3836	398	}
ashleymills	0:714293de3836	399	}
ashleymills	0:714293de3836	400	fp_clamp(c);
ashleymills	0:714293de3836	401	}
ashleymills	0:714293de3836	402
ashleymills	0:714293de3836	403	/* generic PxQ multiplier */
ashleymills	0:714293de3836	404	void fp_mul_comba(fp_int A, fp_int B, fp_int *C)
ashleymills	0:714293de3836	405	{
ashleymills	0:714293de3836	406	int ix, iy, iz, tx, ty, pa;
ashleymills	0:714293de3836	407	fp_digit c0, c1, c2, tmpx, tmpy;
ashleymills	0:714293de3836	408	fp_int tmp, *dst;
ashleymills	0:714293de3836	409
ashleymills	0:714293de3836	410	COMBA_START;
ashleymills	0:714293de3836	411	COMBA_CLEAR;
ashleymills	0:714293de3836	412
ashleymills	0:714293de3836	413	/* get size of output and trim */
ashleymills	0:714293de3836	414	pa = A->used + B->used;
ashleymills	0:714293de3836	415	if (pa >= FP_SIZE) {
ashleymills	0:714293de3836	416	pa = FP_SIZE-1;
ashleymills	0:714293de3836	417	}
ashleymills	0:714293de3836	418
ashleymills	0:714293de3836	419	if (A == C \|\| B == C) {
ashleymills	0:714293de3836	420	fp_zero(&tmp);
ashleymills	0:714293de3836	421	dst = &tmp;
ashleymills	0:714293de3836	422	} else {
ashleymills	0:714293de3836	423	fp_zero(C);
ashleymills	0:714293de3836	424	dst = C;
ashleymills	0:714293de3836	425	}
ashleymills	0:714293de3836	426
ashleymills	0:714293de3836	427	for (ix = 0; ix < pa; ix++) {
ashleymills	0:714293de3836	428	/* get offsets into the two bignums */
ashleymills	0:714293de3836	429	ty = MIN(ix, B->used-1);
ashleymills	0:714293de3836	430	tx = ix - ty;
ashleymills	0:714293de3836	431
ashleymills	0:714293de3836	432	/* setup temp aliases */
ashleymills	0:714293de3836	433	tmpx = A->dp + tx;
ashleymills	0:714293de3836	434	tmpy = B->dp + ty;
ashleymills	0:714293de3836	435
ashleymills	0:714293de3836	436	/* this is the number of times the loop will iterrate, essentially its
ashleymills	0:714293de3836	437	while (tx++ < a->used && ty-- >= 0) { ... }
ashleymills	0:714293de3836	438	*/
ashleymills	0:714293de3836	439	iy = MIN(A->used-tx, ty+1);
ashleymills	0:714293de3836	440
ashleymills	0:714293de3836	441	/* execute loop */
ashleymills	0:714293de3836	442	COMBA_FORWARD;
ashleymills	0:714293de3836	443	for (iz = 0; iz < iy; ++iz) {
ashleymills	0:714293de3836	444	/* TAO change COMBA_ADD back to MULADD */
ashleymills	0:714293de3836	445	MULADD(tmpx++, tmpy--);
ashleymills	0:714293de3836	446	}
ashleymills	0:714293de3836	447
ashleymills	0:714293de3836	448	/* store term */
ashleymills	0:714293de3836	449	COMBA_STORE(dst->dp[ix]);
ashleymills	0:714293de3836	450	}
ashleymills	0:714293de3836	451	COMBA_FINI;
ashleymills	0:714293de3836	452
ashleymills	0:714293de3836	453	dst->used = pa;
ashleymills	0:714293de3836	454	dst->sign = A->sign ^ B->sign;
ashleymills	0:714293de3836	455	fp_clamp(dst);
ashleymills	0:714293de3836	456	fp_copy(dst, C);
ashleymills	0:714293de3836	457	}
ashleymills	0:714293de3836	458
ashleymills	0:714293de3836	459	/* a/b => cb + d == a */
ashleymills	0:714293de3836	460	int fp_div(fp_int a, fp_int b, fp_int c, fp_int d)
ashleymills	0:714293de3836	461	{
ashleymills	0:714293de3836	462	fp_int q, x, y, t1, t2;
ashleymills	0:714293de3836	463	int n, t, i, norm, neg;
ashleymills	0:714293de3836	464
ashleymills	0:714293de3836	465	/* is divisor zero ? */
ashleymills	0:714293de3836	466	if (fp_iszero (b) == 1) {
ashleymills	0:714293de3836	467	return FP_VAL;
ashleymills	0:714293de3836	468	}
ashleymills	0:714293de3836	469
ashleymills	0:714293de3836	470	/* if a < b then q=0, r = a */
ashleymills	0:714293de3836	471	if (fp_cmp_mag (a, b) == FP_LT) {
ashleymills	0:714293de3836	472	if (d != NULL) {
ashleymills	0:714293de3836	473	fp_copy (a, d);
ashleymills	0:714293de3836	474	}
ashleymills	0:714293de3836	475	if (c != NULL) {
ashleymills	0:714293de3836	476	fp_zero (c);
ashleymills	0:714293de3836	477	}
ashleymills	0:714293de3836	478	return FP_OKAY;
ashleymills	0:714293de3836	479	}
ashleymills	0:714293de3836	480
ashleymills	0:714293de3836	481	fp_init(&q);
ashleymills	0:714293de3836	482	q.used = a->used + 2;
ashleymills	0:714293de3836	483
ashleymills	0:714293de3836	484	fp_init(&t1);
ashleymills	0:714293de3836	485	fp_init(&t2);
ashleymills	0:714293de3836	486	fp_init_copy(&x, a);
ashleymills	0:714293de3836	487	fp_init_copy(&y, b);
ashleymills	0:714293de3836	488
ashleymills	0:714293de3836	489	/* fix the sign */
ashleymills	0:714293de3836	490	neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
ashleymills	0:714293de3836	491	x.sign = y.sign = FP_ZPOS;
ashleymills	0:714293de3836	492
ashleymills	0:714293de3836	493	/* normalize both x and y, ensure that y >= b/2, [b == 2*DIGIT_BIT] /
ashleymills	0:714293de3836	494	norm = fp_count_bits(&y) % DIGIT_BIT;
ashleymills	0:714293de3836	495	if (norm < (int)(DIGIT_BIT-1)) {
ashleymills	0:714293de3836	496	norm = (DIGIT_BIT-1) - norm;
ashleymills	0:714293de3836	497	fp_mul_2d (&x, norm, &x);
ashleymills	0:714293de3836	498	fp_mul_2d (&y, norm, &y);
ashleymills	0:714293de3836	499	} else {
ashleymills	0:714293de3836	500	norm = 0;
ashleymills	0:714293de3836	501	}
ashleymills	0:714293de3836	502
ashleymills	0:714293de3836	503	/* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
ashleymills	0:714293de3836	504	n = x.used - 1;
ashleymills	0:714293de3836	505	t = y.used - 1;
ashleymills	0:714293de3836	506
ashleymills	0:714293de3836	507	/* while (x >= ybn-t) do { q[n-t] += 1; x -= yb*{n-t} } /
ashleymills	0:714293de3836	508	fp_lshd (&y, n - t); /* y = yb{n-t} /
ashleymills	0:714293de3836	509
ashleymills	0:714293de3836	510	while (fp_cmp (&x, &y) != FP_LT) {
ashleymills	0:714293de3836	511	++(q.dp[n - t]);
ashleymills	0:714293de3836	512	fp_sub (&x, &y, &x);
ashleymills	0:714293de3836	513	}
ashleymills	0:714293de3836	514
ashleymills	0:714293de3836	515	/* reset y by shifting it back down */
ashleymills	0:714293de3836	516	fp_rshd (&y, n - t);
ashleymills	0:714293de3836	517
ashleymills	0:714293de3836	518	/* step 3. for i from n down to (t + 1) */
ashleymills	0:714293de3836	519	for (i = n; i >= (t + 1); i--) {
ashleymills	0:714293de3836	520	if (i > x.used) {
ashleymills	0:714293de3836	521	continue;
ashleymills	0:714293de3836	522	}
ashleymills	0:714293de3836	523
ashleymills	0:714293de3836	524	/* step 3.1 if xi == yt then set q{i-t-1} to b-1,
ashleymills	0:714293de3836	525	* otherwise set q{i-t-1} to (xib + x{i-1})/yt /
ashleymills	0:714293de3836	526	if (x.dp[i] == y.dp[t]) {
ashleymills	0:714293de3836	527	q.dp[i - t - 1] = (fp_digit) ((((fp_word)1) << DIGIT_BIT) - 1);
ashleymills	0:714293de3836	528	} else {
ashleymills	0:714293de3836	529	fp_word tmp;
ashleymills	0:714293de3836	530	tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
ashleymills	0:714293de3836	531	tmp \|= ((fp_word) x.dp[i - 1]);
ashleymills	0:714293de3836	532	tmp /= ((fp_word)y.dp[t]);
ashleymills	0:714293de3836	533	q.dp[i - t - 1] = (fp_digit) (tmp);
ashleymills	0:714293de3836	534	}
ashleymills	0:714293de3836	535
ashleymills	0:714293de3836	536	/* while (q{i-t-1} * (yt * b + y{t-1})) >
ashleymills	0:714293de3836	537	xi * b*2 + xi-1 b + xi-2
ashleymills	0:714293de3836	538
ashleymills	0:714293de3836	539	do q{i-t-1} -= 1;
ashleymills	0:714293de3836	540	*/
ashleymills	0:714293de3836	541	q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
ashleymills	0:714293de3836	542	do {
ashleymills	0:714293de3836	543	q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
ashleymills	0:714293de3836	544
ashleymills	0:714293de3836	545	/* find left hand */
ashleymills	0:714293de3836	546	fp_zero (&t1);
ashleymills	0:714293de3836	547	t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
ashleymills	0:714293de3836	548	t1.dp[1] = y.dp[t];
ashleymills	0:714293de3836	549	t1.used = 2;
ashleymills	0:714293de3836	550	fp_mul_d (&t1, q.dp[i - t - 1], &t1);
ashleymills	0:714293de3836	551
ashleymills	0:714293de3836	552	/* find right hand */
ashleymills	0:714293de3836	553	t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
ashleymills	0:714293de3836	554	t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
ashleymills	0:714293de3836	555	t2.dp[2] = x.dp[i];
ashleymills	0:714293de3836	556	t2.used = 3;
ashleymills	0:714293de3836	557	} while (fp_cmp_mag(&t1, &t2) == FP_GT);
ashleymills	0:714293de3836	558
ashleymills	0:714293de3836	559	/* step 3.3 x = x - q{i-t-1} * y * b*{i-t-1} /
ashleymills	0:714293de3836	560	fp_mul_d (&y, q.dp[i - t - 1], &t1);
ashleymills	0:714293de3836	561	fp_lshd (&t1, i - t - 1);
ashleymills	0:714293de3836	562	fp_sub (&x, &t1, &x);
ashleymills	0:714293de3836	563
ashleymills	0:714293de3836	564	/* if x < 0 then { x = x + yb{i-t-1}; q{i-t-1} -= 1; } /
ashleymills	0:714293de3836	565	if (x.sign == FP_NEG) {
ashleymills	0:714293de3836	566	fp_copy (&y, &t1);
ashleymills	0:714293de3836	567	fp_lshd (&t1, i - t - 1);
ashleymills	0:714293de3836	568	fp_add (&x, &t1, &x);
ashleymills	0:714293de3836	569	q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
ashleymills	0:714293de3836	570	}
ashleymills	0:714293de3836	571	}
ashleymills	0:714293de3836	572
ashleymills	0:714293de3836	573	/* now q is the quotient and x is the remainder
ashleymills	0:714293de3836	574	* [which we have to normalize]
ashleymills	0:714293de3836	575	*/
ashleymills	0:714293de3836	576
ashleymills	0:714293de3836	577	/* get sign before writing to c */
ashleymills	0:714293de3836	578	x.sign = x.used == 0 ? FP_ZPOS : a->sign;
ashleymills	0:714293de3836	579
ashleymills	0:714293de3836	580	if (c != NULL) {
ashleymills	0:714293de3836	581	fp_clamp (&q);
ashleymills	0:714293de3836	582	fp_copy (&q, c);
ashleymills	0:714293de3836	583	c->sign = neg;
ashleymills	0:714293de3836	584	}
ashleymills	0:714293de3836	585
ashleymills	0:714293de3836	586	if (d != NULL) {
ashleymills	0:714293de3836	587	fp_div_2d (&x, norm, &x, NULL);
ashleymills	0:714293de3836	588
ashleymills	0:714293de3836	589	/* the following is a kludge, essentially we were seeing the right remainder but
ashleymills	0:714293de3836	590	with excess digits that should have been zero
ashleymills	0:714293de3836	591	*/
ashleymills	0:714293de3836	592	for (i = b->used; i < x.used; i++) {
ashleymills	0:714293de3836	593	x.dp[i] = 0;
ashleymills	0:714293de3836	594	}
ashleymills	0:714293de3836	595	fp_clamp(&x);
ashleymills	0:714293de3836	596	fp_copy (&x, d);
ashleymills	0:714293de3836	597	}
ashleymills	0:714293de3836	598
ashleymills	0:714293de3836	599	return FP_OKAY;
ashleymills	0:714293de3836	600	}
ashleymills	0:714293de3836	601
ashleymills	0:714293de3836	602	/* b = a/2 */
ashleymills	0:714293de3836	603	void fp_div_2(fp_int * a, fp_int * b)
ashleymills	0:714293de3836	604	{
ashleymills	0:714293de3836	605	int x, oldused;
ashleymills	0:714293de3836	606
ashleymills	0:714293de3836	607	oldused = b->used;
ashleymills	0:714293de3836	608	b->used = a->used;
ashleymills	0:714293de3836	609	{
ashleymills	0:714293de3836	610	register fp_digit r, rr, tmpa, tmpb;
ashleymills	0:714293de3836	611
ashleymills	0:714293de3836	612	/* source alias */
ashleymills	0:714293de3836	613	tmpa = a->dp + b->used - 1;
ashleymills	0:714293de3836	614
ashleymills	0:714293de3836	615	/* dest alias */
ashleymills	0:714293de3836	616	tmpb = b->dp + b->used - 1;
ashleymills	0:714293de3836	617
ashleymills	0:714293de3836	618	/* carry */
ashleymills	0:714293de3836	619	r = 0;
ashleymills	0:714293de3836	620	for (x = b->used - 1; x >= 0; x--) {
ashleymills	0:714293de3836	621	/* get the carry for the next iteration */
ashleymills	0:714293de3836	622	rr = *tmpa & 1;
ashleymills	0:714293de3836	623
ashleymills	0:714293de3836	624	/* shift the current digit, add in carry and store */
ashleymills	0:714293de3836	625	tmpb-- = (tmpa-- >> 1) \| (r << (DIGIT_BIT - 1));
ashleymills	0:714293de3836	626
ashleymills	0:714293de3836	627	/* forward carry to next iteration */
ashleymills	0:714293de3836	628	r = rr;
ashleymills	0:714293de3836	629	}
ashleymills	0:714293de3836	630
ashleymills	0:714293de3836	631	/* zero excess digits */
ashleymills	0:714293de3836	632	tmpb = b->dp + b->used;
ashleymills	0:714293de3836	633	for (x = b->used; x < oldused; x++) {
ashleymills	0:714293de3836	634	*tmpb++ = 0;
ashleymills	0:714293de3836	635	}
ashleymills	0:714293de3836	636	}
ashleymills	0:714293de3836	637	b->sign = a->sign;
ashleymills	0:714293de3836	638	fp_clamp (b);
ashleymills	0:714293de3836	639	}
ashleymills	0:714293de3836	640
ashleymills	0:714293de3836	641	/* c = a / 2*b /
ashleymills	0:714293de3836	642	void fp_div_2d(fp_int a, int b, fp_int c, fp_int *d)
ashleymills	0:714293de3836	643	{
ashleymills	0:714293de3836	644	fp_digit D, r, rr;
ashleymills	0:714293de3836	645	int x;
ashleymills	0:714293de3836	646	fp_int t;
ashleymills	0:714293de3836	647
ashleymills	0:714293de3836	648	/* if the shift count is <= 0 then we do no work */
ashleymills	0:714293de3836	649	if (b <= 0) {
ashleymills	0:714293de3836	650	fp_copy (a, c);
ashleymills	0:714293de3836	651	if (d != NULL) {
ashleymills	0:714293de3836	652	fp_zero (d);
ashleymills	0:714293de3836	653	}
ashleymills	0:714293de3836	654	return;
ashleymills	0:714293de3836	655	}
ashleymills	0:714293de3836	656
ashleymills	0:714293de3836	657	fp_init(&t);
ashleymills	0:714293de3836	658
ashleymills	0:714293de3836	659	/* get the remainder */
ashleymills	0:714293de3836	660	if (d != NULL) {
ashleymills	0:714293de3836	661	fp_mod_2d (a, b, &t);
ashleymills	0:714293de3836	662	}
ashleymills	0:714293de3836	663
ashleymills	0:714293de3836	664	/* copy */
ashleymills	0:714293de3836	665	fp_copy(a, c);
ashleymills	0:714293de3836	666
ashleymills	0:714293de3836	667	/* shift by as many digits in the bit count */
ashleymills	0:714293de3836	668	if (b >= (int)DIGIT_BIT) {
ashleymills	0:714293de3836	669	fp_rshd (c, b / DIGIT_BIT);
ashleymills	0:714293de3836	670	}
ashleymills	0:714293de3836	671
ashleymills	0:714293de3836	672	/* shift any bit count < DIGIT_BIT */
ashleymills	0:714293de3836	673	D = (fp_digit) (b % DIGIT_BIT);
ashleymills	0:714293de3836	674	if (D != 0) {
ashleymills	0:714293de3836	675	register fp_digit *tmpc, mask, shift;
ashleymills	0:714293de3836	676
ashleymills	0:714293de3836	677	/* mask */
ashleymills	0:714293de3836	678	mask = (((fp_digit)1) << D) - 1;
ashleymills	0:714293de3836	679
ashleymills	0:714293de3836	680	/* shift for lsb */
ashleymills	0:714293de3836	681	shift = DIGIT_BIT - D;
ashleymills	0:714293de3836	682
ashleymills	0:714293de3836	683	/* alias */
ashleymills	0:714293de3836	684	tmpc = c->dp + (c->used - 1);
ashleymills	0:714293de3836	685
ashleymills	0:714293de3836	686	/* carry */
ashleymills	0:714293de3836	687	r = 0;
ashleymills	0:714293de3836	688	for (x = c->used - 1; x >= 0; x--) {
ashleymills	0:714293de3836	689	/* get the lower bits of this word in a temp */
ashleymills	0:714293de3836	690	rr = *tmpc & mask;
ashleymills	0:714293de3836	691
ashleymills	0:714293de3836	692	/* shift the current word and mix in the carry bits from the previous word */
ashleymills	0:714293de3836	693	tmpc = (tmpc >> D) \| (r << shift);
ashleymills	0:714293de3836	694	--tmpc;
ashleymills	0:714293de3836	695
ashleymills	0:714293de3836	696	/* set the carry to the carry bits of the current word found above */
ashleymills	0:714293de3836	697	r = rr;
ashleymills	0:714293de3836	698	}
ashleymills	0:714293de3836	699	}
ashleymills	0:714293de3836	700	fp_clamp (c);
ashleymills	0:714293de3836	701	if (d != NULL) {
ashleymills	0:714293de3836	702	fp_copy (&t, d);
ashleymills	0:714293de3836	703	}
ashleymills	0:714293de3836	704	}
ashleymills	0:714293de3836	705
ashleymills	0:714293de3836	706	/* c = a mod b, 0 <= c < b */
ashleymills	0:714293de3836	707	int fp_mod(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	708	{
ashleymills	0:714293de3836	709	fp_int t;
ashleymills	0:714293de3836	710	int err;
ashleymills	0:714293de3836	711
ashleymills	0:714293de3836	712	fp_zero(&t);
ashleymills	0:714293de3836	713	if ((err = fp_div(a, b, NULL, &t)) != FP_OKAY) {
ashleymills	0:714293de3836	714	return err;
ashleymills	0:714293de3836	715	}
ashleymills	0:714293de3836	716	if (t.sign != b->sign) {
ashleymills	0:714293de3836	717	fp_add(&t, b, c);
ashleymills	0:714293de3836	718	} else {
ashleymills	0:714293de3836	719	fp_copy(&t, c);
ashleymills	0:714293de3836	720	}
ashleymills	0:714293de3836	721	return FP_OKAY;
ashleymills	0:714293de3836	722	}
ashleymills	0:714293de3836	723
ashleymills	0:714293de3836	724	/* c = a mod 2*d /
ashleymills	0:714293de3836	725	void fp_mod_2d(fp_int a, int b, fp_int c)
ashleymills	0:714293de3836	726	{
ashleymills	0:714293de3836	727	int x;
ashleymills	0:714293de3836	728
ashleymills	0:714293de3836	729	/* zero if count less than or equal to zero */
ashleymills	0:714293de3836	730	if (b <= 0) {
ashleymills	0:714293de3836	731	fp_zero(c);
ashleymills	0:714293de3836	732	return;
ashleymills	0:714293de3836	733	}
ashleymills	0:714293de3836	734
ashleymills	0:714293de3836	735	/* get copy of input */
ashleymills	0:714293de3836	736	fp_copy(a, c);
ashleymills	0:714293de3836	737
ashleymills	0:714293de3836	738	/* if 2*d is larger than we just return /
ashleymills	0:714293de3836	739	if (b >= (DIGIT_BIT * a->used)) {
ashleymills	0:714293de3836	740	return;
ashleymills	0:714293de3836	741	}
ashleymills	0:714293de3836	742
ashleymills	0:714293de3836	743	/* zero digits above the last digit of the modulus */
ashleymills	0:714293de3836	744	for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
ashleymills	0:714293de3836	745	c->dp[x] = 0;
ashleymills	0:714293de3836	746	}
ashleymills	0:714293de3836	747	/* clear the digit that is not completely outside/inside the modulus */
ashleymills	0:714293de3836	748	c->dp[b / DIGIT_BIT] &= ~((fp_digit)0) >> (DIGIT_BIT - b);
ashleymills	0:714293de3836	749	fp_clamp (c);
ashleymills	0:714293de3836	750	}
ashleymills	0:714293de3836	751
ashleymills	0:714293de3836	752	static int fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
ashleymills	0:714293de3836	753	{
ashleymills	0:714293de3836	754	fp_int x, y, u, v, A, B, C, D;
ashleymills	0:714293de3836	755	int res;
ashleymills	0:714293de3836	756
ashleymills	0:714293de3836	757	/* b cannot be negative */
ashleymills	0:714293de3836	758	if (b->sign == FP_NEG \|\| fp_iszero(b) == 1) {
ashleymills	0:714293de3836	759	return FP_VAL;
ashleymills	0:714293de3836	760	}
ashleymills	0:714293de3836	761
ashleymills	0:714293de3836	762	/* init temps */
ashleymills	0:714293de3836	763	fp_init(&x); fp_init(&y);
ashleymills	0:714293de3836	764	fp_init(&u); fp_init(&v);
ashleymills	0:714293de3836	765	fp_init(&A); fp_init(&B);
ashleymills	0:714293de3836	766	fp_init(&C); fp_init(&D);
ashleymills	0:714293de3836	767
ashleymills	0:714293de3836	768	/* x = a, y = b */
ashleymills	0:714293de3836	769	if ((res = fp_mod(a, b, &x)) != FP_OKAY) {
ashleymills	0:714293de3836	770	return res;
ashleymills	0:714293de3836	771	}
ashleymills	0:714293de3836	772	fp_copy(b, &y);
ashleymills	0:714293de3836	773
ashleymills	0:714293de3836	774	/* 2. [modified] if x,y are both even then return an error! */
ashleymills	0:714293de3836	775	if (fp_iseven (&x) == 1 && fp_iseven (&y) == 1) {
ashleymills	0:714293de3836	776	return FP_VAL;
ashleymills	0:714293de3836	777	}
ashleymills	0:714293de3836	778
ashleymills	0:714293de3836	779	/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
ashleymills	0:714293de3836	780	fp_copy (&x, &u);
ashleymills	0:714293de3836	781	fp_copy (&y, &v);
ashleymills	0:714293de3836	782	fp_set (&A, 1);
ashleymills	0:714293de3836	783	fp_set (&D, 1);
ashleymills	0:714293de3836	784
ashleymills	0:714293de3836	785	top:
ashleymills	0:714293de3836	786	/* 4. while u is even do */
ashleymills	0:714293de3836	787	while (fp_iseven (&u) == 1) {
ashleymills	0:714293de3836	788	/* 4.1 u = u/2 */
ashleymills	0:714293de3836	789	fp_div_2 (&u, &u);
ashleymills	0:714293de3836	790
ashleymills	0:714293de3836	791	/* 4.2 if A or B is odd then */
ashleymills	0:714293de3836	792	if (fp_isodd (&A) == 1 \|\| fp_isodd (&B) == 1) {
ashleymills	0:714293de3836	793	/* A = (A+y)/2, B = (B-x)/2 */
ashleymills	0:714293de3836	794	fp_add (&A, &y, &A);
ashleymills	0:714293de3836	795	fp_sub (&B, &x, &B);
ashleymills	0:714293de3836	796	}
ashleymills	0:714293de3836	797	/* A = A/2, B = B/2 */
ashleymills	0:714293de3836	798	fp_div_2 (&A, &A);
ashleymills	0:714293de3836	799	fp_div_2 (&B, &B);
ashleymills	0:714293de3836	800	}
ashleymills	0:714293de3836	801
ashleymills	0:714293de3836	802	/* 5. while v is even do */
ashleymills	0:714293de3836	803	while (fp_iseven (&v) == 1) {
ashleymills	0:714293de3836	804	/* 5.1 v = v/2 */
ashleymills	0:714293de3836	805	fp_div_2 (&v, &v);
ashleymills	0:714293de3836	806
ashleymills	0:714293de3836	807	/* 5.2 if C or D is odd then */
ashleymills	0:714293de3836	808	if (fp_isodd (&C) == 1 \|\| fp_isodd (&D) == 1) {
ashleymills	0:714293de3836	809	/* C = (C+y)/2, D = (D-x)/2 */
ashleymills	0:714293de3836	810	fp_add (&C, &y, &C);
ashleymills	0:714293de3836	811	fp_sub (&D, &x, &D);
ashleymills	0:714293de3836	812	}
ashleymills	0:714293de3836	813	/* C = C/2, D = D/2 */
ashleymills	0:714293de3836	814	fp_div_2 (&C, &C);
ashleymills	0:714293de3836	815	fp_div_2 (&D, &D);
ashleymills	0:714293de3836	816	}
ashleymills	0:714293de3836	817
ashleymills	0:714293de3836	818	/* 6. if u >= v then */
ashleymills	0:714293de3836	819	if (fp_cmp (&u, &v) != FP_LT) {
ashleymills	0:714293de3836	820	/* u = u - v, A = A - C, B = B - D */
ashleymills	0:714293de3836	821	fp_sub (&u, &v, &u);
ashleymills	0:714293de3836	822	fp_sub (&A, &C, &A);
ashleymills	0:714293de3836	823	fp_sub (&B, &D, &B);
ashleymills	0:714293de3836	824	} else {
ashleymills	0:714293de3836	825	/* v - v - u, C = C - A, D = D - B */
ashleymills	0:714293de3836	826	fp_sub (&v, &u, &v);
ashleymills	0:714293de3836	827	fp_sub (&C, &A, &C);
ashleymills	0:714293de3836	828	fp_sub (&D, &B, &D);
ashleymills	0:714293de3836	829	}
ashleymills	0:714293de3836	830
ashleymills	0:714293de3836	831	/* if not zero goto step 4 */
ashleymills	0:714293de3836	832	if (fp_iszero (&u) == 0)
ashleymills	0:714293de3836	833	goto top;
ashleymills	0:714293de3836	834
ashleymills	0:714293de3836	835	/* now a = C, b = D, gcd == gv /
ashleymills	0:714293de3836	836
ashleymills	0:714293de3836	837	/* if v != 1 then there is no inverse */
ashleymills	0:714293de3836	838	if (fp_cmp_d (&v, 1) != FP_EQ) {
ashleymills	0:714293de3836	839	return FP_VAL;
ashleymills	0:714293de3836	840	}
ashleymills	0:714293de3836	841
ashleymills	0:714293de3836	842	/* if its too low */
ashleymills	0:714293de3836	843	while (fp_cmp_d(&C, 0) == FP_LT) {
ashleymills	0:714293de3836	844	fp_add(&C, b, &C);
ashleymills	0:714293de3836	845	}
ashleymills	0:714293de3836	846
ashleymills	0:714293de3836	847	/* too big */
ashleymills	0:714293de3836	848	while (fp_cmp_mag(&C, b) != FP_LT) {
ashleymills	0:714293de3836	849	fp_sub(&C, b, &C);
ashleymills	0:714293de3836	850	}
ashleymills	0:714293de3836	851
ashleymills	0:714293de3836	852	/* C is now the inverse */
ashleymills	0:714293de3836	853	fp_copy(&C, c);
ashleymills	0:714293de3836	854	return FP_OKAY;
ashleymills	0:714293de3836	855	}
ashleymills	0:714293de3836	856
ashleymills	0:714293de3836	857	/* c = 1/a (mod b) for odd b only */
ashleymills	0:714293de3836	858	int fp_invmod(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	859	{
ashleymills	0:714293de3836	860	fp_int x, y, u, v, B, D;
ashleymills	0:714293de3836	861	int neg;
ashleymills	0:714293de3836	862
ashleymills	0:714293de3836	863	/* 2. [modified] b must be odd */
ashleymills	0:714293de3836	864	if (fp_iseven (b) == FP_YES) {
ashleymills	0:714293de3836	865	return fp_invmod_slow(a,b,c);
ashleymills	0:714293de3836	866	}
ashleymills	0:714293de3836	867
ashleymills	0:714293de3836	868	/* init all our temps */
ashleymills	0:714293de3836	869	fp_init(&x); fp_init(&y);
ashleymills	0:714293de3836	870	fp_init(&u); fp_init(&v);
ashleymills	0:714293de3836	871	fp_init(&B); fp_init(&D);
ashleymills	0:714293de3836	872
ashleymills	0:714293de3836	873	/* x == modulus, y == value to invert */
ashleymills	0:714293de3836	874	fp_copy(b, &x);
ashleymills	0:714293de3836	875
ashleymills	0:714293de3836	876	/* we need y = \|a\| */
ashleymills	0:714293de3836	877	fp_abs(a, &y);
ashleymills	0:714293de3836	878
ashleymills	0:714293de3836	879	/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
ashleymills	0:714293de3836	880	fp_copy(&x, &u);
ashleymills	0:714293de3836	881	fp_copy(&y, &v);
ashleymills	0:714293de3836	882	fp_set (&D, 1);
ashleymills	0:714293de3836	883
ashleymills	0:714293de3836	884	top:
ashleymills	0:714293de3836	885	/* 4. while u is even do */
ashleymills	0:714293de3836	886	while (fp_iseven (&u) == FP_YES) {
ashleymills	0:714293de3836	887	/* 4.1 u = u/2 */
ashleymills	0:714293de3836	888	fp_div_2 (&u, &u);
ashleymills	0:714293de3836	889
ashleymills	0:714293de3836	890	/* 4.2 if B is odd then */
ashleymills	0:714293de3836	891	if (fp_isodd (&B) == FP_YES) {
ashleymills	0:714293de3836	892	fp_sub (&B, &x, &B);
ashleymills	0:714293de3836	893	}
ashleymills	0:714293de3836	894	/* B = B/2 */
ashleymills	0:714293de3836	895	fp_div_2 (&B, &B);
ashleymills	0:714293de3836	896	}
ashleymills	0:714293de3836	897
ashleymills	0:714293de3836	898	/* 5. while v is even do */
ashleymills	0:714293de3836	899	while (fp_iseven (&v) == FP_YES) {
ashleymills	0:714293de3836	900	/* 5.1 v = v/2 */
ashleymills	0:714293de3836	901	fp_div_2 (&v, &v);
ashleymills	0:714293de3836	902
ashleymills	0:714293de3836	903	/* 5.2 if D is odd then */
ashleymills	0:714293de3836	904	if (fp_isodd (&D) == FP_YES) {
ashleymills	0:714293de3836	905	/* D = (D-x)/2 */
ashleymills	0:714293de3836	906	fp_sub (&D, &x, &D);
ashleymills	0:714293de3836	907	}
ashleymills	0:714293de3836	908	/* D = D/2 */
ashleymills	0:714293de3836	909	fp_div_2 (&D, &D);
ashleymills	0:714293de3836	910	}
ashleymills	0:714293de3836	911
ashleymills	0:714293de3836	912	/* 6. if u >= v then */
ashleymills	0:714293de3836	913	if (fp_cmp (&u, &v) != FP_LT) {
ashleymills	0:714293de3836	914	/* u = u - v, B = B - D */
ashleymills	0:714293de3836	915	fp_sub (&u, &v, &u);
ashleymills	0:714293de3836	916	fp_sub (&B, &D, &B);
ashleymills	0:714293de3836	917	} else {
ashleymills	0:714293de3836	918	/* v - v - u, D = D - B */
ashleymills	0:714293de3836	919	fp_sub (&v, &u, &v);
ashleymills	0:714293de3836	920	fp_sub (&D, &B, &D);
ashleymills	0:714293de3836	921	}
ashleymills	0:714293de3836	922
ashleymills	0:714293de3836	923	/* if not zero goto step 4 */
ashleymills	0:714293de3836	924	if (fp_iszero (&u) == FP_NO) {
ashleymills	0:714293de3836	925	goto top;
ashleymills	0:714293de3836	926	}
ashleymills	0:714293de3836	927
ashleymills	0:714293de3836	928	/* now a = C, b = D, gcd == gv /
ashleymills	0:714293de3836	929
ashleymills	0:714293de3836	930	/* if v != 1 then there is no inverse */
ashleymills	0:714293de3836	931	if (fp_cmp_d (&v, 1) != FP_EQ) {
ashleymills	0:714293de3836	932	return FP_VAL;
ashleymills	0:714293de3836	933	}
ashleymills	0:714293de3836	934
ashleymills	0:714293de3836	935	/* b is now the inverse */
ashleymills	0:714293de3836	936	neg = a->sign;
ashleymills	0:714293de3836	937	while (D.sign == FP_NEG) {
ashleymills	0:714293de3836	938	fp_add (&D, b, &D);
ashleymills	0:714293de3836	939	}
ashleymills	0:714293de3836	940	fp_copy (&D, c);
ashleymills	0:714293de3836	941	c->sign = neg;
ashleymills	0:714293de3836	942	return FP_OKAY;
ashleymills	0:714293de3836	943	}
ashleymills	0:714293de3836	944
ashleymills	0:714293de3836	945	/* d = a * b (mod c) */
ashleymills	0:714293de3836	946	int fp_mulmod(fp_int a, fp_int b, fp_int c, fp_int d)
ashleymills	0:714293de3836	947	{
ashleymills	0:714293de3836	948	fp_int tmp;
ashleymills	0:714293de3836	949	fp_zero(&tmp);
ashleymills	0:714293de3836	950	fp_mul(a, b, &tmp);
ashleymills	0:714293de3836	951	return fp_mod(&tmp, c, d);
ashleymills	0:714293de3836	952	}
ashleymills	0:714293de3836	953
ashleymills	0:714293de3836	954	#ifdef TFM_TIMING_RESISTANT
ashleymills	0:714293de3836	955
ashleymills	0:714293de3836	956	/* timing resistant montgomery ladder based exptmod
ashleymills	0:714293de3836	957
ashleymills	0:714293de3836	958	Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
ashleymills	0:714293de3836	959	*/
ashleymills	0:714293de3836	960	static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills	0:714293de3836	961	{
ashleymills	0:714293de3836	962	fp_int R[2];
ashleymills	0:714293de3836	963	fp_digit buf, mp;
ashleymills	0:714293de3836	964	int err, bitcnt, digidx, y;
ashleymills	0:714293de3836	965
ashleymills	0:714293de3836	966	/* now setup montgomery */
ashleymills	0:714293de3836	967	if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
ashleymills	0:714293de3836	968	return err;
ashleymills	0:714293de3836	969	}
ashleymills	0:714293de3836	970
ashleymills	0:714293de3836	971	fp_init(&R[0]);
ashleymills	0:714293de3836	972	fp_init(&R[1]);
ashleymills	0:714293de3836	973
ashleymills	0:714293de3836	974	/* now we need R mod m */
ashleymills	0:714293de3836	975	fp_montgomery_calc_normalization (&R[0], P);
ashleymills	0:714293de3836	976
ashleymills	0:714293de3836	977	/* now set R[0][1] to G * R mod m */
ashleymills	0:714293de3836	978	if (fp_cmp_mag(P, G) != FP_GT) {
ashleymills	0:714293de3836	979	/* G > P so we reduce it first */
ashleymills	0:714293de3836	980	fp_mod(G, P, &R[1]);
ashleymills	0:714293de3836	981	} else {
ashleymills	0:714293de3836	982	fp_copy(G, &R[1]);
ashleymills	0:714293de3836	983	}
ashleymills	0:714293de3836	984	fp_mulmod (&R[1], &R[0], P, &R[1]);
ashleymills	0:714293de3836	985
ashleymills	0:714293de3836	986	/* for j = t-1 downto 0 do
ashleymills	0:714293de3836	987	r_!k = R0*R1; r_k = r_k^2
ashleymills	0:714293de3836	988	*/
ashleymills	0:714293de3836	989
ashleymills	0:714293de3836	990	/* set initial mode and bit cnt */
ashleymills	0:714293de3836	991	bitcnt = 1;
ashleymills	0:714293de3836	992	buf = 0;
ashleymills	0:714293de3836	993	digidx = X->used - 1;
ashleymills	0:714293de3836	994
ashleymills	0:714293de3836	995	for (;;) {
ashleymills	0:714293de3836	996	/* grab next digit as required */
ashleymills	0:714293de3836	997	if (--bitcnt == 0) {
ashleymills	0:714293de3836	998	/* if digidx == -1 we are out of digits so break */
ashleymills	0:714293de3836	999	if (digidx == -1) {
ashleymills	0:714293de3836	1000	break;
ashleymills	0:714293de3836	1001	}
ashleymills	0:714293de3836	1002	/* read next digit and reset bitcnt */
ashleymills	0:714293de3836	1003	buf = X->dp[digidx--];
ashleymills	0:714293de3836	1004	bitcnt = (int)DIGIT_BIT;
ashleymills	0:714293de3836	1005	}
ashleymills	0:714293de3836	1006
ashleymills	0:714293de3836	1007	/* grab the next msb from the exponent */
ashleymills	0:714293de3836	1008	y = (int)(buf >> (DIGIT_BIT - 1)) & 1;
ashleymills	0:714293de3836	1009	buf <<= (fp_digit)1;
ashleymills	0:714293de3836	1010
ashleymills	0:714293de3836	1011	/* do ops */
ashleymills	0:714293de3836	1012	fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
ashleymills	0:714293de3836	1013	fp_sqr(&R[y], &R[y]); fp_montgomery_reduce(&R[y], P, mp);
ashleymills	0:714293de3836	1014	}
ashleymills	0:714293de3836	1015
ashleymills	0:714293de3836	1016	fp_montgomery_reduce(&R[0], P, mp);
ashleymills	0:714293de3836	1017	fp_copy(&R[0], Y);
ashleymills	0:714293de3836	1018	return FP_OKAY;
ashleymills	0:714293de3836	1019	}
ashleymills	0:714293de3836	1020
ashleymills	0:714293de3836	1021	#else
ashleymills	0:714293de3836	1022
ashleymills	0:714293de3836	1023	/* y = g**x (mod b)
ashleymills	0:714293de3836	1024	* Some restrictions... x must be positive and < b
ashleymills	0:714293de3836	1025	*/
ashleymills	0:714293de3836	1026	static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills	0:714293de3836	1027	{
ashleymills	0:714293de3836	1028	fp_int M[64], res;
ashleymills	0:714293de3836	1029	fp_digit buf, mp;
ashleymills	0:714293de3836	1030	int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
ashleymills	0:714293de3836	1031
ashleymills	0:714293de3836	1032	/* find window size */
ashleymills	0:714293de3836	1033	x = fp_count_bits (X);
ashleymills	0:714293de3836	1034	if (x <= 21) {
ashleymills	0:714293de3836	1035	winsize = 1;
ashleymills	0:714293de3836	1036	} else if (x <= 36) {
ashleymills	0:714293de3836	1037	winsize = 3;
ashleymills	0:714293de3836	1038	} else if (x <= 140) {
ashleymills	0:714293de3836	1039	winsize = 4;
ashleymills	0:714293de3836	1040	} else if (x <= 450) {
ashleymills	0:714293de3836	1041	winsize = 5;
ashleymills	0:714293de3836	1042	} else {
ashleymills	0:714293de3836	1043	winsize = 6;
ashleymills	0:714293de3836	1044	}
ashleymills	0:714293de3836	1045
ashleymills	0:714293de3836	1046	/* init M array */
ashleymills	0:714293de3836	1047	XMEMSET(M, 0, sizeof(M));
ashleymills	0:714293de3836	1048
ashleymills	0:714293de3836	1049	/* now setup montgomery */
ashleymills	0:714293de3836	1050	if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
ashleymills	0:714293de3836	1051	return err;
ashleymills	0:714293de3836	1052	}
ashleymills	0:714293de3836	1053
ashleymills	0:714293de3836	1054	/* setup result */
ashleymills	0:714293de3836	1055	fp_init(&res);
ashleymills	0:714293de3836	1056
ashleymills	0:714293de3836	1057	/* create M table
ashleymills	0:714293de3836	1058	*
ashleymills	0:714293de3836	1059	* The M table contains powers of the input base, e.g. M[x] = G^x mod P
ashleymills	0:714293de3836	1060	*
ashleymills	0:714293de3836	1061	* The first half of the table is not computed though accept for M[0] and M[1]
ashleymills	0:714293de3836	1062	*/
ashleymills	0:714293de3836	1063
ashleymills	0:714293de3836	1064	/* now we need R mod m */
ashleymills	0:714293de3836	1065	fp_montgomery_calc_normalization (&res, P);
ashleymills	0:714293de3836	1066
ashleymills	0:714293de3836	1067	/* now set M[1] to G * R mod m */
ashleymills	0:714293de3836	1068	if (fp_cmp_mag(P, G) != FP_GT) {
ashleymills	0:714293de3836	1069	/* G > P so we reduce it first */
ashleymills	0:714293de3836	1070	fp_mod(G, P, &M[1]);
ashleymills	0:714293de3836	1071	} else {
ashleymills	0:714293de3836	1072	fp_copy(G, &M[1]);
ashleymills	0:714293de3836	1073	}
ashleymills	0:714293de3836	1074	fp_mulmod (&M[1], &res, P, &M[1]);
ashleymills	0:714293de3836	1075
ashleymills	0:714293de3836	1076	/* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
ashleymills	0:714293de3836	1077	fp_copy (&M[1], &M[1 << (winsize - 1)]);
ashleymills	0:714293de3836	1078	for (x = 0; x < (winsize - 1); x++) {
ashleymills	0:714293de3836	1079	fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
ashleymills	0:714293de3836	1080	fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
ashleymills	0:714293de3836	1081	}
ashleymills	0:714293de3836	1082
ashleymills	0:714293de3836	1083	/* create upper table */
ashleymills	0:714293de3836	1084	for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
ashleymills	0:714293de3836	1085	fp_mul(&M[x - 1], &M[1], &M[x]);
ashleymills	0:714293de3836	1086	fp_montgomery_reduce(&M[x], P, mp);
ashleymills	0:714293de3836	1087	}
ashleymills	0:714293de3836	1088
ashleymills	0:714293de3836	1089	/* set initial mode and bit cnt */
ashleymills	0:714293de3836	1090	mode = 0;
ashleymills	0:714293de3836	1091	bitcnt = 1;
ashleymills	0:714293de3836	1092	buf = 0;
ashleymills	0:714293de3836	1093	digidx = X->used - 1;
ashleymills	0:714293de3836	1094	bitcpy = 0;
ashleymills	0:714293de3836	1095	bitbuf = 0;
ashleymills	0:714293de3836	1096
ashleymills	0:714293de3836	1097	for (;;) {
ashleymills	0:714293de3836	1098	/* grab next digit as required */
ashleymills	0:714293de3836	1099	if (--bitcnt == 0) {
ashleymills	0:714293de3836	1100	/* if digidx == -1 we are out of digits so break */
ashleymills	0:714293de3836	1101	if (digidx == -1) {
ashleymills	0:714293de3836	1102	break;
ashleymills	0:714293de3836	1103	}
ashleymills	0:714293de3836	1104	/* read next digit and reset bitcnt */
ashleymills	0:714293de3836	1105	buf = X->dp[digidx--];
ashleymills	0:714293de3836	1106	bitcnt = (int)DIGIT_BIT;
ashleymills	0:714293de3836	1107	}
ashleymills	0:714293de3836	1108
ashleymills	0:714293de3836	1109	/* grab the next msb from the exponent */
ashleymills	0:714293de3836	1110	y = (int)(buf >> (DIGIT_BIT - 1)) & 1;
ashleymills	0:714293de3836	1111	buf <<= (fp_digit)1;
ashleymills	0:714293de3836	1112
ashleymills	0:714293de3836	1113	/* if the bit is zero and mode == 0 then we ignore it
ashleymills	0:714293de3836	1114	* These represent the leading zero bits before the first 1 bit
ashleymills	0:714293de3836	1115	* in the exponent. Technically this opt is not required but it
ashleymills	0:714293de3836	1116	* does lower the # of trivial squaring/reductions used
ashleymills	0:714293de3836	1117	*/
ashleymills	0:714293de3836	1118	if (mode == 0 && y == 0) {
ashleymills	0:714293de3836	1119	continue;
ashleymills	0:714293de3836	1120	}
ashleymills	0:714293de3836	1121
ashleymills	0:714293de3836	1122	/* if the bit is zero and mode == 1 then we square */
ashleymills	0:714293de3836	1123	if (mode == 1 && y == 0) {
ashleymills	0:714293de3836	1124	fp_sqr(&res, &res);
ashleymills	0:714293de3836	1125	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1126	continue;
ashleymills	0:714293de3836	1127	}
ashleymills	0:714293de3836	1128
ashleymills	0:714293de3836	1129	/* else we add it to the window */
ashleymills	0:714293de3836	1130	bitbuf \|= (y << (winsize - ++bitcpy));
ashleymills	0:714293de3836	1131	mode = 2;
ashleymills	0:714293de3836	1132
ashleymills	0:714293de3836	1133	if (bitcpy == winsize) {
ashleymills	0:714293de3836	1134	/* ok window is filled so square as required and multiply */
ashleymills	0:714293de3836	1135	/* square first */
ashleymills	0:714293de3836	1136	for (x = 0; x < winsize; x++) {
ashleymills	0:714293de3836	1137	fp_sqr(&res, &res);
ashleymills	0:714293de3836	1138	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1139	}
ashleymills	0:714293de3836	1140
ashleymills	0:714293de3836	1141	/* then multiply */
ashleymills	0:714293de3836	1142	fp_mul(&res, &M[bitbuf], &res);
ashleymills	0:714293de3836	1143	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1144
ashleymills	0:714293de3836	1145	/* empty window and reset */
ashleymills	0:714293de3836	1146	bitcpy = 0;
ashleymills	0:714293de3836	1147	bitbuf = 0;
ashleymills	0:714293de3836	1148	mode = 1;
ashleymills	0:714293de3836	1149	}
ashleymills	0:714293de3836	1150	}
ashleymills	0:714293de3836	1151
ashleymills	0:714293de3836	1152	/* if bits remain then square/multiply */
ashleymills	0:714293de3836	1153	if (mode == 2 && bitcpy > 0) {
ashleymills	0:714293de3836	1154	/* square then multiply if the bit is set */
ashleymills	0:714293de3836	1155	for (x = 0; x < bitcpy; x++) {
ashleymills	0:714293de3836	1156	fp_sqr(&res, &res);
ashleymills	0:714293de3836	1157	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1158
ashleymills	0:714293de3836	1159	/* get next bit of the window */
ashleymills	0:714293de3836	1160	bitbuf <<= 1;
ashleymills	0:714293de3836	1161	if ((bitbuf & (1 << winsize)) != 0) {
ashleymills	0:714293de3836	1162	/* then multiply */
ashleymills	0:714293de3836	1163	fp_mul(&res, &M[1], &res);
ashleymills	0:714293de3836	1164	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1165	}
ashleymills	0:714293de3836	1166	}
ashleymills	0:714293de3836	1167	}
ashleymills	0:714293de3836	1168
ashleymills	0:714293de3836	1169	/* fixup result if Montgomery reduction is used
ashleymills	0:714293de3836	1170	* recall that any value in a Montgomery system is
ashleymills	0:714293de3836	1171	* actually multiplied by R mod n. So we have
ashleymills	0:714293de3836	1172	* to reduce one more time to cancel out the factor
ashleymills	0:714293de3836	1173	* of R.
ashleymills	0:714293de3836	1174	*/
ashleymills	0:714293de3836	1175	fp_montgomery_reduce(&res, P, mp);
ashleymills	0:714293de3836	1176
ashleymills	0:714293de3836	1177	/* swap res with Y */
ashleymills	0:714293de3836	1178	fp_copy (&res, Y);
ashleymills	0:714293de3836	1179	return FP_OKAY;
ashleymills	0:714293de3836	1180	}
ashleymills	0:714293de3836	1181
ashleymills	0:714293de3836	1182	#endif
ashleymills	0:714293de3836	1183
ashleymills	0:714293de3836	1184	int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills	0:714293de3836	1185	{
ashleymills	0:714293de3836	1186	/* prevent overflows */
ashleymills	0:714293de3836	1187	if (P->used > (FP_SIZE/2)) {
ashleymills	0:714293de3836	1188	return FP_VAL;
ashleymills	0:714293de3836	1189	}
ashleymills	0:714293de3836	1190
ashleymills	0:714293de3836	1191	if (X->sign == FP_NEG) {
ashleymills	0:714293de3836	1192	#ifndef POSITIVE_EXP_ONLY /* reduce stack if assume no negatives */
ashleymills	0:714293de3836	1193	int err;
ashleymills	0:714293de3836	1194	fp_int tmp;
ashleymills	0:714293de3836	1195
ashleymills	0:714293de3836	1196	/* yes, copy G and invmod it */
ashleymills	0:714293de3836	1197	fp_copy(G, &tmp);
ashleymills	0:714293de3836	1198	if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
ashleymills	0:714293de3836	1199	return err;
ashleymills	0:714293de3836	1200	}
ashleymills	0:714293de3836	1201	X->sign = FP_ZPOS;
ashleymills	0:714293de3836	1202	err = _fp_exptmod(&tmp, X, P, Y);
ashleymills	0:714293de3836	1203	if (X != Y) {
ashleymills	0:714293de3836	1204	X->sign = FP_NEG;
ashleymills	0:714293de3836	1205	}
ashleymills	0:714293de3836	1206	return err;
ashleymills	0:714293de3836	1207	#else
ashleymills	0:714293de3836	1208	return FP_VAL;
ashleymills	0:714293de3836	1209	#endif
ashleymills	0:714293de3836	1210	}
ashleymills	0:714293de3836	1211	else {
ashleymills	0:714293de3836	1212	/* Positive exponent so just exptmod */
ashleymills	0:714293de3836	1213	return _fp_exptmod(G, X, P, Y);
ashleymills	0:714293de3836	1214	}
ashleymills	0:714293de3836	1215	}
ashleymills	0:714293de3836	1216
ashleymills	0:714293de3836	1217	/* computes a = 2*b /
ashleymills	0:714293de3836	1218	void fp_2expt(fp_int *a, int b)
ashleymills	0:714293de3836	1219	{
ashleymills	0:714293de3836	1220	int z;
ashleymills	0:714293de3836	1221
ashleymills	0:714293de3836	1222	/* zero a as per default */
ashleymills	0:714293de3836	1223	fp_zero (a);
ashleymills	0:714293de3836	1224
ashleymills	0:714293de3836	1225	if (b < 0) {
ashleymills	0:714293de3836	1226	return;
ashleymills	0:714293de3836	1227	}
ashleymills	0:714293de3836	1228
ashleymills	0:714293de3836	1229	z = b / DIGIT_BIT;
ashleymills	0:714293de3836	1230	if (z >= FP_SIZE) {
ashleymills	0:714293de3836	1231	return;
ashleymills	0:714293de3836	1232	}
ashleymills	0:714293de3836	1233
ashleymills	0:714293de3836	1234	/* set the used count of where the bit will go */
ashleymills	0:714293de3836	1235	a->used = z + 1;
ashleymills	0:714293de3836	1236
ashleymills	0:714293de3836	1237	/* put the single bit in its place */
ashleymills	0:714293de3836	1238	a->dp[z] = ((fp_digit)1) << (b % DIGIT_BIT);
ashleymills	0:714293de3836	1239	}
ashleymills	0:714293de3836	1240
ashleymills	0:714293de3836	1241	/* b = aa /
ashleymills	0:714293de3836	1242	void fp_sqr(fp_int A, fp_int B)
ashleymills	0:714293de3836	1243	{
ashleymills	0:714293de3836	1244	int y = A->used;
ashleymills	0:714293de3836	1245
ashleymills	0:714293de3836	1246	/* call generic if we're out of range */
ashleymills	0:714293de3836	1247	if (y + y > FP_SIZE) {
ashleymills	0:714293de3836	1248	fp_sqr_comba(A, B);
ashleymills	0:714293de3836	1249	return ;
ashleymills	0:714293de3836	1250	}
ashleymills	0:714293de3836	1251
ashleymills	0:714293de3836	1252	#if defined(TFM_SQR3)
ashleymills	0:714293de3836	1253	if (y <= 3) {
ashleymills	0:714293de3836	1254	fp_sqr_comba3(A,B);
ashleymills	0:714293de3836	1255	return;
ashleymills	0:714293de3836	1256	}
ashleymills	0:714293de3836	1257	#endif
ashleymills	0:714293de3836	1258	#if defined(TFM_SQR4)
ashleymills	0:714293de3836	1259	if (y == 4) {
ashleymills	0:714293de3836	1260	fp_sqr_comba4(A,B);
ashleymills	0:714293de3836	1261	return;
ashleymills	0:714293de3836	1262	}
ashleymills	0:714293de3836	1263	#endif
ashleymills	0:714293de3836	1264	#if defined(TFM_SQR6)
ashleymills	0:714293de3836	1265	if (y <= 6) {
ashleymills	0:714293de3836	1266	fp_sqr_comba6(A,B);
ashleymills	0:714293de3836	1267	return;
ashleymills	0:714293de3836	1268	}
ashleymills	0:714293de3836	1269	#endif
ashleymills	0:714293de3836	1270	#if defined(TFM_SQR7)
ashleymills	0:714293de3836	1271	if (y == 7) {
ashleymills	0:714293de3836	1272	fp_sqr_comba7(A,B);
ashleymills	0:714293de3836	1273	return;
ashleymills	0:714293de3836	1274	}
ashleymills	0:714293de3836	1275	#endif
ashleymills	0:714293de3836	1276	#if defined(TFM_SQR8)
ashleymills	0:714293de3836	1277	if (y == 8) {
ashleymills	0:714293de3836	1278	fp_sqr_comba8(A,B);
ashleymills	0:714293de3836	1279	return;
ashleymills	0:714293de3836	1280	}
ashleymills	0:714293de3836	1281	#endif
ashleymills	0:714293de3836	1282	#if defined(TFM_SQR9)
ashleymills	0:714293de3836	1283	if (y == 9) {
ashleymills	0:714293de3836	1284	fp_sqr_comba9(A,B);
ashleymills	0:714293de3836	1285	return;
ashleymills	0:714293de3836	1286	}
ashleymills	0:714293de3836	1287	#endif
ashleymills	0:714293de3836	1288	#if defined(TFM_SQR12)
ashleymills	0:714293de3836	1289	if (y <= 12) {
ashleymills	0:714293de3836	1290	fp_sqr_comba12(A,B);
ashleymills	0:714293de3836	1291	return;
ashleymills	0:714293de3836	1292	}
ashleymills	0:714293de3836	1293	#endif
ashleymills	0:714293de3836	1294	#if defined(TFM_SQR17)
ashleymills	0:714293de3836	1295	if (y <= 17) {
ashleymills	0:714293de3836	1296	fp_sqr_comba17(A,B);
ashleymills	0:714293de3836	1297	return;
ashleymills	0:714293de3836	1298	}
ashleymills	0:714293de3836	1299	#endif
ashleymills	0:714293de3836	1300	#if defined(TFM_SMALL_SET)
ashleymills	0:714293de3836	1301	if (y <= 16) {
ashleymills	0:714293de3836	1302	fp_sqr_comba_small(A,B);
ashleymills	0:714293de3836	1303	return;
ashleymills	0:714293de3836	1304	}
ashleymills	0:714293de3836	1305	#endif
ashleymills	0:714293de3836	1306	#if defined(TFM_SQR20)
ashleymills	0:714293de3836	1307	if (y <= 20) {
ashleymills	0:714293de3836	1308	fp_sqr_comba20(A,B);
ashleymills	0:714293de3836	1309	return;
ashleymills	0:714293de3836	1310	}
ashleymills	0:714293de3836	1311	#endif
ashleymills	0:714293de3836	1312	#if defined(TFM_SQR24)
ashleymills	0:714293de3836	1313	if (y <= 24) {
ashleymills	0:714293de3836	1314	fp_sqr_comba24(A,B);
ashleymills	0:714293de3836	1315	return;
ashleymills	0:714293de3836	1316	}
ashleymills	0:714293de3836	1317	#endif
ashleymills	0:714293de3836	1318	#if defined(TFM_SQR28)
ashleymills	0:714293de3836	1319	if (y <= 28) {
ashleymills	0:714293de3836	1320	fp_sqr_comba28(A,B);
ashleymills	0:714293de3836	1321	return;
ashleymills	0:714293de3836	1322	}
ashleymills	0:714293de3836	1323	#endif
ashleymills	0:714293de3836	1324	#if defined(TFM_SQR32)
ashleymills	0:714293de3836	1325	if (y <= 32) {
ashleymills	0:714293de3836	1326	fp_sqr_comba32(A,B);
ashleymills	0:714293de3836	1327	return;
ashleymills	0:714293de3836	1328	}
ashleymills	0:714293de3836	1329	#endif
ashleymills	0:714293de3836	1330	#if defined(TFM_SQR48)
ashleymills	0:714293de3836	1331	if (y <= 48) {
ashleymills	0:714293de3836	1332	fp_sqr_comba48(A,B);
ashleymills	0:714293de3836	1333	return;
ashleymills	0:714293de3836	1334	}
ashleymills	0:714293de3836	1335	#endif
ashleymills	0:714293de3836	1336	#if defined(TFM_SQR64)
ashleymills	0:714293de3836	1337	if (y <= 64) {
ashleymills	0:714293de3836	1338	fp_sqr_comba64(A,B);
ashleymills	0:714293de3836	1339	return;
ashleymills	0:714293de3836	1340	}
ashleymills	0:714293de3836	1341	#endif
ashleymills	0:714293de3836	1342	fp_sqr_comba(A, B);
ashleymills	0:714293de3836	1343	}
ashleymills	0:714293de3836	1344
ashleymills	0:714293de3836	1345	/* generic comba squarer */
ashleymills	0:714293de3836	1346	void fp_sqr_comba(fp_int A, fp_int B)
ashleymills	0:714293de3836	1347	{
ashleymills	0:714293de3836	1348	int pa, ix, iz;
ashleymills	0:714293de3836	1349	fp_digit c0, c1, c2;
ashleymills	0:714293de3836	1350	fp_int tmp, *dst;
ashleymills	0:714293de3836	1351	#ifdef TFM_ISO
ashleymills	0:714293de3836	1352	fp_word tt;
ashleymills	0:714293de3836	1353	#endif
ashleymills	0:714293de3836	1354
ashleymills	0:714293de3836	1355	/* get size of output and trim */
ashleymills	0:714293de3836	1356	pa = A->used + A->used;
ashleymills	0:714293de3836	1357	if (pa >= FP_SIZE) {
ashleymills	0:714293de3836	1358	pa = FP_SIZE-1;
ashleymills	0:714293de3836	1359	}
ashleymills	0:714293de3836	1360
ashleymills	0:714293de3836	1361	/* number of output digits to produce */
ashleymills	0:714293de3836	1362	COMBA_START;
ashleymills	0:714293de3836	1363	COMBA_CLEAR;
ashleymills	0:714293de3836	1364
ashleymills	0:714293de3836	1365	if (A == B) {
ashleymills	0:714293de3836	1366	fp_zero(&tmp);
ashleymills	0:714293de3836	1367	dst = &tmp;
ashleymills	0:714293de3836	1368	} else {
ashleymills	0:714293de3836	1369	fp_zero(B);
ashleymills	0:714293de3836	1370	dst = B;
ashleymills	0:714293de3836	1371	}
ashleymills	0:714293de3836	1372
ashleymills	0:714293de3836	1373	for (ix = 0; ix < pa; ix++) {
ashleymills	0:714293de3836	1374	int tx, ty, iy;
ashleymills	0:714293de3836	1375	fp_digit tmpy, tmpx;
ashleymills	0:714293de3836	1376
ashleymills	0:714293de3836	1377	/* get offsets into the two bignums */
ashleymills	0:714293de3836	1378	ty = MIN(A->used-1, ix);
ashleymills	0:714293de3836	1379	tx = ix - ty;
ashleymills	0:714293de3836	1380
ashleymills	0:714293de3836	1381	/* setup temp aliases */
ashleymills	0:714293de3836	1382	tmpx = A->dp + tx;
ashleymills	0:714293de3836	1383	tmpy = A->dp + ty;
ashleymills	0:714293de3836	1384
ashleymills	0:714293de3836	1385	/* this is the number of times the loop will iterrate,
ashleymills	0:714293de3836	1386	while (tx++ < a->used && ty-- >= 0) { ... }
ashleymills	0:714293de3836	1387	*/
ashleymills	0:714293de3836	1388	iy = MIN(A->used-tx, ty+1);
ashleymills	0:714293de3836	1389
ashleymills	0:714293de3836	1390	/* now for squaring tx can never equal ty
ashleymills	0:714293de3836	1391	* we halve the distance since they approach
ashleymills	0:714293de3836	1392	* at a rate of 2x and we have to round because
ashleymills	0:714293de3836	1393	* odd cases need to be executed
ashleymills	0:714293de3836	1394	*/
ashleymills	0:714293de3836	1395	iy = MIN(iy, (ty-tx+1)>>1);
ashleymills	0:714293de3836	1396
ashleymills	0:714293de3836	1397	/* forward carries */
ashleymills	0:714293de3836	1398	COMBA_FORWARD;
ashleymills	0:714293de3836	1399
ashleymills	0:714293de3836	1400	/* execute loop */
ashleymills	0:714293de3836	1401	for (iz = 0; iz < iy; iz++) {
ashleymills	0:714293de3836	1402	SQRADD2(tmpx++, tmpy--);
ashleymills	0:714293de3836	1403	}
ashleymills	0:714293de3836	1404
ashleymills	0:714293de3836	1405	/* even columns have the square term in them */
ashleymills	0:714293de3836	1406	if ((ix&1) == 0) {
ashleymills	0:714293de3836	1407	/* TAO change COMBA_ADD back to SQRADD */
ashleymills	0:714293de3836	1408	SQRADD(A->dp[ix>>1], A->dp[ix>>1]);
ashleymills	0:714293de3836	1409	}
ashleymills	0:714293de3836	1410
ashleymills	0:714293de3836	1411	/* store it */
ashleymills	0:714293de3836	1412	COMBA_STORE(dst->dp[ix]);
ashleymills	0:714293de3836	1413	}
ashleymills	0:714293de3836	1414
ashleymills	0:714293de3836	1415	COMBA_FINI;
ashleymills	0:714293de3836	1416
ashleymills	0:714293de3836	1417	/* setup dest */
ashleymills	0:714293de3836	1418	dst->used = pa;
ashleymills	0:714293de3836	1419	fp_clamp (dst);
ashleymills	0:714293de3836	1420	if (dst != B) {
ashleymills	0:714293de3836	1421	fp_copy(dst, B);
ashleymills	0:714293de3836	1422	}
ashleymills	0:714293de3836	1423	}
ashleymills	0:714293de3836	1424
ashleymills	0:714293de3836	1425	int fp_cmp(fp_int a, fp_int b)
ashleymills	0:714293de3836	1426	{
ashleymills	0:714293de3836	1427	if (a->sign == FP_NEG && b->sign == FP_ZPOS) {
ashleymills	0:714293de3836	1428	return FP_LT;
ashleymills	0:714293de3836	1429	} else if (a->sign == FP_ZPOS && b->sign == FP_NEG) {
ashleymills	0:714293de3836	1430	return FP_GT;
ashleymills	0:714293de3836	1431	} else {
ashleymills	0:714293de3836	1432	/* compare digits */
ashleymills	0:714293de3836	1433	if (a->sign == FP_NEG) {
ashleymills	0:714293de3836	1434	/* if negative compare opposite direction */
ashleymills	0:714293de3836	1435	return fp_cmp_mag(b, a);
ashleymills	0:714293de3836	1436	} else {
ashleymills	0:714293de3836	1437	return fp_cmp_mag(a, b);
ashleymills	0:714293de3836	1438	}
ashleymills	0:714293de3836	1439	}
ashleymills	0:714293de3836	1440	}
ashleymills	0:714293de3836	1441
ashleymills	0:714293de3836	1442	/* compare against a single digit */
ashleymills	0:714293de3836	1443	int fp_cmp_d(fp_int *a, fp_digit b)
ashleymills	0:714293de3836	1444	{
ashleymills	0:714293de3836	1445	/* compare based on sign */
ashleymills	0:714293de3836	1446	if ((b && a->used == 0) \|\| a->sign == FP_NEG) {
ashleymills	0:714293de3836	1447	return FP_LT;
ashleymills	0:714293de3836	1448	}
ashleymills	0:714293de3836	1449
ashleymills	0:714293de3836	1450	/* compare based on magnitude */
ashleymills	0:714293de3836	1451	if (a->used > 1) {
ashleymills	0:714293de3836	1452	return FP_GT;
ashleymills	0:714293de3836	1453	}
ashleymills	0:714293de3836	1454
ashleymills	0:714293de3836	1455	/* compare the only digit of a to b */
ashleymills	0:714293de3836	1456	if (a->dp[0] > b) {
ashleymills	0:714293de3836	1457	return FP_GT;
ashleymills	0:714293de3836	1458	} else if (a->dp[0] < b) {
ashleymills	0:714293de3836	1459	return FP_LT;
ashleymills	0:714293de3836	1460	} else {
ashleymills	0:714293de3836	1461	return FP_EQ;
ashleymills	0:714293de3836	1462	}
ashleymills	0:714293de3836	1463
ashleymills	0:714293de3836	1464	}
ashleymills	0:714293de3836	1465
ashleymills	0:714293de3836	1466	int fp_cmp_mag(fp_int a, fp_int b)
ashleymills	0:714293de3836	1467	{
ashleymills	0:714293de3836	1468	int x;
ashleymills	0:714293de3836	1469
ashleymills	0:714293de3836	1470	if (a->used > b->used) {
ashleymills	0:714293de3836	1471	return FP_GT;
ashleymills	0:714293de3836	1472	} else if (a->used < b->used) {
ashleymills	0:714293de3836	1473	return FP_LT;
ashleymills	0:714293de3836	1474	} else {
ashleymills	0:714293de3836	1475	for (x = a->used - 1; x >= 0; x--) {
ashleymills	0:714293de3836	1476	if (a->dp[x] > b->dp[x]) {
ashleymills	0:714293de3836	1477	return FP_GT;
ashleymills	0:714293de3836	1478	} else if (a->dp[x] < b->dp[x]) {
ashleymills	0:714293de3836	1479	return FP_LT;
ashleymills	0:714293de3836	1480	}
ashleymills	0:714293de3836	1481	}
ashleymills	0:714293de3836	1482	}
ashleymills	0:714293de3836	1483	return FP_EQ;
ashleymills	0:714293de3836	1484	}
ashleymills	0:714293de3836	1485
ashleymills	0:714293de3836	1486	/* setups the montgomery reduction */
ashleymills	0:714293de3836	1487	int fp_montgomery_setup(fp_int a, fp_digit rho)
ashleymills	0:714293de3836	1488	{
ashleymills	0:714293de3836	1489	fp_digit x, b;
ashleymills	0:714293de3836	1490
ashleymills	0:714293de3836	1491	/* fast inversion mod 2**k
ashleymills	0:714293de3836	1492	*
ashleymills	0:714293de3836	1493	* Based on the fact that
ashleymills	0:714293de3836	1494	*
ashleymills	0:714293de3836	1495	* XA = 1 (mod 2n) => (X(2-XA)) A = 1 (mod 22n)
ashleymills	0:714293de3836	1496	* => 2XA - XXA*A = 1
ashleymills	0:714293de3836	1497	* => 2*(1) - (1) = 1
ashleymills	0:714293de3836	1498	*/
ashleymills	0:714293de3836	1499	b = a->dp[0];
ashleymills	0:714293de3836	1500
ashleymills	0:714293de3836	1501	if ((b & 1) == 0) {
ashleymills	0:714293de3836	1502	return FP_VAL;
ashleymills	0:714293de3836	1503	}
ashleymills	0:714293de3836	1504
ashleymills	0:714293de3836	1505	x = (((b + 2) & 4) << 1) + b; /* here xa==1 mod 24 /
ashleymills	0:714293de3836	1506	x = 2 - b x; /* here xa==1 mod 28 /
ashleymills	0:714293de3836	1507	x = 2 - b x; /* here xa==1 mod 216 /
ashleymills	0:714293de3836	1508	x = 2 - b x; /* here xa==1 mod 232 /
ashleymills	0:714293de3836	1509	#ifdef FP_64BIT
ashleymills	0:714293de3836	1510	x = 2 - b x; /* here xa==1 mod 264 /
ashleymills	0:714293de3836	1511	#endif
ashleymills	0:714293de3836	1512
ashleymills	0:714293de3836	1513	/* rho = -1/m mod b */
ashleymills	0:714293de3836	1514	*rho = (fp_digit) (((fp_word) 1 << ((fp_word) DIGIT_BIT)) - ((fp_word)x));
ashleymills	0:714293de3836	1515
ashleymills	0:714293de3836	1516	return FP_OKAY;
ashleymills	0:714293de3836	1517	}
ashleymills	0:714293de3836	1518
ashleymills	0:714293de3836	1519	/* computes a = B**n mod b without division or multiplication useful for
ashleymills	0:714293de3836	1520	* normalizing numbers in a Montgomery system.
ashleymills	0:714293de3836	1521	*/
ashleymills	0:714293de3836	1522	void fp_montgomery_calc_normalization(fp_int a, fp_int b)
ashleymills	0:714293de3836	1523	{
ashleymills	0:714293de3836	1524	int x, bits;
ashleymills	0:714293de3836	1525
ashleymills	0:714293de3836	1526	/* how many bits of last digit does b use */
ashleymills	0:714293de3836	1527	bits = fp_count_bits (b) % DIGIT_BIT;
ashleymills	0:714293de3836	1528	if (!bits) bits = DIGIT_BIT;
ashleymills	0:714293de3836	1529
ashleymills	0:714293de3836	1530	/* compute A = B^(n-1) * 2^(bits-1) */
ashleymills	0:714293de3836	1531	if (b->used > 1) {
ashleymills	0:714293de3836	1532	fp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1);
ashleymills	0:714293de3836	1533	} else {
ashleymills	0:714293de3836	1534	fp_set(a, 1);
ashleymills	0:714293de3836	1535	bits = 1;
ashleymills	0:714293de3836	1536	}
ashleymills	0:714293de3836	1537
ashleymills	0:714293de3836	1538	/* now compute C = A * B mod b */
ashleymills	0:714293de3836	1539	for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
ashleymills	0:714293de3836	1540	fp_mul_2 (a, a);
ashleymills	0:714293de3836	1541	if (fp_cmp_mag (a, b) != FP_LT) {
ashleymills	0:714293de3836	1542	s_fp_sub (a, b, a);
ashleymills	0:714293de3836	1543	}
ashleymills	0:714293de3836	1544	}
ashleymills	0:714293de3836	1545	}
ashleymills	0:714293de3836	1546
ashleymills	0:714293de3836	1547
ashleymills	0:714293de3836	1548	#ifdef TFM_SMALL_MONT_SET
ashleymills	0:714293de3836	1549	#include "fp_mont_small.i"
ashleymills	0:714293de3836	1550	#endif
ashleymills	0:714293de3836	1551
ashleymills	0:714293de3836	1552	/* computes x/R == x (mod N) via Montgomery Reduction */
ashleymills	0:714293de3836	1553	void fp_montgomery_reduce(fp_int a, fp_int m, fp_digit mp)
ashleymills	0:714293de3836	1554	{
ashleymills	0:714293de3836	1555	fp_digit c[FP_SIZE], _c, tmpm, mu = 0;
ashleymills	0:714293de3836	1556	int oldused, x, y, pa;
ashleymills	0:714293de3836	1557
ashleymills	0:714293de3836	1558	/* bail if too large */
ashleymills	0:714293de3836	1559	if (m->used > (FP_SIZE/2)) {
ashleymills	0:714293de3836	1560	(void)mu; /* shut up compiler */
ashleymills	0:714293de3836	1561	return;
ashleymills	0:714293de3836	1562	}
ashleymills	0:714293de3836	1563
ashleymills	0:714293de3836	1564	#ifdef TFM_SMALL_MONT_SET
ashleymills	0:714293de3836	1565	if (m->used <= 16) {
ashleymills	0:714293de3836	1566	fp_montgomery_reduce_small(a, m, mp);
ashleymills	0:714293de3836	1567	return;
ashleymills	0:714293de3836	1568	}
ashleymills	0:714293de3836	1569	#endif
ashleymills	0:714293de3836	1570
ashleymills	0:714293de3836	1571
ashleymills	0:714293de3836	1572	/* now zero the buff */
ashleymills	0:714293de3836	1573	XMEMSET(c, 0, sizeof c);
ashleymills	0:714293de3836	1574	pa = m->used;
ashleymills	0:714293de3836	1575
ashleymills	0:714293de3836	1576	/* copy the input */
ashleymills	0:714293de3836	1577	oldused = a->used;
ashleymills	0:714293de3836	1578	for (x = 0; x < oldused; x++) {
ashleymills	0:714293de3836	1579	c[x] = a->dp[x];
ashleymills	0:714293de3836	1580	}
ashleymills	0:714293de3836	1581	MONT_START;
ashleymills	0:714293de3836	1582
ashleymills	0:714293de3836	1583	for (x = 0; x < pa; x++) {
ashleymills	0:714293de3836	1584	fp_digit cy = 0;
ashleymills	0:714293de3836	1585	/* get Mu for this round */
ashleymills	0:714293de3836	1586	LOOP_START;
ashleymills	0:714293de3836	1587	_c = c + x;
ashleymills	0:714293de3836	1588	tmpm = m->dp;
ashleymills	0:714293de3836	1589	y = 0;
ashleymills	0:714293de3836	1590	#if (defined(TFM_SSE2) \|\| defined(TFM_X86_64))
ashleymills	0:714293de3836	1591	for (; y < (pa & ~7); y += 8) {
ashleymills	0:714293de3836	1592	INNERMUL8;
ashleymills	0:714293de3836	1593	_c += 8;
ashleymills	0:714293de3836	1594	tmpm += 8;
ashleymills	0:714293de3836	1595	}
ashleymills	0:714293de3836	1596	#endif
ashleymills	0:714293de3836	1597
ashleymills	0:714293de3836	1598	for (; y < pa; y++) {
ashleymills	0:714293de3836	1599	INNERMUL;
ashleymills	0:714293de3836	1600	++_c;
ashleymills	0:714293de3836	1601	}
ashleymills	0:714293de3836	1602	LOOP_END;
ashleymills	0:714293de3836	1603	while (cy) {
ashleymills	0:714293de3836	1604	PROPCARRY;
ashleymills	0:714293de3836	1605	++_c;
ashleymills	0:714293de3836	1606	}
ashleymills	0:714293de3836	1607	}
ashleymills	0:714293de3836	1608
ashleymills	0:714293de3836	1609	/* now copy out */
ashleymills	0:714293de3836	1610	_c = c + pa;
ashleymills	0:714293de3836	1611	tmpm = a->dp;
ashleymills	0:714293de3836	1612	for (x = 0; x < pa+1; x++) {
ashleymills	0:714293de3836	1613	tmpm++ = _c++;
ashleymills	0:714293de3836	1614	}
ashleymills	0:714293de3836	1615
ashleymills	0:714293de3836	1616	for (; x < oldused; x++) {
ashleymills	0:714293de3836	1617	*tmpm++ = 0;
ashleymills	0:714293de3836	1618	}
ashleymills	0:714293de3836	1619
ashleymills	0:714293de3836	1620	MONT_FINI;
ashleymills	0:714293de3836	1621
ashleymills	0:714293de3836	1622	a->used = pa+1;
ashleymills	0:714293de3836	1623	fp_clamp(a);
ashleymills	0:714293de3836	1624
ashleymills	0:714293de3836	1625	/* if A >= m then A = A - m */
ashleymills	0:714293de3836	1626	if (fp_cmp_mag (a, m) != FP_LT) {
ashleymills	0:714293de3836	1627	s_fp_sub (a, m, a);
ashleymills	0:714293de3836	1628	}
ashleymills	0:714293de3836	1629	}
ashleymills	0:714293de3836	1630
ashleymills	0:714293de3836	1631	void fp_read_unsigned_bin(fp_int a, unsigned char b, int c)
ashleymills	0:714293de3836	1632	{
ashleymills	0:714293de3836	1633	/* zero the int */
ashleymills	0:714293de3836	1634	fp_zero (a);
ashleymills	0:714293de3836	1635
ashleymills	0:714293de3836	1636	/* If we know the endianness of this architecture, and we're using
ashleymills	0:714293de3836	1637	32-bit fp_digits, we can optimize this */
ashleymills	0:714293de3836	1638	#if (defined(LITTLE_ENDIAN_ORDER) \|\| defined(BIG_ENDIAN_ORDER)) && !defined(FP_64BIT)
ashleymills	0:714293de3836	1639	/* But not for both simultaneously */
ashleymills	0:714293de3836	1640	#if defined(LITTLE_ENDIAN_ORDER) && defined(BIG_ENDIAN_ORDER)
ashleymills	0:714293de3836	1641	#error Both LITTLE_ENDIAN_ORDER and BIG_ENDIAN_ORDER defined.
ashleymills	0:714293de3836	1642	#endif
ashleymills	0:714293de3836	1643	{
ashleymills	0:714293de3836	1644	unsigned char pd = (unsigned char )a->dp;
ashleymills	0:714293de3836	1645
ashleymills	0:714293de3836	1646	if ((unsigned)c > (FP_SIZE * sizeof(fp_digit))) {
ashleymills	0:714293de3836	1647	int excess = c - (FP_SIZE * sizeof(fp_digit));
ashleymills	0:714293de3836	1648	c -= excess;
ashleymills	0:714293de3836	1649	b += excess;
ashleymills	0:714293de3836	1650	}
ashleymills	0:714293de3836	1651	a->used = (c + sizeof(fp_digit) - 1)/sizeof(fp_digit);
ashleymills	0:714293de3836	1652	/* read the bytes in */
ashleymills	0:714293de3836	1653	#ifdef BIG_ENDIAN_ORDER
ashleymills	0:714293de3836	1654	{
ashleymills	0:714293de3836	1655	/* Use Duff's device to unroll the loop. */
ashleymills	0:714293de3836	1656	int idx = (c - 1) & ~3;
ashleymills	0:714293de3836	1657	switch (c % 4) {
ashleymills	0:714293de3836	1658	case 0: do { pd[idx+0] = *b++;
ashleymills	0:714293de3836	1659	case 3: pd[idx+1] = *b++;
ashleymills	0:714293de3836	1660	case 2: pd[idx+2] = *b++;
ashleymills	0:714293de3836	1661	case 1: pd[idx+3] = *b++;
ashleymills	0:714293de3836	1662	idx -= 4;
ashleymills	0:714293de3836	1663	} while ((c -= 4) > 0);
ashleymills	0:714293de3836	1664	}
ashleymills	0:714293de3836	1665	}
ashleymills	0:714293de3836	1666	#else
ashleymills	0:714293de3836	1667	for (c -= 1; c >= 0; c -= 1) {
ashleymills	0:714293de3836	1668	pd[c] = *b++;
ashleymills	0:714293de3836	1669	}
ashleymills	0:714293de3836	1670	#endif
ashleymills	0:714293de3836	1671	}
ashleymills	0:714293de3836	1672	#else
ashleymills	0:714293de3836	1673	/* read the bytes in */
ashleymills	0:714293de3836	1674	for (; c > 0; c--) {
ashleymills	0:714293de3836	1675	fp_mul_2d (a, 8, a);
ashleymills	0:714293de3836	1676	a->dp[0] \|= *b++;
ashleymills	0:714293de3836	1677	a->used += 1;
ashleymills	0:714293de3836	1678	}
ashleymills	0:714293de3836	1679	#endif
ashleymills	0:714293de3836	1680	fp_clamp (a);
ashleymills	0:714293de3836	1681	}
ashleymills	0:714293de3836	1682
ashleymills	0:714293de3836	1683	void fp_to_unsigned_bin(fp_int a, unsigned char b)
ashleymills	0:714293de3836	1684	{
ashleymills	0:714293de3836	1685	int x;
ashleymills	0:714293de3836	1686	fp_int t;
ashleymills	0:714293de3836	1687
ashleymills	0:714293de3836	1688	fp_init_copy(&t, a);
ashleymills	0:714293de3836	1689
ashleymills	0:714293de3836	1690	x = 0;
ashleymills	0:714293de3836	1691	while (fp_iszero (&t) == FP_NO) {
ashleymills	0:714293de3836	1692	b[x++] = (unsigned char) (t.dp[0] & 255);
ashleymills	0:714293de3836	1693	fp_div_2d (&t, 8, &t, NULL);
ashleymills	0:714293de3836	1694	}
ashleymills	0:714293de3836	1695	fp_reverse (b, x);
ashleymills	0:714293de3836	1696	}
ashleymills	0:714293de3836	1697
ashleymills	0:714293de3836	1698	int fp_unsigned_bin_size(fp_int *a)
ashleymills	0:714293de3836	1699	{
ashleymills	0:714293de3836	1700	int size = fp_count_bits (a);
ashleymills	0:714293de3836	1701	return (size / 8 + ((size & 7) != 0 ? 1 : 0));
ashleymills	0:714293de3836	1702	}
ashleymills	0:714293de3836	1703
ashleymills	0:714293de3836	1704	void fp_set(fp_int *a, fp_digit b)
ashleymills	0:714293de3836	1705	{
ashleymills	0:714293de3836	1706	fp_zero(a);
ashleymills	0:714293de3836	1707	a->dp[0] = b;
ashleymills	0:714293de3836	1708	a->used = a->dp[0] ? 1 : 0;
ashleymills	0:714293de3836	1709	}
ashleymills	0:714293de3836	1710
ashleymills	0:714293de3836	1711	int fp_count_bits (fp_int * a)
ashleymills	0:714293de3836	1712	{
ashleymills	0:714293de3836	1713	int r;
ashleymills	0:714293de3836	1714	fp_digit q;
ashleymills	0:714293de3836	1715
ashleymills	0:714293de3836	1716	/* shortcut */
ashleymills	0:714293de3836	1717	if (a->used == 0) {
ashleymills	0:714293de3836	1718	return 0;
ashleymills	0:714293de3836	1719	}
ashleymills	0:714293de3836	1720
ashleymills	0:714293de3836	1721	/* get number of digits and add that */
ashleymills	0:714293de3836	1722	r = (a->used - 1) * DIGIT_BIT;
ashleymills	0:714293de3836	1723
ashleymills	0:714293de3836	1724	/* take the last digit and count the bits in it */
ashleymills	0:714293de3836	1725	q = a->dp[a->used - 1];
ashleymills	0:714293de3836	1726	while (q > ((fp_digit) 0)) {
ashleymills	0:714293de3836	1727	++r;
ashleymills	0:714293de3836	1728	q >>= ((fp_digit) 1);
ashleymills	0:714293de3836	1729	}
ashleymills	0:714293de3836	1730	return r;
ashleymills	0:714293de3836	1731	}
ashleymills	0:714293de3836	1732
ashleymills	0:714293de3836	1733	void fp_lshd(fp_int *a, int x)
ashleymills	0:714293de3836	1734	{
ashleymills	0:714293de3836	1735	int y;
ashleymills	0:714293de3836	1736
ashleymills	0:714293de3836	1737	/* move up and truncate as required */
ashleymills	0:714293de3836	1738	y = MIN(a->used + x - 1, (int)(FP_SIZE-1));
ashleymills	0:714293de3836	1739
ashleymills	0:714293de3836	1740	/* store new size */
ashleymills	0:714293de3836	1741	a->used = y + 1;
ashleymills	0:714293de3836	1742
ashleymills	0:714293de3836	1743	/* move digits */
ashleymills	0:714293de3836	1744	for (; y >= x; y--) {
ashleymills	0:714293de3836	1745	a->dp[y] = a->dp[y-x];
ashleymills	0:714293de3836	1746	}
ashleymills	0:714293de3836	1747
ashleymills	0:714293de3836	1748	/* zero lower digits */
ashleymills	0:714293de3836	1749	for (; y >= 0; y--) {
ashleymills	0:714293de3836	1750	a->dp[y] = 0;
ashleymills	0:714293de3836	1751	}
ashleymills	0:714293de3836	1752
ashleymills	0:714293de3836	1753	/* clamp digits */
ashleymills	0:714293de3836	1754	fp_clamp(a);
ashleymills	0:714293de3836	1755	}
ashleymills	0:714293de3836	1756
ashleymills	0:714293de3836	1757	void fp_rshd(fp_int *a, int x)
ashleymills	0:714293de3836	1758	{
ashleymills	0:714293de3836	1759	int y;
ashleymills	0:714293de3836	1760
ashleymills	0:714293de3836	1761	/* too many digits just zero and return */
ashleymills	0:714293de3836	1762	if (x >= a->used) {
ashleymills	0:714293de3836	1763	fp_zero(a);
ashleymills	0:714293de3836	1764	return;
ashleymills	0:714293de3836	1765	}
ashleymills	0:714293de3836	1766
ashleymills	0:714293de3836	1767	/* shift */
ashleymills	0:714293de3836	1768	for (y = 0; y < a->used - x; y++) {
ashleymills	0:714293de3836	1769	a->dp[y] = a->dp[y+x];
ashleymills	0:714293de3836	1770	}
ashleymills	0:714293de3836	1771
ashleymills	0:714293de3836	1772	/* zero rest */
ashleymills	0:714293de3836	1773	for (; y < a->used; y++) {
ashleymills	0:714293de3836	1774	a->dp[y] = 0;
ashleymills	0:714293de3836	1775	}
ashleymills	0:714293de3836	1776
ashleymills	0:714293de3836	1777	/* decrement count */
ashleymills	0:714293de3836	1778	a->used -= x;
ashleymills	0:714293de3836	1779	fp_clamp(a);
ashleymills	0:714293de3836	1780	}
ashleymills	0:714293de3836	1781
ashleymills	0:714293de3836	1782	/* reverse an array, used for radix code */
ashleymills	0:714293de3836	1783	void fp_reverse (unsigned char *s, int len)
ashleymills	0:714293de3836	1784	{
ashleymills	0:714293de3836	1785	int ix, iy;
ashleymills	0:714293de3836	1786	unsigned char t;
ashleymills	0:714293de3836	1787
ashleymills	0:714293de3836	1788	ix = 0;
ashleymills	0:714293de3836	1789	iy = len - 1;
ashleymills	0:714293de3836	1790	while (ix < iy) {
ashleymills	0:714293de3836	1791	t = s[ix];
ashleymills	0:714293de3836	1792	s[ix] = s[iy];
ashleymills	0:714293de3836	1793	s[iy] = t;
ashleymills	0:714293de3836	1794	++ix;
ashleymills	0:714293de3836	1795	--iy;
ashleymills	0:714293de3836	1796	}
ashleymills	0:714293de3836	1797	}
ashleymills	0:714293de3836	1798
ashleymills	0:714293de3836	1799
ashleymills	0:714293de3836	1800	/* c = a - b */
ashleymills	0:714293de3836	1801	void fp_sub_d(fp_int a, fp_digit b, fp_int c)
ashleymills	0:714293de3836	1802	{
ashleymills	0:714293de3836	1803	fp_int tmp;
ashleymills	0:714293de3836	1804	fp_set(&tmp, b);
ashleymills	0:714293de3836	1805	fp_sub(a, &tmp, c);
ashleymills	0:714293de3836	1806	}
ashleymills	0:714293de3836	1807
ashleymills	0:714293de3836	1808
ashleymills	0:714293de3836	1809	/* CyaSSL callers from normal lib */
ashleymills	0:714293de3836	1810
ashleymills	0:714293de3836	1811	/* init a new mp_int */
ashleymills	0:714293de3836	1812	int mp_init (mp_int * a)
ashleymills	0:714293de3836	1813	{
ashleymills	0:714293de3836	1814	if (a)
ashleymills	0:714293de3836	1815	fp_init(a);
ashleymills	0:714293de3836	1816	return MP_OKAY;
ashleymills	0:714293de3836	1817	}
ashleymills	0:714293de3836	1818
ashleymills	0:714293de3836	1819	/* clear one (frees) */
ashleymills	0:714293de3836	1820	void mp_clear (mp_int * a)
ashleymills	0:714293de3836	1821	{
ashleymills	0:714293de3836	1822	fp_zero(a);
ashleymills	0:714293de3836	1823	}
ashleymills	0:714293de3836	1824
ashleymills	0:714293de3836	1825	/* handle up to 6 inits */
ashleymills	0:714293de3836	1826	int mp_init_multi(mp_int* a, mp_int* b, mp_int* c, mp_int* d, mp_int* e, mp_int* f)
ashleymills	0:714293de3836	1827	{
ashleymills	0:714293de3836	1828	if (a)
ashleymills	0:714293de3836	1829	fp_init(a);
ashleymills	0:714293de3836	1830	if (b)
ashleymills	0:714293de3836	1831	fp_init(b);
ashleymills	0:714293de3836	1832	if (c)
ashleymills	0:714293de3836	1833	fp_init(c);
ashleymills	0:714293de3836	1834	if (d)
ashleymills	0:714293de3836	1835	fp_init(d);
ashleymills	0:714293de3836	1836	if (e)
ashleymills	0:714293de3836	1837	fp_init(e);
ashleymills	0:714293de3836	1838	if (f)
ashleymills	0:714293de3836	1839	fp_init(f);
ashleymills	0:714293de3836	1840
ashleymills	0:714293de3836	1841	return MP_OKAY;
ashleymills	0:714293de3836	1842	}
ashleymills	0:714293de3836	1843
ashleymills	0:714293de3836	1844	/* high level addition (handles signs) */
ashleymills	0:714293de3836	1845	int mp_add (mp_int * a, mp_int * b, mp_int * c)
ashleymills	0:714293de3836	1846	{
ashleymills	0:714293de3836	1847	fp_add(a, b, c);
ashleymills	0:714293de3836	1848	return MP_OKAY;
ashleymills	0:714293de3836	1849	}
ashleymills	0:714293de3836	1850
ashleymills	0:714293de3836	1851	/* high level subtraction (handles signs) */
ashleymills	0:714293de3836	1852	int mp_sub (mp_int * a, mp_int * b, mp_int * c)
ashleymills	0:714293de3836	1853	{
ashleymills	0:714293de3836	1854	fp_sub(a, b, c);
ashleymills	0:714293de3836	1855	return MP_OKAY;
ashleymills	0:714293de3836	1856	}
ashleymills	0:714293de3836	1857
ashleymills	0:714293de3836	1858	/* high level multiplication (handles sign) */
ashleymills	0:714293de3836	1859	int mp_mul (mp_int * a, mp_int * b, mp_int * c)
ashleymills	0:714293de3836	1860	{
ashleymills	0:714293de3836	1861	fp_mul(a, b, c);
ashleymills	0:714293de3836	1862	return MP_OKAY;
ashleymills	0:714293de3836	1863	}
ashleymills	0:714293de3836	1864
ashleymills	0:714293de3836	1865	/* d = a * b (mod c) */
ashleymills	0:714293de3836	1866	int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
ashleymills	0:714293de3836	1867	{
ashleymills	0:714293de3836	1868	return fp_mulmod(a, b, c, d);
ashleymills	0:714293de3836	1869	}
ashleymills	0:714293de3836	1870
ashleymills	0:714293de3836	1871	/* c = a mod b, 0 <= c < b */
ashleymills	0:714293de3836	1872	int mp_mod (mp_int * a, mp_int * b, mp_int * c)
ashleymills	0:714293de3836	1873	{
ashleymills	0:714293de3836	1874	return fp_mod (a, b, c);
ashleymills	0:714293de3836	1875	}
ashleymills	0:714293de3836	1876
ashleymills	0:714293de3836	1877	/* hac 14.61, pp608 */
ashleymills	0:714293de3836	1878	int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
ashleymills	0:714293de3836	1879	{
ashleymills	0:714293de3836	1880	return fp_invmod(a, b, c);
ashleymills	0:714293de3836	1881	}
ashleymills	0:714293de3836	1882
ashleymills	0:714293de3836	1883	/* this is a shell function that calls either the normal or Montgomery
ashleymills	0:714293de3836	1884	* exptmod functions. Originally the call to the montgomery code was
ashleymills	0:714293de3836	1885	* embedded in the normal function but that wasted alot of stack space
ashleymills	0:714293de3836	1886	* for nothing (since 99% of the time the Montgomery code would be called)
ashleymills	0:714293de3836	1887	*/
ashleymills	0:714293de3836	1888	int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
ashleymills	0:714293de3836	1889	{
ashleymills	0:714293de3836	1890	return fp_exptmod(G, X, P, Y);
ashleymills	0:714293de3836	1891	}
ashleymills	0:714293de3836	1892
ashleymills	0:714293de3836	1893	/* compare two ints (signed)*/
ashleymills	0:714293de3836	1894	int mp_cmp (mp_int * a, mp_int * b)
ashleymills	0:714293de3836	1895	{
ashleymills	0:714293de3836	1896	return fp_cmp(a, b);
ashleymills	0:714293de3836	1897	}
ashleymills	0:714293de3836	1898
ashleymills	0:714293de3836	1899	/* compare a digit */
ashleymills	0:714293de3836	1900	int mp_cmp_d(mp_int * a, mp_digit b)
ashleymills	0:714293de3836	1901	{
ashleymills	0:714293de3836	1902	return fp_cmp_d(a, b);
ashleymills	0:714293de3836	1903	}
ashleymills	0:714293de3836	1904
ashleymills	0:714293de3836	1905	/* get the size for an unsigned equivalent */
ashleymills	0:714293de3836	1906	int mp_unsigned_bin_size (mp_int * a)
ashleymills	0:714293de3836	1907	{
ashleymills	0:714293de3836	1908	return fp_unsigned_bin_size(a);
ashleymills	0:714293de3836	1909	}
ashleymills	0:714293de3836	1910
ashleymills	0:714293de3836	1911	/* store in unsigned [big endian] format */
ashleymills	0:714293de3836	1912	int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
ashleymills	0:714293de3836	1913	{
ashleymills	0:714293de3836	1914	fp_to_unsigned_bin(a,b);
ashleymills	0:714293de3836	1915	return MP_OKAY;
ashleymills	0:714293de3836	1916	}
ashleymills	0:714293de3836	1917
ashleymills	0:714293de3836	1918	/* reads a unsigned char array, assumes the msb is stored first [big endian] */
ashleymills	0:714293de3836	1919	int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
ashleymills	0:714293de3836	1920	{
ashleymills	0:714293de3836	1921	fp_read_unsigned_bin(a, (unsigned char *)b, c);
ashleymills	0:714293de3836	1922	return MP_OKAY;
ashleymills	0:714293de3836	1923	}
ashleymills	0:714293de3836	1924
ashleymills	0:714293de3836	1925
ashleymills	0:714293de3836	1926	int mp_sub_d(fp_int a, fp_digit b, fp_int c)
ashleymills	0:714293de3836	1927	{
ashleymills	0:714293de3836	1928	fp_sub_d(a, b, c);
ashleymills	0:714293de3836	1929	return MP_OKAY;
ashleymills	0:714293de3836	1930	}
ashleymills	0:714293de3836	1931
ashleymills	0:714293de3836	1932
ashleymills	0:714293de3836	1933	/* fast math conversion */
ashleymills	0:714293de3836	1934	int mp_copy(fp_int* a, fp_int* b)
ashleymills	0:714293de3836	1935	{
ashleymills	0:714293de3836	1936	fp_copy(a, b);
ashleymills	0:714293de3836	1937	return MP_OKAY;
ashleymills	0:714293de3836	1938	}
ashleymills	0:714293de3836	1939
ashleymills	0:714293de3836	1940
ashleymills	0:714293de3836	1941	/* fast math conversion */
ashleymills	0:714293de3836	1942	int mp_isodd(mp_int* a)
ashleymills	0:714293de3836	1943	{
ashleymills	0:714293de3836	1944	return fp_isodd(a);
ashleymills	0:714293de3836	1945	}
ashleymills	0:714293de3836	1946
ashleymills	0:714293de3836	1947
ashleymills	0:714293de3836	1948	/* fast math conversion */
ashleymills	0:714293de3836	1949	int mp_iszero(mp_int* a)
ashleymills	0:714293de3836	1950	{
ashleymills	0:714293de3836	1951	return fp_iszero(a);
ashleymills	0:714293de3836	1952	}
ashleymills	0:714293de3836	1953
ashleymills	0:714293de3836	1954
ashleymills	0:714293de3836	1955	/* fast math conversion */
ashleymills	0:714293de3836	1956	int mp_count_bits (mp_int* a)
ashleymills	0:714293de3836	1957	{
ashleymills	0:714293de3836	1958	return fp_count_bits(a);
ashleymills	0:714293de3836	1959	}
ashleymills	0:714293de3836	1960
ashleymills	0:714293de3836	1961
ashleymills	0:714293de3836	1962	/* fast math wrappers */
ashleymills	0:714293de3836	1963	int mp_set_int(fp_int *a, fp_digit b)
ashleymills	0:714293de3836	1964	{
ashleymills	0:714293de3836	1965	fp_set(a, b);
ashleymills	0:714293de3836	1966	return MP_OKAY;
ashleymills	0:714293de3836	1967	}
ashleymills	0:714293de3836	1968
ashleymills	0:714293de3836	1969
ashleymills	0:714293de3836	1970	#if defined(CYASSL_KEY_GEN) \|\| defined (HAVE_ECC)
ashleymills	0:714293de3836	1971
ashleymills	0:714293de3836	1972	/* c = a * a (mod b) */
ashleymills	0:714293de3836	1973	int fp_sqrmod(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	1974	{
ashleymills	0:714293de3836	1975	fp_int tmp;
ashleymills	0:714293de3836	1976	fp_zero(&tmp);
ashleymills	0:714293de3836	1977	fp_sqr(a, &tmp);
ashleymills	0:714293de3836	1978	return fp_mod(&tmp, b, c);
ashleymills	0:714293de3836	1979	}
ashleymills	0:714293de3836	1980
ashleymills	0:714293de3836	1981	/* fast math conversion */
ashleymills	0:714293de3836	1982	int mp_sqrmod(mp_int a, mp_int b, mp_int *c)
ashleymills	0:714293de3836	1983	{
ashleymills	0:714293de3836	1984	return fp_sqrmod(a, b, c);
ashleymills	0:714293de3836	1985	}
ashleymills	0:714293de3836	1986
ashleymills	0:714293de3836	1987	/* fast math conversion */
ashleymills	0:714293de3836	1988	int mp_montgomery_calc_normalization(mp_int a, mp_int b)
ashleymills	0:714293de3836	1989	{
ashleymills	0:714293de3836	1990	fp_montgomery_calc_normalization(a, b);
ashleymills	0:714293de3836	1991	return MP_OKAY;
ashleymills	0:714293de3836	1992	}
ashleymills	0:714293de3836	1993
ashleymills	0:714293de3836	1994	#endif /* CYASSL_KEYGEN \|\| HAVE_ECC */
ashleymills	0:714293de3836	1995
ashleymills	0:714293de3836	1996
ashleymills	0:714293de3836	1997	#ifdef CYASSL_KEY_GEN
ashleymills	0:714293de3836	1998
ashleymills	0:714293de3836	1999	void fp_gcd(fp_int a, fp_int b, fp_int *c);
ashleymills	0:714293de3836	2000	void fp_lcm(fp_int a, fp_int b, fp_int *c);
ashleymills	0:714293de3836	2001	int fp_isprime(fp_int *a);
ashleymills	0:714293de3836	2002	int fp_cnt_lsb(fp_int *a);
ashleymills	0:714293de3836	2003
ashleymills	0:714293de3836	2004	int mp_gcd(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	2005	{
ashleymills	0:714293de3836	2006	fp_gcd(a, b, c);
ashleymills	0:714293de3836	2007	return MP_OKAY;
ashleymills	0:714293de3836	2008	}
ashleymills	0:714293de3836	2009
ashleymills	0:714293de3836	2010
ashleymills	0:714293de3836	2011	int mp_lcm(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	2012	{
ashleymills	0:714293de3836	2013	fp_lcm(a, b, c);
ashleymills	0:714293de3836	2014	return MP_OKAY;
ashleymills	0:714293de3836	2015	}
ashleymills	0:714293de3836	2016
ashleymills	0:714293de3836	2017
ashleymills	0:714293de3836	2018	int mp_prime_is_prime(mp_int* a, int t, int* result)
ashleymills	0:714293de3836	2019	{
ashleymills	0:714293de3836	2020	(void)t;
ashleymills	0:714293de3836	2021	*result = fp_isprime(a);
ashleymills	0:714293de3836	2022	return MP_OKAY;
ashleymills	0:714293de3836	2023	}
ashleymills	0:714293de3836	2024
ashleymills	0:714293de3836	2025
ashleymills	0:714293de3836	2026
ashleymills	0:714293de3836	2027	static int s_is_power_of_two(fp_digit b, int *p)
ashleymills	0:714293de3836	2028	{
ashleymills	0:714293de3836	2029	int x;
ashleymills	0:714293de3836	2030
ashleymills	0:714293de3836	2031	/* fast return if no power of two */
ashleymills	0:714293de3836	2032	if ((b==0) \|\| (b & (b-1))) {
ashleymills	0:714293de3836	2033	return 0;
ashleymills	0:714293de3836	2034	}
ashleymills	0:714293de3836	2035
ashleymills	0:714293de3836	2036	for (x = 0; x < DIGIT_BIT; x++) {
ashleymills	0:714293de3836	2037	if (b == (((fp_digit)1)<<x)) {
ashleymills	0:714293de3836	2038	*p = x;
ashleymills	0:714293de3836	2039	return 1;
ashleymills	0:714293de3836	2040	}
ashleymills	0:714293de3836	2041	}
ashleymills	0:714293de3836	2042	return 0;
ashleymills	0:714293de3836	2043	}
ashleymills	0:714293de3836	2044
ashleymills	0:714293de3836	2045	/* a/b => cb + d == a */
ashleymills	0:714293de3836	2046	static int fp_div_d(fp_int a, fp_digit b, fp_int c, fp_digit *d)
ashleymills	0:714293de3836	2047	{
ashleymills	0:714293de3836	2048	fp_int q;
ashleymills	0:714293de3836	2049	fp_word w;
ashleymills	0:714293de3836	2050	fp_digit t;
ashleymills	0:714293de3836	2051	int ix;
ashleymills	0:714293de3836	2052
ashleymills	0:714293de3836	2053	/* cannot divide by zero */
ashleymills	0:714293de3836	2054	if (b == 0) {
ashleymills	0:714293de3836	2055	return FP_VAL;
ashleymills	0:714293de3836	2056	}
ashleymills	0:714293de3836	2057
ashleymills	0:714293de3836	2058	/* quick outs */
ashleymills	0:714293de3836	2059	if (b == 1 \|\| fp_iszero(a) == 1) {
ashleymills	0:714293de3836	2060	if (d != NULL) {
ashleymills	0:714293de3836	2061	*d = 0;
ashleymills	0:714293de3836	2062	}
ashleymills	0:714293de3836	2063	if (c != NULL) {
ashleymills	0:714293de3836	2064	fp_copy(a, c);
ashleymills	0:714293de3836	2065	}
ashleymills	0:714293de3836	2066	return FP_OKAY;
ashleymills	0:714293de3836	2067	}
ashleymills	0:714293de3836	2068
ashleymills	0:714293de3836	2069	/* power of two ? */
ashleymills	0:714293de3836	2070	if (s_is_power_of_two(b, &ix) == 1) {
ashleymills	0:714293de3836	2071	if (d != NULL) {
ashleymills	0:714293de3836	2072	*d = a->dp[0] & ((((fp_digit)1)<<ix) - 1);
ashleymills	0:714293de3836	2073	}
ashleymills	0:714293de3836	2074	if (c != NULL) {
ashleymills	0:714293de3836	2075	fp_div_2d(a, ix, c, NULL);
ashleymills	0:714293de3836	2076	}
ashleymills	0:714293de3836	2077	return FP_OKAY;
ashleymills	0:714293de3836	2078	}
ashleymills	0:714293de3836	2079
ashleymills	0:714293de3836	2080	/* no easy answer [c'est la vie]. Just division */
ashleymills	0:714293de3836	2081	fp_init(&q);
ashleymills	0:714293de3836	2082
ashleymills	0:714293de3836	2083	q.used = a->used;
ashleymills	0:714293de3836	2084	q.sign = a->sign;
ashleymills	0:714293de3836	2085	w = 0;
ashleymills	0:714293de3836	2086	for (ix = a->used - 1; ix >= 0; ix--) {
ashleymills	0:714293de3836	2087	w = (w << ((fp_word)DIGIT_BIT)) \| ((fp_word)a->dp[ix]);
ashleymills	0:714293de3836	2088
ashleymills	0:714293de3836	2089	if (w >= b) {
ashleymills	0:714293de3836	2090	t = (fp_digit)(w / b);
ashleymills	0:714293de3836	2091	w -= ((fp_word)t) * ((fp_word)b);
ashleymills	0:714293de3836	2092	} else {
ashleymills	0:714293de3836	2093	t = 0;
ashleymills	0:714293de3836	2094	}
ashleymills	0:714293de3836	2095	q.dp[ix] = (fp_digit)t;
ashleymills	0:714293de3836	2096	}
ashleymills	0:714293de3836	2097
ashleymills	0:714293de3836	2098	if (d != NULL) {
ashleymills	0:714293de3836	2099	*d = (fp_digit)w;
ashleymills	0:714293de3836	2100	}
ashleymills	0:714293de3836	2101
ashleymills	0:714293de3836	2102	if (c != NULL) {
ashleymills	0:714293de3836	2103	fp_clamp(&q);
ashleymills	0:714293de3836	2104	fp_copy(&q, c);
ashleymills	0:714293de3836	2105	}
ashleymills	0:714293de3836	2106
ashleymills	0:714293de3836	2107	return FP_OKAY;
ashleymills	0:714293de3836	2108	}
ashleymills	0:714293de3836	2109
ashleymills	0:714293de3836	2110
ashleymills	0:714293de3836	2111	/* c = a mod b, 0 <= c < b */
ashleymills	0:714293de3836	2112	static int fp_mod_d(fp_int a, fp_digit b, fp_digit c)
ashleymills	0:714293de3836	2113	{
ashleymills	0:714293de3836	2114	return fp_div_d(a, b, NULL, c);
ashleymills	0:714293de3836	2115	}
ashleymills	0:714293de3836	2116
ashleymills	0:714293de3836	2117
ashleymills	0:714293de3836	2118	/* Miller-Rabin test of "a" to the base of "b" as described in
ashleymills	0:714293de3836	2119	* HAC pp. 139 Algorithm 4.24
ashleymills	0:714293de3836	2120	*
ashleymills	0:714293de3836	2121	* Sets result to 0 if definitely composite or 1 if probably prime.
ashleymills	0:714293de3836	2122	* Randomly the chance of error is no more than 1/4 and often
ashleymills	0:714293de3836	2123	* very much lower.
ashleymills	0:714293de3836	2124	*/
ashleymills	0:714293de3836	2125	static void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result)
ashleymills	0:714293de3836	2126	{
ashleymills	0:714293de3836	2127	fp_int n1, y, r;
ashleymills	0:714293de3836	2128	int s, j;
ashleymills	0:714293de3836	2129
ashleymills	0:714293de3836	2130	/* default */
ashleymills	0:714293de3836	2131	*result = FP_NO;
ashleymills	0:714293de3836	2132
ashleymills	0:714293de3836	2133	/* ensure b > 1 */
ashleymills	0:714293de3836	2134	if (fp_cmp_d(b, 1) != FP_GT) {
ashleymills	0:714293de3836	2135	return;
ashleymills	0:714293de3836	2136	}
ashleymills	0:714293de3836	2137
ashleymills	0:714293de3836	2138	/* get n1 = a - 1 */
ashleymills	0:714293de3836	2139	fp_init_copy(&n1, a);
ashleymills	0:714293de3836	2140	fp_sub_d(&n1, 1, &n1);
ashleymills	0:714293de3836	2141
ashleymills	0:714293de3836	2142	/* set 2*s r = n1 */
ashleymills	0:714293de3836	2143	fp_init_copy(&r, &n1);
ashleymills	0:714293de3836	2144
ashleymills	0:714293de3836	2145	/* count the number of least significant bits
ashleymills	0:714293de3836	2146	* which are zero
ashleymills	0:714293de3836	2147	*/
ashleymills	0:714293de3836	2148	s = fp_cnt_lsb(&r);
ashleymills	0:714293de3836	2149
ashleymills	0:714293de3836	2150	/* now divide n - 1 by 2*s /
ashleymills	0:714293de3836	2151	fp_div_2d (&r, s, &r, NULL);
ashleymills	0:714293de3836	2152
ashleymills	0:714293de3836	2153	/* compute y = b*r mod a /
ashleymills	0:714293de3836	2154	fp_init(&y);
ashleymills	0:714293de3836	2155	fp_exptmod(b, &r, a, &y);
ashleymills	0:714293de3836	2156
ashleymills	0:714293de3836	2157	/* if y != 1 and y != n1 do */
ashleymills	0:714293de3836	2158	if (fp_cmp_d (&y, 1) != FP_EQ && fp_cmp (&y, &n1) != FP_EQ) {
ashleymills	0:714293de3836	2159	j = 1;
ashleymills	0:714293de3836	2160	/* while j <= s-1 and y != n1 */
ashleymills	0:714293de3836	2161	while ((j <= (s - 1)) && fp_cmp (&y, &n1) != FP_EQ) {
ashleymills	0:714293de3836	2162	fp_sqrmod (&y, a, &y);
ashleymills	0:714293de3836	2163
ashleymills	0:714293de3836	2164	/* if y == 1 then composite */
ashleymills	0:714293de3836	2165	if (fp_cmp_d (&y, 1) == FP_EQ) {
ashleymills	0:714293de3836	2166	return;
ashleymills	0:714293de3836	2167	}
ashleymills	0:714293de3836	2168	++j;
ashleymills	0:714293de3836	2169	}
ashleymills	0:714293de3836	2170
ashleymills	0:714293de3836	2171	/* if y != n1 then composite */
ashleymills	0:714293de3836	2172	if (fp_cmp (&y, &n1) != FP_EQ) {
ashleymills	0:714293de3836	2173	return;
ashleymills	0:714293de3836	2174	}
ashleymills	0:714293de3836	2175	}
ashleymills	0:714293de3836	2176
ashleymills	0:714293de3836	2177	/* probably prime now */
ashleymills	0:714293de3836	2178	*result = FP_YES;
ashleymills	0:714293de3836	2179	}
ashleymills	0:714293de3836	2180
ashleymills	0:714293de3836	2181
ashleymills	0:714293de3836	2182	/* a few primes */
ashleymills	0:714293de3836	2183	static const fp_digit primes[256] = {
ashleymills	0:714293de3836	2184	0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
ashleymills	0:714293de3836	2185	0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
ashleymills	0:714293de3836	2186	0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
ashleymills	0:714293de3836	2187	0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
ashleymills	0:714293de3836	2188	0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
ashleymills	0:714293de3836	2189	0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
ashleymills	0:714293de3836	2190	0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
ashleymills	0:714293de3836	2191	0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
ashleymills	0:714293de3836	2192
ashleymills	0:714293de3836	2193	0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
ashleymills	0:714293de3836	2194	0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
ashleymills	0:714293de3836	2195	0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
ashleymills	0:714293de3836	2196	0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
ashleymills	0:714293de3836	2197	0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
ashleymills	0:714293de3836	2198	0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
ashleymills	0:714293de3836	2199	0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
ashleymills	0:714293de3836	2200	0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
ashleymills	0:714293de3836	2201
ashleymills	0:714293de3836	2202	0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
ashleymills	0:714293de3836	2203	0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
ashleymills	0:714293de3836	2204	0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
ashleymills	0:714293de3836	2205	0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
ashleymills	0:714293de3836	2206	0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
ashleymills	0:714293de3836	2207	0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
ashleymills	0:714293de3836	2208	0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
ashleymills	0:714293de3836	2209	0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
ashleymills	0:714293de3836	2210
ashleymills	0:714293de3836	2211	0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
ashleymills	0:714293de3836	2212	0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
ashleymills	0:714293de3836	2213	0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
ashleymills	0:714293de3836	2214	0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
ashleymills	0:714293de3836	2215	0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
ashleymills	0:714293de3836	2216	0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
ashleymills	0:714293de3836	2217	0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
ashleymills	0:714293de3836	2218	0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
ashleymills	0:714293de3836	2219	};
ashleymills	0:714293de3836	2220
ashleymills	0:714293de3836	2221	int fp_isprime(fp_int *a)
ashleymills	0:714293de3836	2222	{
ashleymills	0:714293de3836	2223	fp_int b;
ashleymills	0:714293de3836	2224	fp_digit d = 0;
ashleymills	0:714293de3836	2225	int r, res;
ashleymills	0:714293de3836	2226
ashleymills	0:714293de3836	2227	/* do trial division */
ashleymills	0:714293de3836	2228	for (r = 0; r < 256; r++) {
ashleymills	0:714293de3836	2229	fp_mod_d(a, primes[r], &d);
ashleymills	0:714293de3836	2230	if (d == 0) {
ashleymills	0:714293de3836	2231	return FP_NO;
ashleymills	0:714293de3836	2232	}
ashleymills	0:714293de3836	2233	}
ashleymills	0:714293de3836	2234
ashleymills	0:714293de3836	2235	/* now do 8 miller rabins */
ashleymills	0:714293de3836	2236	fp_init(&b);
ashleymills	0:714293de3836	2237	for (r = 0; r < 8; r++) {
ashleymills	0:714293de3836	2238	fp_set(&b, primes[r]);
ashleymills	0:714293de3836	2239	fp_prime_miller_rabin(a, &b, &res);
ashleymills	0:714293de3836	2240	if (res == FP_NO) {
ashleymills	0:714293de3836	2241	return FP_NO;
ashleymills	0:714293de3836	2242	}
ashleymills	0:714293de3836	2243	}
ashleymills	0:714293de3836	2244	return FP_YES;
ashleymills	0:714293de3836	2245	}
ashleymills	0:714293de3836	2246
ashleymills	0:714293de3836	2247
ashleymills	0:714293de3836	2248	/* c = [a, b] */
ashleymills	0:714293de3836	2249	void fp_lcm(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	2250	{
ashleymills	0:714293de3836	2251	fp_int t1, t2;
ashleymills	0:714293de3836	2252
ashleymills	0:714293de3836	2253	fp_init(&t1);
ashleymills	0:714293de3836	2254	fp_init(&t2);
ashleymills	0:714293de3836	2255	fp_gcd(a, b, &t1);
ashleymills	0:714293de3836	2256	if (fp_cmp_mag(a, b) == FP_GT) {
ashleymills	0:714293de3836	2257	fp_div(a, &t1, &t2, NULL);
ashleymills	0:714293de3836	2258	fp_mul(b, &t2, c);
ashleymills	0:714293de3836	2259	} else {
ashleymills	0:714293de3836	2260	fp_div(b, &t1, &t2, NULL);
ashleymills	0:714293de3836	2261	fp_mul(a, &t2, c);
ashleymills	0:714293de3836	2262	}
ashleymills	0:714293de3836	2263	}
ashleymills	0:714293de3836	2264
ashleymills	0:714293de3836	2265
ashleymills	0:714293de3836	2266	static const int lnz[16] = {
ashleymills	0:714293de3836	2267	4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
ashleymills	0:714293de3836	2268	};
ashleymills	0:714293de3836	2269
ashleymills	0:714293de3836	2270	/* Counts the number of lsbs which are zero before the first zero bit */
ashleymills	0:714293de3836	2271	int fp_cnt_lsb(fp_int *a)
ashleymills	0:714293de3836	2272	{
ashleymills	0:714293de3836	2273	int x;
ashleymills	0:714293de3836	2274	fp_digit q, qq;
ashleymills	0:714293de3836	2275
ashleymills	0:714293de3836	2276	/* easy out */
ashleymills	0:714293de3836	2277	if (fp_iszero(a) == 1) {
ashleymills	0:714293de3836	2278	return 0;
ashleymills	0:714293de3836	2279	}
ashleymills	0:714293de3836	2280
ashleymills	0:714293de3836	2281	/* scan lower digits until non-zero */
ashleymills	0:714293de3836	2282	for (x = 0; x < a->used && a->dp[x] == 0; x++);
ashleymills	0:714293de3836	2283	q = a->dp[x];
ashleymills	0:714293de3836	2284	x *= DIGIT_BIT;
ashleymills	0:714293de3836	2285
ashleymills	0:714293de3836	2286	/* now scan this digit until a 1 is found */
ashleymills	0:714293de3836	2287	if ((q & 1) == 0) {
ashleymills	0:714293de3836	2288	do {
ashleymills	0:714293de3836	2289	qq = q & 15;
ashleymills	0:714293de3836	2290	x += lnz[qq];
ashleymills	0:714293de3836	2291	q >>= 4;
ashleymills	0:714293de3836	2292	} while (qq == 0);
ashleymills	0:714293de3836	2293	}
ashleymills	0:714293de3836	2294	return x;
ashleymills	0:714293de3836	2295	}
ashleymills	0:714293de3836	2296
ashleymills	0:714293de3836	2297
ashleymills	0:714293de3836	2298	/* c = (a, b) */
ashleymills	0:714293de3836	2299	void fp_gcd(fp_int a, fp_int b, fp_int *c)
ashleymills	0:714293de3836	2300	{
ashleymills	0:714293de3836	2301	fp_int u, v, r;
ashleymills	0:714293de3836	2302
ashleymills	0:714293de3836	2303	/* either zero than gcd is the largest */
ashleymills	0:714293de3836	2304	if (fp_iszero (a) == 1 && fp_iszero (b) == 0) {
ashleymills	0:714293de3836	2305	fp_abs (b, c);
ashleymills	0:714293de3836	2306	return;
ashleymills	0:714293de3836	2307	}
ashleymills	0:714293de3836	2308	if (fp_iszero (a) == 0 && fp_iszero (b) == 1) {
ashleymills	0:714293de3836	2309	fp_abs (a, c);
ashleymills	0:714293de3836	2310	return;
ashleymills	0:714293de3836	2311	}
ashleymills	0:714293de3836	2312
ashleymills	0:714293de3836	2313	/* optimized. At this point if a == 0 then
ashleymills	0:714293de3836	2314	* b must equal zero too
ashleymills	0:714293de3836	2315	*/
ashleymills	0:714293de3836	2316	if (fp_iszero (a) == 1) {
ashleymills	0:714293de3836	2317	fp_zero(c);
ashleymills	0:714293de3836	2318	return;
ashleymills	0:714293de3836	2319	}
ashleymills	0:714293de3836	2320
ashleymills	0:714293de3836	2321	/* sort inputs */
ashleymills	0:714293de3836	2322	if (fp_cmp_mag(a, b) != FP_LT) {
ashleymills	0:714293de3836	2323	fp_init_copy(&u, a);
ashleymills	0:714293de3836	2324	fp_init_copy(&v, b);
ashleymills	0:714293de3836	2325	} else {
ashleymills	0:714293de3836	2326	fp_init_copy(&u, b);
ashleymills	0:714293de3836	2327	fp_init_copy(&v, a);
ashleymills	0:714293de3836	2328	}
ashleymills	0:714293de3836	2329
ashleymills	0:714293de3836	2330	fp_zero(&r);
ashleymills	0:714293de3836	2331	while (fp_iszero(&v) == FP_NO) {
ashleymills	0:714293de3836	2332	fp_mod(&u, &v, &r);
ashleymills	0:714293de3836	2333	fp_copy(&v, &u);
ashleymills	0:714293de3836	2334	fp_copy(&r, &v);
ashleymills	0:714293de3836	2335	}
ashleymills	0:714293de3836	2336	fp_copy(&u, c);
ashleymills	0:714293de3836	2337	}
ashleymills	0:714293de3836	2338
ashleymills	0:714293de3836	2339	#endif /* CYASSL_KEY_GEN */
ashleymills	0:714293de3836	2340
ashleymills	0:714293de3836	2341
ashleymills	0:714293de3836	2342	#if defined(HAVE_ECC) \|\| !defined(NO_PWDBASED)
ashleymills	0:714293de3836	2343	/* c = a + b */
ashleymills	0:714293de3836	2344	void fp_add_d(fp_int a, fp_digit b, fp_int c)
ashleymills	0:714293de3836	2345	{
ashleymills	0:714293de3836	2346	fp_int tmp;
ashleymills	0:714293de3836	2347	fp_set(&tmp, b);
ashleymills	0:714293de3836	2348	fp_add(a,&tmp,c);
ashleymills	0:714293de3836	2349	}
ashleymills	0:714293de3836	2350
ashleymills	0:714293de3836	2351	/* external compatibility */
ashleymills	0:714293de3836	2352	int mp_add_d(fp_int a, fp_digit b, fp_int c)
ashleymills	0:714293de3836	2353	{
ashleymills	0:714293de3836	2354	fp_add_d(a, b, c);
ashleymills	0:714293de3836	2355	return MP_OKAY;
ashleymills	0:714293de3836	2356	}
ashleymills	0:714293de3836	2357
ashleymills	0:714293de3836	2358	#endif /* HAVE_ECC \|\| !NO_PWDBASED */
ashleymills	0:714293de3836	2359
ashleymills	0:714293de3836	2360
ashleymills	0:714293de3836	2361	#ifdef HAVE_ECC
ashleymills	0:714293de3836	2362
ashleymills	0:714293de3836	2363	/* chars used in radix conversions */
ashleymills	0:714293de3836	2364	const char *fp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
ashleymills	0:714293de3836	2365
ashleymills	0:714293de3836	2366	static int fp_read_radix(fp_int a, const char str, int radix)
ashleymills	0:714293de3836	2367	{
ashleymills	0:714293de3836	2368	int y, neg;
ashleymills	0:714293de3836	2369	char ch;
ashleymills	0:714293de3836	2370
ashleymills	0:714293de3836	2371	/* make sure the radix is ok */
ashleymills	0:714293de3836	2372	if (radix < 2 \|\| radix > 64) {
ashleymills	0:714293de3836	2373	return FP_VAL;
ashleymills	0:714293de3836	2374	}
ashleymills	0:714293de3836	2375
ashleymills	0:714293de3836	2376	/* if the leading digit is a
ashleymills	0:714293de3836	2377	* minus set the sign to negative.
ashleymills	0:714293de3836	2378	*/
ashleymills	0:714293de3836	2379	if (*str == '-') {
ashleymills	0:714293de3836	2380	++str;
ashleymills	0:714293de3836	2381	neg = FP_NEG;
ashleymills	0:714293de3836	2382	} else {
ashleymills	0:714293de3836	2383	neg = FP_ZPOS;
ashleymills	0:714293de3836	2384	}
ashleymills	0:714293de3836	2385
ashleymills	0:714293de3836	2386	/* set the integer to the default of zero */
ashleymills	0:714293de3836	2387	fp_zero (a);
ashleymills	0:714293de3836	2388
ashleymills	0:714293de3836	2389	/* process each digit of the string */
ashleymills	0:714293de3836	2390	while (*str) {
ashleymills	0:714293de3836	2391	/* if the radix < 36 the conversion is case insensitive
ashleymills	0:714293de3836	2392	* this allows numbers like 1AB and 1ab to represent the same value
ashleymills	0:714293de3836	2393	* [e.g. in hex]
ashleymills	0:714293de3836	2394	*/
ashleymills	0:714293de3836	2395	ch = (char) ((radix < 36) ? XTOUPPER(str) : str);
ashleymills	0:714293de3836	2396	for (y = 0; y < 64; y++) {
ashleymills	0:714293de3836	2397	if (ch == fp_s_rmap[y]) {
ashleymills	0:714293de3836	2398	break;
ashleymills	0:714293de3836	2399	}
ashleymills	0:714293de3836	2400	}
ashleymills	0:714293de3836	2401
ashleymills	0:714293de3836	2402	/* if the char was found in the map
ashleymills	0:714293de3836	2403	* and is less than the given radix add it
ashleymills	0:714293de3836	2404	* to the number, otherwise exit the loop.
ashleymills	0:714293de3836	2405	*/
ashleymills	0:714293de3836	2406	if (y < radix) {
ashleymills	0:714293de3836	2407	fp_mul_d (a, (fp_digit) radix, a);
ashleymills	0:714293de3836	2408	fp_add_d (a, (fp_digit) y, a);
ashleymills	0:714293de3836	2409	} else {
ashleymills	0:714293de3836	2410	break;
ashleymills	0:714293de3836	2411	}
ashleymills	0:714293de3836	2412	++str;
ashleymills	0:714293de3836	2413	}
ashleymills	0:714293de3836	2414
ashleymills	0:714293de3836	2415	/* set the sign only if a != 0 */
ashleymills	0:714293de3836	2416	if (fp_iszero(a) != FP_YES) {
ashleymills	0:714293de3836	2417	a->sign = neg;
ashleymills	0:714293de3836	2418	}
ashleymills	0:714293de3836	2419	return FP_OKAY;
ashleymills	0:714293de3836	2420	}
ashleymills	0:714293de3836	2421
ashleymills	0:714293de3836	2422	/* fast math conversion */
ashleymills	0:714293de3836	2423	int mp_read_radix(mp_int a, const char str, int radix)
ashleymills	0:714293de3836	2424	{
ashleymills	0:714293de3836	2425	return fp_read_radix(a, str, radix);
ashleymills	0:714293de3836	2426	}
ashleymills	0:714293de3836	2427
ashleymills	0:714293de3836	2428	/* fast math conversion */
ashleymills	0:714293de3836	2429	int mp_set(fp_int *a, fp_digit b)
ashleymills	0:714293de3836	2430	{
ashleymills	0:714293de3836	2431	fp_set(a,b);
ashleymills	0:714293de3836	2432	return MP_OKAY;
ashleymills	0:714293de3836	2433	}
ashleymills	0:714293de3836	2434
ashleymills	0:714293de3836	2435	/* fast math conversion */
ashleymills	0:714293de3836	2436	int mp_sqr(fp_int A, fp_int B)
ashleymills	0:714293de3836	2437	{
ashleymills	0:714293de3836	2438	fp_sqr(A, B);
ashleymills	0:714293de3836	2439	return MP_OKAY;
ashleymills	0:714293de3836	2440	}
ashleymills	0:714293de3836	2441
ashleymills	0:714293de3836	2442	/* fast math conversion */
ashleymills	0:714293de3836	2443	int mp_montgomery_reduce(fp_int a, fp_int m, fp_digit mp)
ashleymills	0:714293de3836	2444	{
ashleymills	0:714293de3836	2445	fp_montgomery_reduce(a, m, mp);
ashleymills	0:714293de3836	2446	return MP_OKAY;
ashleymills	0:714293de3836	2447	}
ashleymills	0:714293de3836	2448
ashleymills	0:714293de3836	2449
ashleymills	0:714293de3836	2450	/* fast math conversion */
ashleymills	0:714293de3836	2451	int mp_montgomery_setup(fp_int a, fp_digit rho)
ashleymills	0:714293de3836	2452	{
ashleymills	0:714293de3836	2453	return fp_montgomery_setup(a, rho);
ashleymills	0:714293de3836	2454	}
ashleymills	0:714293de3836	2455
ashleymills	0:714293de3836	2456	int mp_div_2(fp_int * a, fp_int * b)
ashleymills	0:714293de3836	2457	{
ashleymills	0:714293de3836	2458	fp_div_2(a, b);
ashleymills	0:714293de3836	2459	return MP_OKAY;
ashleymills	0:714293de3836	2460	}
ashleymills	0:714293de3836	2461
ashleymills	0:714293de3836	2462
ashleymills	0:714293de3836	2463	int mp_init_copy(fp_int * a, fp_int * b)
ashleymills	0:714293de3836	2464	{
ashleymills	0:714293de3836	2465	fp_init_copy(a, b);
ashleymills	0:714293de3836	2466	return MP_OKAY;
ashleymills	0:714293de3836	2467	}
ashleymills	0:714293de3836	2468
ashleymills	0:714293de3836	2469
ashleymills	0:714293de3836	2470
ashleymills	0:714293de3836	2471	#endif /* HAVE_ECC */
ashleymills	0:714293de3836	2472
ashleymills	0:714293de3836	2473	#endif /* USE_FAST_MATH */

Repository toolbox

Export to desktop IDE

Repository details

Type:	Library
Created:	05 Sep 2013
Imports:	18
Forks:	0
Commits:	2
Dependents:	2
Dependencies:	0
Followers:	4

Developers

Ashley Mills

Francois Berder

The code in this repository is Apache licensed.

ctaocrypt/src/tfm.c@0:714293de3836, 2013-09-05 (annotated)

Who changed what in which revision?

Repository toolbox

Repository details

Developers

Important Information for this Arm website

Access Warning