# Copyright 2015 Brian Smith.
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND BRIAN SMITH AND THE AUTHORS DISCLAIM
# ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES
# OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL BRIAN SMITH OR THE AUTHORS
# BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY
# DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN
# AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

# Usage: python lib/gfp_generate.py --outdir <dir>

modp_template = """
{
  { %(limbs)s
  },

  %(bits)d / 8 / sizeof(Limb), /* count */
  %(bits)d / 8, /* bytes */
  %(bits)d, /* bits */

  gfp_%(curve_name)s_mul_mont_%(name)s,

  %(k)s, /* k */

  /* rr */
  { %(rr)s
  },
}
"""

nistz_declarations_template = """
extern void ecp_nistz%(bits)d_add(Limb *r, const Limb *a, const Limb *b);
extern void ecp_nistz%(bits)d_mul_mont(Limb *r, const Limb *a, const Limb *b);
extern void ecp_nistz%(bits)d_neg(Limb *r, const Limb *s);
extern void ecp_nistz%(bits)d_sqr_mont(Limb *r, const Limb *s);
extern void ecp_nistz%(bits)d_sub(Limb *r, const Limb *a, const Limb *b);

extern void ecp_nistz%(bits)d_point_double(JacobianPoint *r,
                                           const JacobianPoint *s);
extern void ecp_nistz%(bits)d_point_add(JacobianPoint *r,
                                        const JacobianPoint *a,
                                        const JacobianPoint *b);
extern void ecp_nistz%(bits)d_point_add_affine(JacobianPoint *r,
                                          const JacobianPoint *a,
                                          const XY00Point *b);

#define gfp_%(name)s_add_mod_q ecp_nistz%(bits)d_add
#define gfp_%(name)s_mul_mont_q ecp_nistz%(bits)d_mul_mont
#define gfp_%(name)s_neg_mod_q ecp_nistz%(bits)d_neg
#define gfp_%(name)s_sqr_mont_q ecp_nistz%(bits)d_sqr_mont
#define gfp_%(name)s_sub_mod_q ecp_nistz%(bits)d_sub

#define gfp_%(name)s_point_double ecp_nistz%(bits)d_point_double
#define gfp_%(name)s_point_add_jacobian ecp_nistz%(bits)d_point_add
#define gfp_%(name)s_point_add_xy00 ecp_nistz%(bits)d_point_add_affine

static void gfp_%(name)s_mul_mont_n(Limb *r, const Limb *a, const Limb *b);
"""

nistz_definitions_template = """
static void gfp_%(name)s_mul_mont_n(Limb *r, const Limb *a, const Limb *b) {
  GFp_limbarray_mul_mont(r, gfp_%(name)s.n.count, a, b, gfp_%(name)s.n.count,
                         &gfp_%(name)s.n);
}
"""

curve_template = """
/* Copyright 2015 Brian Smith.
 *
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
 * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
 * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */

/* This entire file was generated by gfp_generate_curve.py. */

#include "gfp_internal.h"

%(declarations)s

static const GFp_Curve gfp_%(name)s = {
  /* q */
  %(q)s,

  /* n */
  %(n)s,

  /* one (Montgomery) */
  { { %(one)s
  } },

  /* three (Montgomery) */
  { { %(three)s
  } },

  /* b (Montgomery) */
  { { %(b)s
  } },

  /* G (Montgomery) */
  { { { %(Gx)s
    } },
    { { %(Gy)s
    } }
  },

  /* q_minus_n */
  { %(q_minus_n)s
  },

  gfp_%(name)s_add_mod_q,
  gfp_%(name)s_neg_mod_q,
  gfp_%(name)s_sqr_mont_q,
  gfp_%(name)s_sub_mod_q,

  gfp_%(name)s_point_double,
  gfp_%(name)s_point_add_jacobian,
  gfp_%(name)s_point_add_xy00,
};

const GFp_Curve *GFp_%(name)s(void) {
  return &gfp_%(name)s;
}

%(definitions)s
"""


