/** This code is in the public domain 
 If you don't want to link libstdc++ with your binary, compile like this:
  gcc -c treeverse.c
  and link with gcc the object treeverse.o.

  -DDEBUG: print debuging infos on what treeverse is doing
  -DMAIN: define a main routine to test the code
   gcc -o treeverse -DMAIN -DDEBUG treeverse.cc

*/

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

typedef signed long long int llint;

#ifndef _STACK_SIZE
#define _STACK_SIZE 200
#endif

#define mystack(type,len) \
typedef struct {\
        type array[len];\
        int pos;\
} stack;\
static type *stack_look(stack *s) { return &s->array[s->pos]; }\
static type *stack_pop(stack *s) { s->pos--; return &s->array[s->pos+1]; }\
static void stack_push(stack *s, type *t) { s->pos++; s->array[s->pos] = *t; }\
static int stack_empty(stack *s) { return s->pos == -1; }


/** Compute the factorial of n
  @param n a long long integer
 */
llint fact(llint n) {
    assert(n >= 0);
    if (n <= 1)
        return 1;
    return fact(n-1)*n;
}
/** Return  /     \ i.e the binomial law.
            | s+t |
            |  t  |
            \     /
  @param s number of snapshots
  @param t number of reexecution
 */
unsigned int beta(unsigned int s,unsigned int t)
{
    if (s > 0 && t > 0) {
       unsigned int ret = 1,i;
       for (i = 1; i <= s; ++i) {
            ret *= t+i;
            ret /= i;
       }
       return ret;
    }
    return fact(s+t)/(fact(s)*fact(t));
}

/** Find the biggest t so that beta(s,t) < m
  @param s the number of checkpoints that can be used
  @param l the length of the evolution to reverse */
static unsigned int find_max_t_under (unsigned int s, unsigned int l) {
    llint c;
    unsigned int t = 0;
    static unsigned int last = 0;


    if (l <= s)
        return 0;
    /* Try to begin with the last result... */
    if (last > 0) 
    {
        t = last;
        if (beta (s, t) > l) 
        {
            do {
                t--;
                c = beta (s, t);
            } while (t > 0 && c > l);
            return t+1;
        }
    }
    else {

        do {
            t++;
            c = beta (s, t);
        } while (c < l);
        return t-1;
    }
    return t;
}  

unsigned int choose_ckp_pos (unsigned int s, unsigned int m) {
    unsigned int t = find_max_t_under (s, m);

    return beta (s, t);
}

struct arg 
{
    // member variables
    int s, m, start, state, ckp_pos;
};
mystack(struct arg, 100)
static struct arg* arg(int a, int b, int c, int state) {
        static struct arg t = { 0, 0, 0, 0, 0 };
        t.s = a; t.m = b; t.start = c;
        return &t;
}

enum { rev=1, finish, forw, pop, look, takeshot };


#ifdef DEBUG
#define debug(args...) printf(args)
#else
#define debug(args...)
#endif

/** A non-reentrant implementation of the treeverse
 checkpointing scheme, this function use a global state to
 drive a little automat, the first
 call initialize it with the value in s and m.
 out is an output parameter when there is forward sweeps
 to do. 
 @param s the maximum number of snapshots to do
 @param m the length of the evolution */
int _treeverse(int *s, int *m, int *out) {
    int start;
    static stack args = { .pos = -1 };


    if (stack_empty(&args))
        stack_push(&args, arg(*s,*m,0,0));

    while (! stack_empty(&args)) {
        struct arg *last = stack_look(&args);

        s = &last->s;
        m = &last->m;
        start = last->start;
        int *state = &last->state;
        int *ckp_pos = &last->ckp_pos;

        switch (*state) {
            case 0:
                if (*m == 1) {
                    *state = 10;
                    debug("Reverse %d\n", start+*m);
                    return rev;
                }
                if (*s > 1) {
                    *ckp_pos = choose_ckp_pos (*s,*m);
                    if (*ckp_pos != 0) {
                        debug("Forward until %d\n", start+*ckp_pos);
                        *out = *ckp_pos;
                        *state = 1;
                        return forw;
            case 1:
                        debug("Takeshot %d\n", start+*ckp_pos);
                        *state = 2;
                        return takeshot;
            case 2:
                        *state = 3;
                        stack_push(&args, arg(*s-1, *m-*ckp_pos, start+*ckp_pos, 0));
                        continue;
            case 3:
                        *state = 4;
                        debug("Remove shot %d\n", start+*ckp_pos);
                        return pop;
            case 4:
                        *state = 5;
                        debug("Look at shot %d\n", start);
                        return look;
            case 5:
                        *state = 10;
                        stack_push(&args, arg(*s, *ckp_pos, start, 0));
                        continue;
                    } else
                        debug ("ca arrive\n");
                } else {
                    static int j;

                    for (j = *m; j > 0; --j) 
                    {
                        if (j != *m) 
                        {
                            *state = 6;
                            debug("Look at shot %d\n", start);
                            return look;
                        }
            case 6:
                        if (j != 1) 
                        {
                            *state = 7;
                            debug("Forward until %d\n", start+j-1);
                            *out = j-1;
                            return forw;
                        }
            case 7:
                        *state = 8;
                        debug("Reverse %d\n", start+j);
                        return rev;
            case 8:
                        ;;
                    }
                }
            case 10:
                ;;
        }
        stack_pop(&args);
    }
    return finish;
}


#ifdef MAIN
int main (int argc, char **argv) {
    int s, l;
    assert(argc > 2);
    s = atoi(argv[1]);
    l = atoi(argv[2]);
    int ret, out;
    int etat[100];
    int where = 1;
    int last[100] = {1}, j = 0;

    memset (etat, 0, sizeof(etat));

    while ((ret = _treeverse(&s,&l,&out)) != finish) {
        switch(ret) {
            case forw:
                {
                    int i;
                    for (i = 0; i < out; ++i) {
                    #ifdef DEBUG
                        etat[where]++;
                    #endif    
                        where++;
                    }
                    break;
                }
            case rev:
                //reverse
                break;
            case takeshot:
                ++j;
                last[j] = where;
                //PUSH
                break;
            case pop:
                //POP
                where = last[j];
                --j;
                break;
            case look:
                //LOOK
                where = last[j];
                break;
        }
    }
    #ifdef DEBUG
    {
      int u;
      for (u = 0; u < 60; ++u)
        printf("%d %d\n", u, etat[u]);
    }
    #endif    
    printf ("max t under (s=%d,l=%d) = %d\n", s, l, find_max_t_under(s,l)+1);
    return 0;
}

#endif
#ifdef DUMPENUM
int main() {
  char *name[] = { "REVERSE", "FINISH", "FORWARD", "POP",
  "LOOK", "TAKESHOT" };
  int value[] = { rev, finish, forw, pop, look, takeshot };

  for (int i = 0; i < sizeof(name)/sizeof(name[0]); ++i)
      printf("      %s = %d\n", name[i], value[i]);
}
#endif