#!c

int cohendiv(int x, int y){
 assert(x>0 && y>0);

  int q=0; int r=x;

  while (r>=y) {
    int a=1; int b=y;
    while (true){
      //x==q*y+r, b==y*a, , r>=0, r>=2*y*a
      if(!(r >= 2*b)) break;
      a = 2*a; b = 2*b;
    }
    r=r-b; q=q+a;
  }

  //x == q*y+r, 
  //-x <= -1, r - y <= -1, -r <= 0
  //q == x % y, r == x % y
  return q;
}

#!c

int divbin(int A, int B){
  /*
  a division algorithm, by Kaldewaij
  returns A//B
  */
  assert(A > 0 && B > 0);
  int q = 0; int r = A; int b = B;
  while(1){
    if (!(r>=b)) break;
    assert(q == 0 && r == A  &&  A >= 0  &&  B > 0);
    b=2*b;
  }

  while(1){
    if (!(b!=B)) break;
    //(A == q*b + r && r >= 0);
    q = 2*q; b = b/2;
    if (r >= b) {
      q = q + 1; r = r - b;
    }
  }
  //q==A/B, r==A%B
  return q;
}

#!c

int hard(int A, int B) {
  //hardware integer division program, by Manna
  //returns q==A//B

  assert(A >= 0);
  assert(B >= 1);

  int r,d,p,q;

  r=A;
  d=B;
  p=1;
  q=0;

  while(1){
    if (!(r >= d)) break;
    // A >= 0 && B > 0 && q == 0 && r == A && d == B*p
    d = 2 * d;
    p  = 2 * p;
  }

  while(1){
    if (!(p!=1)) break;
    // A == q*B+r
    // d == B*p
    d=d/2;
    p=p/2;

    if(r>=d){
      r=r-d;
      q=q+p;
    }
  }

  // r == A % B
  // q == A / B
  return q;
}

#!c 
int mannadiv (int x1, int x2){
  //Division algorithm from
  //"Z. Manna, Mathematical Theory of Computation, McGraw-Hill, 1974"
  //return x1 // x2

  int y1,y2,y3;
  y1 = 0; y2 = 0; y3 = x1;
  while(1){
    if(!(y3 != 0)) break;
    // y1* x2 + y2 + y3 == x1
    if (y2 + 1 == x2) {
      y1 = y1 + 1;
      y2 = 0;
      y3 = y3 - 1;
    }
    else {
      y2 = y2 + 1;
      y3 = y3 - 1;
    }
  }

  // y1 == x1 / x2
  return y1;
}

#!c
int sqrt1(int n){
  /* computing the floor of the square root of a natural number */
  assert(n >= 0);

  int a,s,t;
  a=0;
  s=1;
  t=1;

  while(1){

    if(!(s <= n)) break;

    assert(a*a <= n);
    assert(t == 2*a + 1);
    assert(s == (a + 1)*(a + 1));

    a=a+1;
    t=t+2;
    s=s+t;
  }

  // assert(a==(int)floor(sqrt(n)));
  return a;
}

#!c

int dijkstra(int n){
  /* program computing the floor of the square root, by Dijkstra */
  assert(n>=0);

  int p,q,r,h;

  p=0;
  q=1;
  r=n;
  h= 0;

  while(1){
    if (!(q<=n)) break;

    q=4*q;
    assert(n >= 0 && p == 0 && r==n );
  }
  //q == 4^n


  while(1){

    if (!(q!=1)) break;

    assert(r < 2*p + q);
    assert(p*p + r*q == n*q);

    //h^2*p - 4*h*n*q + 4*h*q*r + 4*n*p*q - p*q^2 - 4*p*q*r == 0
    //h^2*n - h^2*r - 4*h*n*p + 4*h*p*r + 4*n^2*q - n*q^2 - 8*n*q*r + q^2*r + 4*q*r^2 == 0
    //h^3 - 12*h*n*q - h*q^2 + 12*h*q*r + 16*n*p*q - 4*p*q^2 - 16*p*q*r == 0

    q=q/4;
    h=p+q;
    p=p/2;
    if (r>=h){
      p=p+q;
      r=r-h;
    } 
  }

  assert(p == (int)floor(sqrt(n)));
  return p;
}

#!c

int freire1 (float a){
  //algorithm for finding the closest integer to square root,
  //from  www.pedrofreire.com/crea2_en.htm?

  float x;
  int r;

  x=a/2.0;
  r=0;

  while(1){
    if (!(x>r)) break;
    x=x-r;
    r=r+1;
  }

  assert(r==(int)round(sqrt(a)));

  return r;
}

