void curve25519_zero(double h[12])
{
  h[0] = 0.0;
  h[1] = 0.0;
  h[2] = 0.0;
  h[3] = 0.0;
  h[4] = 0.0;
  h[5] = 0.0;
  h[6] = 0.0;
  h[7] = 0.0;
  h[8] = 0.0;
  h[9] = 0.0;
  h[10] = 0.0;
  h[11] = 0.0;
}

void curve25519_one(double h[12])
{
  h[0] = 1.0;
  h[1] = 0.0;
  h[2] = 0.0;
  h[3] = 0.0;
  h[4] = 0.0;
  h[5] = 0.0;
  h[6] = 0.0;
  h[7] = 0.0;
  h[8] = 0.0;
  h[9] = 0.0;
  h[10] = 0.0;
  h[11] = 0.0;
}

void curve25519_copy(double h[12],const double f[12])
{
  h[0] = f[0];
  h[1] = f[1];
  h[2] = f[2];
  h[3] = f[3];
  h[4] = f[4];
  h[5] = f[5];
  h[6] = f[6];
  h[7] = f[7];
  h[8] = f[8];
  h[9] = f[9];
  h[10] = f[10];
  h[11] = f[11];
}

void curve25519_mulcarry(double h[12],const double f[12],const double g[12])
{
  double t[23];
  curve25519_mul(t,f,g);
  curve25519_carry(h,t);
}

void curve25519_sqcarry(double h[12],const double f[12])
{
  double t[23];
  curve25519_sq(t,f);
  curve25519_carry(h,t);
}

#define mult curve25519_mulcarry
#define square curve25519_sqcarry

void curve25519_recip(double out[12],const double z[12])
{
  double z2[12];
  double z9[12];
  double z11[12];
  double z2_5_0[12];
  double z2_10_0[12];
  double z2_20_0[12];
  double z2_50_0[12];
  double z2_100_0[12];
  double t0[12];
  double t1[12];
  int i;

  /* 2 */ square(z2,z);
  /* 4 */ square(t1,z2);
  /* 8 */ square(t0,t1);
  /* 9 */ mult(z9,t0,z);
  /* 11 */ mult(z11,z9,z2);
  /* 22 */ square(t0,z11);
  /* 2^5 - 2^0 = 31 */ mult(z2_5_0,t0,z9);

  /* 2^6 - 2^1 */ square(t0,z2_5_0);
  /* 2^7 - 2^2 */ square(t1,t0);
  /* 2^8 - 2^3 */ square(t0,t1);
  /* 2^9 - 2^4 */ square(t1,t0);
  /* 2^10 - 2^5 */ square(t0,t1);
  /* 2^10 - 2^0 */ mult(z2_10_0,t0,z2_5_0);

  /* 2^11 - 2^1 */ square(t0,z2_10_0);
  /* 2^12 - 2^2 */ square(t1,t0);
  /* 2^20 - 2^10 */ for (i = 2;i < 10;i += 2) { square(t0,t1); square(t1,t0); }
  /* 2^20 - 2^0 */ mult(z2_20_0,t1,z2_10_0);

  /* 2^21 - 2^1 */ square(t0,z2_20_0);
  /* 2^22 - 2^2 */ square(t1,t0);
  /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { square(t0,t1); square(t1,t0); }
  /* 2^40 - 2^0 */ mult(t0,t1,z2_20_0);

  /* 2^41 - 2^1 */ square(t1,t0);
  /* 2^42 - 2^2 */ square(t0,t1);
  /* 2^50 - 2^10 */ for (i = 2;i < 10;i += 2) { square(t1,t0); square(t0,t1); }
  /* 2^50 - 2^0 */ mult(z2_50_0,t0,z2_10_0);

  /* 2^51 - 2^1 */ square(t0,z2_50_0);
  /* 2^52 - 2^2 */ square(t1,t0);
  /* 2^100 - 2^50 */ for (i = 2;i < 50;i += 2) { square(t0,t1); square(t1,t0); }
  /* 2^100 - 2^0 */ mult(z2_100_0,t1,z2_50_0);

  /* 2^101 - 2^1 */ square(t1,z2_100_0);
  /* 2^102 - 2^2 */ square(t0,t1);
  /* 2^200 - 2^100 */ for (i = 2;i < 100;i += 2) { square(t1,t0); square(t0,t1); }
  /* 2^200 - 2^0 */ mult(t1,t0,z2_100_0);

  /* 2^201 - 2^1 */ square(t0,t1);
  /* 2^202 - 2^2 */ square(t1,t0);
  /* 2^250 - 2^50 */ for (i = 2;i < 50;i += 2) { square(t0,t1); square(t1,t0); }
  /* 2^250 - 2^0 */ mult(t0,t1,z2_50_0);

  /* 2^251 - 2^1 */ square(t1,t0);
  /* 2^252 - 2^2 */ square(t0,t1);
  /* 2^253 - 2^3 */ square(t1,t0);
  /* 2^254 - 2^4 */ square(t0,t1);
  /* 2^255 - 2^5 */ square(t1,t0);
  /* 2^255 - 21 */ mult(out,t1,z11);
}