import math
import random
import sys

def fill_in_template(template, value_map, indent):
    # In order to enhance readability, the definitions of the templates put the
    # '"""' delimiters on separate lines from the first/last lines of the
    # template, but that results in unwanted blank lines at the beginning and
    # end. |splitlines| removes the last blank line but not the first one.
    lines = template.splitlines()
    if lines[0] == "":
        del lines[0]
    unindented = "\n".join(lines) % value_map

    # Indent all non-blank lines with |indent| spaces, except for the first
    # line. The first line isn't indented because the embedding template
    # provides its indention.
    unindented_lines = unindented.splitlines()
    indented_lines = unindented_lines[0:1] + \
                     [(' ' * indent) + l if l != "" else l
                     for l in unindented_lines[1:]]

    return '\n'.join(indented_lines)

def to_montgomery(x, p):
    return (x * 2**p.bit_length()) % p

# Given a number |x|, return a generator of a sequence |a| such that
#
#     x == a[0] + a[1]*2**limb_bits + ...
#
def little_endian_limbs(x, limb_bits):
    if x < 0:
        raise ValueError("x must be positive");
    if x == 0:
        yield [0]
        return

    while x != 0:
        x, digit = divmod(x, 2**limb_bits)
        yield digit

def format_hex_or_decimal(x, limb_bits):
    """format_hex_or_decimal formats a positive integer in decimal if its
       magnitude, when interpreted as a two's complement signed integer,
       is small, and formats it in fixed-width hexadecimal otherwise. This
       inconsistent formatting makes the special structure of the large
       values used in high-performance curves easier to see."""
    if x < 0:
        raise ValueError("x must be positive")
    if x <= 0x100:
        return "%d" % x
    if x >= 2**limb_bits - 0x100:
        return "(Limb)(%d)" % (x - (2**limb_bits)) # two's complement to negative
    return ("0x%0" + ("%d" % (limb_bits / 4)) + "x") % x

def group_list(x, group_len):
    """pair_list returns a list with the items in x grouped into sub-arrays,
       |group_len| at a time. For example, |group_list([1,2,3,4,5,6], 2)|
       returns |[[1, 2], [3, 4], [5, 6]]|."""
    return [x[i:i+group_len] for i in range(0, len(x), group_len)]

