http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/BIG.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/BIG.go b/go/amcl-go/BIG.go new file mode 100644 index 0000000..a1c5184 --- /dev/null +++ b/go/amcl-go/BIG.go @@ -0,0 +1,956 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +/* AMCL BIG number class */ + +package amcl + +import "strconv" + +//import "fmt" + +type BIG struct { + w [NLEN]int64 +} + +func NewBIG() *BIG { + b := new(BIG) + for i := 0; i < NLEN; i++ { + b.w[i] = 0 + } + return b +} + +func NewBIGint(x int) *BIG { + b := new(BIG) + b.w[0] = int64(x) + for i := 1; i < NLEN; i++ { + b.w[i] = 0 + } + return b +} + +func NewBIGcopy(x *BIG) *BIG { + b := new(BIG) + for i := 0; i < NLEN; i++ { + b.w[i] = x.w[i] + } + return b +} + +func NewBIGdcopy(x *DBIG) *BIG { + b := new(BIG) + for i := 0; i < NLEN; i++ { + b.w[i] = x.w[i] + } + return b +} + +func NewBIGints(x [NLEN]int64) *BIG { + b := new(BIG) + for i := 0; i < NLEN; i++ { + b.w[i] = x[i] + } + return b +} + +func (r *BIG) get(i int) int64 { + return r.w[i] +} + +func (r *BIG) set(i int, x int64) { + r.w[i] = x +} + +func (r *BIG) xortop(x int64) { + r.w[NLEN-1] ^= x +} + +func (r *BIG) ortop(x int64) { + r.w[NLEN-1] |= x +} + +/* test for zero */ +func (r *BIG) iszilch() bool { + for i := 0; i < NLEN; i++ { + if r.w[i] != 0 { + return false + } + } + return true +} + +/* set to zero */ +func (r *BIG) zero() { + for i := 0; i < NLEN; i++ { + r.w[i] = 0 + } +} + +/* Test for equal to one */ +func (r *BIG) isunity() bool { + for i := 1; i < NLEN; i++ { + if r.w[i] != 0 { + return false + } + } + if r.w[0] != 1 { + return false + } + return true +} + +/* set to one */ +func (r *BIG) one() { + r.w[0] = 1 + for i := 1; i < NLEN; i++ { + r.w[i] = 0 + } +} + +/* Copy from another BIG */ +func (r *BIG) copy(x *BIG) { + for i := 0; i < NLEN; i++ { + r.w[i] = x.w[i] + } +} + +/* Copy from another DBIG */ +func (r *BIG) dcopy(x *DBIG) { + for i := 0; i < NLEN; i++ { + r.w[i] = x.w[i] + } +} + +/* calculate Field Excess */ +func EXCESS(a *BIG) int64 { + return ((a.w[NLEN-1] & OMASK) >> (MODBITS % BASEBITS)) +} + +/* normalise BIG - force all digits < 2^BASEBITS */ +func (r *BIG) norm() int64 { + var carry int64 = 0 + for i := 0; i < NLEN-1; i++ { + d := r.w[i] + carry + r.w[i] = d & MASK + carry = d >> BASEBITS + } + r.w[NLEN-1] = (r.w[NLEN-1] + carry) + + return (r.w[NLEN-1] >> ((8 * MODBYTES) % BASEBITS)) +} + +/* Conditional swap of two bigs depending on d using XOR - no branches */ +func (r *BIG) cswap(b *BIG, d int32) { + var c = int64(d) + c = ^(c - 1) + + for i := 0; i < NLEN; i++ { + t := c & (r.w[i] ^ b.w[i]) + r.w[i] ^= t + b.w[i] ^= t + } +} + +func (r *BIG) cmove(g *BIG, d int32) { + var b = int64(-d) + + for i := 0; i < NLEN; i++ { + r.w[i] ^= (r.w[i] ^ g.w[i]) & b + } +} + +/* Shift right by less than a word */ +func (r *BIG) fshr(k uint) int64 { + w := r.w[0] & ((int64(1) << k) - 1) /* shifted out part */ + for i := 0; i < NLEN-1; i++ { + r.w[i] = (r.w[i] >> k) | ((r.w[i+1] << (BASEBITS - k)) & MASK) + } + r.w[NLEN-1] = r.w[NLEN-1] >> k + return w +} + +/* general shift right */ +func (r *BIG) shr(k uint) { + n := (k % BASEBITS) + m := int(k / BASEBITS) + for i := 0; i < NLEN-m-1; i++ { + r.w[i] = (r.w[m+i] >> n) | ((r.w[m+i+1] << (BASEBITS - n)) & MASK) + } + r.w[NLEN-m-1] = r.w[NLEN-1] >> n + for i := NLEN - m; i < NLEN; i++ { + r.w[i] = 0 + } +} + +/* Shift right by less than a word */ +func (r *BIG) fshl(k uint) int64 { + r.w[NLEN-1] = (r.w[NLEN-1] << k) | (r.w[NLEN-2] >> (BASEBITS - k)) + for i := NLEN - 2; i > 0; i-- { + r.w[i] = ((r.w[i] << k) & MASK) | (r.w[i-1] >> (BASEBITS - k)) + } + r.w[0] = (r.w[0] << k) & MASK + return (r.w[NLEN-1] >> ((8 * MODBYTES) % BASEBITS)) /* return excess - only used in ff.c */ +} + +/* general shift left */ +func (r *BIG) shl(k uint) { + n := k % BASEBITS + m := int(k / BASEBITS) + + r.w[NLEN-1] = (r.w[NLEN-1-m] << n) | (r.w[NLEN-m-2] >> (BASEBITS - n)) + for i := NLEN - 2; i > m; i-- { + r.w[i] = ((r.w[i-m] << n) & MASK) | (r.w[i-m-1] >> (BASEBITS - n)) + } + r.w[m] = (r.w[0] << n) & MASK + for i := 0; i < m; i++ { + r.w[i] = 0 + } +} + +/* return number of bits */ +func (r *BIG) nbits() int { + k := NLEN - 1 + r.norm() + for k >= 0 && r.w[k] == 0 { + k-- + } + if k < 0 { + return 0 + } + bts := int(BASEBITS) * k + c := r.w[k] + for c != 0 { + c /= 2 + bts++ + } + return bts +} + +/* Convert to Hex String */ +func (r *BIG) toString() string { + s := "" + len := r.nbits() + + if len%4 == 0 { + len /= 4 + } else { + len /= 4 + len++ + + } + MB := int(MODBYTES * 2) + if len < MB { + len = MB + } + + for i := len - 1; i >= 0; i-- { + b := NewBIGcopy(r) + + b.shr(uint(i * 4)) + s += strconv.FormatInt(b.w[0]&15, 16) + } + return s +} + +func (r *BIG) add(x *BIG) { + for i := 0; i < NLEN; i++ { + r.w[i] = r.w[i] + x.w[i] + } +} + +/* return this+x */ +func (r *BIG) plus(x *BIG) *BIG { + s := new(BIG) + for i := 0; i < NLEN; i++ { + s.w[i] = r.w[i] + x.w[i] + } + return s +} + +/* this+=x, where x is int */ +func (r *BIG) inc(x int) { + r.norm() + r.w[0] += int64(x) +} + +/* return this-x */ +func (r *BIG) minus(x *BIG) *BIG { + d := new(BIG) + for i := 0; i < NLEN; i++ { + d.w[i] = r.w[i] - x.w[i] + } + return d +} + +/* this-=x */ +func (r *BIG) sub(x *BIG) { + for i := 0; i < NLEN; i++ { + r.w[i] = r.w[i] - x.w[i] + } +} + +/* reverse subtract this=x-this */ +func (r *BIG) rsub(x *BIG) { + for i := 0; i < NLEN; i++ { + r.w[i] = x.w[i] - r.w[i] + } +} + +/* this-=x, where x is int */ +func (r *BIG) dec(x int) { + r.norm() + r.w[0] -= int64(x) +} + +/* this*=x, where x is small int<NEXCESS */ +func (r *BIG) imul(c int) { + for i := 0; i < NLEN; i++ { + r.w[i] *= int64(c) + } +} + +/* convert this BIG to byte array */ +func (r *BIG) tobytearray(b []byte, n int) { + r.norm() + c := NewBIGcopy(r) + + for i := int(MODBYTES) - 1; i >= 0; i-- { + b[i+n] = byte(c.w[0]) + c.fshr(8) + } +} + +/* convert from byte array to BIG */ +func frombytearray(b []byte, n int) *BIG { + m := NewBIG() + for i := 0; i < int(MODBYTES); i++ { + m.fshl(8) + m.w[0] += int64(b[i+n] & 0xff) + } + return m +} + +func (r *BIG) toBytes(b []byte) { + r.tobytearray(b, 0) +} + +func fromBytes(b []byte) *BIG { + return frombytearray(b, 0) +} + +/* set this[i]+=x*y+c, and return high part */ + +func (r *BIG) muladd(a int64, b int64, c int64, i int) int64 { + x0 := a & HMASK + x1 := (a >> HBITS) + y0 := b & HMASK + y1 := (b >> HBITS) + bot := x0 * y0 + top := x1 * y1 + mid := x0*y1 + x1*y0 + x0 = mid & HMASK + x1 = (mid >> HBITS) + bot += x0 << HBITS + bot += c + bot += r.w[i] + top += x1 + carry := bot >> BASEBITS + bot &= MASK + top += carry + r.w[i] = bot + return top +} + +/* this*=x, where x is >NEXCESS */ +func (r *BIG) pmul(c int) int64 { + var carry int64 = 0 + r.norm() + for i := 0; i < NLEN; i++ { + ak := r.w[i] + r.w[i] = 0 + carry = r.muladd(ak, int64(c), carry, i) + } + return carry +} + +/* this*=c and catch overflow in DBIG */ +func (r *BIG) pxmul(c int) *DBIG { + m := NewDBIG() + var carry int64 = 0 + for j := 0; j < NLEN; j++ { + carry = m.muladd(r.w[j], int64(c), carry, j) + } + m.w[NLEN] = carry + return m +} + +/* divide by 3 */ +func (r *BIG) div3() int { + var carry int64 = 0 + r.norm() + base := (int64(1) << BASEBITS) + for i := NLEN - 1; i >= 0; i-- { + ak := (carry*base + r.w[i]) + r.w[i] = ak / 3 + carry = ak % 3 + } + return int(carry) +} + +/* return a*b where result fits in a BIG */ +func smul(a *BIG, b *BIG) *BIG { + var carry int64 + c := NewBIG() + for i := 0; i < NLEN; i++ { + carry = 0 + for j := 0; j < NLEN; j++ { + if i+j < NLEN { + carry = c.muladd(a.w[i], b.w[j], carry, i+j) + } + } + } + return c +} + +/* Compare a and b, return 0 if a==b, -1 if a<b, +1 if a>b. Inputs must be normalised */ +func comp(a *BIG, b *BIG) int { + for i := NLEN - 1; i >= 0; i-- { + if a.w[i] == b.w[i] { + continue + } + if a.w[i] > b.w[i] { + return 1 + } else { + return -1 + } + } + return 0 +} + +/* return parity */ +func (r *BIG) parity() int { + return int(r.w[0] % 2) +} + +/* return n-th bit */ +func (r *BIG) bit(n int) int { + if (r.w[n/int(BASEBITS)] & (int64(1) << (uint(n) % BASEBITS))) > 0 { + return 1 + } + return 0 +} + +/* return n last bits */ +func (r *BIG) lastbits(n int) int { + msk := (1 << uint(n)) - 1 + r.norm() + return (int(r.w[0])) & msk +} + +/* set x = x mod 2^m */ +func (r *BIG) mod2m(m uint) { + wd := int(m / BASEBITS) + bt := m % BASEBITS + msk := (int64(1) << bt) - 1 + r.w[wd] &= msk + for i := wd + 1; i < NLEN; i++ { + r.w[i] = 0 + } +} + +/* Arazi and Qi inversion mod 256 */ +func invmod256(a int) int { + var t1 int = 0 + c := (a >> 1) & 1 + t1 += c + t1 &= 1 + t1 = 2 - t1 + t1 <<= 1 + U := t1 + 1 + + // i=2 + b := a & 3 + t1 = U * b + t1 >>= 2 + c = (a >> 2) & 3 + t2 := (U * c) & 3 + t1 += t2 + t1 *= U + t1 &= 3 + t1 = 4 - t1 + t1 <<= 2 + U += t1 + + // i=4 + b = a & 15 + t1 = U * b + t1 >>= 4 + c = (a >> 4) & 15 + t2 = (U * c) & 15 + t1 += t2 + t1 *= U + t1 &= 15 + t1 = 16 - t1 + t1 <<= 4 + U += t1 + + return U +} + +/* a=1/a mod 2^256. This is very fast! */ +func (r *BIG) invmod2m() { + U := NewBIG() + b := NewBIG() + c := NewBIG() + + U.inc(invmod256(r.lastbits(8))) + + for i := 8; i < 256; i <<= 1 { + ui := uint(i) + b.copy(r) + b.mod2m(ui) + t1 := smul(U, b) + t1.shr(ui) + c.copy(r) + c.shr(ui) + c.mod2m(ui) + + t2 := smul(U, c) + t2.mod2m(ui) + t1.add(t2) + b = smul(t1, U) + t1.copy(b) + t1.mod2m(ui) + + t2.one() + t2.shl(ui) + t1.rsub(t2) + t1.norm() + t1.shl(ui) + U.add(t1) + } + r.copy(U) +} + +/* reduce this mod m */ +func (r *BIG) mod(m *BIG) { + r.norm() + if comp(r, m) < 0 { + return + } + + m.fshl(1) + k := 1 + + for comp(r, m) >= 0 { + m.fshl(1) + k++ + } + + for k > 0 { + m.fshr(1) + if comp(r, m) >= 0 { + r.sub(m) + r.norm() + } + k-- + } +} + +/* divide this by m */ +func (r *BIG) div(m *BIG) { + k := 0 + r.norm() + e := NewBIGint(1) + b := NewBIGcopy(r) + r.zero() + + for comp(b, m) >= 0 { + e.fshl(1) + m.fshl(1) + k++ + } + + for k > 0 { + m.fshr(1) + e.fshr(1) + if comp(b, m) >= 0 { + r.add(e) + r.norm() + b.sub(m) + b.norm() + } + k-- + } +} + +/* get 8*MODBYTES size random number */ +func random(rng *RAND) *BIG { + m := NewBIG() + var j int = 0 + var r byte = 0 + /* generate random BIG */ + for i := 0; i < 8*int(MODBYTES); i++ { + if j == 0 { + r = rng.GetByte() + } else { + r >>= 1 + } + + b := int64(r & 1) + m.shl(1) + m.w[0] += b // m.inc(b) + j++ + j &= 7 + } + return m +} + +/* Create random BIG in portable way, one bit at a time */ +func randomnum(q *BIG, rng *RAND) *BIG { + d := NewDBIG() + var j int = 0 + var r byte = 0 + for i := 0; i < 2*int(MODBITS); i++ { + if j == 0 { + r = rng.GetByte() + } else { + r >>= 1 + } + + b := int64(r & 1) + d.shl(1) + d.w[0] += b // m.inc(b); + j++ + j &= 7 + } + m := d.mod(q) + return m +} + +/* return NAF value as +/- 1, 3 or 5. x and x3 should be normed. +nbs is number of bits processed, and nzs is number of trailing 0s detected */ +func nafbits(x *BIG, x3 *BIG, i int) [3]int { + var n [3]int + var j int + nb := x3.bit(i) - x.bit(i) + + n[1] = 1 + n[0] = 0 + if nb == 0 { + n[0] = 0 + return n + } + if i == 0 { + n[0] = nb + return n + } + if nb > 0 { + n[0] = 1 + } else { + n[0] = (-1) + } + + for j = i - 1; j > 0; j-- { + n[1]++ + n[0] *= 2 + nb = x3.bit(j) - x.bit(j) + if nb > 0 { + n[0] += 1 + } + if nb < 0 { + n[0] -= 1 + } + if n[0] > 5 || n[0] < -5 { + break + } + } + + if n[0]%2 != 0 && j != 0 { /* backtrack */ + if nb > 0 { + n[0] = (n[0] - 1) / 2 + } + if nb < 0 { + n[0] = (n[0] + 1) / 2 + } + n[1]-- + } + for n[0]%2 == 0 { /* remove trailing zeros */ + n[0] /= 2 + n[2]++ + n[1]-- + } + return n +} + +/* return a*b as DBIG */ +func mul(a *BIG, b *BIG) *DBIG { + c := NewDBIG() + var carry int64 + a.norm() + b.norm() + + for i := 0; i < NLEN; i++ { + carry = 0 + for j := 0; j < NLEN; j++ { + carry = c.muladd(a.w[i], b.w[j], carry, i+j) + } + c.w[NLEN+i] = carry + } + + return c +} + +/* return a^2 as DBIG */ +func sqr(a *BIG) *DBIG { + c := NewDBIG() + var carry int64 + a.norm() + for i := 0; i < NLEN; i++ { + carry = 0 + for j := i + 1; j < NLEN; j++ { + carry = c.muladd(2*a.w[i], a.w[j], carry, i+j) + } + c.w[NLEN+i] = carry + } + + for i := 0; i < NLEN; i++ { + c.w[2*i+1] += c.muladd(a.w[i], a.w[i], 0, 2*i) + } + c.norm() + return c +} + +/* reduce a DBIG to a BIG using the appropriate form of the modulus */ +func mod(d *DBIG) *BIG { + var b *BIG + if MODTYPE == PSEUDO_MERSENNE { + t := d.split(MODBITS) + b = NewBIGdcopy(d) + + v := t.pmul(int(MConst)) + tw := t.w[NLEN-1] + t.w[NLEN-1] &= TMASK + t.w[0] += (MConst * ((tw >> TBITS) + (v << (BASEBITS - TBITS)))) + + b.add(t) + } + if MODTYPE == MONTGOMERY_FRIENDLY { + for i := 0; i < NLEN; i++ { + d.w[NLEN+i] += d.muladd(d.w[i], MConst-1, d.w[i], NLEN+i-1) + } + b = NewBIG() + + for i := 0; i < NLEN; i++ { + b.w[i] = d.w[NLEN+i] + } + } + + if MODTYPE == NOT_SPECIAL { + md := NewBIGints(Modulus) + var carry, m int64 + for i := 0; i < NLEN; i++ { + if MConst == -1 { + m = (-d.w[i]) & MASK + } else { + if MConst == 1 { + m = d.w[i] + } else { + m = (MConst * d.w[i]) & MASK + } + } + + carry = 0 + for j := 0; j < NLEN; j++ { + carry = d.muladd(m, md.w[j], carry, i+j) + } + d.w[NLEN+i] += carry + } + + b = NewBIG() + for i := 0; i < NLEN; i++ { + b.w[i] = d.w[NLEN+i] + } + + } + b.norm() + return b +} + +/* return a*b mod m */ +func modmul(a, b, m *BIG) *BIG { + a.mod(m) + b.mod(m) + d := mul(a, b) + return d.mod(m) +} + +/* return a^2 mod m */ +func modsqr(a, m *BIG) *BIG { + a.mod(m) + d := sqr(a) + return d.mod(m) +} + +/* return -a mod m */ +func modneg(a, m *BIG) *BIG { + a.mod(m) + return m.minus(a) +} + +/* return this^e mod m */ +func (r *BIG) powmod(e *BIG, m *BIG) *BIG { + r.norm() + e.norm() + a := NewBIGint(1) + z := NewBIGcopy(e) + s := NewBIGcopy(r) + for true { + bt := z.parity() + z.fshr(1) + if bt == 1 { + a = modmul(a, s, m) + } + if z.iszilch() { + break + } + s = modsqr(s, m) + } + return a +} + +/* Jacobi Symbol (this/p). Returns 0, 1 or -1 */ +func (r *BIG) jacobi(p *BIG) int { + m := 0 + t := NewBIGint(0) + x := NewBIGint(0) + n := NewBIGint(0) + zilch := NewBIGint(0) + one := NewBIGint(1) + if p.parity() == 0 || comp(r, zilch) == 0 || comp(p, one) <= 0 { + return 0 + } + r.norm() + x.copy(r) + n.copy(p) + x.mod(p) + + for comp(n, one) > 0 { + if comp(x, zilch) == 0 { + return 0 + } + n8 := n.lastbits(3) + k := 0 + for x.parity() == 0 { + k++ + x.shr(1) + } + if k%2 == 1 { + m += (n8*n8 - 1) / 8 + } + m += (n8 - 1) * (x.lastbits(2) - 1) / 4 + t.copy(n) + t.mod(x) + n.copy(x) + x.copy(t) + m %= 2 + + } + if m == 0 { + return 1 + } + return -1 +} + +/* this=1/this mod p. Binary method */ +func (r *BIG) invmodp(p *BIG) { + r.mod(p) + u := NewBIGcopy(r) + + v := NewBIGcopy(p) + x1 := NewBIGint(1) + x2 := NewBIGint(0) + t := NewBIGint(0) + one := NewBIGint(1) + for comp(u, one) != 0 && comp(v, one) != 0 { + for u.parity() == 0 { + u.shr(1) + if x1.parity() != 0 { + x1.add(p) + x1.norm() + } + x1.shr(1) + } + for v.parity() == 0 { + v.shr(1) + if x2.parity() != 0 { + x2.add(p) + x2.norm() + } + x2.shr(1) + } + if comp(u, v) >= 0 { + u.sub(v) + u.norm() + if comp(x1, x2) >= 0 { + x1.sub(x2) + } else { + t.copy(p) + t.sub(x2) + x1.add(t) + } + x1.norm() + } else { + v.sub(u) + v.norm() + if comp(x2, x1) >= 0 { + x2.sub(x1) + } else { + t.copy(p) + t.sub(x1) + x2.add(t) + } + x2.norm() + } + } + if comp(u, one) == 0 { + r.copy(x1) + } else { + r.copy(x2) + } +} + +/* +func main() { + a := NewBIGint(3) + m := NewBIGints(Modulus) + + fmt.Printf("Modulus= "+m.toString()) + fmt.Printf("\n") + + + e := NewBIGcopy(m); + e.dec(7); e.norm(); + fmt.Printf("Exponent= "+e.toString()) + fmt.Printf("\n") + a=a.powmod(e,m); + fmt.Printf("Result= "+a.toString()) +} +*/
http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/DBIG.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/DBIG.go b/go/amcl-go/DBIG.go new file mode 100644 index 0000000..98314b6 --- /dev/null +++ b/go/amcl-go/DBIG.go @@ -0,0 +1,260 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +/* AMCL double length DBIG number class */ + +package amcl + +import "strconv" + +type DBIG struct { + w [2 * NLEN]int64 +} + +func NewDBIG() *DBIG { + b := new(DBIG) + for i := 0; i < DNLEN; i++ { + b.w[i] = 0 + } + return b +} + +func NewDBIGcopy(x *DBIG) *DBIG { + b := new(DBIG) + for i := 0; i < DNLEN; i++ { + b.w[i] = x.w[i] + } + return b +} + +func NewDBIGscopy(x *BIG) *DBIG { + b := new(DBIG) + for i := 0; i < NLEN-1; i++ { + b.w[i] = x.w[i] + } + b.w[NLEN-1] = x.get(NLEN-1) & MASK /* top word normalized */ + b.w[NLEN] = x.get(NLEN-1) >> BASEBITS + + for i := NLEN + 1; i < DNLEN; i++ { + b.w[i] = 0 + } + return b +} + +/* set this[i]+=x*y+c, and return high part */ + +func (r *DBIG) muladd(a int64, b int64, c int64, i int) int64 { + x0 := a & HMASK + x1 := (a >> HBITS) + y0 := b & HMASK + y1 := (b >> HBITS) + bot := x0 * y0 + top := x1 * y1 + mid := x0*y1 + x1*y0 + x0 = mid & HMASK + x1 = (mid >> HBITS) + bot += x0 << HBITS + bot += c + bot += r.w[i] + top += x1 + carry := bot >> BASEBITS + bot &= MASK + top += carry + r.w[i] = bot + return top +} + +/* normalise this */ +func (r *DBIG) norm() { + var carry int64 = 0 + for i := 0; i < DNLEN-1; i++ { + d := r.w[i] + carry + r.w[i] = d & MASK + carry = d >> BASEBITS + } + r.w[DNLEN-1] = (r.w[DNLEN-1] + carry) +} + +/* split DBIG at position n, return higher half, keep lower half */ +func (r *DBIG) split(n uint) *BIG { + t := NewBIG() + m := n % BASEBITS + carry := r.w[DNLEN-1] << (BASEBITS - m) + + for i := DNLEN - 2; i >= NLEN-1; i-- { + nw := (r.w[i] >> m) | carry + carry = (r.w[i] << (BASEBITS - m)) & MASK + t.set(i-NLEN+1, nw) + } + r.w[NLEN-1] &= ((int64(1) << m) - 1) + return t +} + +/* Compare a and b, return 0 if a==b, -1 if a<b, +1 if a>b. Inputs must be normalised */ +func dcomp(a *DBIG, b *DBIG) int { + for i := DNLEN - 1; i >= 0; i-- { + if a.w[i] == b.w[i] { + continue + } + if a.w[i] > b.w[i] { + return 1 + } else { + return -1 + } + } + return 0 +} + +func (r *DBIG) add(x *DBIG) { + for i := 0; i < DNLEN; i++ { + r.w[i] = r.w[i] + x.w[i] + } +} + +/* this-=x */ +func (r *DBIG) sub(x *DBIG) { + for i := 0; i < DNLEN; i++ { + r.w[i] = r.w[i] - x.w[i] + } +} + +/* general shift left */ +func (r *DBIG) shl(k uint) { + n := k % BASEBITS + m := int(k / BASEBITS) + + r.w[DNLEN-1] = (r.w[DNLEN-1-m] << n) | (r.w[DNLEN-m-2] >> (BASEBITS - n)) + for i := DNLEN - 2; i > m; i-- { + r.w[i] = ((r.w[i-m] << n) & MASK) | (r.w[i-m-1] >> (BASEBITS - n)) + } + r.w[m] = (r.w[0] << n) & MASK + for i := 0; i < m; i++ { + r.w[i] = 0 + } +} + +/* general shift right */ +func (r *DBIG) shr(k uint) { + n := (k % BASEBITS) + m := int(k / BASEBITS) + for i := 0; i < DNLEN-m-1; i++ { + r.w[i] = (r.w[m+i] >> n) | ((r.w[m+i+1] << (BASEBITS - n)) & MASK) + } + r.w[DNLEN-m-1] = r.w[DNLEN-1] >> n + for i := DNLEN - m; i < DNLEN; i++ { + r.w[i] = 0 + } +} + +/* reduces this DBIG mod a BIG, and returns the BIG */ +func (r *DBIG) mod(c *BIG) *BIG { + r.norm() + m := NewDBIGscopy(c) + + if dcomp(r, m) < 0 { + return NewBIGdcopy(r) + } + + m.shl(1) + k := 1 + + for dcomp(r, m) >= 0 { + m.shl(1) + k++ + } + + for k > 0 { + m.shr(1) + if dcomp(r, m) >= 0 { + r.sub(m) + r.norm() + } + k-- + } + return NewBIGdcopy(r) +} + +/* return this/c */ +func (r *DBIG) div(c *BIG) *BIG { + k := 0 + m := NewDBIGscopy(c) + a := NewBIGint(0) + e := NewBIGint(1) + r.norm() + + for dcomp(r, m) >= 0 { + e.fshl(1) + m.shl(1) + k++ + } + + for k > 0 { + m.shr(1) + e.shr(1) + if dcomp(r, m) > 0 { + a.add(e) + a.norm() + r.sub(m) + r.norm() + } + k-- + } + return a +} + +/* Convert to Hex String */ +func (r *DBIG) toString() string { + s := "" + len := r.nbits() + + if len%4 == 0 { + len /= 4 + } else { + len /= 4 + len++ + + } + + for i := len - 1; i >= 0; i-- { + b := NewDBIGcopy(r) + + b.shr(uint(i * 4)) + s += strconv.FormatInt(b.w[0]&15, 16) + } + return s +} + +/* return number of bits */ +func (r *DBIG) nbits() int { + k := DNLEN - 1 + r.norm() + for k >= 0 && r.w[k] == 0 { + k-- + } + if k < 0 { + return 0 + } + bts := int(BASEBITS) * k + c := r.w[k] + for c != 0 { + c /= 2 + bts++ + } + return bts +} http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/ECDH.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/ECDH.go b/go/amcl-go/ECDH.go new file mode 100644 index 0000000..20718eb --- /dev/null +++ b/go/amcl-go/ECDH.go @@ -0,0 +1,657 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +/* Elliptic Curve API high-level functions */ + +package amcl + +import "fmt" + +const ECDH_INVALID_PUBLIC_KEY int = -2 +const ECDH_ERROR int = -3 +const ECDH_INVALID int = -4 +const ECDH_EFS int = int(MODBYTES) +const ECDH_EGS int = int(MODBYTES) +const ECDH_EAS int = 16 +const ECDH_EBS int = 16 + +/* Convert Integer to n-byte array */ +func inttoBytes(n int, len int) []byte { + var b []byte + var i int + for i = 0; i < len; i++ { + b = append(b, 0) + } + i = len + for n > 0 && i > 0 { + i-- + b[i] = byte(n & 0xff) + n /= 256 + } + return b +} + +/* Key Derivation Functions */ +/* Input octet Z */ +/* Output key of length olen */ +func KDF1(Z []byte, olen int) []byte { + /* NOTE: the parameter olen is the length of the output K in bytes */ + H := NewHASH() + hlen := 32 + var K []byte + k := 0 + + for i := 0; i < olen; i++ { + K = append(K, 0) + } + + cthreshold := olen / hlen + if olen%hlen != 0 { + cthreshold++ + } + + for counter := 0; counter < cthreshold; counter++ { + H.Process_array(Z) + if counter > 0 { + H.Process_num(int32(counter)) + } + B := H.Hash() + if k+hlen > olen { + for i := 0; i < olen%hlen; i++ { + K[k] = B[i] + k++ + } + } else { + for i := 0; i < hlen; i++ { + K[k] = B[i] + k++ + } + } + } + return K +} + +func KDF2(Z []byte, P []byte, olen int) []byte { + /* NOTE: the parameter olen is the length of the output k in bytes */ + H := NewHASH() + hlen := 32 + var K []byte + + k := 0 + + for i := 0; i < olen; i++ { + K = append(K, 0) + } + + cthreshold := olen / hlen + if olen%hlen != 0 { + cthreshold++ + } + + for counter := 1; counter <= cthreshold; counter++ { + H.Process_array(Z) + H.Process_num(int32(counter)) + H.Process_array(P) + B := H.Hash() + if k+hlen > olen { + for i := 0; i < olen%hlen; i++ { + K[k] = B[i] + k++ + } + } else { + for i := 0; i < hlen; i++ { + K[k] = B[i] + k++ + } + } + } + return K +} + +/* Password based Key Derivation Function */ +/* Input password p, salt s, and repeat count */ +/* Output key of length olen */ +func PBKDF2(Pass []byte, Salt []byte, rep int, olen int) []byte { + d := olen / 32 + if olen%32 != 0 { + d++ + } + var F [ECDH_EFS]byte + var U [ECDH_EFS]byte + + var S []byte + + //byte[] S=new byte[Salt.length+4]; + + var K []byte + //byte[] K=new byte[d*EFS]; + //opt:=0 + + for i := 1; i <= d; i++ { + for j := 0; j < len(Salt); j++ { + S = append(S, Salt[j]) + } + N := inttoBytes(i, 4) + for j := 0; j < 4; j++ { + S = append(S, N[j]) + } + + HMAC(S, Pass, F[:]) + + for j := 0; j < ECDH_EFS; j++ { + U[j] = F[j] + } + for j := 2; j <= rep; j++ { + HMAC(U[:], Pass, U[:]) + for k := 0; k < ECDH_EFS; k++ { + F[k] ^= U[k] + } + } + for j := 0; j < ECDH_EFS; j++ { + K = append(K, F[j]) + } + } + var key []byte + for i := 0; i < olen; i++ { + key = append(key, K[i]) + } + return key +} + +/* Calculate HMAC of m using key k. HMAC is tag of length olen */ +func HMAC(M []byte, K []byte, tag []byte) int { + /* Input is from an octet m * + * olen is requested output length in bytes. k is the key * + * The output is the calculated tag */ + var B [32]byte + var K0 [64]byte + olen := len(tag) + + b := len(K0) + if olen < 4 || olen > 32 { + return 0 + } + + for i := 0; i < b; i++ { + K0[i] = 0 + } + + H := NewHASH() + + if len(K) > b { + H.Process_array(K) + B = H.Hash() + for i := 0; i < 32; i++ { + K0[i] = B[i] + } + } else { + for i := 0; i < len(K); i++ { + K0[i] = K[i] + } + } + + for i := 0; i < b; i++ { + K0[i] ^= 0x36 + } + H.Process_array(K0[:]) + H.Process_array(M) + B = H.Hash() + + for i := 0; i < b; i++ { + K0[i] ^= 0x6a + } + H.Process_array(K0[:]) + H.Process_array(B[:]) + B = H.Hash() + + for i := 0; i < olen; i++ { + tag[i] = B[i] + } + + return 1 +} + +/* AES encryption/decryption. Encrypt byte array M using key K and returns ciphertext */ +func AES_CBC_IV0_ENCRYPT(K []byte, M []byte) []byte { /* AES CBC encryption, with Null IV and key K */ + /* Input is from an octet string M, output is to an octet string C */ + /* Input is padded as necessary to make up a full final block */ + a := NewAES() + fin := false + + var buff [16]byte + var C []byte + + a.Init(aes_CBC, K, nil) + + ipt := 0 //opt:=0 + var i int + for true { + for i = 0; i < 16; i++ { + if ipt < len(M) { + buff[i] = M[ipt] + ipt++ + } else { + fin = true + break + } + } + if fin { + break + } + a.Encrypt(buff[:]) + for i = 0; i < 16; i++ { + C = append(C, buff[i]) + } + } + + /* last block, filled up to i-th index */ + + padlen := 16 - i + for j := i; j < 16; j++ { + buff[j] = byte(padlen) + } + + a.Encrypt(buff[:]) + + for i = 0; i < 16; i++ { + C = append(C, buff[i]) + } + a.End() + return C +} + +/* returns plaintext if all consistent, else returns null string */ +func AES_CBC_IV0_DECRYPT(K []byte, C []byte) []byte { /* padding is removed */ + a := NewAES() + var buff [16]byte + var MM []byte + var M []byte + + var i int + ipt := 0 + opt := 0 + + a.Init(aes_CBC, K, nil) + + if len(C) == 0 { + return nil + } + ch := C[ipt] + ipt++ + + fin := false + + for true { + for i = 0; i < 16; i++ { + buff[i] = ch + if ipt >= len(C) { + fin = true + break + } else { + ch = C[ipt] + ipt++ + } + } + a.Decrypt(buff[:]) + if fin { + break + } + for i = 0; i < 16; i++ { + MM = append(MM, buff[i]) + opt++ + } + } + + a.End() + bad := false + padlen := int(buff[15]) + if i != 15 || padlen < 1 || padlen > 16 { + bad = true + } + if padlen >= 2 && padlen <= 16 { + for i = 16 - padlen; i < 16; i++ { + if buff[i] != byte(padlen) { + bad = true + } + } + } + + if !bad { + for i = 0; i < 16-padlen; i++ { + MM = append(MM, buff[i]) + opt++ + } + } + + if bad { + return nil + } + + for i = 0; i < opt; i++ { + M = append(M, MM[i]) + } + + return M +} + +/* Calculate a public/private EC GF(p) key pair W,S where W=S.G mod EC(p), + * where S is the secret key and W is the public key + * and G is fixed generator. + * If RNG is NULL then the private key is provided externally in S + * otherwise it is generated randomly internally */ +func ECDH_KEY_PAIR_GENERATE(RNG *RAND, S []byte, W []byte) int { + res := 0 + var T [ECDH_EFS]byte + var s *BIG + var G *ECP + + gx := NewBIGints(CURVE_Gx) + if CURVETYPE != MONTGOMERY { + gy := NewBIGints(CURVE_Gy) + G = NewECPbigs(gx, gy) + } else { + G = NewECPbig(gx) + } + + r := NewBIGints(CURVE_Order) + + if RNG == nil { + s = fromBytes(S) + } else { + s = randomnum(r, RNG) + + s.toBytes(T[:]) + for i := 0; i < ECDH_EGS; i++ { + S[i] = T[i] + } + } + + WP := G.mul(s) + + WP.toBytes(W) + + return res +} + +/* validate public key. Set full=true for fuller check */ +func ECDH_PUBLIC_KEY_VALIDATE(full bool, W []byte) int { + WP := ECP_fromBytes(W) + res := 0 + + r := NewBIGints(CURVE_Order) + + if WP.is_infinity() { + res = ECDH_INVALID_PUBLIC_KEY + } + if res == 0 && full { + WP = WP.mul(r) + if !WP.is_infinity() { + res = ECDH_INVALID_PUBLIC_KEY + } + } + return res +} + +/* IEEE-1363 Diffie-Hellman online calculation Z=S.WD */ +func ECPSVDP_DH(S []byte, WD []byte, Z []byte) int { + res := 0 + var T [ECDH_EFS]byte + + s := fromBytes(S) + + W := ECP_fromBytes(WD) + if W.is_infinity() { + res = ECDH_ERROR + } + + if res == 0 { + r := NewBIGints(CURVE_Order) + s.mod(r) + W = W.mul(s) + if W.is_infinity() { + res = ECDH_ERROR + } else { + W.getX().toBytes(T[:]) + for i := 0; i < ECDH_EFS; i++ { + Z[i] = T[i] + } + } + } + return res +} + +/* IEEE ECDSA Signature, C and D are signature on F using private key S */ +func ECPSP_DSA(RNG *RAND, S []byte, F []byte, C []byte, D []byte) int { + var T [ECDH_EFS]byte + + H := NewHASH() + H.Process_array(F) + B := H.Hash() + + gx := NewBIGints(CURVE_Gx) + gy := NewBIGints(CURVE_Gy) + + G := NewECPbigs(gx, gy) + r := NewBIGints(CURVE_Order) + + s := fromBytes(S) + f := fromBytes(B[:]) + + c := NewBIGint(0) + d := NewBIGint(0) + V := NewECP() + + for d.iszilch() { + u := randomnum(r, RNG) + + V.copy(G) + V = V.mul(u) + vx := V.getX() + c.copy(vx) + c.mod(r) + if c.iszilch() { + continue + } + u.invmodp(r) + d.copy(modmul(s, c, r)) + d.add(f) + d.copy(modmul(u, d, r)) + } + + c.toBytes(T[:]) + for i := 0; i < ECDH_EFS; i++ { + C[i] = T[i] + } + d.toBytes(T[:]) + for i := 0; i < ECDH_EFS; i++ { + D[i] = T[i] + } + return 0 +} + +/* IEEE1363 ECDSA Signature Verification. Signature C and D on F is verified using public key W */ +func ECPVP_DSA(W []byte, F []byte, C []byte, D []byte) int { + res := 0 + + H := NewHASH() + H.Process_array(F) + B := H.Hash() + + gx := NewBIGints(CURVE_Gx) + gy := NewBIGints(CURVE_Gy) + + G := NewECPbigs(gx, gy) + r := NewBIGints(CURVE_Order) + + c := fromBytes(C) + d := fromBytes(D) + f := fromBytes(B[:]) + + if c.iszilch() || comp(c, r) >= 0 || d.iszilch() || comp(d, r) >= 0 { + res = ECDH_INVALID + } + + if res == 0 { + d.invmodp(r) + f.copy(modmul(f, d, r)) + h2 := modmul(c, d, r) + + WP := ECP_fromBytes(W) + if WP.is_infinity() { + res = ECDH_ERROR + } else { + P := NewECP() + P.copy(WP) + + P = P.mul2(h2, G, f) + + if P.is_infinity() { + res = ECDH_INVALID + } else { + d = P.getX() + d.mod(r) + + if comp(d, c) != 0 { + res = ECDH_INVALID + } + } + } + } + + return res +} + +/* IEEE1363 ECIES encryption. Encryption of plaintext M uses public key W and produces ciphertext V,C,T */ +func ECIES_ENCRYPT(P1 []byte, P2 []byte, RNG *RAND, W []byte, M []byte, V []byte, T []byte) []byte { + var Z [ECDH_EFS]byte + var VZ [3*ECDH_EFS + 1]byte + var K1 [ECDH_EAS]byte + var K2 [ECDH_EAS]byte + var U [ECDH_EGS]byte + + if ECDH_KEY_PAIR_GENERATE(RNG, U[:], V) != 0 { + return nil + } + if ECPSVDP_DH(U[:], W, Z[:]) != 0 { + return nil + } + + for i := 0; i < 2*ECDH_EFS+1; i++ { + VZ[i] = V[i] + } + for i := 0; i < ECDH_EFS; i++ { + VZ[2*ECDH_EFS+1+i] = Z[i] + } + + K := KDF2(VZ[:], P1, ECDH_EFS) + + for i := 0; i < ECDH_EAS; i++ { + K1[i] = K[i] + K2[i] = K[ECDH_EAS+i] + } + + C := AES_CBC_IV0_ENCRYPT(K1[:], M) + + L2 := inttoBytes(len(P2), 8) + + var AC []byte + + for i := 0; i < len(C); i++ { + AC = append(AC, C[i]) + } + for i := 0; i < len(P2); i++ { + AC = append(AC, P2[i]) + } + for i := 0; i < 8; i++ { + AC = append(AC, L2[i]) + } + + HMAC(AC, K2[:], T) + + return C +} + +/* IEEE1363 ECIES decryption. Decryption of ciphertext V,C,T using private key U outputs plaintext M */ +func ECIES_DECRYPT(P1 []byte, P2 []byte, V []byte, C []byte, T []byte, U []byte) []byte { + var Z [ECDH_EFS]byte + var VZ [3*ECDH_EFS + 1]byte + var K1 [ECDH_EAS]byte + var K2 [ECDH_EAS]byte + + var TAG []byte = T[:] + + if ECPSVDP_DH(U, V, Z[:]) != 0 { + return nil + } + + for i := 0; i < 2*ECDH_EFS+1; i++ { + VZ[i] = V[i] + } + for i := 0; i < ECDH_EFS; i++ { + VZ[2*ECDH_EFS+1+i] = Z[i] + } + + K := KDF2(VZ[:], P1, ECDH_EFS) + + for i := 0; i < ECDH_EAS; i++ { + K1[i] = K[i] + K2[i] = K[ECDH_EAS+i] + } + + M := AES_CBC_IV0_DECRYPT(K1[:], C) + + if M == nil { + return nil + } + + L2 := inttoBytes(len(P2), 8) + + var AC []byte + + for i := 0; i < len(C); i++ { + AC = append(AC, C[i]) + } + for i := 0; i < len(P2); i++ { + AC = append(AC, P2[i]) + } + for i := 0; i < 8; i++ { + AC = append(AC, L2[i]) + } + + HMAC(AC, K2[:], TAG) + + same := true + for i := 0; i < len(T); i++ { + if T[i] != TAG[i] { + same = false + } + } + if !same { + return nil + } + + return M +} + +func ECDH_printBinary(array []byte) { + for i := 0; i < len(array); i++ { + fmt.Printf("%02x", array[i]) + } + fmt.Printf("\n") +} http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/ECP.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/ECP.go b/go/amcl-go/ECP.go new file mode 100644 index 0000000..3ed1d04 --- /dev/null +++ b/go/amcl-go/ECP.go @@ -0,0 +1,1076 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +package amcl + +//import "fmt" + +/* Elliptic Curve Point Structure */ + +type ECP struct { + x *FP + y *FP + z *FP + INF bool +} + +/* Constructors */ +func NewECP() *ECP { + E := new(ECP) + E.x = NewFPint(0) + E.y = NewFPint(0) + E.z = NewFPint(0) + E.INF = true + return E +} + +/* set (x,y) from two BIGs */ +func NewECPbigs(ix *BIG, iy *BIG) *ECP { + E := new(ECP) + E.x = NewFPbig(ix) + E.y = NewFPbig(iy) + E.z = NewFPint(1) + rhs := RHS(E.x) + + if CURVETYPE == MONTGOMERY { + if rhs.jacobi() == 1 { + E.INF = false + } else { + E.inf() + } + } else { + y2 := NewFPcopy(E.y) + y2.sqr() + if y2.equals(rhs) { + E.INF = false + } else { + E.inf() + } + } + return E +} + +/* set (x,y) from BIG and a bit */ +func NewECPbigint(ix *BIG, s int) *ECP { + E := new(ECP) + E.x = NewFPbig(ix) + E.y = NewFPint(0) + rhs := RHS(E.x) + E.z = NewFPint(1) + if rhs.jacobi() == 1 { + ny := rhs.sqrt() + if ny.redc().parity() != s { + ny.neg() + } + E.y.copy(ny) + E.INF = false + } else { + E.inf() + } + return E +} + +/* set from x - calculate y from curve equation */ +func NewECPbig(ix *BIG) *ECP { + E := new(ECP) + E.x = NewFPbig(ix) + E.y = NewFPint(0) + rhs := RHS(E.x) + E.z = NewFPint(1) + if rhs.jacobi() == 1 { + if CURVETYPE != MONTGOMERY { + E.y.copy(rhs.sqrt()) + } + E.INF = false + } else { + E.INF = true + } + return E +} + +/* test for O point-at-infinity */ +func (E *ECP) is_infinity() bool { + if CURVETYPE == EDWARDS { + E.x.reduce() + E.y.reduce() + E.z.reduce() + return (E.x.iszilch() && E.y.equals(E.z)) + } else { + return E.INF + } +} + +/* Conditional swap of P and Q dependant on d */ +func (E *ECP) cswap(Q *ECP, d int32) { + E.x.cswap(Q.x, d) + if CURVETYPE != MONTGOMERY { + E.y.cswap(Q.y, d) + } + E.z.cswap(Q.z, d) + if CURVETYPE != EDWARDS { + bd := true + if d == 0 { + bd = false + } + bd = bd && (E.INF != Q.INF) + E.INF = (bd != E.INF) + Q.INF = (bd != Q.INF) + } +} + +/* Conditional move of Q to P dependant on d */ +func (E *ECP) cmove(Q *ECP, d int32) { + E.x.cmove(Q.x, d) + if CURVETYPE != MONTGOMERY { + E.y.cmove(Q.y, d) + } + E.z.cmove(Q.z, d) + if CURVETYPE != EDWARDS { + bd := true + if d == 0 { + bd = false + } + E.INF = (E.INF != ((E.INF != Q.INF) && bd)) + } +} + +/* return 1 if b==c, no branching */ +func teq(b int32, c int32) int32 { + x := b ^ c + x -= 1 // if x=0, x now -1 + return ((x >> 31) & 1) +} + +/* this=P */ +func (E *ECP) copy(P *ECP) { + E.x.copy(P.x) + if CURVETYPE != MONTGOMERY { + E.y.copy(P.y) + } + E.z.copy(P.z) + E.INF = P.INF +} + +/* this=-this */ +func (E *ECP) neg() { + if E.is_infinity() { + return + } + if CURVETYPE == WEIERSTRASS { + E.y.neg() + E.y.reduce() + } + if CURVETYPE == EDWARDS { + E.x.neg() + E.x.reduce() + } + return +} + +/* Constant time select from pre-computed table */ +func (E *ECP) selector(W []*ECP, b int32) { + MP := NewECP() + m := b >> 31 + babs := (b ^ m) - m + + babs = (babs - 1) / 2 + + E.cmove(W[0], teq(babs, 0)) // conditional move + E.cmove(W[1], teq(babs, 1)) + E.cmove(W[2], teq(babs, 2)) + E.cmove(W[3], teq(babs, 3)) + E.cmove(W[4], teq(babs, 4)) + E.cmove(W[5], teq(babs, 5)) + E.cmove(W[6], teq(babs, 6)) + E.cmove(W[7], teq(babs, 7)) + + MP.copy(E) + MP.neg() + E.cmove(MP, (m & 1)) +} + +/* set this=O */ +func (E *ECP) inf() { + E.INF = true + E.x.zero() + E.y.one() + E.z.one() +} + +/* Test P == Q */ +func (E *ECP) equals(Q *ECP) bool { + if E.is_infinity() && Q.is_infinity() { + return true + } + if E.is_infinity() || Q.is_infinity() { + return false + } + if CURVETYPE == WEIERSTRASS { + zs2 := NewFPcopy(E.z) + zs2.sqr() + zo2 := NewFPcopy(Q.z) + zo2.sqr() + zs3 := NewFPcopy(zs2) + zs3.mul(E.z) + zo3 := NewFPcopy(zo2) + zo3.mul(Q.z) + zs2.mul(Q.x) + zo2.mul(E.x) + if !zs2.equals(zo2) { + return false + } + zs3.mul(Q.y) + zo3.mul(E.y) + if !zs3.equals(zo3) { + return false + } + } else { + a := NewFPint(0) + b := NewFPint(0) + a.copy(E.x) + a.mul(Q.z) + a.reduce() + b.copy(Q.x) + b.mul(E.z) + b.reduce() + if !a.equals(b) { + return false + } + if CURVETYPE == EDWARDS { + a.copy(E.y) + a.mul(Q.z) + a.reduce() + b.copy(Q.y) + b.mul(E.z) + b.reduce() + if !a.equals(b) { + return false + } + } + } + return true +} + +/* Calculate RHS of curve equation */ +func RHS(x *FP) *FP { + x.norm() + r := NewFPcopy(x) + r.sqr() + + if CURVETYPE == WEIERSTRASS { // x^3+Ax+B + b := NewFPbig(NewBIGints(CURVE_B)) + r.mul(x) + if CURVE_A == -3 { + cx := NewFPcopy(x) + cx.imul(3) + cx.neg() + cx.norm() + r.add(cx) + } + r.add(b) + } + if CURVETYPE == EDWARDS { // (Ax^2-1)/(Bx^2-1) + b := NewFPbig(NewBIGints(CURVE_B)) + + one := NewFPint(1) + b.mul(r) + b.sub(one) + if CURVE_A == -1 { + r.neg() + } + r.sub(one) + b.inverse() + r.mul(b) + } + if CURVETYPE == MONTGOMERY { // x^3+Ax^2+x + x3 := NewFPint(0) + x3.copy(r) + x3.mul(x) + r.imul(CURVE_A) + r.add(x3) + r.add(x) + } + r.reduce() + return r +} + +/* set to affine - from (x,y,z) to (x,y) */ +func (E *ECP) affine() { + if E.is_infinity() { + return + } + one := NewFPint(1) + if E.z.equals(one) { + return + } + E.z.inverse() + if CURVETYPE == WEIERSTRASS { + z2 := NewFPcopy(E.z) + z2.sqr() + E.x.mul(z2) + E.x.reduce() + E.y.mul(z2) + E.y.mul(E.z) + E.y.reduce() + } + if CURVETYPE == EDWARDS { + E.x.mul(E.z) + E.x.reduce() + E.y.mul(E.z) + E.y.reduce() + } + if CURVETYPE == MONTGOMERY { + E.x.mul(E.z) + E.x.reduce() + } + E.z.one() +} + +/* extract x as a BIG */ +func (E *ECP) getX() *BIG { + E.affine() + return E.x.redc() +} + +/* extract y as a BIG */ +func (E *ECP) getY() *BIG { + E.affine() + return E.y.redc() +} + +/* get sign of Y */ +func (E *ECP) getS() int { + E.affine() + y := E.getY() + return y.parity() +} + +/* extract x as an FP */ +func (E *ECP) getx() *FP { + return E.x +} + +/* extract y as an FP */ +func (E *ECP) gety() *FP { + return E.y +} + +/* extract z as an FP */ +func (E *ECP) getz() *FP { + return E.z +} + +/* convert to byte array */ +func (E *ECP) toBytes(b []byte) { + var t [int(MODBYTES)]byte + MB := int(MODBYTES) + if CURVETYPE != MONTGOMERY { + b[0] = 0x04 + } else { + b[0] = 0x02 + } + + E.affine() + E.x.redc().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i+1] = t[i] + } + if CURVETYPE != MONTGOMERY { + E.y.redc().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i+MB+1] = t[i] + } + } +} + +/* convert from byte array to point */ +func ECP_fromBytes(b []byte) *ECP { + var t [int(MODBYTES)]byte + MB := int(MODBYTES) + p := NewBIGints(Modulus) + + for i := 0; i < MB; i++ { + t[i] = b[i+1] + } + px := fromBytes(t[:]) + if comp(px, p) >= 0 { + return NewECP() + } + + if b[0] == 0x04 { + for i := 0; i < MB; i++ { + t[i] = b[i+MB+1] + } + py := fromBytes(t[:]) + if comp(py, p) >= 0 { + return NewECP() + } + return NewECPbigs(px, py) + } else { + return NewECPbig(px) + } +} + +/* convert to hex string */ +func (E *ECP) toString() string { + if E.is_infinity() { + return "infinity" + } + E.affine() + if CURVETYPE == MONTGOMERY { + return "(" + E.x.redc().toString() + ")" + } else { + return "(" + E.x.redc().toString() + "," + E.y.redc().toString() + ")" + } +} + +/* this*=2 */ +func (E *ECP) dbl() { + if CURVETYPE == WEIERSTRASS { + if E.INF { + return + } + if E.y.iszilch() { + E.inf() + return + } + + w1 := NewFPcopy(E.x) + w6 := NewFPcopy(E.z) + w2 := NewFPint(0) + w3 := NewFPcopy(E.x) + w8 := NewFPcopy(E.x) + + if CURVE_A == -3 { + w6.sqr() + w1.copy(w6) + w1.neg() + w3.add(w1) + + w8.add(w6) + + w3.mul(w8) + w8.copy(w3) + w8.imul(3) + } else { + w1.sqr() + w8.copy(w1) + w8.imul(3) + } + + w2.copy(E.y) + w2.sqr() + w3.copy(E.x) + w3.mul(w2) + w3.imul(4) + w1.copy(w3) + w1.neg() + // w1.norm(); + + E.x.copy(w8) + E.x.sqr() + E.x.add(w1) + E.x.add(w1) + // x.reduce(); + E.x.norm() + + E.z.mul(E.y) + E.z.add(E.z) + + w2.add(w2) + w2.sqr() + w2.add(w2) + w3.sub(E.x) + E.y.copy(w8) + E.y.mul(w3) + // w2.norm(); + E.y.sub(w2) + // y.reduce(); + // z.reduce(); + E.y.norm() + E.z.norm() + + } + if CURVETYPE == EDWARDS { + C := NewFPcopy(E.x) + D := NewFPcopy(E.y) + H := NewFPcopy(E.z) + J := NewFPint(0) + + E.x.mul(E.y) + E.x.add(E.x) + C.sqr() + D.sqr() + if CURVE_A == -1 { + C.neg() + } + E.y.copy(C) + E.y.add(D) + // y.norm(); + H.sqr() + H.add(H) + E.z.copy(E.y) + J.copy(E.y) + J.sub(H) + E.x.mul(J) + C.sub(D) + E.y.mul(C) + E.z.mul(J) + + E.x.norm() + E.y.norm() + E.z.norm() + } + if CURVETYPE == MONTGOMERY { + A := NewFPcopy(E.x) + B := NewFPcopy(E.x) + AA := NewFPint(0) + BB := NewFPint(0) + C := NewFPint(0) + + if E.INF { + return + } + + A.add(E.z) + AA.copy(A) + AA.sqr() + B.sub(E.z) + BB.copy(B) + BB.sqr() + C.copy(AA) + C.sub(BB) + // C.norm(); + + E.x.copy(AA) + E.x.mul(BB) + + A.copy(C) + A.imul((CURVE_A + 2) / 4) + + BB.add(A) + E.z.copy(BB) + E.z.mul(C) + // x.reduce(); + // z.reduce(); + E.x.norm() + E.z.norm() + } + return +} + +/* this+=Q */ +func (E *ECP) add(Q *ECP) { + if CURVETYPE == WEIERSTRASS { + if E.INF { + E.copy(Q) + return + } + if Q.INF { + return + } + + aff := false + + one := NewFPint(1) + if Q.z.equals(one) { + aff = true + } + + var A, C *FP + B := NewFPcopy(E.z) + D := NewFPcopy(E.z) + if !aff { + A = NewFPcopy(Q.z) + C = NewFPcopy(Q.z) + + A.sqr() + B.sqr() + C.mul(A) + D.mul(B) + + A.mul(E.x) + C.mul(E.y) + } else { + A = NewFPcopy(E.x) + C = NewFPcopy(E.y) + + B.sqr() + D.mul(B) + } + + B.mul(Q.x) + B.sub(A) + D.mul(Q.y) + D.sub(C) + + if B.iszilch() { + if D.iszilch() { + E.dbl() + return + } else { + E.INF = true + return + } + } + + if !aff { + E.z.mul(Q.z) + } + E.z.mul(B) + + e := NewFPcopy(B) + e.sqr() + B.mul(e) + A.mul(e) + + e.copy(A) + e.add(A) + e.add(B) + E.x.copy(D) + E.x.sqr() + E.x.sub(e) + + A.sub(E.x) + E.y.copy(A) + E.y.mul(D) + C.mul(B) + E.y.sub(C) + + // x.reduce(); + // y.reduce(); + // z.reduce(); + E.x.norm() + E.y.norm() + E.z.norm() + } + if CURVETYPE == EDWARDS { + b := NewFPbig(NewBIGints(CURVE_B)) + A := NewFPcopy(E.z) + B := NewFPint(0) + C := NewFPcopy(E.x) + D := NewFPcopy(E.y) + EE := NewFPint(0) + F := NewFPint(0) + G := NewFPint(0) + //H:=NewFPint(0) + //I:=NewFPint(0) + + A.mul(Q.z) + B.copy(A) + B.sqr() + C.mul(Q.x) + D.mul(Q.y) + + EE.copy(C) + EE.mul(D) + EE.mul(b) + F.copy(B) + F.sub(EE) + G.copy(B) + G.add(EE) + C.add(D) + + if CURVE_A == 1 { + EE.copy(D) + D.sub(C) + } + + B.copy(E.x) + B.add(E.y) + D.copy(Q.x) + D.add(Q.y) + B.mul(D) + B.sub(C) + B.mul(F) + E.x.copy(A) + E.x.mul(B) + + if CURVE_A == 1 { + C.copy(EE) + C.mul(G) + } + if CURVE_A == -1 { + C.mul(G) + } + E.y.copy(A) + E.y.mul(C) + E.z.copy(F) + E.z.mul(G) + // x.reduce(); y.reduce(); z.reduce(); + E.x.norm() + E.y.norm() + E.z.norm() + } + return +} + +/* Differential Add for Montgomery curves. this+=Q where W is this-Q and is affine. */ +func (E *ECP) dadd(Q *ECP, W *ECP) { + A := NewFPcopy(E.x) + B := NewFPcopy(E.x) + C := NewFPcopy(Q.x) + D := NewFPcopy(Q.x) + DA := NewFPint(0) + CB := NewFPint(0) + + A.add(E.z) + B.sub(E.z) + + C.add(Q.z) + D.sub(Q.z) + + DA.copy(D) + DA.mul(A) + CB.copy(C) + CB.mul(B) + + A.copy(DA) + A.add(CB) + A.sqr() + B.copy(DA) + B.sub(CB) + B.sqr() + + E.x.copy(A) + E.z.copy(W.x) + E.z.mul(B) + + if E.z.iszilch() { + E.inf() + } else { + E.INF = false + } + + // x.reduce(); + E.x.norm() +} + +/* this-=Q */ +func (E *ECP) sub(Q *ECP) { + Q.neg() + E.add(Q) + Q.neg() +} + +func multiaffine(m int, P []*ECP) { + t1 := NewFPint(0) + t2 := NewFPint(0) + + var work []*FP + + for i := 0; i < m; i++ { + work = append(work, NewFPint(0)) + } + + work[0].one() + work[1].copy(P[0].z) + + for i := 2; i < m; i++ { + work[i].copy(work[i-1]) + work[i].mul(P[i-1].z) + } + + t1.copy(work[m-1]) + t1.mul(P[m-1].z) + t1.inverse() + t2.copy(P[m-1].z) + work[m-1].mul(t1) + + for i := m - 2; ; i-- { + if i == 0 { + work[0].copy(t1) + work[0].mul(t2) + break + } + work[i].mul(t2) + work[i].mul(t1) + t2.mul(P[i].z) + } + /* now work[] contains inverses of all Z coordinates */ + + for i := 0; i < m; i++ { + P[i].z.one() + t1.copy(work[i]) + t1.sqr() + P[i].x.mul(t1) + t1.mul(work[i]) + P[i].y.mul(t1) + } +} + +/* constant time multiply by small integer of length bts - use ladder */ +func (E *ECP) pinmul(e int32, bts int32) *ECP { + if CURVETYPE == MONTGOMERY { + return E.mul(NewBIGint(int(e))) + } else { + P := NewECP() + R0 := NewECP() + R1 := NewECP() + R1.copy(E) + + for i := bts - 1; i >= 0; i-- { + b := (e >> uint32(i)) & 1 + P.copy(R1) + P.add(R0) + R0.cswap(R1, b) + R1.copy(P) + R0.dbl() + R0.cswap(R1, b) + } + P.copy(R0) + P.affine() + return P + } +} + +/* return e.this */ + +func (E *ECP) mul(e *BIG) *ECP { + if e.iszilch() || E.is_infinity() { + return NewECP() + } + P := NewECP() + if CURVETYPE == MONTGOMERY { + /* use Ladder */ + D := NewECP() + R0 := NewECP() + R0.copy(E) + R1 := NewECP() + R1.copy(E) + R1.dbl() + D.copy(E) + D.affine() + nb := e.nbits() + for i := nb - 2; i >= 0; i-- { + b := int32(e.bit(i)) + P.copy(R1) + P.dadd(R0, D) + R0.cswap(R1, b) + R1.copy(P) + R0.dbl() + R0.cswap(R1, b) + } + P.copy(R0) + } else { + // fixed size windows + mt := NewBIG() + t := NewBIG() + Q := NewECP() + C := NewECP() + + var W []*ECP + var w [1 + (NLEN*int(BASEBITS)+3)/4]int8 + + E.affine() + + Q.copy(E) + Q.dbl() + + W = append(W, NewECP()) + W[0].copy(E) + + for i := 1; i < 8; i++ { + W = append(W, NewECP()) + W[i].copy(W[i-1]) + W[i].add(Q) + } + + // convert the table to affine + if CURVETYPE == WEIERSTRASS { + multiaffine(8, W[:]) + } + + // make exponent odd - add 2P if even, P if odd + t.copy(e) + s := int32(t.parity()) + t.inc(1) + t.norm() + ns := int32(t.parity()) + mt.copy(t) + mt.inc(1) + mt.norm() + t.cmove(mt, s) + Q.cmove(E, ns) + C.copy(Q) + + nb := 1 + (t.nbits()+3)/4 + + // convert exponent to signed 4-bit window + for i := 0; i < nb; i++ { + w[i] = int8(t.lastbits(5) - 16) + t.dec(int(w[i])) + t.norm() + t.fshr(4) + } + w[nb] = int8(t.lastbits(5)) + + P.copy(W[(int(w[nb])-1)/2]) + for i := nb - 1; i >= 0; i-- { + Q.selector(W, int32(w[i])) + P.dbl() + P.dbl() + P.dbl() + P.dbl() + P.add(Q) + } + P.sub(C) /* apply correction */ + } + P.affine() + return P +} + +/* Return e.this+f.Q */ + +func (E *ECP) mul2(e *BIG, Q *ECP, f *BIG) *ECP { + te := NewBIG() + tf := NewBIG() + mt := NewBIG() + S := NewECP() + T := NewECP() + C := NewECP() + var W []*ECP + //ECP[] W=new ECP[8]; + var w [1 + (NLEN*int(BASEBITS)+1)/2]int8 + + E.affine() + Q.affine() + + te.copy(e) + tf.copy(f) + + // precompute table + for i := 0; i < 8; i++ { + W = append(W, NewECP()) + } + W[1].copy(E) + W[1].sub(Q) + W[2].copy(E) + W[2].add(Q) + S.copy(Q) + S.dbl() + W[0].copy(W[1]) + W[0].sub(S) + W[3].copy(W[2]) + W[3].add(S) + T.copy(E) + T.dbl() + W[5].copy(W[1]) + W[5].add(T) + W[6].copy(W[2]) + W[6].add(T) + W[4].copy(W[5]) + W[4].sub(S) + W[7].copy(W[6]) + W[7].add(S) + + // convert the table to affine + if CURVETYPE == WEIERSTRASS { + multiaffine(8, W) + } + + // if multiplier is odd, add 2, else add 1 to multiplier, and add 2P or P to correction + + s := int32(te.parity()) + te.inc(1) + te.norm() + ns := int32(te.parity()) + mt.copy(te) + mt.inc(1) + mt.norm() + te.cmove(mt, s) + T.cmove(E, ns) + C.copy(T) + + s = int32(tf.parity()) + tf.inc(1) + tf.norm() + ns = int32(tf.parity()) + mt.copy(tf) + mt.inc(1) + mt.norm() + tf.cmove(mt, s) + S.cmove(Q, ns) + C.add(S) + + mt.copy(te) + mt.add(tf) + mt.norm() + nb := 1 + (mt.nbits()+1)/2 + + // convert exponent to signed 2-bit window + for i := 0; i < nb; i++ { + a := (te.lastbits(3) - 4) + te.dec(int(a)) + te.norm() + te.fshr(2) + b := (tf.lastbits(3) - 4) + tf.dec(int(b)) + tf.norm() + tf.fshr(2) + w[i] = int8(4*a + b) + } + w[nb] = int8(4*te.lastbits(3) + tf.lastbits(3)) + S.copy(W[(w[nb]-1)/2]) + + for i := nb - 1; i >= 0; i-- { + T.selector(W, int32(w[i])) + S.dbl() + S.dbl() + S.add(T) + } + S.sub(C) /* apply correction */ + S.affine() + return S +} + +/* +func main() { + Gx:=NewBIGints(CURVE_Gx); + var Gy *BIG + var P *ECP + + if CURVETYPE!=MONTGOMERY {Gy=NewBIGints(CURVE_Gy)} + r:=NewBIGints(CURVE_Order) + + //r.dec(7); + + fmt.Printf("Gx= "+Gx.toString()) + fmt.Printf("\n") + + if CURVETYPE!=MONTGOMERY { + fmt.Printf("Gy= "+Gy.toString()) + fmt.Printf("\n") + } + + if CURVETYPE!=MONTGOMERY { + P=NewECPbigs(Gx,Gy) + } else {P=NewECPbig(Gx)} + + fmt.Printf("P= "+P.toString()); + fmt.Printf("\n") + + R:=P.mul(r); + //for (int i=0;i<10000;i++) + // R=P.mul(r); + + fmt.Printf("R= "+R.toString()) + fmt.Printf("\n") +} +*/ http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/ECP2.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/ECP2.go b/go/amcl-go/ECP2.go new file mode 100644 index 0000000..6770378 --- /dev/null +++ b/go/amcl-go/ECP2.go @@ -0,0 +1,672 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +/* AMCL Weierstrass elliptic curve functions over FP2 */ + +package amcl + +//import "fmt" + +type ECP2 struct { + x *FP2 + y *FP2 + z *FP2 + INF bool +} + +func NewECP2() *ECP2 { + E := new(ECP2) + E.x = NewFP2int(0) + E.y = NewFP2int(1) + E.z = NewFP2int(1) + E.INF = true + return E +} + +/* Test this=O? */ +func (E *ECP2) is_infinity() bool { + return E.INF +} + +/* copy this=P */ +func (E *ECP2) copy(P *ECP2) { + E.x.copy(P.x) + E.y.copy(P.y) + E.z.copy(P.z) + E.INF = P.INF +} + +/* set this=O */ +func (E *ECP2) inf() { + E.INF = true + E.x.zero() + E.y.zero() + E.z.zero() +} + +/* set this=-this */ +func (E *ECP2) neg() { + if E.is_infinity() { + return + } + E.y.neg() + E.y.reduce() +} + +/* Conditional move of Q to P dependant on d */ +func (E *ECP2) cmove(Q *ECP2, d int32) { + E.x.cmove(Q.x, d) + E.y.cmove(Q.y, d) + E.z.cmove(Q.z, d) + + var bd bool + if d == 0 { + bd = false + } else { + bd = true + } + E.INF = (E.INF != (E.INF != Q.INF) && bd) +} + +/* Constant time select from pre-computed table */ +func (E *ECP2) selector(W []*ECP2, b int32) { + MP := NewECP2() + m := b >> 31 + babs := (b ^ m) - m + + babs = (babs - 1) / 2 + + E.cmove(W[0], teq(babs, 0)) // conditional move + E.cmove(W[1], teq(babs, 1)) + E.cmove(W[2], teq(babs, 2)) + E.cmove(W[3], teq(babs, 3)) + E.cmove(W[4], teq(babs, 4)) + E.cmove(W[5], teq(babs, 5)) + E.cmove(W[6], teq(babs, 6)) + E.cmove(W[7], teq(babs, 7)) + + MP.copy(E) + MP.neg() + E.cmove(MP, (m & 1)) +} + +/* Test if P == Q */ +func (E *ECP2) equals(Q *ECP2) bool { + if E.is_infinity() && Q.is_infinity() { + return true + } + if E.is_infinity() || Q.is_infinity() { + return false + } + + zs2 := NewFP2copy(E.z) + zs2.sqr() + zo2 := NewFP2copy(Q.z) + zo2.sqr() + zs3 := NewFP2copy(zs2) + zs3.mul(E.z) + zo3 := NewFP2copy(zo2) + zo3.mul(Q.z) + zs2.mul(Q.x) + zo2.mul(E.x) + if !zs2.equals(zo2) { + return false + } + zs3.mul(Q.y) + zo3.mul(E.y) + if !zs3.equals(zo3) { + return false + } + + return true +} + +/* set to Affine - (x,y,z) to (x,y) */ +func (E *ECP2) affine() { + if E.is_infinity() { + return + } + one := NewFP2int(1) + if E.z.equals(one) { + return + } + E.z.inverse() + + z2 := NewFP2copy(E.z) + z2.sqr() + E.x.mul(z2) + E.x.reduce() + E.y.mul(z2) + E.y.mul(E.z) + E.y.reduce() + E.z.copy(one) +} + +/* extract affine x as FP2 */ +func (E *ECP2) getX() *FP2 { + E.affine() + return E.x +} + +/* extract affine y as FP2 */ +func (E *ECP2) getY() *FP2 { + E.affine() + return E.y +} + +/* extract projective x */ +func (E *ECP2) getx() *FP2 { + return E.x +} + +/* extract projective y */ +func (E *ECP2) gety() *FP2 { + return E.y +} + +/* extract projective z */ +func (E *ECP2) getz() *FP2 { + return E.z +} + +/* convert to byte array */ +func (E *ECP2) toBytes(b []byte) { + var t [int(MODBYTES)]byte + MB := int(MODBYTES) + + E.affine() + E.x.getA().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i] = t[i] + } + E.x.getB().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i+MB] = t[i] + } + + E.y.getA().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i+2*MB] = t[i] + } + E.y.getB().toBytes(t[:]) + for i := 0; i < MB; i++ { + b[i+3*MB] = t[i] + } +} + +/* convert from byte array to point */ +func ECP2_fromBytes(b []byte) *ECP2 { + var t [int(MODBYTES)]byte + MB := int(MODBYTES) + + for i := 0; i < MB; i++ { + t[i] = b[i] + } + ra := fromBytes(t[:]) + for i := 0; i < MB; i++ { + t[i] = b[i+MB] + } + rb := fromBytes(t[:]) + rx := NewFP2bigs(ra, rb) + + for i := 0; i < MB; i++ { + t[i] = b[i+2*MB] + } + ra = fromBytes(t[:]) + for i := 0; i < MB; i++ { + t[i] = b[i+3*MB] + } + rb = fromBytes(t[:]) + ry := NewFP2bigs(ra, rb) + + return NewECP2fp2s(rx, ry) +} + +/* convert this to hex string */ +func (E *ECP2) toString() string { + if E.is_infinity() { + return "infinity" + } + E.affine() + return "(" + E.x.toString() + "," + E.y.toString() + ")" +} + +/* Calculate RHS of twisted curve equation x^3+B/i */ +func RHS2(x *FP2) *FP2 { + x.norm() + r := NewFP2copy(x) + r.sqr() + b := NewFP2big(NewBIGints(CURVE_B)) + b.div_ip() + r.mul(x) + r.add(b) + + r.reduce() + return r +} + +/* construct this from (x,y) - but set to O if not on curve */ +func NewECP2fp2s(ix *FP2, iy *FP2) *ECP2 { + E := new(ECP2) + E.x = NewFP2copy(ix) + E.y = NewFP2copy(iy) + E.z = NewFP2int(1) + rhs := RHS2(E.x) + y2 := NewFP2copy(E.y) + y2.sqr() + if y2.equals(rhs) { + E.INF = false + } else { + E.x.zero() + E.INF = true + } + return E +} + +/* construct this from x - but set to O if not on curve */ +func NewECP2fp2(ix *FP2) *ECP2 { + E := new(ECP2) + E.x = NewFP2copy(ix) + E.y = NewFP2int(1) + E.z = NewFP2int(1) + rhs := RHS2(E.x) + if rhs.sqrt() { + E.y.copy(rhs) + E.INF = false + } else { + E.x.zero() + E.INF = true + } + return E +} + +/* this+=this */ +func (E *ECP2) dbl() int { + if E.INF { + return -1 + } + if E.y.iszilch() { + E.inf() + return -1 + } + + w1 := NewFP2copy(E.x) + w2 := NewFP2int(0) + w3 := NewFP2copy(E.x) + w8 := NewFP2copy(E.x) + + w1.sqr() + w8.copy(w1) + w8.imul(3) + + w2.copy(E.y) + w2.sqr() + w3.copy(E.x) + w3.mul(w2) + w3.imul(4) + w1.copy(w3) + w1.neg() + // w1.norm(); + + E.x.copy(w8) + E.x.sqr() + E.x.add(w1) + E.x.add(w1) + E.x.norm() + + E.z.mul(E.y) + E.z.add(E.z) + + w2.add(w2) + w2.sqr() + w2.add(w2) + w3.sub(E.x) + E.y.copy(w8) + E.y.mul(w3) + // w2.norm(); + E.y.sub(w2) + + E.y.norm() + E.z.norm() + + return 1 +} + +/* this+=Q - return 0 for add, 1 for double, -1 for O */ +func (E *ECP2) add(Q *ECP2) int { + if E.INF { + E.copy(Q) + return -1 + } + if Q.INF { + return -1 + } + + aff := false + + if Q.z.isunity() { + aff = true + } + + var A, C *FP2 + B := NewFP2copy(E.z) + D := NewFP2copy(E.z) + if !aff { + A = NewFP2copy(Q.z) + C = NewFP2copy(Q.z) + + A.sqr() + B.sqr() + C.mul(A) + D.mul(B) + + A.mul(E.x) + C.mul(E.y) + } else { + A = NewFP2copy(E.x) + C = NewFP2copy(E.y) + + B.sqr() + D.mul(B) + } + + B.mul(Q.x) + B.sub(A) + D.mul(Q.y) + D.sub(C) + + if B.iszilch() { + if D.iszilch() { + E.dbl() + return 1 + } else { + E.INF = true + return -1 + } + } + + if !aff { + E.z.mul(Q.z) + } + E.z.mul(B) + + e := NewFP2copy(B) + e.sqr() + B.mul(e) + A.mul(e) + + e.copy(A) + e.add(A) + e.add(B) + E.x.copy(D) + E.x.sqr() + E.x.sub(e) + + A.sub(E.x) + E.y.copy(A) + E.y.mul(D) + C.mul(B) + E.y.sub(C) + + E.x.norm() + E.y.norm() + E.z.norm() + + return 0 +} + +/* set this-=Q */ +func (E *ECP2) sub(Q *ECP2) int { + Q.neg() + D := E.add(Q) + Q.neg() + return D +} + +/* set this*=q, where q is Modulus, using Frobenius */ +func (E *ECP2) frob(X *FP2) { + if E.INF { + return + } + X2 := NewFP2copy(X) + X2.sqr() + E.x.conj() + E.y.conj() + E.z.conj() + E.z.reduce() + E.x.mul(X2) + E.y.mul(X2) + E.y.mul(X) +} + +/* normalises m-array of ECP2 points. Requires work vector of m FP2s */ + +func multiaffine2(m int, P []*ECP2) { + t1 := NewFP2int(0) + t2 := NewFP2int(0) + + var work []*FP2 + + for i := 0; i < m; i++ { + work = append(work, NewFP2int(0)) + } + + work[0].one() + work[1].copy(P[0].z) + + for i := 2; i < m; i++ { + work[i].copy(work[i-1]) + work[i].mul(P[i-1].z) + } + + t1.copy(work[m-1]) + t1.mul(P[m-1].z) + + t1.inverse() + + t2.copy(P[m-1].z) + work[m-1].mul(t1) + + for i := m - 2; ; i-- { + if i == 0 { + work[0].copy(t1) + work[0].mul(t2) + break + } + work[i].mul(t2) + work[i].mul(t1) + t2.mul(P[i].z) + } + /* now work[] contains inverses of all Z coordinates */ + + for i := 0; i < m; i++ { + P[i].z.one() + t1.copy(work[i]) + t1.sqr() + P[i].x.mul(t1) + t1.mul(work[i]) + P[i].y.mul(t1) + } +} + +/* P*=e */ +func (E *ECP2) mul(e *BIG) *ECP2 { + /* fixed size windows */ + mt := NewBIG() + t := NewBIG() + P := NewECP2() + Q := NewECP2() + C := NewECP2() + + if E.is_infinity() { + return NewECP2() + } + + var W []*ECP2 + var w [1 + (NLEN*int(BASEBITS)+3)/4]int8 + + E.affine() + + /* precompute table */ + Q.copy(E) + Q.dbl() + + W = append(W, NewECP2()) + W[0].copy(E) + + for i := 1; i < 8; i++ { + W = append(W, NewECP2()) + W[i].copy(W[i-1]) + W[i].add(Q) + } + + /* convert the table to affine */ + + multiaffine2(8, W[:]) + + /* make exponent odd - add 2P if even, P if odd */ + t.copy(e) + s := int32(t.parity()) + t.inc(1) + t.norm() + ns := int32(t.parity()) + mt.copy(t) + mt.inc(1) + mt.norm() + t.cmove(mt, s) + Q.cmove(E, ns) + C.copy(Q) + + nb := 1 + (t.nbits()+3)/4 + /* convert exponent to signed 4-bit window */ + for i := 0; i < nb; i++ { + w[i] = int8(t.lastbits(5) - 16) + t.dec(int(w[i])) + t.norm() + t.fshr(4) + } + w[nb] = int8(t.lastbits(5)) + + P.copy(W[(w[nb]-1)/2]) + for i := nb - 1; i >= 0; i-- { + Q.selector(W, int32(w[i])) + P.dbl() + P.dbl() + P.dbl() + P.dbl() + P.add(Q) + } + P.sub(C) + P.affine() + return P +} + +/* P=u0.Q0+u1*Q1+u2*Q2+u3*Q3 */ +func mul4(Q []*ECP2, u []*BIG) *ECP2 { + var a [4]int8 + T := NewECP2() + C := NewECP2() + P := NewECP2() + + var W []*ECP2 + + mt := NewBIG() + var t []*BIG + + var w [NLEN*int(BASEBITS) + 1]int8 + + for i := 0; i < 4; i++ { + t = append(t, NewBIGcopy(u[i])) + Q[i].affine() + } + + /* precompute table */ + + W = append(W, NewECP2()) + W[0].copy(Q[0]) + W[0].sub(Q[1]) + W = append(W, NewECP2()) + W[1].copy(W[0]) + W = append(W, NewECP2()) + W[2].copy(W[0]) + W = append(W, NewECP2()) + W[3].copy(W[0]) + W = append(W, NewECP2()) + W[4].copy(Q[0]) + W[4].add(Q[1]) + W = append(W, NewECP2()) + W[5].copy(W[4]) + W = append(W, NewECP2()) + W[6].copy(W[4]) + W = append(W, NewECP2()) + W[7].copy(W[4]) + + T.copy(Q[2]) + T.sub(Q[3]) + W[1].sub(T) + W[2].add(T) + W[5].sub(T) + W[6].add(T) + T.copy(Q[2]) + T.add(Q[3]) + W[0].sub(T) + W[3].add(T) + W[4].sub(T) + W[7].add(T) + + multiaffine2(8, W[:]) + + /* if multiplier is even add 1 to multiplier, and add P to correction */ + mt.zero() + C.inf() + for i := 0; i < 4; i++ { + if t[i].parity() == 0 { + t[i].inc(1) + t[i].norm() + C.add(Q[i]) + } + mt.add(t[i]) + mt.norm() + } + + nb := 1 + mt.nbits() + + /* convert exponent to signed 1-bit window */ + for j := 0; j < nb; j++ { + for i := 0; i < 4; i++ { + a[i] = int8(t[i].lastbits(2) - 2) + t[i].dec(int(a[i])) + t[i].norm() + t[i].fshr(1) + } + w[j] = (8*a[0] + 4*a[1] + 2*a[2] + a[3]) + } + w[nb] = int8(8*t[0].lastbits(2) + 4*t[1].lastbits(2) + 2*t[2].lastbits(2) + t[3].lastbits(2)) + + P.copy(W[(w[nb]-1)/2]) + for i := nb - 1; i >= 0; i-- { + T.selector(W, int32(w[i])) + P.dbl() + P.add(T) + } + P.sub(C) /* apply correction */ + + P.affine() + return P +} http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/FF.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/FF.go b/go/amcl-go/FF.go new file mode 100644 index 0000000..9e6e68c --- /dev/null +++ b/go/amcl-go/FF.go @@ -0,0 +1,926 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +package amcl + +//import "fmt" + +const P_MBITS uint = MODBYTES * 8 +const P_MB uint = (P_MBITS % BASEBITS) +const P_OMASK int64 = (int64(-1) << (P_MBITS % BASEBITS)) +const P_FEXCESS int64 = (int64(1) << (BASEBITS*uint(NLEN) - P_MBITS)) +const P_TBITS uint = (P_MBITS % BASEBITS) + +type FF struct { + length int + v []*BIG +} + +func (F *FF) P_EXCESS() int64 { + return ((F.v[F.length-1].get(NLEN-1) & P_OMASK) >> (P_MB)) +} + +/* Constructors */ +func NewFFint(n int) *FF { + F := new(FF) + for i := 0; i < n; i++ { + F.v = append(F.v, NewBIG()) + } + F.length = n + return F +} + +func NewFFints(x [][5]int64, n int) *FF { + F := new(FF) + for i := 0; i < n; i++ { + F.v = append(F.v, NewBIGints(x[i])) + } + F.length = n + return F +} + +/* set to zero */ +func (F *FF) zero() { + for i := 0; i < F.length; i++ { + F.v[i].zero() + } +} + +func (F *FF) getlen() int { + return F.length +} + +/* set to integer */ +func (F *FF) set(m int) { + F.zero() + F.v[0].set(0, int64(m)) +} + +/* copy from FF b */ +func (F *FF) copy(b *FF) { + for i := 0; i < F.length; i++ { + F.v[i].copy(b.v[i]) + } +} + +/* x=y<<n */ +func (F *FF) dsucopy(b *FF) { + for i := 0; i < b.length; i++ { + F.v[b.length+i].copy(b.v[i]) + F.v[i].zero() + } +} + +/* x=y */ +func (F *FF) dscopy(b *FF) { + for i := 0; i < b.length; i++ { + F.v[i].copy(b.v[i]) + F.v[b.length+i].zero() + } +} + +/* x=y>>n */ +func (F *FF) sducopy(b *FF) { + for i := 0; i < F.length; i++ { + F.v[i].copy(b.v[F.length+i]) + } +} + +func (F *FF) one() { + F.v[0].one() + for i := 1; i < F.length; i++ { + F.v[i].zero() + } +} + +/* test equals 0 */ +func (F *FF) iszilch() bool { + for i := 0; i < F.length; i++ { + if !F.v[i].iszilch() { + return false + } + } + return true +} + +/* shift right by 256-bit words */ +func (F *FF) shrw(n int) { + for i := 0; i < n; i++ { + F.v[i].copy(F.v[i+n]) + F.v[i+n].zero() + } +} + +/* shift left by 256-bit words */ +func (F *FF) shlw(n int) { + for i := 0; i < n; i++ { + F.v[n+i].copy(F.v[i]) + F.v[i].zero() + } +} + +/* extract last bit */ +func (F *FF) parity() int { + return F.v[0].parity() +} + +func (F *FF) lastbits(m int) int { + return F.v[0].lastbits(m) +} + +/* compare x and y - must be normalised, and of same length */ +func ff_comp(a *FF, b *FF) int { + for i := a.length - 1; i >= 0; i-- { + j := comp(a.v[i], b.v[i]) + if j != 0 { + return j + } + } + return 0 +} + +/* recursive add */ +func (F *FF) radd(vp int, x *FF, xp int, y *FF, yp int, n int) { + for i := 0; i < n; i++ { + F.v[vp+i].copy(x.v[xp+i]) + F.v[vp+i].add(y.v[yp+i]) + } +} + +/* recursive inc */ +func (F *FF) rinc(vp int, y *FF, yp int, n int) { + for i := 0; i < n; i++ { + F.v[vp+i].add(y.v[yp+i]) + } +} + +/* recursive sub */ +func (F *FF) rsub(vp int, x *FF, xp int, y *FF, yp int, n int) { + for i := 0; i < n; i++ { + F.v[vp+i].copy(x.v[xp+i]) + F.v[vp+i].sub(y.v[yp+i]) + } +} + +/* recursive dec */ +func (F *FF) rdec(vp int, y *FF, yp int, n int) { + for i := 0; i < n; i++ { + F.v[vp+i].sub(y.v[yp+i]) + } +} + +/* simple add */ +func (F *FF) add(b *FF) { + for i := 0; i < F.length; i++ { + F.v[i].add(b.v[i]) + } +} + +/* simple sub */ +func (F *FF) sub(b *FF) { + for i := 0; i < F.length; i++ { + F.v[i].sub(b.v[i]) + } +} + +/* reverse sub */ +func (F *FF) revsub(b *FF) { + for i := 0; i < F.length; i++ { + F.v[i].rsub(b.v[i]) + } +} + +/* normalise - but hold any overflow in top part unless n<0 */ +func (F *FF) rnorm(vp int, n int) { + trunc := false + var carry int64 + if n < 0 { /* -v n signals to do truncation */ + n = -n + trunc = true + } + for i := 0; i < n-1; i++ { + carry = F.v[vp+i].norm() + F.v[vp+i].xortop(carry << P_TBITS) + F.v[vp+i+1].inc(int(carry)) + } + carry = F.v[vp+n-1].norm() + if trunc { + F.v[vp+n-1].xortop(carry << P_TBITS) + } + +} + +func (F *FF) norm() { + F.rnorm(0, F.length) +} + +/* increment/decrement by a small integer */ +func (F *FF) inc(m int) { + F.v[0].inc(m) + F.norm() +} + +func (F *FF) dec(m int) { + F.v[0].dec(m) + F.norm() +} + +/* shift left by one bit */ +func (F *FF) shl() { + var delay_carry int = 0 + for i := 0; i < F.length-1; i++ { + carry := F.v[i].fshl(1) + F.v[i].inc(delay_carry) + F.v[i].xortop(carry << P_TBITS) + delay_carry = int(carry) + } + F.v[F.length-1].fshl(1) + F.v[F.length-1].inc(delay_carry) +} + +/* shift right by one bit */ + +func (F *FF) shr() { + for i := F.length - 1; i > 0; i-- { + carry := F.v[i].fshr(1) + F.v[i-1].ortop(carry << P_TBITS) + } + F.v[0].fshr(1) +} + +/* Convert to Hex String */ +func (F *FF) toString() string { + F.norm() + s := "" + for i := F.length - 1; i >= 0; i-- { + s += F.v[i].toString() + } + return s +} + +/* Convert FFs to/from byte arrays */ +func (F *FF) toBytes(b []byte) { + for i := 0; i < F.length; i++ { + F.v[i].tobytearray(b, (F.length-i-1)*int(MODBYTES)) + } +} + +func ff_fromBytes(x *FF, b []byte) { + for i := 0; i < x.length; i++ { + x.v[i] = frombytearray(b, (x.length-i-1)*int(MODBYTES)) + } +} + +/* in-place swapping using xor - side channel resistant - lengths must be the same */ +func ff_cswap(a *FF, b *FF, d int32) { + for i := 0; i < a.length; i++ { + a.v[i].cswap(b.v[i], d) + } +} + +/* z=x*y, t is workspace */ +func (F *FF) karmul(vp int, x *FF, xp int, y *FF, yp int, t *FF, tp int, n int) { + if n == 1 { + d := mul(x.v[xp], y.v[yp]) + F.v[vp+1] = d.split(8 * MODBYTES) + F.v[vp].dcopy(d) + return + } + nd2 := n / 2 + F.radd(vp, x, xp, x, xp+nd2, nd2) + F.radd(vp+nd2, y, yp, y, yp+nd2, nd2) + t.karmul(tp, F, vp, F, vp+nd2, t, tp+n, nd2) + F.karmul(vp, x, xp, y, yp, t, tp+n, nd2) + F.karmul(vp+n, x, xp+nd2, y, yp+nd2, t, tp+n, nd2) + t.rdec(tp, F, vp, n) + t.rdec(tp, F, vp+n, n) + F.rinc(vp+nd2, t, tp, n) + F.rnorm(vp, 2*n) +} + +func (F *FF) karsqr(vp int, x *FF, xp int, t *FF, tp int, n int) { + if n == 1 { + d := sqr(x.v[xp]) + F.v[vp+1].copy(d.split(8 * MODBYTES)) + F.v[vp].dcopy(d) + return + } + + nd2 := n / 2 + F.karsqr(vp, x, xp, t, tp+n, nd2) + F.karsqr(vp+n, x, xp+nd2, t, tp+n, nd2) + t.karmul(tp, x, xp, x, xp+nd2, t, tp+n, nd2) + F.rinc(vp+nd2, t, tp, n) + F.rinc(vp+nd2, t, tp, n) + F.rnorm(vp+nd2, n) +} + +/* Calculates Least Significant bottom half of x*y */ +func (F *FF) karmul_lower(vp int, x *FF, xp int, y *FF, yp int, t *FF, tp int, n int) { + if n == 1 { /* only calculate bottom half of product */ + F.v[vp].copy(smul(x.v[xp], y.v[yp])) + return + } + nd2 := n / 2 + + F.karmul(vp, x, xp, y, yp, t, tp+n, nd2) + t.karmul_lower(tp, x, xp+nd2, y, yp, t, tp+n, nd2) + F.rinc(vp+nd2, t, tp, nd2) + t.karmul_lower(tp, x, xp, y, yp+nd2, t, tp+n, nd2) + F.rinc(vp+nd2, t, tp, nd2) + F.rnorm(vp+nd2, -nd2) /* truncate it */ +} + +/* Calculates Most Significant upper half of x*y, given lower part */ +func (F *FF) karmul_upper(x *FF, y *FF, t *FF, n int) { + nd2 := n / 2 + F.radd(n, x, 0, x, nd2, nd2) + F.radd(n+nd2, y, 0, y, nd2, nd2) + + t.karmul(0, F, n+nd2, F, n, t, n, nd2) /* t = (a0+a1)(b0+b1) */ + F.karmul(n, x, nd2, y, nd2, t, n, nd2) /* z[n]= a1*b1 */ + /* z[0-nd2]=l(a0b0) z[nd2-n]= h(a0b0)+l(t)-l(a0b0)-l(a1b1) */ + t.rdec(0, F, n, n) /* t=t-a1b1 */ + F.rinc(nd2, F, 0, nd2) /* z[nd2-n]+=l(a0b0) = h(a0b0)+l(t)-l(a1b1) */ + F.rdec(nd2, t, 0, nd2) /* z[nd2-n]=h(a0b0)+l(t)-l(a1b1)-l(t-a1b1)=h(a0b0) */ + F.rnorm(0, -n) /* a0b0 now in z - truncate it */ + t.rdec(0, F, 0, n) /* (a0+a1)(b0+b1) - a0b0 */ + F.rinc(nd2, t, 0, n) + + F.rnorm(nd2, n) +} + +/* z=x*y. Assumes x and y are of same length. */ +func ff_mul(x *FF, y *FF) *FF { + n := x.length + z := NewFFint(2 * n) + t := NewFFint(2 * n) + z.karmul(0, x, 0, y, 0, t, 0, n) + return z +} + +/* return low part of product this*y */ +func (F *FF) lmul(y *FF) { + n := F.length + t := NewFFint(2 * n) + x := NewFFint(n) + x.copy(F) + F.karmul_lower(0, x, 0, y, 0, t, 0, n) +} + +/* Set b=b mod c */ +func (F *FF) mod(c *FF) { + var k int = 1 + + F.norm() + if ff_comp(F, c) < 0 { + return + } + + c.shl() + for ff_comp(F, c) >= 0 { + c.shl() + k++ + } + + for k > 0 { + c.shr() + if ff_comp(F, c) >= 0 { + F.sub(c) + F.norm() + } + k-- + } +} + +/* z=x^2 */ +func ff_sqr(x *FF) *FF { + n := x.length + z := NewFFint(2 * n) + t := NewFFint(2 * n) + z.karsqr(0, x, 0, t, 0, n) + return z +} + +/* return This mod modulus, N is modulus, ND is Montgomery Constant */ +func (F *FF) reduce(N *FF, ND *FF) *FF { /* fast karatsuba Montgomery reduction */ + n := N.length + t := NewFFint(2 * n) + r := NewFFint(n) + m := NewFFint(n) + + r.sducopy(F) + m.karmul_lower(0, F, 0, ND, 0, t, 0, n) + F.karmul_upper(N, m, t, n) + m.sducopy(F) + + r.add(N) + r.sub(m) + r.norm() + + return r + +} + +/* Set r=this mod b */ +/* this is of length - 2*n */ +/* r,b is of length - n */ +func (F *FF) dmod(b *FF) *FF { + n := b.length + m := NewFFint(2 * n) + x := NewFFint(2 * n) + r := NewFFint(n) + + x.copy(F) + x.norm() + m.dsucopy(b) + k := 256 * n + + for k > 0 { + m.shr() + + if ff_comp(x, m) >= 0 { + x.sub(m) + x.norm() + } + k-- + } + + r.copy(x) + r.mod(b) + return r +} + +/* Set return=1/this mod p. Binary method - a<p on entry */ + +func (F *FF) invmodp(p *FF) { + n := p.length + + u := NewFFint(n) + v := NewFFint(n) + x1 := NewFFint(n) + x2 := NewFFint(n) + t := NewFFint(n) + one := NewFFint(n) + + one.one() + u.copy(F) + v.copy(p) + x1.copy(one) + x2.zero() + + // reduce n in here as well! + for ff_comp(u, one) != 0 && ff_comp(v, one) != 0 { + for u.parity() == 0 { + u.shr() + if x1.parity() != 0 { + x1.add(p) + x1.norm() + } + x1.shr() + } + for v.parity() == 0 { + v.shr() + if x2.parity() != 0 { + x2.add(p) + x2.norm() + } + x2.shr() + } + if ff_comp(u, v) >= 0 { + u.sub(v) + u.norm() + if ff_comp(x1, x2) >= 0 { + x1.sub(x2) + } else { + t.copy(p) + t.sub(x2) + x1.add(t) + } + x1.norm() + } else { + v.sub(u) + v.norm() + if ff_comp(x2, x1) >= 0 { + x2.sub(x1) + } else { + t.copy(p) + t.sub(x1) + x2.add(t) + } + x2.norm() + } + } + if ff_comp(u, one) == 0 { + F.copy(x1) + } else { + F.copy(x2) + } +} + +/* nresidue mod m */ +func (F *FF) nres(m *FF) { + n := m.length + d := NewFFint(2 * n) + d.dsucopy(F) + F.copy(d.dmod(m)) +} + +func (F *FF) redc(m *FF, ND *FF) { + n := m.length + d := NewFFint(2 * n) + F.mod(m) + d.dscopy(F) + F.copy(d.reduce(m, ND)) + F.mod(m) +} + +func (F *FF) mod2m(m int) { + for i := m; i < F.length; i++ { + F.v[i].zero() + } +} + +/* U=1/a mod 2^m - Arazi & Qi */ +func (F *FF) invmod2m() *FF { + n := F.length + + b := NewFFint(n) + c := NewFFint(n) + U := NewFFint(n) + + U.zero() + U.v[0].copy(F.v[0]) + U.v[0].invmod2m() + + for i := 1; i < n; i <<= 1 { + b.copy(F) + b.mod2m(i) + t := ff_mul(U, b) + t.shrw(i) + b.copy(t) + c.copy(F) + c.shrw(i) + c.mod2m(i) + c.lmul(U) + c.mod2m(i) + + b.add(c) + b.norm() + b.lmul(U) + b.mod2m(i) + + c.one() + c.shlw(i) + b.revsub(c) + b.norm() + b.shlw(i) + U.add(b) + } + U.norm() + return U +} + +func (F *FF) random(rng *RAND) { + n := F.length + for i := 0; i < n; i++ { + F.v[i].copy(random(rng)) + } + /* make sure top bit is 1 */ + for F.v[n-1].nbits() < int(MODBYTES*8) { + F.v[n-1].copy(random(rng)) + } +} + +/* generate random x less than p */ +func (F *FF) randomnum(p *FF, rng *RAND) { + n := F.length + d := NewFFint(2 * n) + + for i := 0; i < 2*n; i++ { + d.v[i].copy(random(rng)) + } + F.copy(d.dmod(p)) +} + +/* this*=y mod p */ +func (F *FF) modmul(y *FF, p *FF, nd *FF) { + ex := F.P_EXCESS() + ey := y.P_EXCESS() + if (ex+1)*(ey+1)+1 >= P_FEXCESS { + F.mod(p) + } + d := ff_mul(F, y) + F.copy(d.reduce(p, nd)) +} + +/* this*=y mod p */ +func (F *FF) modsqr(p *FF, nd *FF) { + ex := F.P_EXCESS() + if (ex+1)*(ex+1)+1 >= P_FEXCESS { + F.mod(p) + } + d := ff_sqr(F) + F.copy(d.reduce(p, nd)) +} + +/* this=this^e mod p using side-channel resistant Montgomery Ladder, for large e */ +func (F *FF) skpow(e *FF, p *FF) { + n := p.length + R0 := NewFFint(n) + R1 := NewFFint(n) + ND := p.invmod2m() + + F.mod(p) + R0.one() + R1.copy(F) + R0.nres(p) + R1.nres(p) + + for i := int(8*MODBYTES)*n - 1; i >= 0; i-- { + b := int32(e.v[i/256].bit(i % 256)) + F.copy(R0) + F.modmul(R1, p, ND) + + ff_cswap(R0, R1, b) + R0.modsqr(p, ND) + + R1.copy(F) + ff_cswap(R0, R1, b) + } + F.copy(R0) + F.redc(p, ND) +} + +/* this =this^e mod p using side-channel resistant Montgomery Ladder, for short e */ +func (F *FF) skpows(e *BIG, p *FF) { + n := p.length + R0 := NewFFint(n) + R1 := NewFFint(n) + ND := p.invmod2m() + + F.mod(p) + R0.one() + R1.copy(F) + R0.nres(p) + R1.nres(p) + + for i := int(8*MODBYTES) - 1; i >= 0; i-- { + b := int32(e.bit(i)) + F.copy(R0) + F.modmul(R1, p, ND) + + ff_cswap(R0, R1, b) + R0.modsqr(p, ND) + + R1.copy(F) + ff_cswap(R0, R1, b) + } + F.copy(R0) + F.redc(p, ND) +} + +/* raise to an integer power - right-to-left method */ +func (F *FF) power(e int, p *FF) { + n := p.length + w := NewFFint(n) + ND := p.invmod2m() + f := true + + w.copy(F) + w.nres(p) + + if e == 2 { + F.copy(w) + F.modsqr(p, ND) + } else { + for true { + if e%2 == 1 { + if f { + F.copy(w) + } else { + F.modmul(w, p, ND) + } + f = false + } + e >>= 1 + if e == 0 { + break + } + w.modsqr(p, ND) + } + } + F.redc(p, ND) +} + +/* this=this^e mod p, faster but not side channel resistant */ +func (F *FF) pow(e *FF, p *FF) { + n := p.length + w := NewFFint(n) + ND := p.invmod2m() + + w.copy(F) + F.one() + F.nres(p) + w.nres(p) + for i := int(8*MODBYTES)*n - 1; i >= 0; i-- { + F.modsqr(p, ND) + b := e.v[i/256].bit(i % 256) + if b == 1 { + F.modmul(w, p, ND) + } + } + F.redc(p, ND) +} + +/* double exponentiation r=x^e.y^f mod p */ +func (F *FF) pow2(e *BIG, y *FF, f *BIG, p *FF) { + n := p.length + xn := NewFFint(n) + yn := NewFFint(n) + xy := NewFFint(n) + ND := p.invmod2m() + + xn.copy(F) + yn.copy(y) + xn.nres(p) + yn.nres(p) + xy.copy(xn) + xy.modmul(yn, p, ND) + F.one() + F.nres(p) + + for i := int(8*MODBYTES) - 1; i >= 0; i-- { + eb := e.bit(i) + fb := f.bit(i) + F.modsqr(p, ND) + if eb == 1 { + if fb == 1 { + F.modmul(xy, p, ND) + } else { + F.modmul(xn, p, ND) + } + } else { + if fb == 1 { + F.modmul(yn, p, ND) + } + } + } + F.redc(p, ND) +} + +func igcd(x int, y int) int { /* integer GCD, returns GCD of x and y */ + var r int + if y == 0 { + return x + } + for true { + r = x % y + if r == 0 { + break + } + x = y + y = r + } + return y +} + +/* quick and dirty check for common factor with n */ +func (F *FF) cfactor(s int) bool { + n := F.length + + x := NewFFint(n) + y := NewFFint(n) + + y.set(s) + x.copy(F) + x.norm() + + x.sub(y) + x.norm() + + for !x.iszilch() && x.parity() == 0 { + x.shr() + } + + for ff_comp(x, y) > 0 { + x.sub(y) + x.norm() + for !x.iszilch() && x.parity() == 0 { + x.shr() + } + } + + g := int(x.v[0].get(0)) + r := igcd(s, g) + if r > 1 { + return true + } + return false +} + +/* Miller-Rabin test for primality. Slow. */ +func prime(p *FF, rng *RAND) bool { + s := 0 + n := p.length + d := NewFFint(n) + x := NewFFint(n) + unity := NewFFint(n) + nm1 := NewFFint(n) + + sf := 4849845 /* 3*5*.. *19 */ + p.norm() + + if p.cfactor(sf) { + return false + } + unity.one() + nm1.copy(p) + nm1.sub(unity) + nm1.norm() + d.copy(nm1) + + for d.parity() == 0 { + d.shr() + s++ + } + if s == 0 { + return false + } + for i := 0; i < 10; i++ { + x.randomnum(p, rng) + x.pow(d, p) + if ff_comp(x, unity) == 0 || ff_comp(x, nm1) == 0 { + continue + } + loop := false + for j := 1; j < s; j++ { + x.power(2, p) + if ff_comp(x, unity) == 0 { + return false + } + if ff_comp(x, nm1) == 0 { + loop = true + break + } + } + if loop { + continue + } + return false + } + return true +} + +/* +func main() { + + var P = [4][5]int64 {{0xAD19A781670957,0x76A79C00965796,0xDEFCC5FC9A9717,0xF02F2940E20E9,0xBF59E34F},{0x6894F31844C908,0x8DADA70E82C79F,0xFD29F3836046F6,0x8C1D874D314DD0,0x46D077B},{0x3C515217813331,0x56680FD1CE935B,0xE55C53EEA8838E,0x92C2F7E14A4A95,0xD945E5B1},{0xACF673E919F5EF,0x6723E7E7DAB446,0x6B6FA69B36EB1B,0xF7D13920ECA300,0xB5FC2165}} + + fmt.Printf("Testing FF\n") + var raw [100]byte + rng:=NewRAND() + + rng.Clean() + for i:=0;i<100;i++ { + raw[i]=byte(i) + } + + rng.Seed(100,raw[:]) + + n:=4 + + x:=NewFFint(n) + x.set(3) + + p:=NewFFints(P[:],n) + + if prime(p,rng) {fmt.Printf("p is a prime\n"); fmt.Printf("\n")} + + e:=NewFFint(n) + e.copy(p) + e.dec(1); e.norm() + + fmt.Printf("e= "+e.toString()) + fmt.Printf("\n") + x.skpow(e,p) + fmt.Printf("x= "+x.toString()) + fmt.Printf("\n") +} +*/ http://git-wip-us.apache.org/repos/asf/incubator-milagro-crypto/blob/85fabaa6/go/amcl-go/FP.go ---------------------------------------------------------------------- diff --git a/go/amcl-go/FP.go b/go/amcl-go/FP.go new file mode 100644 index 0000000..c8a4d62 --- /dev/null +++ b/go/amcl-go/FP.go @@ -0,0 +1,288 @@ +/* +Licensed to the Apache Software Foundation (ASF) under one +or more contributor license agreements. See the NOTICE file +distributed with this work for additional information +regarding copyright ownership. The ASF licenses this file +to you under the Apache License, Version 2.0 (the +"License"); you may not use this file except in compliance +with the License. You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, +software distributed under the License is distributed on an +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +KIND, either express or implied. See the License for the +specific language governing permissions and limitations +under the License. +*/ + +/* Finite Field arithmetic */ +/* CLINT mod p functions */ + +package amcl + +//import "fmt" + +var p BIG = BIG{w: [NLEN]int64(Modulus)} + +type FP struct { + x *BIG +} + +/* Constructors */ +func NewFPint(a int) *FP { + F := new(FP) + F.x = NewBIGint(a) + F.nres() + return F +} + +func NewFPbig(a *BIG) *FP { + F := new(FP) + F.x = NewBIGcopy(a) + F.nres() + return F +} + +func NewFPcopy(a *FP) *FP { + F := new(FP) + F.x = NewBIGcopy(a.x) + return F +} + +func (F *FP) toString() string { + return F.redc().toString() +} + +/* convert to Montgomery n-residue form */ +func (F *FP) nres() { + if MODTYPE != PSEUDO_MERSENNE { + d := NewDBIGscopy(F.x) + d.shl(uint(NLEN) * BASEBITS) + F.x.copy(d.mod(&p)) + } +} + +/* convert back to regular form */ +func (F *FP) redc() *BIG { + if MODTYPE != PSEUDO_MERSENNE { + d := NewDBIGscopy(F.x) + return mod(d) + } else { + r := NewBIGcopy(F.x) + return r + } +} + +/* reduce this mod Modulus */ +func (F *FP) reduce() { + F.x.mod(&p) +} + +/* test this=0? */ +func (F *FP) iszilch() bool { + F.reduce() + return F.x.iszilch() +} + +/* copy from FP b */ +func (F *FP) copy(b *FP) { + F.x.copy(b.x) +} + +/* set this=0 */ +func (F *FP) zero() { + F.x.zero() +} + +/* set this=1 */ +func (F *FP) one() { + F.x.one() + F.nres() +} + +/* normalise this */ +func (F *FP) norm() { + F.x.norm() +} + +/* swap FPs depending on d */ +func (F *FP) cswap(b *FP, d int32) { + F.x.cswap(b.x, d) +} + +/* copy FPs depending on d */ +func (F *FP) cmove(b *FP, d int32) { + F.x.cmove(b.x, d) +} + +/* this*=b mod Modulus */ +func (F *FP) mul(b *FP) { + ea := EXCESS(F.x) + eb := EXCESS(b.x) + + if (ea+1)*(eb+1)+1 >= FEXCESS { + F.reduce() + } + + d := mul(F.x, b.x) + F.x.copy(mod(d)) +} + +/* this = -this mod Modulus */ +func (F *FP) neg() { + m := NewBIGcopy(&p) + + F.norm() + + ov := EXCESS(F.x) + sb := uint(1) + for ov != 0 { + sb++ + ov >>= 1 + } + + m.fshl(sb) + F.x.rsub(m) + + if EXCESS(F.x) >= FEXCESS { + F.reduce() + } +} + +/* this*=c mod Modulus, where c is a small int */ +func (F *FP) imul(c int) { + F.norm() + s := false + if c < 0 { + c = -c + s = true + } + afx := (EXCESS(F.x)+1)*(int64(c)+1) + 1 + if c < NEXCESS && afx < FEXCESS { + F.x.imul(c) + } else { + if afx < FEXCESS { + F.x.pmul(c) + } else { + d := F.x.pxmul(c) + F.x.copy(d.mod(&p)) + } + } + if s { + F.neg() + } + F.norm() +} + +/* this*=this mod Modulus */ +func (F *FP) sqr() { + ea := EXCESS(F.x) + if (ea+1)*(ea+1)+1 >= FEXCESS { + F.reduce() + } + + d := sqr(F.x) + + F.x.copy(mod(d)) +} + +/* this+=b */ +func (F *FP) add(b *FP) { + F.x.add(b.x) + if EXCESS(F.x)+2 >= FEXCESS { + F.reduce() + } +} + +/* this-=b */ +func (F *FP) sub(b *FP) { + n := NewFPcopy(b) + n.neg() + F.add(n) +} + +/* this/=2 mod Modulus */ +func (F *FP) div2() { + F.x.norm() + if F.x.parity() == 0 { + F.x.fshr(1) + } else { + F.x.add(&p) + F.x.norm() + F.x.fshr(1) + } +} + +/* this=1/this mod Modulus */ +func (F *FP) inverse() { + r := F.redc() + r.invmodp(&p) + F.x.copy(r) + F.nres() +} + +/* return TRUE if this==a */ +func (F *FP) equals(a *FP) bool { + a.reduce() + F.reduce() + if comp(a.x, F.x) == 0 { + return true + } + return false +} + +/* return this^e mod Modulus */ +func (F *FP) pow(e *BIG) *FP { + r := NewFPint(1) + e.norm() + F.x.norm() + m := NewFPcopy(F) + for true { + bt := e.parity() + e.fshr(1) + if bt == 1 { + r.mul(m) + } + if e.iszilch() { + break + } + m.sqr() + } + r.x.mod(&p) + return r +} + +/* return sqrt(this) mod Modulus */ +func (F *FP) sqrt() *FP { + F.reduce() + b := NewBIGcopy(&p) + if MOD8 == 5 { + b.dec(5) + b.norm() + b.shr(3) + i := NewFPcopy(F) + i.x.shl(1) + v := i.pow(b) + i.mul(v) + i.mul(v) + i.x.dec(1) + r := NewFPcopy(F) + r.mul(v) + r.mul(i) + r.reduce() + return r + } else { + b.inc(1) + b.norm() + b.shr(2) + return F.pow(b) + } +} + +/* return jacobi symbol (this/Modulus) */ +func (F *FP) jacobi() int { + w := F.redc() + return w.jacobi(&p) +}