void curve25519(unsigned char ek[32],
  const unsigned char e[32],
  const unsigned char k[32])
{
  int pos;
  int b;
  double x1[12];
  double xzm[24]; /* x[m], z[m] */
  double xzm1[24]; /* x[m+1], z[m+1] */
  double xzmb[24]; /* x[m+b], z[m+b] */
  double xzm1b[24]; /* x[m+1-b], z[m+1-b] */
  double a0[24]; /* x[m+b]+z[m+b], x[m+b]-z[m+b] */
  double b0[24]; /* (x[m+b]+z[m+b])^2, (x[m+b]-z[m+b])^2 */
  double s[12]; /* (x[m+b]+z[m+b])^2 - (x[m+b]-z[m+b])^2 */
  double t[12]; /* 121665 s */
  double u[12]; /* t + (x[m+b]+z[m+b])^2 */
  double a1[24]; /* x[m+1-b]+z[m+1-b], x[m+1-b]-z[m+1-b] */
  double b1[24]; /* (x[m+b]-z[m+b])(x[m+1-b]+z[m+1-b]), (x[m+b]+z[m+b])(x[m+1-b]-z[m+1-b]) */
  double c1[24]; /* +, - */
  double r[12]; /* (-)^2 */
  double xznb[24]; /* x[2m+2b], z[2m+2b] */
  double xzn1b[24]; /* x[2m+1], z[2m+1] */
  double recip[12];
  double result[12];

  curve25519_todouble(x1,k);
  curve25519_one(xzm);
  curve25519_zero(xzm + 12);
  curve25519_copy(xzm1,x1);
  curve25519_one(xzm1 + 12);

  for (pos = 254;pos >= 0;--pos) {
    b = e[pos / 8] >> (pos & 7);
    b &= 1;
    curve25519_select(xzmb,xzm1b,xzm,xzm1,b);
    curve25519_add(a0,xzmb,xzmb + 12);
    curve25519_sub(a0 + 12,xzmb,xzmb + 12);
    curve25519_add(a1,xzm1b,xzm1b + 12);
    curve25519_sub(a1 + 12,xzm1b,xzm1b + 12);
    curve25519_sqcarry(b0,a0);
    curve25519_sqcarry(b0 + 12,a0 + 12);
    curve25519_mulcarry(b1,a1,a0 + 12);
    curve25519_mulcarry(b1 + 12,a1 + 12,a0);
    curve25519_add(c1,b1,b1 + 12);
    curve25519_sub(c1 + 12,b1,b1 + 12);
    curve25519_sqcarry(r,c1 + 12);
    curve25519_sub(s,b0,b0 + 12);
    curve25519_121665(t,s);
    curve25519_add(u,t,b0);
    curve25519_mulcarry(xznb,b0,b0 + 12);
    curve25519_mulcarry(xznb + 12,s,u);
    curve25519_sqcarry(xzn1b,c1);
    curve25519_mulcarry(xzn1b + 12,r,x1);
    curve25519_select(xzm,xzm1,xznb,xzn1b,b);
  }

  curve25519_recip(recip,xzm + 12);
  curve25519_mulcarry(result,recip,xzm);
  curve25519_freeze(ek,result);
}

