#include <Arduino.h>
#include <avr/io.h>
#include <stdlib.h>

// Define constants and function prototypes for Ed25519
#define add Assembly_Ed_add
#define sub Assembly_Ed_sub
#define modulo Assembly_Ed_mod

typedef struct
{
  unsigned char Ed[32];
}
field_element;

extern "C"
{
  void sub(field_element *r, const field_element *x, const field_element *y);
  void add(field_element *r, const field_element *x, const field_element *y);
  void modulo(field_element *r, unsigned char *C);
  char bignum_sub_prime(unsigned char* r, const unsigned char* a);
  char Assembly_Ed_square(unsigned char* r, const unsigned char* a);
  char Assembly_Ed_mul(unsigned char* r, const unsigned char* a, const unsigned char* b);
  unsigned char scalar_sub_halforder(unsigned char* r, const unsigned char* a);
  char prime_sub_bignum(unsigned char* r, const unsigned char* a);
}

void finalize_result(field_element *r);
void conditional_swap(field_element *r, const field_element *x, unsigned char b);
void mul(field_element *r, const field_element *x, const field_element *y);
void square(field_element *r, const field_element *x);
void Mul_Inverse(field_element *r, const field_element *x);
void cswap(field_element *p_coord, char b);
void ed25519_ladder(field_element *xr, field_element *yr, field_element *zr, unsigned char scalar[32]);

// d = (-121665/121666) mod prime
static const field_element d = {{0xA3, 0x78, 0x59, 0x13, 0xCA, 0x4D, 0xEB, 0x75, 0xAB, 0xD8, 0x41, 0x41, 0x4D, 0x0A, 0x70, 0x00, 0x98, 0xE8, 0x79, 0x77, 0x79, 0x40, 0xC7, 0x8C, 0x73, 0xFE, 0x6F, 0x2B, 0xEE, 0x6C, 0x03, 0x52}};

// a = (-1) mod prime
static const field_element a = {{0xEC, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x7F}};

//Order = (2^255)-19 ============> Half_Order = Order/2
static const field_element Half_Order = {{0xF6,	0xE9,	0x7A,	0x2E,	0x8D,	0x31,	0x09,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08}};