#!c

int freire2(float a){
  //algorithm for finding the closest integer to the cubic root,
  //from www.pedrofreire.com/crea2_en.htm? */ 
  //a is a float

  float x,s;
  int r;

  x=a;
  r=1;
  s=3.25;

  while(1){

    if (!(x-s > 0.0)) break;

    //(4*r*r*r - 6*r*r + 3*r) + (4*x - 4*a) == 1
    //4*s - 12*r*r == 1

    x = x - s;
    s = s + 6 * r + 3;
    r = r + 1;
  }

  //r == (int)round(pow(a,(1.0/3.0)))

  return r;
}

#!c

int cohencu(int a){
  /* printing consecutive cube, by Cohen */

  int n,x,y,z;
  n=0; x=0; y=1; z=6;

  while(1){

    if (!(n<=a)) break;
    //z == 6*n + 6);
    //y == 3*n*n + 3*n + 1);
    //x == n*n*n);

    n=n+1;
    x=x+y;
    y=y+z;
    z=z+6;
  }
  return x;
}

#!c
int euclidex1 (int x, int y){

  int a,b,p,q,r,s;

  a=x;
  b=y;
  p=1;
  q=0;
  r=0;
  s=1;

  while(1){ 
    if (!(a!=b)) break;
    //1 == p*s - r*q, a == y*r + x*p, b == x*q + y*s
    if (a>b){
      a = a-b; 
      p = p-q; 
      r = r-s;
    }
    else {
      b = b-a; 
      q = q-p; 
      s = s-r;}
  }

  return a;
}

#!c
int euclidex2 (int x, int y){
  /* extended Euclid's algorithm */

  int a,b,p,q,r,s;

  a=x; b=y; p=1; q=0; r=0; s=1;

  while(1){

    if (!(b!=0)) break;

    int c=a,k=0;

    while(1){

      if (!(c>=b)) break;

      //a == k*b+c);
      //a == y*r+x*p);
      //b == x*q+y*s);

      c=c-b;
      k=k+1;
    }

    a=b;
    b=c;

    int temp;
    temp=p;
    p=q;
    q=temp-q*k;
    temp=r;
    r=s;
    s=temp-s*k;
  }

  //mygcd(x,y)==y*r+x*p );
  return a;
}

#!c
/* def xgcd(a,b): */
/* if b == 0: */
/* return [1,0,a] */
/* else: */
/* x,y,d = xgcd(b, a%b) */
/* return [y, x - (a//b)*y, d] # Note that a//b is floor(a/b) in Python. */

#!c

int euclidex3(int x, int y){
  /* extended Euclid's algorithm */

  int a,b,p,q,r,s;

  a=x; b=y;  p=1;  q=0;  r=0;   s=1;

  //a==y*r+x*p
  //b==x*q+y*s

  while(1){ 

    if(!(b!=0)) break;

    int c=a,k=0;

    while(1){
      if (!(c>=b)) break;
      int d=1,v=b;

      while(1){

    if (!(c>= 2*v)) break;
        //a == y*r+x*p, b == x*q+y*s, a == k*b+c, v == b*d

        d = 2*d;
        v = 2*v;

      }
      c=c-v;
      k=k+d;

      //a==y*r+x*p, b==x*q+y*s, a==k*b+c

    }

    a=b;
    b=c;
    int temp;
    temp=p;
    p=q;
    q=temp-q*k;
    temp=r;
    r=s;
    s=temp-s*k;
  }

  //mygcd(x,y)==y*r+x*p );

  return a;
}

#!c
int lcm1(int a, int b){
  /* algorithm for computing simultaneously the GCD and the LCM, 
     by Sankaranarayanan */

  int x,y,u,v;

  x=a;
  y=b;
  u=b;
  v=0;

  while(1){
    if (!(x!=y)) break;

    //x*u + y*v == a*b

    while(1){
      if (!(x>y)) break;

      x=x-y;
      v=v+u;
    }

    while(1){
      if (!(x<y)) break;

      y=y-x;
      u=u+v;
    }
  }
  //x==gcd(a,b)

  int r = u+v; 
  return r; //lcm
}

#!c


int lcm2 (int a, int b){
  /* algorithm for computing simultaneously the GCD and the LCM, by Dijkstra */

  int x,y,u,v;

  x=a;
  y=b;
  u=b;
  v=a;

  while(1){ 
    if (!(x!=y)) break;

    //x*u + y*v == 2*a*b
    if (x>y){
      x=x-y;
      v=v+u;
    }
    else {
      y=y-x;
      u=u+v;
    }

  }

  //x==gcd(a,b)
  int r = (u+v)/2;
  return r; //lcm
}