const double curve25519_constants[34] = {
  170141183460469231731687303715884105728.0 /* 2^127 */
, 0.0000000000000000000000000000000000000058774717541114375398436826861112283890933277838604376075437585313920862972736358642578125 /* 2^-127 */
, 0.000000000000000000000000000000000000000000000000000000000000000000000000000328174405093588895764681370786415183702408013848578692630900731866406488520172228202917285131946952609300797169953687216685813759979613974358814143528206841438077390193939208984375 /* 19 2^-255 */
, 2.0
, 28334198897217871282176.0 /* 2^22 * (2^52 + 2^51) */
, 59421121885698253195157962752.0 /* 2^43 * (2^52 + 2^51) */
, 124615124604835863084731911901282304.0 /* 2^64 * (2^52 + 2^51) */
, 261336857795280739939871698507597986398208.0 /* 2^85 * (2^52 + 2^51) */
, 1096126227998177188652763624537212264741949407232.0 /* 2^107 * (2^52 + 2^51) */
, 2298743311298833287537520540725463775428108683275403264.0 /* 2^128 * (2^52 + 2^51) */
, 4820814132776970826625886277023487807566608981348378505904128.0 /* 2^149 * (2^52 + 2^51) */
, 10109980000181489923000130657632361502613929158452714680413853843456.0 /* 2^170 * (2^52 + 2^51) */
, 42404329554681223909999140017830044379859613525014854994918548831022874624.0 /* 2^192 * (2^52 + 2^51) */
, 88928324534258838085302516486672313231311348223211953182303424518077283563470848.0 /* 2^213 * (2^52 + 2^51) */
, 186496213653669990808268343055057815037671056549005394040173991334934811379700015824896.0 /* 2^234 * (2^52 + 2^51) */
, 391110907456221328563541572174600606921881931583859760122138966276041209554560647587212296192.0 /* 2^255 * (2^52 + 2^51) */
, 6755399441055744.0 /* 2^0 * (2^52 + 2^51) */
, 1901475900342344102245054808064.0 /* 2^48 * (2^52 + 2^51) */
, 2090694862362245919518973588060783891185664.0 /* 2^88 * (2^52 + 2^51) */
, 647038720011615355072008362088471136167291466140973739546486645981184.0 /* 2^176 * (2^52 + 2^51) */
, 711426596274070704682420131893378505850490785785695625458427396144618268507766784.0 /* 2^216 * (2^52 + 2^51) */
, 121665.0
, 6755399449444333.0 /* 2^0 * (2^52 + 2^51) + 2^23 - 19 */
, 442721875361206894592.0 /* 2^16 * (2^52 + 2^51) + 2^44 - 2^23 */
, 7427640272605752204627804160.0 /* 2^40 * (2^52 + 2^51) + 2^65 - 2^44 */
, 124615124682207078646580031663374336.0 /* 2^64 * (2^52 + 2^51) + 2^86 - 2^65 */
, 8166777130620999410295262025251276455936.0 /* 2^80 * (2^52 + 2^51) + 2^108 - 2^86 */
, 137015779180336557904918721567173964800259522560.0 /* 2^104 * (2^52 + 2^51) + 2^129 - 2^108 */
, 2298743312726080299678746579906822818128388956121726976.0 /* 2^128 * (2^52 + 2^51) + 2^150 - 2^129 */
, 150650444642434264338055416678249473103434705427350293577728.0 /* 2^144 * (2^52 + 2^51) + 2^171 - 2^150 */
, 2527495012599572958368040938390492740922277718870265303681684471808.0 /* 2^168 * (2^52 + 2^51) + 2^193 - 2^171 */
, 42404329581009284272934965918948023629198375189405796640344937716531593216.0 /* 2^192 * (2^52 + 2^51) + 2^214 - 2^193 */
, 2779010196909533136417296800389882207187762333976165590536224320212295324008448.0 /* 2^208 * (2^52 + 2^51) + 2^235 - 2^214 */
, 46624053471313487106754409151039646785714819847682721233992002490370053548687778381824.0 /* 2^232 * (2^52 + 2^51) + 2^255 - 2^235 */
} ;