#include <stdint.h>
#include <string.h>
#include <math.h>
#include "fips202.h"
#include "params.h"
#include "sign.h"
#include "poly.h"
#include "polz.h"
#include "gaussian.h"
#include "pdata.h"
#include "fips202.h"

typedef struct {
  poly pol0;
  polx pol1;
  polx pol2;
} polvec;

static void commit(poly *t1,
                   polz *t2,
                   polvec *r,
                   const int32_t m[3],
                   const poly *a1,
                   const poly *a2,
                   const polx *a3)
{
  int i;
  uint32_t mx[6];
  polx g;
  poly *f = &g.vec[0];

  poly_binary(&r->pol0);
  poly_binary(&r->pol1.vec[0]);
  poly_binary(&r->pol2.vec[0]);
  polx_frompoly(&r->pol1,&r->pol1.vec[0]);
  polx_frompoly(&r->pol2,&r->pol2.vec[0]);

  polx_ntt(&r->pol2);
  polx_pointwise(&g,a3,&r->pol2);
  polx_invntt(&g);
  polx_add(&g,&g,&r->pol1);

  mx[0] = ((((((int64_t)m[2] << 31) + m[1]) % P0) << 31) + m[0]) % P0;
  mx[1] = ((((((int64_t)m[2] << 31) + m[1]) % P1) << 31) + m[0]) % P1;
  mx[2] = ((((((int64_t)m[2] << 31) + m[1]) % P2) << 31) + m[0]) % P2;
  mx[3] = ((((((int64_t)m[2] << 31) + m[1]) % P3) << 31) + m[0]) % P3;
  mx[4] = ((((((int64_t)m[2] << 31) + m[1]) % P4) << 31) + m[0]) % P4;
  mx[5] = ((((((int64_t)m[2] << 31) + m[1]) % P5) << 31) + m[0]) % P5;

  for(i=0;i<6;i++)
    g.vec[i].coeffs[0] += mx[i];

  polx_reduce(&g);
  polz_frompolx(t2,&g);

  poly_ntt(&r->pol0,PDATA0);
  polx_ntt(&r->pol1);

  poly_pointwise(t1,a1,&r->pol1.vec[0],PDATA0);
  poly_pointwise(f,a2,&r->pol2.vec[0],PDATA0);
  poly_add(t1,t1,f);
  poly_scale(t1,t1,PDATA0,(((1ULL << 32) % P0) << 32) % P0);
  poly_add(t1,t1,&r->pol0);
}

static void challenge(poly *c,
                      unsigned char mu[32],
                      const poly *w,
                      const poly *wprime,
                      const poly *w1,
                      const poly *w5,
                      const polz *ws,
                      const polz *x,
                      const polz *x1,
                      const polz *x5)
{
  unsigned char inbuf[32+4*4*N+4*10*N] = {0};
  unsigned char outbuf[SHAKE128_RATE];
  int i,j,pos;
  keccak_state state;
  uint32_t signs;

  memcpy(inbuf,mu,32);
  memcpy(inbuf+32,(unsigned char *)w,4*N);
  memcpy(inbuf+32+4*N,(unsigned char *)wprime,4*N);
  memcpy(inbuf+32+8*N,(unsigned char *)w1,4*N);
  memcpy(inbuf+32+12*N,(unsigned char *)w5,4*N);
  polz_pack(inbuf+32+16*N,ws);
  polz_pack(inbuf+32+26*N,x);
  polz_pack(inbuf+32+36*N,x1);
  polz_pack(inbuf+32+46*N,x5);

  shake128_absorb(&state,inbuf,sizeof(inbuf));
  shake128_squeezeblocks(outbuf,1,&state);

  signs = 0;
  for(i=0;i<4;i++)
    signs |= (uint32_t)outbuf[i] << 8*i;

  pos = 4;

  for(i=0;i<N;i++)
    c->coeffs[i] = 0;

  for(i=4070;i<4096;i++) {
    do {
      j  = outbuf[pos++];
      j |= ((int)outbuf[pos++] & 0xF) << 8;
    } while(j > i);

    c->coeffs[i] = c->coeffs[j];
    c->coeffs[j] = (signs & 1) ? -1 : 1;
    signs >>= 1;
  }
}