int main()
{
  Serial.begin(500000);
  Serial.println();

  //scalar = 5
  unsigned char scalar[32] = {0x05};

  //Order = (2^255)-19 = 7237005577332262213973186563042994240857116359379907606001950938285454250989 = 0x2aee5456Fdc8959dEcB68627Cb2D7BCDaE08e8e1
  //Half_Order = 3618502788666131106986593281521497120428558179689953803000975469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c09318d2e7ae9f6

  //scalar = Half_Order - 100 000 000 000 = 3618502788666131106986593281521497120428558179689953803000975469042727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c093175e60401f6
  //unsigned char scalar[32] = {0xF6,	0x01,	0x04,	0xE6,	0x75,	0x31,	0x09,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 = 3618502788666131106986593281521497120428558179689953803000975468142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c0930a459d5d9f6
  //unsigned char scalar[32] = {0xF6,	0xD9,	0xD5,	0x59,	0xA4,	0x30,	0x09,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 = 3618502788666131106986593281521497120428558179689953803000974469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c05a40e89b469f6
  //unsigned char scalar[32] = {0xF6,	0x69,	0xB4,	0x89,	0x0E,	0xA4,	0x05,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 000  = 3618502788666131106986593281521497120428558179689953802999975469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b1e287ad98716e9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x16,	0x87,	0xD9,	0x7A,	0x28,	0x1E,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179689953802000975469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce34f63f83c74fdae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0xDA,	0x4F,	0xC7,	0x83,	0x3F,	0xF6,	0x34,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179689952803000975469142727125494 = 0x80000000000000000000000000000000a6f7cef517afaa9103a649f8d7ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0x8D,	0x9F,	0x64,	0x3A,	0x10,	0xA9,	0xFA,	0x7A,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179688953803000975469142727125494 = 0x80000000000000000000000000000000a6f7cef4e40a02e8c38b150467ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0x46,	0x50,	0xB1,	0x38,	0x8C,	0x2E,	0xA0,	0x40,	0x4E,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order - 1000 000 000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558178689953803000975469142727125494 = 0x80000000000000000000000000000000a6f7ce2b24f319ae59443a2ee7ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0xEE,	0xA2,	0x43,	0x94,	0xE5,	0x9A,	0x31,	0x4F,	0xB2,	0xE2,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //---------------------------------------------------------------------------------------------------------

  //scalar = Half_Order + 100 000 000 000 = 3618502788666131106986593281521497120428558179689953803000975469242727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c0931a476f1d1f6
  //unsigned char scalar[32] = {0xF6,	0xD1,	0xF1,	0x76,	0xA4,	0x31,	0x09,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 = 3618502788666131106986593281521497120428558179689953803000975470142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c093276031ff9f6
  //unsigned char scalar[32] = {0xF6,	0xF9,	0x1F,	0x03,	0x76,	0x32,	0x09,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00, 0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 = 3618502788666131106986593281521497120428558179689953803000976469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b2c0cbf0bd34169f6
  //unsigned char scalar[32] = {0xF6,	0x69,	0x41,	0xD3,	0x0B,	0xBF,	0x0C,	0x2C,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00, 0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 000  = 3618502788666131106986593281521497120428558179689953803001975469142727125494 = 0x80000000000000000000000000000000a6f7cef517bce6b39e9e840d5dee9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0xDE,	0xD5,	0x40,	0xE8,	0xE9,	0x39,	0x6B,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179689953804000975469142727125494 = 0x80000000000000000000000000000000a6f7cef517bcea161d2df530d1ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x1A,	0x0D,	0x53,	0xDF,	0xD2,	0x61,	0xA1,	0xCE,	0x7B,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179689954803000975469142727125494 = 0x80000000000000000000000000000000a6f7cef517ca22d47d7fe7acf7ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0xCF,	0x7A,	0xFE,	0xD7,	0x47,	0x2D,	0xA2,	0x7C,	0x51,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558179690953803000975469142727125494 = 0x80000000000000000000000000000000a6f7cef54b6fca7cbd9b1ca167ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0x16,	0xCA,	0xB1,	0xD9,	0xCB,	0xA7,	0xFC,	0xB6,	0x54,	0xEF,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //scalar = Half_Order + 1000 000 000 000 000 000 000 000 000 000 = 3618502788666131106986593281521497120428558180689953803000975469142727125494 = 0x80000000000000000000000000000000a6f7cfbf0a86b3b727e1f776e7ae9f6
  //unsigned char scalar[32] = {0xF6,	0xE9,	0x7A,	0x6E,	0x77,	0x1F,	0x7E,	0x72,	0x3B,	0x6B,	0xA8,	0xF0,	0xFB,	0x7C,	0x6F,	0x0A,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x00,	0x08};

  //-----------------------------------------------------------------------------------------------------------------------------

  // Gx coordinate of base point = 15112221349535400772501151409588531511454012693041857206046113283949847762202
  field_element Gx[32] = {{0x1A,	0xD5,	0x25,	0x8F,	0x60,	0x2D,	0x56,	0xC9,	0xB2,	0xA7,	0x25,	0x95,	0x60,	0xC7,	0x2C,	0x69,	0x5C,	0xDC,	0xD6,	0xFD,	0x31,	0xE2,	0xA4,	0xC0,	0xFE,	0x53,	0x6E,	0xCD,	0xD3,	0x36,	0x69,	0x21}};
  //Gy coordinate of base point yeilds from -------> ((u-1) * pow(u+1, prime-2, prime)) % prime, where u = 9
  field_element Gy[32] = {{0x58,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66,	0x66}};
  field_element Gz[32] = {{0x01}};

  // Initialize timer and start measuring
  TCCR1A = 0;           // Reset Timer 1 control register A
  TCCR1B = bit(CS10);   // Set clock source to no prescaler (start Timer 1)
  TCNT1 = 0;            // Reset Timer 1 counter

  //scalar[0] &= 248;
  //scalar[31] &= 127;
  //scalar[31] |= 64;

  ed25519_ladder(Gx, Gy, Gz, scalar);

  // Stop Timer 1 and read the number of cycles
  unsigned int cycles = TCNT1;  // Read the timer value
  Serial.print("Number of cycles = ");
  Serial.println(cycles - 1);   // Subtract 1 for the overhead of starting/stopping the timer

  Serial.print("Result = ");
  for (int i = 0; i < 32; i++)
  {
    Serial.print(Gx->Ed[i], HEX);
  }
  Serial.println();
}

