#include <stdio.h>

double powers[12] = {
  1.0 /* 2^0 */
, 4194304.0 /* 2^22 */
, 8796093022208.0 /* 2^43 */
, 18446744073709551616.0 /* 2^64 */
, 38685626227668133590597632.0 /* 2^85 */
, 162259276829213363391578010288128.0 /* 2^107 */
, 340282366920938463463374607431768211456.0 /* 2^128 */
, 713623846352979940529142984724747568191373312.0 /* 2^149 */
, 1496577676626844588240573268701473812127674924007424.0 /* 2^170 */
, 6277101735386680763835789423207666416102355444464034512896.0 /* 2^192 */
, 13164036458569648337239753460458804039861886925068638906788872192.0 /* 2^213 */
, 27606985387162255149739023449108101809804435888681546220650096895197184.0 /* 2^234 */
} ;

char s[32];
char t[32];
char u[32];

double f[12];
double g[12];
double h[23];

extern void curve25519_todouble(double [12],const char [32]);
extern void curve25519_mul(double [23],const double [12],const double [12]);
extern void curve25519_carry(double [12],const double [23]);
extern void curve25519_freeze(char [32],const double [12]);

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 */
} ;

int main(int argc,char **argv)
{
  int i;
  int loop;

  for (loop = 0;loop < 100;++loop) {
    for (i = 0;i < 32;++i) s[i] = random(); s[31] &= 127;
    printf("f=0"); for (i = 0;i < 32;++i) printf("%+d*256^%d",255 & (int) s[i],i); printf(";\n");
    curve25519_todouble(f,s);
    for (i = 0;i < 32;++i) s[i] = random(); s[31] &= 127;
    printf("g=0"); for (i = 0;i < 32;++i) printf("%+d*256^%d",255 & (int) s[i],i); printf(";\n");
    curve25519_todouble(g,s);
    curve25519_mul(h,f,g);
    curve25519_carry(f,h);
    curve25519_freeze(s,f);
    printf("h=0"); for (i = 0;i < 32;++i) printf("%+d*256^%d",255 & (int) s[i],i); printf(";\n");
    printf("h-((f*g) %% (2^255-19))\n");
  }

  return 0;
}