int sign(signature *sig,
         const poly *a1,
         const poly *a2,
         const polx *a3,
         const polx a[2],
         const polx b[2],
         const polx s[4],
         const int32_t id[3])
{
  int tries = 0;
  double d;
  unsigned char mu[32] = {0};
  polvec r,rprime;
  polx v[2],sprime[2];
  polx y,yprime,y1,y5;
  polx ys[6],g;
  poly w,wprime,w1,w5;
  poly *f0 = &g.vec[0];
  poly *f1 = &g.vec[1];
  polz x,x1,x5;
  polz ws;
  int32_t m[3];
  uint64_t tmp;

  /* m */
  m[1] = id[0] << 9;
  m[2] = (id[0] >> 22) | (id[1] << 9);
  tmp = ((uint64_t)id[1] >> 9) | ((uint64_t)id[2] << 22);
  tmp *= 65;
  m[0] = tmp;
  m[1] += tmp >> 31;
  m[2] += m[1] >> 31;
  m[0] &= (1U << 31) - 1;
  m[1] &= (1U << 31) - 1;
  m[2] &= (1U << 18) - 1;

  /* r, t, r', t' */
  commit(&sig->t1,&sig->t2,&r,id,a1,a2,a3);
  commit(&sig->tprime1,&sig->tprime2,&rprime,m,a1,a2,a3);
  polz_topolx(&v[0],&sig->t2);
  polx_ntt(&v[0]);
  polx_reduce(&v[0]);
  polz_topolx(&v[1],&sig->tprime2);
  polx_ntt(&v[1]);
  polx_reduce(&v[1]);

  /* mu */

  /* v, s' */
  polx_add(&v[0],&b[0],&v[0]);
  polx_add(&v[1],&b[1],&v[1]);

  polx_pointwise(&sprime[0],&r.pol1,&s[2]);
  polx_pointwise_add(&sprime[0],&rprime.pol1,&s[3],&sprime[0]);
  m[0] = 0; m[1] = 0; m[2] = -4;
  polx_scale(&sprime[0],&sprime[0],m);
  polx_pointwise(&sprime[1],&r.pol2,&s[2]);
  polx_pointwise_add(&sprime[1],&rprime.pol2,&s[3],&sprime[1]);
  polx_scale(&sprime[1],&sprime[1],m);

start:
  tries++;

  /* y, w */
  poly_gaussian16(&sig->z[0]);
  poly_gaussian16(&sig->z[1]);
  poly_gaussian16(&sig->z[2]);
  polx_frompoly(&y,&sig->z[2]);
  *f0 = sig->z[1];
  poly_ntt(f0,PDATA0);
  polx_ntt(&y);
  poly_pointwise(&w,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&y.vec[0],PDATA0);
  poly_add(&w,&w,f0);
  poly_invntt(&w,PDATA0);
  poly_add(&w,&w,&sig->z[0]);
  poly_reduce(&w,PDATA0);
  poly_caddp(&w,PDATA0);

  /* yprime, wprime */
  poly_gaussian16(&sig->zprime[0]);
  poly_gaussian16(&sig->zprime[1]);
  poly_gaussian16(&sig->zprime[2]);
  polx_frompoly(&yprime,&sig->zprime[2]);
  *f0 = sig->zprime[1];
  poly_ntt(f0,PDATA0);
  polx_ntt(&yprime);
  poly_pointwise(&wprime,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&yprime.vec[0],PDATA0);
  poly_add(&wprime,&wprime,f0);
  poly_invntt(&wprime,PDATA0);
  poly_add(&wprime,&wprime,&sig->zprime[0]);
  poly_reduce(&wprime,PDATA0);
  poly_caddp(&wprime,PDATA0);

  /* y1, w1 */
  poly_gaussian16(&sig->z1[0]);
  poly_gaussian16(&sig->z1[1]);
  poly_gaussian16(&sig->z1[2]);
  polx_frompoly(&y1,&sig->z1[2]);
  *f0 = sig->z1[1];
  poly_ntt(f0,PDATA0);
  polx_ntt(&y1);
  poly_pointwise(&w1,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&y1.vec[0],PDATA0);
  poly_add(&w1,&w1,f0);
  poly_invntt(&w1,PDATA0);
  poly_add(&w1,&w1,&sig->z1[0]);
  //poly_sigmam1(&w1,&w1);
  poly_reduce(&w1,PDATA0);
  poly_caddp(&w1,PDATA0);

  /* y5, w5 */
  poly_gaussian16(&sig->z5[0]);
  poly_gaussian16(&sig->z5[1]);
  poly_gaussian16(&sig->z5[2]);
  polx_frompoly(&y5,&sig->z5[2]);
  *f0 = sig->z5[1];
  poly_ntt(f0,PDATA0);
  polx_ntt(&y5);
  poly_pointwise(&w5,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&y5.vec[0],PDATA0);
  poly_add(&w5,&w5,f0);
  poly_invntt(&w5,PDATA0);
  poly_add(&w5,&w5,&sig->z5[0]);
  //poly_sigma5(&w5,&w5);
  poly_reduce(&w5,PDATA0);
  poly_caddp(&w5,PDATA0);

  /* x, x1, x5 */
  polx_pointwise(&ys[0],a3,&y);
  polx_pointwise(&ys[1],a3,&yprime);
  polx_pointwise(&ys[2],a3,&y1);
  polx_pointwise(&ys[3],a3,&y5);
  polx_invntt(&ys[0]);
  polx_invntt(&ys[1]);
  polx_invntt(&ys[2]);
  polx_invntt(&ys[3]);
  polx_add_poly(&ys[0],&ys[0],&sig->z[1]);
  polx_add_poly(&ys[1],&ys[1],&sig->zprime[1]);
  polx_add_poly(&ys[2],&ys[2],&sig->z1[1]);
  polx_add_poly(&ys[3],&ys[3],&sig->z5[1]);
  polx_sigmam1(&ys[2],&ys[2]);
  polx_sigma5(&ys[3],&ys[3]);
  polx_sub(&ys[2],&ys[0],&ys[2]);
  polx_sub(&ys[3],&ys[0],&ys[3]);
  m[0] = 0;
  m[1] = 0;
  m[2] = 1U << 10;
  polx_scale(&ys[0],&ys[0],m);
  polx_sub(&ys[1],&ys[0],&ys[1]);
  polz_frompolx(&x,&ys[1]);
  polz_frompolx(&x1,&ys[2]);
  polz_frompolx(&x5,&ys[3]);

  /* ys */
  polz_gaussian64(&sig->zs[0]);
  polz_cmodq(&sig->zs[0]); // FIXME
  polz_topolx(&ys[0],&sig->zs[0]);
  polx_ntt(&ys[0]);
  polz_gaussian64(&sig->zs[1]);
  polz_cmodq(&sig->zs[1]); // FIXME
  polz_topolx(&ys[1],&sig->zs[1]);
  polx_ntt(&ys[1]);
  polz_gaussian64(&sig->zs[2]);
  polz_cmodq(&sig->zs[2]); // FIXME
  polz_topolx(&ys[2],&sig->zs[2]);
  polx_ntt(&ys[2]);
  polz_gaussian64(&sig->zs[3]);
  polz_cmodq(&sig->zs[3]); // FIXME
  polz_topolx(&ys[3],&sig->zs[3]);
  polx_ntt(&ys[3]);
  polz_gaussian70(&sig->zs[4]);
  polz_cmodq(&sig->zs[4]); // FIXME
  polz_topolx(&ys[4],&sig->zs[4]);
  polz_gaussian70(&sig->zs[5]);
  polz_cmodq(&sig->zs[5]); // FIXME
  polz_topolx(&ys[5],&sig->zs[5]);
  polx_ntt(&ys[5]);

  /* ws */
  polx_pointwise(&g,&a[0],&ys[0]);
  polx_pointwise_add(&g,&a[1],&ys[1],&g);
  polx_pointwise_add(&g,&v[0],&ys[2],&g);
  polx_reduce(&g);
  polx_pointwise_add(&g,&v[1],&ys[3],&g);
  polx_pointwise_add(&g,a3,&ys[5],&g);
  polx_reduce(&g);
  polx_invntt(&g);
  polx_add(&g,&g,&ys[4]);
  polz_frompolx(&ws,&g);

  /* c */
  challenge(&sig->c,mu,&w,&wprime,&w1,&w5,&ws,&x,&x1,&x5);

  /* z, z', z1, z5 */
  *f0 = sig->c;
  poly_ntt(f0,PDATA0);
  poly_pointwise(f1,f0,&r.pol0,PDATA0);
  poly_invntt(f1,PDATA0);
  d  = -poly_sprod(&sig->z[0],&sig->z[0]);
  poly_add(&sig->z[0],&sig->z[0],f1);
  poly_reduce(&sig->z[0],PDATA0);
  d += poly_sprod(&sig->z[0],&sig->z[0]);
  poly_pointwise(f1,f0,&r.pol1.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z[1],&sig->z[1]);
  poly_add(&sig->z[1],&sig->z[1],f1);
  poly_reduce(&sig->z[1],PDATA0);
  d += poly_sprod(&sig->z[1],&sig->z[1]);
  poly_pointwise(f1,f0,&r.pol2.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z[2],&sig->z[2]);
  poly_add(&sig->z[2],&sig->z[2],f1);
  poly_reduce(&sig->z[2],PDATA0);
  d += poly_sprod(&sig->z[2],&sig->z[2]);

  poly_pointwise(f1,f0,&rprime.pol0,PDATA0);
  poly_invntt(f1,PDATA0);
  d  = -poly_sprod(&sig->zprime[0],&sig->zprime[0]);
  poly_add(&sig->zprime[0],&sig->zprime[0],f1);
  poly_reduce(&sig->zprime[0],PDATA0);
  d += poly_sprod(&sig->zprime[0],&sig->zprime[0]);
  poly_pointwise(f1,f0,&rprime.pol1.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->zprime[1],&sig->zprime[1]);
  poly_add(&sig->zprime[1],&sig->zprime[1],f1);
  poly_reduce(&sig->zprime[1],PDATA0);
  d += poly_sprod(&sig->zprime[1],&sig->zprime[1]);
  poly_pointwise(f1,f0,&rprime.pol2.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->zprime[2],&sig->zprime[2]);
  poly_add(&sig->zprime[2],&sig->zprime[2],f1);
  poly_reduce(&sig->zprime[2],PDATA0);
  d += poly_sprod(&sig->zprime[2],&sig->zprime[2]);

  poly_sigmam1(f0,&sig->c);
  poly_ntt(f0,PDATA0);
  poly_pointwise(f1,f0,&r.pol0,PDATA0);
  poly_invntt(f1,PDATA0);
  d  = -poly_sprod(&sig->z1[0],&sig->z1[0]);
  poly_add(&sig->z1[0],&sig->z1[0],f1);
  poly_reduce(&sig->z1[0],PDATA0);
  d += poly_sprod(&sig->z1[0],&sig->z1[0]);
  poly_pointwise(f1,f0,&r.pol1.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z1[1],&sig->z1[1]);
  poly_add(&sig->z1[1],&sig->z1[1],f1);
  poly_reduce(&sig->z1[1],PDATA0);
  d += poly_sprod(&sig->z1[1],&sig->z1[1]);
  poly_pointwise(f1,f0,&r.pol2.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z1[2],&sig->z1[2]);
  poly_add(&sig->z1[2],&sig->z1[2],f1);
  poly_reduce(&sig->z1[2],PDATA0);
  d += poly_sprod(&sig->z1[2],&sig->z1[2]);

  poly_sigma5inv(f0,&sig->c);
  poly_ntt(f0,PDATA0);
  poly_pointwise(f1,f0,&r.pol0,PDATA0);
  poly_invntt(f1,PDATA0);
  d  = -poly_sprod(&sig->z5[0],&sig->z5[0]);
  poly_add(&sig->z5[0],&sig->z5[0],f1);
  poly_reduce(&sig->z5[0],PDATA0);
  d += poly_sprod(&sig->z5[0],&sig->z5[0]);
  poly_pointwise(f1,f0,&r.pol1.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z5[1],&sig->z5[1]);
  poly_add(&sig->z5[1],&sig->z5[1],f1);
  poly_reduce(&sig->z5[1],PDATA0);
  d += poly_sprod(&sig->z5[1],&sig->z5[1]);
  poly_pointwise(f1,f0,&r.pol2.vec[0],PDATA0);
  poly_invntt(f1,PDATA0);
  d -= poly_sprod(&sig->z5[2],&sig->z5[2]);
  poly_add(&sig->z5[2],&sig->z5[2],f1);
  poly_reduce(&sig->z5[2],PDATA0);
  d += poly_sprod(&sig->z5[2],&sig->z5[2]);

  d *= 1/ldexp(2*1.55*1.55,32);
  d += log(3);
  if(!BerExp(d))
    goto start;

  /* zs */
  *f0 = sig->c;
  polx_frompoly(&g,f0);
  polx_ntt(&g);
  polx_pointwise(&ys[0],&g,&s[0]);
  polx_pointwise(&ys[1],&g,&s[1]);
  polx_pointwise(&ys[2],&g,&s[2]);
  polx_pointwise(&ys[3],&g,&s[3]);
  polx_pointwise(&ys[4],&g,&sprime[0]);
  polx_pointwise(&ys[5],&g,&sprime[1]);
  polx_invntt(&ys[0]);
  polx_invntt(&ys[1]);
  polx_invntt(&ys[2]);
  polx_invntt(&ys[3]);
  polx_invntt(&ys[4]);
  polx_invntt(&ys[5]);
  polz_frompolx(&x,&ys[0]);
  polz_cmodq(&x);
  double d0 = 1/ldexp(2*1.55*1.55,128);
  d  = polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[0])*d0;
  polz_add(&sig->zs[0],&sig->zs[0],&x);
  polz_frompolx(&x,&ys[1]);
  polz_cmodq(&x);
  d += polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[1])*d0;
  polz_add(&sig->zs[1],&sig->zs[1],&x);
  polz_frompolx(&x,&ys[2]);
  polz_cmodq(&x);
  d += polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[2])*d0;
  polz_add(&sig->zs[2],&sig->zs[2],&x);
  polz_frompolx(&x,&ys[3]);
  polz_cmodq(&x);
  d += polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[3])*d0;
  polz_add(&sig->zs[3],&sig->zs[3],&x);

  //d *= 1/ldexp(2*1.55*1.55,128);
  d += log(3);
  if(!BerExp(d))
    goto start;

  polz_frompolx(&x,&ys[4]);
  polz_cmodq(&x);
  d0 = 1/ldexp(2*1.55*1.55,140);
  d  = polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[4])*d0;
  polz_add(&sig->zs[4],&sig->zs[4],&x);
  polz_frompolx(&x,&ys[5]);
  polz_cmodq(&x);
  d  = polz_sprod(&x,&x)*d0;
  d += 2*polz_sprod(&x,&sig->zs[5])*d0;
  polz_add(&sig->zs[5],&sig->zs[5],&x);

  //d *= 1/ldexp(2*1.55*1.55,140);
  d += log(3);
  if(!BerExp(d))
    goto start;

  return tries;
}