void ed25519_ladder(field_element *xr, field_element *yr, field_element *zr, unsigned char scalar[32])
{
  field_element p_coord[6];  // We need 6 working arrays for two points for Ed25519 ladder
  unsigned char Current_bit, swap = 0;
  signed char j = 6;
  p_coord[0] = *xr;  // X1
  p_coord[1] = *yr;  // Y1
  p_coord[2] = *zr;  // Z1

  p_coord[3] = {{0x00}}; // X0 = 0
  p_coord[4] = {{0x01}}; // Y0 = 1
  p_coord[5] = {{0x01}}; // Z0 = 1

  // Check if scalar >= Half_Order
  field_element subResult;
  unsigned char borrow;
  cli();  // Disable interrupts
  borrow = scalar_sub_halforder(subResult.Ed, scalar);
  sei();  // Enable interrupts

  // If borrow == 0, scalar is greater than or equal to Half_Order

  // Check if scalar >= Half_Order
  if (borrow == 0)
  {
    // The scalar was >= Half_Order, and the subtraction was already done in rt. So, update the scalar with the result from rt
    //(which now holds scalar - Half_Order)
    for (int i = 0; i < 32; i++)
    {
      scalar[i] = subResult.Ed[i];
    }
    // Recursively call the ladder with the updated scalar

    ed25519_ladder(xr, yr, zr, scalar);
    // After laddering, reflect the X affine coordinate by subtracting from prime
    cli();
    prime_sub_bignum(xr->Ed, xr->Ed);
    sei();

    // Immediate return to avoid redundant calculations
    return;  // <-- This line prevents the function from continuing after the recursive call
  }

  for (signed char i = 31; i >= 0; i--)
  {
    while (j >= 0)
    {
      Current_bit = 1 & (scalar[i] >> j);
      swap ^= Current_bit;
      cswap(p_coord, swap);
      swap = Current_bit;
      ladderstep(p_coord);
      j -= 1;
    }
    j = 7;
  }
  cswap(p_coord, swap);
  *xr = p_coord[3];  // X0
  *yr = p_coord[4];  // Y0
  *zr = p_coord[5];  // Z0

  // Affine coordinate calculation (inserted inside the function)
  Mul_Inverse(zr, zr);        // Invert Z coordinate
  mul(xr, xr, zr);       // Multiply X by the inverse of Z to get X affine
  finalize_result(xr);            // Freeze the X affine coordinate
}

void cswap(field_element *p_coord, char b)
{
  field_element temp_store;

  conditional_swap(&temp_store, p_coord, b);
  conditional_swap(p_coord, p_coord + 3, b);
  conditional_swap(p_coord + 3, &temp_store, b);

  conditional_swap(&temp_store, p_coord + 1, b);
  conditional_swap(p_coord + 1, p_coord + 4, b);
  conditional_swap(p_coord + 4, &temp_store, b);

  conditional_swap(&temp_store, p_coord + 2, b);
  conditional_swap(p_coord + 2, p_coord + 5, b);
  conditional_swap(p_coord + 5, &temp_store, b);
}

void conditional_swap(field_element *r, const field_element *x, unsigned char b)
{
  unsigned long mask = b;
  mask = -mask;
  for (unsigned char i = 0; i < 32; i++)
  {
    r->Ed[i] ^= mask & (x->Ed[i] ^ r->Ed[i]);
  }
}

void finalize_result(field_element *r)
{
  unsigned char bwr;
  field_element tempResult;
  bwr = bignum_sub_prime(tempResult.Ed, r->Ed);
  conditional_swap(r, &tempResult, 1 - bwr);
}

void mul(field_element *r, const field_element *x, const field_element *y)
{
  unsigned char t[64];
  cli();  // Disable interrupts
  Assembly_Ed_mul(t, x->Ed, y->Ed);
  sei();  // Enable interrupts
  modulo(r, t);
}

void square(field_element *r, const field_element *x)
{
  unsigned char t[64];
  cli();  // Disable interrupts
  Assembly_Ed_square(t, x->Ed);
  sei();  // Enable interrupts
  modulo(r, t);
}

