+/******************************************************************************/
+
+void bigint_reduce(bigint_t *a, const bigint_t *r){
+ uint8_t rfbs = GET_FBS(r);
+ if (r->length_W == 0 || a->length_W == 0) {
+ return;
+ }
+
+ if (bigint_length_b(a) + 3 > bigint_length_b(r)) {
+ if ((r->length_W * sizeof(bigint_word_t) <= 4) && (a->length_W * sizeof(bigint_word_t) <= 4)) {
+ uint32_t p = 0, q = 0;
+ memcpy(&p, a->wordv, a->length_W * sizeof(bigint_word_t));
+ memcpy(&q, r->wordv, r->length_W * sizeof(bigint_word_t));
+ p %= q;
+ memcpy(a->wordv, &p, a->length_W * sizeof(bigint_word_t));
+ a->length_W = r->length_W;
+ bigint_adjust(a);
+ return;
+ }
+ unsigned shift;
+ while (a->length_W > r->length_W) {
+ shift = (a->length_W - r->length_W) * CHAR_BIT * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
+ bigint_sub_u_bitscale(a, r, shift);
+ }
+ while ((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)) {
+ shift = GET_FBS(a) - rfbs - 1;
+ bigint_sub_u_bitscale(a, r, shift);
+ }
+ }
+ while (bigint_cmp_u(a, r) >= 0) {
+ bigint_sub_u(a, a, r);
+ }
+ bigint_adjust(a);
+}
+
+/******************************************************************************/
+
+/* calculate dest = a**exp % r */
+/* using square&multiply */
+void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
+ if (a->length_W == 0) {
+ bigint_set_zero(dest);
+ return;
+ }
+
+ if(exp->length_W == 0){
+ dest->info = 0;
+ dest->length_W = 1;
+ dest->wordv[0] = 1;
+ return;
+ }
+
+ bigint_t res, base;
+ bigint_word_t t;
+ bigint_length_t i;
+ uint8_t j;
+
+ ALLOC_BIGINT_WORDS(base_w, MAX(a->length_W, r->length_W));
+ ALLOC_BIGINT_WORDS(res_w, r->length_W * 2);
+
+ res.wordv = res_w;
+ base.wordv = base_w;
+ bigint_copy(&base, a);
+ bigint_reduce(&base, r);
+ res.wordv[0] = 1;
+ res.length_W = 1;
+ res.info = 0;
+ bigint_adjust(&res);
+ bigint_copy(&res, &base);
+ uint8_t flag = 0;
+ t = exp->wordv[exp->length_W - 1];
+ for (i = exp->length_W; i > 0; --i) {
+ t = exp->wordv[i - 1];
+ for (j = BIGINT_WORD_SIZE; j > 0; --j) {
+ if (!flag) {
+ if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
+ flag = 1;
+ }
+ } else {
+ bigint_square(&res, &res);
+ bigint_reduce(&res, r);
+ if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
+ bigint_mul_u(&res, &res, &base);
+ bigint_reduce(&res, r);
+ }
+ }
+ t <<= 1;
+ }
+ }
+
+ SET_POS(&res);
+ bigint_copy(dest, &res);
+ FREE(res_w);
+ FREE(base_w);
+}
+
+/******************************************************************************/
+/* gcd <-- gcd(x,y) a*x+b*y=gcd */
+void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
+ bigint_length_t i = 0;
+ if(x->length_W == 0 || y->length_W == 0){
+ if(gcd){
+ bigint_set_zero(gcd);
+ }
+ if(a){
+ bigint_set_zero(a);
+ }
+ if(b){
+ bigint_set_zero(b);
+ }
+ return;
+ }
+ if(x->length_W == 1 && x->wordv[0] == 1){
+ if(gcd){
+ gcd->length_W = 1;
+ gcd->wordv[0] = 1;
+ gcd->info = 0;
+ }
+ if(a){
+ a->length_W = 1;
+ a->wordv[0] = 1;
+ SET_POS(a);
+ bigint_adjust(a);
+ }
+ if(b){
+ bigint_set_zero(b);
+ }
+ return;
+ }
+ if(y->length_W == 1 && y->wordv[0] == 1){
+ if(gcd){
+ gcd->length_W = 1;
+ gcd->wordv[0] = 1;
+ gcd->info = 0;
+ }
+ if(b){
+ b->length_W = 1;
+ b->wordv[0] = 1;
+ SET_POS(b);
+ bigint_adjust(b);
+ }
+ if(a){
+ bigint_set_zero(a);
+ }
+ return;
+ }
+
+ while(x->wordv[i] == 0 && y->wordv[i] == 0){
+ ++i;
+ }
+
+ ALLOC_BIGINT_WORDS(g_w, i + 2);
+ ALLOC_BIGINT_WORDS(x_w, x->length_W - i);
+ ALLOC_BIGINT_WORDS(y_w, y->length_W - i);
+ ALLOC_BIGINT_WORDS(u_w, x->length_W - i);
+ ALLOC_BIGINT_WORDS(v_w, y->length_W - i);
+ ALLOC_BIGINT_WORDS(a_w, y->length_W + 2);
+ ALLOC_BIGINT_WORDS(c_w, y->length_W + 2);
+ ALLOC_BIGINT_WORDS(b_w, x->length_W + 2);
+ ALLOC_BIGINT_WORDS(d_w, x->length_W + 2);
+
+ bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
+
+ g.wordv = g_w;
+ x_.wordv = x_w;
+ y_.wordv = y_w;
+ memset(g_w, 0, i * sizeof(bigint_word_t));
+ g_w[i] = 1;
+ g.length_W = i + 1;
+ g.info = 0;
+ x_.info = y_.info = 0;
+ x_.length_W = x->length_W - i;
+ y_.length_W = y->length_W - i;
+ memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
+ memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
+
+ for(i = 0; (x_.wordv[0] & ((bigint_word_t)1 << i)) == 0 && (y_.wordv[0] & ((bigint_word_t)1 << i)) == 0; ++i)
+ ;
+
+ bigint_adjust(&x_);
+ bigint_adjust(&y_);
+
+ if(i){
+ bigint_shiftleft_bits(&g, i);
+ bigint_shiftright_bits(&x_, i);
+ bigint_shiftright_bits(&y_, i);
+ }
+
+ u.wordv = u_w;
+ v.wordv = v_w;
+ a_.wordv = a_w;
+ b_.wordv = b_w;
+ c_.wordv = c_w;
+ d_.wordv = d_w;
+
+ bigint_copy(&u, &x_);
+ bigint_copy(&v, &y_);
+ a_.wordv[0] = 1;
+ a_.length_W = 1;
+ a_.info = 0;
+ d_.wordv[0] = 1;
+ d_.length_W = 1;
+ d_.info = 0;
+ bigint_set_zero(&b_);
+ bigint_set_zero(&c_);
+ do {
+ while ((u.wordv[0] & 1) == 0) {
+ bigint_shiftright_1bit(&u);
+ if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
+ bigint_add_s(&a_, &a_, &y_);
+ bigint_sub_s(&b_, &b_, &x_);
+ }
+ bigint_shiftright_1bit(&a_);
+ bigint_shiftright_1bit(&b_);
+ }
+ while ((v.wordv[0] & 1) == 0) {
+ bigint_shiftright_1bit(&v);
+ if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
+ bigint_add_s(&c_, &c_, &y_);
+ bigint_sub_s(&d_, &d_, &x_);
+ }
+ bigint_shiftright_1bit(&c_);
+ bigint_shiftright_1bit(&d_);
+ }
+ if(bigint_cmp_u(&u, &v) >= 0){
+ bigint_sub_u(&u, &u, &v);
+ bigint_sub_s(&a_, &a_, &c_);
+ bigint_sub_s(&b_, &b_, &d_);
+ }else{
+ bigint_sub_u(&v, &v, &u);
+ bigint_sub_s(&c_, &c_, &a_);
+ bigint_sub_s(&d_, &d_, &b_);
+ }
+ } while(u.length_W);
+ if(gcd){
+ bigint_mul_s(gcd, &v, &g);
+ }
+ if(a){
+ bigint_copy(a, &c_);
+ }
+ if(b){
+ bigint_copy(b, &d_);
+ }
+
+ FREE(d_w);
+ FREE(b_w);
+ FREE(c_w);
+ FREE(a_w);
+ FREE(v_w);
+ FREE(u_w);
+ FREE(y_w);
+ FREE(x_w);
+ FREE(g_w);
+}
+
+/******************************************************************************/
+
+void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
+ bigint_gcdext(NULL, dest, NULL, a, m);
+ while(dest->info&BIGINT_NEG_MASK){
+ bigint_add_s(dest, dest, m);
+ }
+}
+
+/******************************************************************************/
+
+void bigint_changeendianess(bigint_t *a){
+ uint8_t t, *p, *q;
+ p = (uint8_t*)a->wordv;
+ q = p + a->length_W * sizeof(bigint_word_t) - 1;
+ while(p < q){
+ t = *p;
+ *p = *q;
+ *q = t;
+ ++p; --q;
+ }
+}
+
+/******************************************************************************/
+
+void bigint_mul_word_u(bigint_t *a, bigint_word_t b){
+ bigint_wordplus_t c0 = 0, c1 = 0;
+ bigint_length_t i;
+
+ if(b == 0){
+ bigint_set_zero(a);
+ return;
+ }
+
+ for(i = 0; i < a->length_W; ++i){
+ c1 = ((bigint_wordplus_t)(a->wordv[i])) * ((bigint_wordplus_t)b);
+ c1 += c0;
+ a->wordv[i] = (bigint_word_t)c1;
+ c0 = c1 >> BIGINT_WORD_SIZE;
+ }
+ if(c0){
+ a->wordv[a->length_W] = (bigint_word_t)c0;
+ a->length_W += 1;
+ }
+ bigint_adjust(a);
+}
+
+/******************************************************************************/
+
+void bigint_clip(bigint_t *dest, bigint_length_t length_W){
+ if(dest->length_W > length_W){
+ dest->length_W = length_W;
+ }
+ bigint_adjust(dest);
+}
+
+/******************************************************************************/
+/*
+ * m_ = m * m'[0]
+ * dest = (a * b) % m (?)
+ */
+
+void bigint_mont_mul(bigint_t *dest, const bigint_t *a, const bigint_t *b, const bigint_t *m, const bigint_t *m_){
+ const bigint_length_t s = MAX(MAX(a->length_W, b->length_W), m->length_W);
+ bigint_t u, t;
+
+ bigint_length_t i;
+
+ if (a->length_W == 0 || b->length_W == 0) {
+ bigint_set_zero(dest);
+ return;
+ }
+ ALLOC_BIGINT_WORDS(u_w, s + 2);
+ ALLOC_BIGINT_WORDS(t_w, s + 2);
+ u.wordv = u_w;
+ u.info = 0;
+ u.length_W = 0;
+ t.wordv = t_w;
+ for (i = 0; i < a->length_W; ++i) {
+ bigint_copy(&t, b);
+ bigint_mul_word_u(&t, a->wordv[i]);
+ bigint_add_u(&u, &u, &t);
+ bigint_copy(&t, m_);
+ if (u.length_W != 0) {
+ bigint_mul_word_u(&t, u.wordv[0]);
+ bigint_add_u(&u, &u, &t);
+ }
+ bigint_shiftright_1word(&u);
+ }
+ for (; i < s; ++i) {
+ bigint_copy(&t, m_);
+ if (u.length_W != 0) {
+ bigint_mul_word_u(&t, u.wordv[0]);
+ bigint_add_u(&u, &u, &t);
+ }
+ bigint_shiftright_1word(&u);
+ }
+ bigint_reduce(&u, m);
+ bigint_copy(dest, &u);
+ FREE(t_w);
+ FREE(u_w);
+}
+
+/******************************************************************************/
+
+void bigint_mont_red(bigint_t *dest, const bigint_t *a, const bigint_t *m, const bigint_t *m_){
+ bigint_t u, t;
+ bigint_length_t i, s = MAX(a->length_W, MAX(m->length_W, m_->length_W));
+
+ if (a->length_W == 0) {
+ bigint_set_zero(dest);
+ return;
+ }
+
+ ALLOC_BIGINT_WORDS(u_w, s + 2);
+ ALLOC_BIGINT_WORDS(t_w, s + 2);
+ t.wordv = t_w;
+ u.wordv = u_w;
+ bigint_copy(&u, a);
+ for (i = 0; i < m->length_W; ++i) {
+ bigint_copy(&t, m_);
+ if (u.length_W != 0) {
+ bigint_mul_word_u(&t, u.wordv[0]);
+ bigint_add_u(&u, &u, &t);
+ }
+ bigint_shiftright_1word(&u);
+ }
+ bigint_reduce(&u, m);
+ bigint_copy(dest, &u);
+ FREE(t_w);
+ FREE(u_w);
+}
+
+/******************************************************************************/
+/*
+ * m_ = m * (- m0^-1 (mod 2^W))
+ */
+void bigint_mont_gen_m_(bigint_t* dest, const bigint_t* m){
+ bigint_word_t x_w[2], m_w_0[1];
+ bigint_t x, m_0;
+ if (m->length_W == 0) {
+ bigint_set_zero(dest);
+ return;
+ }
+ if ((m->wordv[0] & 1) == 0) {
+#if DEBUG
+ printf_P(PSTR("ERROR: m must not be even, m = "));
+ bigint_print_hex(m);
+ putchar('\n');
+ uart0_flush();
+#endif
+ return;
+ }
+ x.wordv = x_w;
+ x.info = 0;
+ x.length_W = 2;
+ x_w[0] = 0;
+ x_w[1] = 1;
+ m_0.wordv = m_w_0;
+ m_0.info = 0;
+ m_0.length_W = 1;
+ m_0.wordv[0] = m->wordv[0];
+ bigint_adjust(&x);
+ bigint_adjust(&m_0);
+ bigint_inverse(dest, &m_0, &x);
+ bigint_sub_s(&x, &x, dest);
+ bigint_copy(dest, m);
+ bigint_mul_word_u(dest, x.wordv[0]);
+
+}
+
+/******************************************************************************/
+
+/*
+ * dest = a * R mod m
+ */
+void bigint_mont_trans(bigint_t *dest, const bigint_t *a, const bigint_t *m){
+ bigint_t t;
+
+ ALLOC_BIGINT_WORDS(t_w, a->length_W + m->length_W);
+ t.wordv = t_w;
+ memset(t_w, 0, m->length_W * sizeof(bigint_word_t));
+ memcpy(&t_w[m->length_W], a->wordv, a->length_W * sizeof(bigint_word_t));
+ t.info = a->info;
+ t.length_W = a->length_W + m->length_W;
+ bigint_reduce(&t, m);
+ bigint_copy(dest, &t);
+ FREE(t_w);
+}
+
+/******************************************************************************/
+
+/* calculate dest = a**exp % r */
+/* using square&multiply */
+void bigint_expmod_u_mont_accel(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r, const bigint_t *m_){
+ if(r->length_W == 0) {
+ return;
+ }
+
+ bigint_t res, ax;
+ bigint_word_t t;
+ bigint_length_t i;
+ uint8_t j;
+
+ if (exp->length_W == 0) {
+ dest->length_W = 1;
+ dest->info = 0;
+ dest->wordv[0] = 1;
+ return;
+ }
+
+ ALLOC_BIGINT_WORDS(res_w, r->length_W * 2);
+ ALLOC_BIGINT_WORDS(ax_w, MAX(r->length_W, a->length_W));
+
+ res.wordv = res_w;
+ ax.wordv = ax_w;
+
+ res.wordv[0] = 1;
+ res.length_W = 1;
+ res.info = 0;
+
+ bigint_copy(&ax, a);
+ bigint_reduce(&ax, r);
+
+ bigint_mont_trans(&ax, &ax, r);
+ bigint_mont_trans(&res, &res, r);
+
+ uint8_t flag = 0;
+ t = exp->wordv[exp->length_W - 1];
+ for (i = exp->length_W; i > 0; --i) {
+ t = exp->wordv[i - 1];
+ for(j = BIGINT_WORD_SIZE; j > 0; --j){
+ if (!flag) {
+ if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
+ flag = 1;
+ }
+ }
+ if (flag) {
+ bigint_square(&res, &res);
+ bigint_mont_red(&res, &res, r, m_);
+ if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
+ bigint_mont_mul(&res, &res, &ax, r, m_);
+ }
+ }
+ t <<= 1;
+ }
+ }
+ SET_POS(&res);
+ bigint_mont_red(dest, &res, r, m_);
+ FREE(ax_w);
+ FREE(res_w);
+}
+
+/******************************************************************************/
+
+void bigint_expmod_u_mont_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
+ if(r->length_W == 0) {
+ return;
+ }
+ if(a->length_W == 0) {
+ bigint_set_zero(dest);
+ return;
+ }
+ bigint_t m_;
+ bigint_word_t m_w_[r->length_W + 1];
+ m_.wordv = m_w_;
+ bigint_mont_gen_m_(&m_, r);
+ bigint_expmod_u_mont_accel(dest, a, exp, r,&m_);
+}
+
+/******************************************************************************/
+
+void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
+#if 0
+ printf("\nDBG: expmod_u (a ** e %% m) <%s %s %d>\n\ta: ", __FILE__, __func__, __LINE__);
+ bigint_print_hex(a);
+ printf("\n\te: ");
+ bigint_print_hex(exp);
+ printf("\n\tm: ");
+ bigint_print_hex(r);
+#endif
+ if (r->wordv[0] & 1) {
+ bigint_expmod_u_mont_sam(dest, a, exp, r);
+ } else {
+ bigint_expmod_u_sam(dest, a, exp, r);
+ }
+}
+
+
+
+
+
+
+
+
+
+