int verify(const signature *sig,
           const poly *a1,
           const poly *a2,
           const polx *a3,
           const polx a[2],
           const polx b[2],
           const polz *u)
{
  int i;
  double d;
  unsigned char mu[32] = {0};
  polx t2,tprime2,v[2];
  polx z,zprime,z1,z5;
  poly w,wprime,w1,w5;
  polz ws,x,x1,x5;
  polx g,h;
  poly *f0 = &g.vec[0];
  poly *f1 = &g.vec[1];
  int32_t s[3];

  d  = poly_sprod(&sig->z[0],&sig->z[0]);
  d += poly_sprod(&sig->z[1],&sig->z[1]);
  d += poly_sprod(&sig->z[2],&sig->z[2]);
  d += poly_sprod(&sig->zprime[0],&sig->zprime[0]);
  d += poly_sprod(&sig->zprime[1],&sig->zprime[1]);
  d += poly_sprod(&sig->zprime[2],&sig->zprime[2]);
  d += poly_sprod(&sig->z1[0],&sig->z1[0]);
  d += poly_sprod(&sig->z1[1],&sig->z1[1]);
  d += poly_sprod(&sig->z1[2],&sig->z1[2]);
  d += poly_sprod(&sig->z5[0],&sig->z5[0]);
  d += poly_sprod(&sig->z5[1],&sig->z5[1]);
  d += poly_sprod(&sig->z5[2],&sig->z5[2]);
  if(sqrt(d) > 1.55*65536*sqrt(24*N))
    return 1;

  d  = polz_sprod(&sig->zs[0],&sig->zs[0]);
  d += polz_sprod(&sig->zs[1],&sig->zs[1]);
  d += polz_sprod(&sig->zs[2],&sig->zs[2]);
  d += polz_sprod(&sig->zs[3],&sig->zs[3]);
  if(sqrt(d) > ldexp(1.55,64)*sqrt(8*N))
    return 1;

  d  = polz_sprod(&sig->zs[4],&sig->zs[4]);
  d += polz_sprod(&sig->zs[5],&sig->zs[5]);
  if(sqrt(d) > ldexp(1.55,70)*sqrt(4*N))
    return 1;

  /* v */
  polz_topolx(&t2,&sig->t2);
  polx_ntt(&t2);
  polx_reduce(&t2);
  polx_add(&v[0],&b[0],&t2);
  polz_topolx(&tprime2,&sig->tprime2);
  polx_ntt(&tprime2);
  polx_reduce(&tprime2);
  polx_add(&v[1],&b[1],&tprime2);

  /* w, w', w1, w5 */
  *f0 = sig->z[1];
  polx_frompoly(&z,&sig->z[2]);
  poly_ntt(f0,PDATA0);
  polx_ntt(&z);
  poly_pointwise(&w,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&z.vec[0],PDATA0);
  poly_add(&w,&w,f0);

  *f1 = sig->c;
  poly_ntt(f1,PDATA0);
  poly_pointwise(f0,f1,&sig->t1,PDATA0);
  poly_sub(&w,&w,f0);
  poly_reduce(&w,PDATA0);
  poly_invntt(&w,PDATA0);
  poly_add(&w,&w,&sig->z[0]);
  poly_reduce(&w,PDATA0);
  poly_caddp(&w,PDATA0);

  *f0 = sig->zprime[1];
  polx_frompoly(&zprime,&sig->zprime[2]);
  poly_ntt(f0,PDATA0);
  polx_ntt(&zprime);
  poly_pointwise(&wprime,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&zprime.vec[0],PDATA0);
  poly_add(&wprime,&wprime,f0);

  poly_pointwise(f0,f1,&sig->tprime1,PDATA0);
  poly_sub(&wprime,&wprime,f0);
  poly_reduce(&wprime,PDATA0);
  poly_invntt(&wprime,PDATA0);
  poly_add(&wprime,&wprime,&sig->zprime[0]);
  poly_reduce(&wprime,PDATA0);
  poly_caddp(&wprime,PDATA0);

  *f0 = sig->z1[1];
  polx_frompoly(&z1,&sig->z1[2]);
  poly_ntt(f0,PDATA0);
  polx_ntt(&z1);
  poly_pointwise(&w1,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&z1.vec[0],PDATA0);
  poly_add(&w1,&w1,f0);

  poly_sigmam1(f0,&sig->c);
  poly_ntt(f0,PDATA0);
  poly_pointwise(f0,f0,&sig->t1,PDATA0);
  poly_sub(&w1,&w1,f0);
  poly_reduce(&w1,PDATA0);
  poly_invntt(&w1,PDATA0);
  poly_add(&w1,&w1,&sig->z1[0]);
  //poly_sigmam1(&w1,&w1);
  poly_reduce(&w1,PDATA0);
  poly_caddp(&w1,PDATA0);

  *f0 = sig->z5[1];
  polx_frompoly(&z5,&sig->z5[2]);
  poly_ntt(f0,PDATA0);
  polx_ntt(&z5);
  poly_pointwise(&w5,a1,f0,PDATA0);
  poly_pointwise(f0,a2,&z5.vec[0],PDATA0);
  poly_add(&w5,&w5,f0);

  poly_sigma5inv(f0,&sig->c);
  poly_ntt(f0,PDATA0);
  poly_pointwise(f0,f0,&sig->t1,PDATA0);
  poly_sub(&w5,&w5,f0);
  poly_reduce(&w5,PDATA0);
  poly_invntt(&w5,PDATA0);
  poly_add(&w5,&w5,&sig->z5[0]);
  //poly_sigma5(&w5,&w5);
  poly_reduce(&w5,PDATA0);
  poly_caddp(&w5,PDATA0);

  /* x, x1, x5 */
  polx_pointwise(&z,a3,&z);
  polx_pointwise(&zprime,a3,&zprime);
  polx_frompoly(&g,&sig->c);
  polx_ntt(&g);
  polx_pointwise(&h,&g,&t2);
  polx_sub(&z,&z,&h);
  polx_pointwise(&h,&g,&tprime2);
  polx_sub(&zprime,&zprime,&h);
  polx_reduce(&z);
  polx_invntt(&z);
  polx_add_poly(&z,&z,&sig->z[1]);
  polx_reduce(&zprime);
  polx_invntt(&zprime);
  polx_add_poly(&zprime,&zprime,&sig->zprime[1]);

  polx_pointwise(&z1,a3,&z1);
  poly_sigmam1(f0,&sig->c);
  polx_frompoly(&g,f0);
  polx_ntt(&g);
  polx_pointwise(&g,&t2,&g);
  polx_sub(&z1,&z1,&g);
  polx_reduce(&z1);
  polx_invntt(&z1);
  polx_add_poly(&z1,&z1,&sig->z1[1]);
  polx_sigmam1(&z1,&z1);
  polx_sub(&z1,&z,&z1);
  polz_frompolx(&x1,&z1);

  polx_pointwise(&z5,a3,&z5);
  poly_sigma5inv(f0,&sig->c);
  polx_frompoly(&g,f0);
  polx_ntt(&g);
  polx_pointwise(&g,&t2,&g);
  polx_sub(&z5,&z5,&g);
  polx_reduce(&z5);
  polx_invntt(&z5);
  polx_add_poly(&z5,&z5,&sig->z5[1]);
  polx_sigma5(&z5,&z5);
  polx_sub(&z5,&z,&z5);
  polz_frompolx(&x5,&z5);

  s[0] = 0; s[1] = 0; s[2] = 1U << 10;
  polx_scale(&z,&z,s);
  polx_sub(&z,&z,&zprime);
  polz_frompolx(&x,&z);

  /* ws */
  polz_topolx(&g,&sig->zs[0]);
  polx_ntt(&g);
  polx_pointwise(&g,&a[0],&g);
  polz_topolx(&h,&sig->zs[1]);
  polx_ntt(&h);
  polx_pointwise(&h,&a[1],&h);
  polx_add(&g,&g,&h);
  polz_topolx(&h,&sig->zs[2]);
  polx_ntt(&h);
  polx_pointwise(&h,&v[0],&h);
  polx_add(&g,&g,&h);
  polx_reduce(&g);
  polz_topolx(&h,&sig->zs[3]);
  polx_ntt(&h);
  polx_pointwise(&h,&v[1],&h);
  polx_add(&g,&g,&h);
  polz_topolx(&h,&sig->zs[5]);
  polx_ntt(&h);
  polx_pointwise(&h,a3,&h);
  polx_add(&g,&g,&h);
  polx_reduce(&g);
  polx_invntt(&g);
  polz_topolx(&h,&sig->zs[4]);
  polx_add(&g,&g,&h);

  polx_frompoly(&h,&sig->c);
  polx_ntt(&h);
  polz_topolx(&v[0],u);
  polx_ntt(&v[0]);
  polx_pointwise(&h,&h,&v[0]);
  polx_invntt(&h);
  polx_sub(&g,&g,&h);
  polx_reduce(&g);
  polz_frompolx(&ws,&g);

  challenge(f0,mu,&w,&wprime,&w1,&w5,&ws,&x,&x1,&x5);

  for(i=0;i<N;i++)
    if(f0->coeffs[i] != sig->c.coeffs[i])
      return 1;

  return 0;
}