#!c

int prodbin (int a, int b){
  /* algorithm for computing the product of two natural numbers */
  /* shift_add */

  assert(b>=0);

  int x,y,z;

  x = a;
  y = b;
  z = 0;

  while(1){ 
    if (!(y!=0)) break;
    // z+x*y==a*b

    if (myodd(y)){
      z = z + x;
      y = y - 1;
    }
    x = 2 * x;
    y = y / 2;

  }

  // z == a*b
  return z; 
}


## Another product of two natural numbers (prod4br)

int prod4br(int x, int y){
  /* algorithm for computing the product of two natural numbers */

  int a,b,p,q;

  a = x;
  b = y;
  p = 1;
  q = 0;

  while(1){ 
    if (!(a!=0 && b!=0)) break;

    //q+a*b*p==x*y)

    if (a % 2 ==0 && b % 2 ==0 ){
      a =a/2;
      b = b/2;
      p = 4*p;
    }
    else if (a % 2 ==1 && b % 2 ==0 ){
      a =a-1;
      q = q+b*p;
    }
    else if (a % 2 ==0 && b % 2 ==1 ){
      b =b-1;
      q = q+a*p;
    }
    else {
      a =a-1;
      b =b-1;
      q = q+(a+b+1) * p;
    }

  }

  //q == x*y
  return q; 

}

#!c
int fermat1(int A, int R){
  // program computing a divisor for factorisation, by Knuth 4.5.4 Alg C ?

  // TODO: check for infinite loop


  //A >= 1
  //(R-1)*(R-1) < A
  //A <= R*R
  //A%2 ==1

  int u,v,r;
  u=2*R+1;
  v=1;
  r=R*R-A;

  // 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1
  while(1){
    if (!(r!=0)) break;

    //4*(A+r) == u*u-v*v-2*u+2*v

    while(1){
      if (!(r>0)) break;
      r=r-v;
      v=v+2;
    }

    while(1){
      if (!(r<0)) break;
      r=r+u;
      u=u+2;
    }

  }

  //u != v
  int o = (u-v)/2;
  return o;
}

#!c

int fermat2(int A, int R){
  /* program computing a divisor for factorisation, by Bressoud */

  //A >= 1
  //(R-1)*(R-1) < A
  //A <= R*R
  //A%2 ==1

  int u,v,r;

  u=2*R+1;
  v=1;
  r=R*R-A;

  //// 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1
  while(1){
    if (!(r!=0)) break;

    //4*(A + r) == u*u - v*v - 2*u + 2*v

    if (r>0){
      r=r-v;
      v=v+2;
    }
    else{
      r=r+u;
      u=u+2;
    }
  }

  //u!=v
  int o = (u-v)/2;
  return o;
}

#!c

int knuth(int n, int a){
  //algorithm searching for a divisor for factorization, by Knuth

  int r,k,q,d,s,t;

  d=a;
  r= n % d;
  t = 0;

  if (d <= 2){
    printf("#E: d - 2 <= 0.  Use a diff val for d\n");
    return 0;
  }
  k=n % (d-2);
  q=4*(n/(d-2) - n/d);
  s=sqrt(n);

  while(1){
    //assert(d*d*q - 2*q*d - 4*r*d + 4*k*d  + 8*r == 8*n);
    //assert(k*t == t*t);
    //assert(d*d*q - 2*d*q - 4*d*r + 4*d*t + 4*a*k - 4*a*t - 8*n + 8*r == 0);
    //assert(d*k - d*t - a*k + a*t == 0); 
    if (!((s>=d)&&(r!=0))) break;

    if (2*r-k+q<0){
      t=r;
      r=2*r-k+q+d+2;
      k=t;
      q=q+4;
      d=d+2;
    } 
    else if ((2*r-k+q>=0)&&(2*r-k+q<d+2)){
      t=r;
      r=2*r-k+q;
      k=t;
      d=d+2;
    }
    else if ((2*r-k+q>=0)&&(2*r-k+q>=d+2)&&(2*r-k+q<2*d+4)){
      t=r;
      r=2*r-k+q-d-2;
      k=t;
      q=q-4;
      d=d+2;
    }
    else {/* ((2*r-k+q>=0)&&(2*r-k+q>=2*d+4)) */
      t=r;
      r=2*r-k+q-2*d-4;
      k=t;
      q=q-8;
      d=d+2;
    }

  }
  return d;
}

#!c