void ladderstep(field_element *p_coord)
{
  field_element A, B, C, D, E, F, G, H, J;
  field_element temp1, temp2;

  // Define variables for the current state and the next state
  field_element *X0 = p_coord + 3; //X0
  field_element *Y0 = p_coord + 4; //Y0
  field_element *Z0 = p_coord + 5; //Z0

  field_element *X1 = p_coord;     //X1
  field_element *Y1 = p_coord + 1; //Y1
  field_element *Z1 = p_coord + 2; //Z1

  // 1. A = Z0 * Z1
  mul(&A, Z0, Z1);  // Common Z coordinate component.
  // 2. B = A^2
  square(&B, &A);  // B = A^2
  // 3. C = X0 * X1
  mul(&C, X0, X1);  // C = X0 * X1
  // 4. D = Y0 * Y1
  mul(&D, Y0, Y1);  // D = Y0 * Y1
  // 5. E = d * C * D
  mul(&E, &C, &D);  // E = C * D
  mul(&E, &E, &d);  // E = d * C * D
  // 6. F = B - E
  sub(&F, &B, &E);  // F = B - E
  // 7. G = B + E
  add(&G, &B, &E);  // G = B + E
  // 8. (X0 + Y0)
  add(&temp1, X0, Y0);  // temp1 = X0 + Y0
  // 9. (X1 + Y1)
  add(&temp2, X1, Y1);  // temp2 = X1 + Y1
  // 10. (X0 + Y0) * (X1 + Y1)
  mul(&temp1, &temp1, &temp2);  // temp1 = (X0 + Y0) * (X1 + Y1)
  // 11. (X0 + Y0) * (X1 + Y1) - C - D
  sub(&temp1, &temp1, &C);  // temp1 = (X0 + Y0) * (X1 + Y1) - C
  sub(&temp1, &temp1, &D);  // temp1 = (X0 + Y0) * (X1 + Y1) - C - D

  // 12. X2 = A * F * ( (X0 + Y0) * (X1 + Y1) - C - D )
  mul(&temp1, &F, &temp1);  // temp1 = F * ((X0 + Y0) * (X1 + Y1) - C - D)
  mul(X1, &A, &temp1);  // X2 = A * temp1
  // 13. Y2 = A * G * (D - a * C)
  mul(&temp1, &a, &C);  // temp1 = a * C
  sub(&temp1, &D, &temp1); // temp1 = D - a * C
  mul(&temp1, &G, &temp1);  // temp1 = G * (D - a * C)
  mul(Y1, &A, &temp1);  // Y2 = A * temp1
  // 14. Z2 = F * G
  mul(Z1, &F, &G);  // Z2 = F * G
  // Doubling in Projective Twisted Coordinates for (X0, Y0, Z0)

  // 15. B = (X0 + Y0)^2
  add(&temp1, X0, Y0);  // temp1 = X0 + Y0
  square(&B, &temp1);  // B = (X0 + Y0)^2
  // 16. C = X0^2
  square(&C, X0);  // C = X0^2
  // 17. D = Y0^2
  square(&D, Y0);  // D = Y0^2
  // 18. E = a * C
  mul(&E, &a, &C);  // E = a * C
  // 19. F = E + D
  add(&F, &E, &D);  // F = E + D
  // 20. H = Z0^2
  square(&H, Z0);  // H = Z0^2
  // 21. J = F - 2 * H
  add(&temp1, &H, &H);  // temp1 = 2 * H
  sub(&J, &F, &temp1);  // J = F - 2 * H
  // 22. X0 = (B - C - D) * J
  sub(&temp1, &B, &C);  // temp1 = B - C
  sub(&temp1, &temp1, &D);  // temp1 = B - C - D
  mul(X0, &temp1, &J);  // X0 = temp1 * J
  // 23. Y0 = F * (E - D)
  sub(&temp1, &E, &D);  // temp1 = E - D
  mul(Y0, &F, &temp1);  // Y0 = F * temp1
  // 24. Z0 = F * J
  mul(Z0, &F, &J);  // Z0 = F * J
}

void Mul_Inverse(field_element *r, const field_element *x)
{
  field_element z2;
  field_element z11;
  field_element z2_10_0;
  field_element z2_50_0;
  field_element z2_100_0;
  field_element t0;
  field_element t1;
  unsigned char i;

  square(&z2, x);
  square(&t1, &z2);
  square(&t0, &t1);
  mul(&z2_10_0, &t0, x);
  mul(&z11, &z2_10_0, &z2);
  square(&t0, &z11);
  mul(&z2_10_0, &t0, &z2_10_0);
  square(&t0, &z2_10_0);
  square(&t1, &t0);
  square(&t0, &t1);
  square(&t1, &t0);
  square(&t0, &t1);
  mul(&z2_10_0, &t0, &z2_10_0);
  square(&t0, &z2_10_0);
  square(&t1, &t0);
  for (i = 2; i < 10; i += 2) {square(&t0, &t1); square(&t1, &t0);}
  mul(&z2_50_0, &t1, &z2_10_0);
  square(&t0, &z2_50_0);
  square(&t1, &t0);
  for (i = 2; i < 20; i += 2) {square(&t0, &t1); square(&t1, &t0);}
  mul(&t0, &t1, &z2_50_0);
  square(&t1, &t0);
  square(&t0, &t1);
  for (i = 2; i < 10; i += 2) {square(&t1, &t0); square(&t0, &t1);}
  mul(&z2_50_0, &t0, &z2_10_0);
  square(&t0, &z2_50_0);
  square(&t1, &t0);
  for (i = 2; i < 50; i += 2) {square(&t0, &t1); square(&t1, &t0);}
  mul(&z2_100_0, &t1, &z2_50_0);
  square(&t1, &z2_100_0);
  square(&t0, &t1);
  for (i = 2; i < 100; i += 2) {square(&t1, &t0); square(&t0, &t1);}
  mul(&t1, &t0, &z2_100_0);
  square(&t0, &t1);
  square(&t1, &t0);
  for (i = 2; i < 50; i += 2) {square(&t0, &t1); square(&t1, &t0);}
  mul(&t0, &t1, &z2_50_0);
  square(&t1, &t0);
  square(&t0, &t1);
  square(&t1, &t0);
  square(&t0, &t1);
  square(&t1, &t0);
  mul(r, &t1, &z11);
}