def format_big_int(x, limbs, limb_bits, indent):
    parts = list(little_endian_limbs(x, limb_bits))
    parts += (limbs - len(parts)) * [0]
    grouped_parts = group_list(parts, 128 // limb_bits)
    formatted_parts = [", ".join([format_hex_or_decimal(x, limb_bits)
                                  for x in group])
                       for group in grouped_parts]
    first = formatted_parts[0]
    rest = [' ' * (indent + 2) + x for x in formatted_parts[1:]]
    return ",\n".join([first] + rest)

# http://rosettacode.org/wiki/Modular_inverse#Python
def modinv(a, m):
    def extended_gcd(aa, bb):
        last_rem, rem = abs(aa), abs(bb)
        x, last_x, y, last_y = 0, 1, 1, 0
        while rem:
            last_rem, (quotient, rem) = rem, divmod(last_rem, rem)
            x, last_x = last_x - quotient * x, x
            y, last_y = last_y - quotient * y, y
        return (last_rem,
                last_x * (-1 if aa < 0 else 1),
                last_y * (-1 if bb < 0 else 1))

    g, x, y = extended_gcd(a, m)
    if g != 1:
        raise ValueError
    return x % m

def format_curve_name(g):
    return "secp%dr1" % g["q"].bit_length()

def format_modp(g, name, limb_count, limb_bits, indent):
    p = g[name]
    values = {
        "bits" : p.bit_length(),
        "limbs" : format_big_int(p, limb_count, limb_bits, 2),
        "name" : name,
        "curve_name" : format_curve_name(g),
        "k": format_hex_or_decimal(modinv(-p, 2**limb_bits), limb_bits),
        "rr": format_big_int(to_montgomery(2**p.bit_length(), p), limb_count,
                             limb_bits, 2)
    }
    return fill_in_template(modp_template, values, indent)

def format_prime_curve(g, arch):
    if g["a"] != -3:
        raise ValueError("Only curves where a == -3 are supported.")
    if g["cofactor"] != 1:
        raise ValueError("Only curves with cofactor 1 are supported.")
    if g["q"] != g["spec_q"]:
        raise ValueError("Polynomial representation of q doesn't match the "
                         "literal version given in the specification.")
    if g["n"] != g["spec_n"]:
        raise ValueError("Polynomial representation of n doesn't match the "
                         "literal version given in the specification.")
    limb_count = (g["q"].bit_length() + arch["bits"] - 1) // arch["bits"]
    limb_bits = arch["bits"]
    values = {
        "bits": g["q"].bit_length(),
        "name": format_curve_name(g),
        "q" : format_modp(g, "q", limb_count, limb_bits, 2),
        "n" : format_modp(g, "n", limb_count, limb_bits, 2),
        "one" : format_big_int(to_montgomery(1, g["q"]), limb_count, limb_bits,
                               4),
        "three" : format_big_int(to_montgomery(3, g["q"]), limb_count,
                                 limb_bits, 4),
        "b" : format_big_int(to_montgomery(g["b"], g["q"]), limb_count,
                             limb_bits, 4),
        "Gx" : format_big_int(to_montgomery(g["Gx"], g["q"]), limb_count,
                              limb_bits, 6),
        "Gy" : format_big_int(to_montgomery(g["Gy"], g["q"]), limb_count,
                              limb_bits, 6),
        "q_minus_n" : format_big_int(g["q"] - g["n"], limb_count, limb_bits,
                                     2),
    }

    declarations_template = nistz_declarations_template
    definitions_template = nistz_definitions_template

    values["declarations"] = fill_in_template(declarations_template, values, 0)
    values["definitions"] = fill_in_template(definitions_template, values, 0)

    return fill_in_template(curve_template, values, 0) + "\n"


# The curve parameters are from
# http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf. |q| and |n| are
# given in the polynomial representation so that their special structures are
# more evident. |spec_q| and |spec_n| are the literal representations of the
# same values given in the specs, and are provided only so that
# |generate_prime_curve_code| can verify that |q| and |n| are correct.

secp256r1 = {
    "q": 2**256 - 2**224 + 2**192 + 2**96 - 1,
    "spec_q" : 115792089210356248762697446949407573530086143415290314195533631308867097853951,
    "n": 2**256 - 2**224 + 2**192 - 2**128 + 0xbce6faada7179e84f3b9cac2fc632551,
    "spec_n" : 115792089210356248762697446949407573529996955224135760342422259061068512044369,
    "a": -3,
    "b": 0x5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b,
    "Gx": 0x6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296,
    "Gy": 0x4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5,
    "cofactor": 1,
}

x86 = {
    "name" : "x86",
    "bits" : 32,
}
x64 = {
    "name" : "x64",
    "bits" : 64,
}
armv4 = {
    "name" : "armv4",
    "bits" : 32,
}
armv8 = {
    "name" : "armv8",
    "bits" : 64,
}
architectures = {
    x86["name"] : x86,
    x64["name"] : x64,
    armv4["name"] : armv4,
    armv8["name"] : armv8,
}

import argparse
import os

def generate_prime_curve_file(g, arch, out_dir):
    code = format_prime_curve(g, arch)
    out_path = os.path.join(out_dir, "gfp_%s.c" % format_curve_name(g))
    with open(out_path, "wb") as f:
        f.write(code)

arg_parser = argparse.ArgumentParser(description='Generate GFp code.')
arg_parser.add_argument("--outdir", help="The output directory", required=True)
arg_parser.add_argument("--arch", help="The architecture to generate code for",
                        choices = architectures.keys(), required=True)
args = arg_parser.parse_args()

arch = architectures[args.arch]
generate_prime_curve_file(secp256r1, arch, args.outdir)