void petter(k){  
  //geo
  int x = 0;
  int y = 0;

  while(1){
    if (!(y != k)) break;

    x = x + pow(y,5);
    y = y + 1;
  }
}

#!c
/*from Petter"s master thesis
  http://www2.cs.tum.edu/~petter/da/da.pdf

IMPORTANT: don't try with very large power (k) 
which makes the number *very* large and becomes problematic for computers.
e.g. 15*8 = -1732076671 on my Macbook Pro
This is a _computer_dependent_problem, not mathematics problem
*/

int geo_pw(const int a, const int z, const int k){
  //computes (z-1) * sum(z^i)[k=0..k-1]

  int i = 0; int s = 0;
  for(; i<k ; ++i) s = s + a*pow(z,i);
  return s;
}

#!c

int geo1(const int z, const int k){
  /* computes x=(z-1)* sum(z^k)[k=0..k-1] , y = z^k
     returns 1+x-y == 0

     The loop is same as geo2 but the result is different.
  */

  assert(k>0);

  int x = 1; int y = z; int c = 1;

  while(1){
    if (!(c < k)) break;

    c = c + 1;
    x = x*z + 1;
    y = y*z;

  }

  x = x *(z - 1);

  //y == pow(z,k)
  //geo_pw(1,z,k)*(z-1) == x
  //1+x-y == 0 //documented inv

  return x;
}

#!c

int geo2(const int z, const int k){
  //computes x = sum(z^k)[k=0..k-1], y = z^(k-1)
  assert (k>0);

  int x = 1; int y = 1; int c = 1;

  while(1){
    if (!(c < k)) break;

    c = c + 1;
    x = x*z + 1;
    y = y*z;
  }

  //y == pow(z,(k-1))
  //geo_pw(1,z,k) == x);

  return x;
}

#!c

int geo3(const int z, const int a, const int k){
  //computes x = sum(a*z^k)[k=0..k-1], y = z^(k-1)
  assert (k>0);

  int x = a; int y = 1;  int c = 1;
  while(1){
    if (!(c < k)) break;

    c = c + 1;
    x = x*z + a;
    y = y*z;

  }

  //y == pow(z,(k-1))
  //geo_pw(a,z,k) == x

  return x;
}

#!c

/*
  # // yielding x=y
  # //todo: too easy, not worth doing
  # def ps1 ( String []
  # y = 0
  # x = 0
  # while( args !=null){
  # y=y +1
  # x=x +1
*/


int ps_pw(const int k, const int z){
  // computes sum(i^z)[k=1..k]
  int i = 1; int s = 0;
  for(; i <= k; ++i) s = s + pow(i,z);
  return s;
}

#!c


int ps1(const int k){
  int y = 0; int x = 0;int c = 0;
  while(1){
    if (!(c < k)) break;

    c = c + 1;
    y = y + 1;
    x = x + 1;
  }

  //y == 1*k
  //x == ps_pw(k,0)
  return x;
}

#!c
int ps2(int k){

  int y = 0; int x = 0; int c = 0;

  while(1){
    //assert((y*y) - 2*x + y == 0);
    if (!(c < k)) break;

    c = c + 1;
    y=y +1;
    x=x+y;

  }

  //y == 1*k
  //x == ps_pw(k,1)
  return x;
}

#!c

int ps3(int k){
  int y = 0;
  int x = 0;
  int c = 0;

  while(1){
    if (!(c < k)) break;

    c = c + 1;
    y=y +1;
    x=y*y+x;

  }

  //y == 1*k
  //x == ps_pw(k,2)
  return x;
}

#!c


int ps4 (int k){
  int y = 0;
  int x = 0;
  int c = 0;

  while(1){
    //assert(4*x-(y*y*y*y)-2*(y*y*y)-(y*y) == 0);

    if (!(c < k)) break;

    c = c +1 ;
    y=y +1;
    x=y*y*y+x;

  }

  //y == 1*k
  //x == ps_pw(k,3)
  return x;
}

#!c


int ps5 (int k){
  int y = 0;
  int x = 0;
  int c = 0;

  while(1){
    if (!(c < k)) break;

    c = c +1 ;
    y=y +1;
    x=y*y*y*y+x;
  }

  //y == 1*k
  //x == ps_pw(k,4)
  return x;
}

#!c

int ps6 (int k){
  int y = 0;
  int x = 0;
  int c = 0;

  while(1){
    if (!(c < k)) break;

    c = c +1 ;
    y=y +1;
    x=y*y*y*y*y+x;
  }

  //y == 1*k
  //x == ps_pw(k,5)
  return x;
}