void
Mopti_shladd (L_Func * fn)
{
  L_Cb *cb1;
  L_Cb *cb2;

  L_alloc_danger_ext = L_create_alloc_pool ("L_Danger_Ext",
					    sizeof (struct L_Danger_Ext), 64);

  MOD ("Flow analysis DOM AD AE LV due to shladd");
  L_do_flow_analysis (fn, DOMINATOR_CB | AVAILABLE_DEFINITION |
		      AVAILABLE_EXPRESSION | LIVE_VARIABLE);

  for (cb1 = fn->first_cb; cb1 != NULL; cb1 = cb1->next_cb)
    while (M_local_shift_add_merge (cb1));

  L_compute_danger_info (fn);
  MOD ("Flow analysis DOM AD AE LV due to shladd");
  L_do_flow_analysis (fn, DOMINATOR_CB | AVAILABLE_DEFINITION |
		      AVAILABLE_EXPRESSION | LIVE_VARIABLE);

  for (cb1 = fn->first_cb; cb1 != NULL; cb1 = cb1->next_cb)
    {
      for (cb2 = fn->first_cb; cb2 != NULL; cb2 = cb2->next_cb)
	{
	  if (cb1->id == cb2->id)
	    continue;
	  while (M_global_shift_add_merge (cb1, cb2));
	}			/* for cb2 */
    }				/* for cb1 */
  L_delete_all_danger_ext (fn);

  L_free_alloc_pool (L_alloc_danger_ext);
  L_alloc_danger_ext = NULL;
}				/* Mopti_shladd */

int solve(int i, int j, int turn)
{
  int &ref = dp[i][j][turn];
  if(ref != -1) return ref;
  if(i == 0 && j == n-1) return ref = 0;
	
  ref = 0;
  int left = MOD(i-1,n);
  if(left != j)
  ref = (turn ? pizza_area[left] : 0) + solve(left,j,!turn);
	
  int right = MOD(j+1,n);
  if(right != i)
  {
	if(turn)
		ref = max(ref, pizza_area[right] + solve(i,right,!turn));
	else
	{
		if(ref == 0)
			ref = solve(i,right,!turn);
		else
			ref = min(ref, solve(i,right,!turn));
	}
  }
  return ref;
}

/* ========================================================================= */
unsigned int adler32(unsigned int adler, const char *buf, unsigned int len)
{
    unsigned long s1 = adler & 0xffff;
    unsigned long s2 = (adler >> 16) & 0xffff;
    int k;

    if (buf == 0) return 1L;

    while (len > 0) {
        k = len < NMAX ? (int)len : NMAX;
        len -= k;
        while (k >= 16) {
            DO16(buf);
            buf += 16;
            k -= 16;
        }
        if (k != 0) do {
            s1 += *buf++;
            s2 += s1;
        } while (--k);
        MOD(s1);
        MOD(s2);
    }
    return (s2 << 16) | s1;
}

/* note that 3rd argument must evaluate to non-NULL - returns hash value
 * of string_to_lookup to add_data() to allow it to skip doing a hash
 */
struct hashnode *
node_lookup(struct hashtable *h, const char *string_to_lookup, unsigned int *rhv)
{
	struct hashnode *r = NULL;
	unsigned int hashval = hash_djb2(string_to_lookup);
	int idx = MOD(hashval, h->maxp);
	struct hashnode *chain;

	*rhv = hashval;

	if (idx < h->p)
		idx = MOD(hashval, (2*h->maxp));

	chain = h->buckets[idx]->next;

	while (chain != h->sentinels[idx])
	{
		/* checking the equivalence of hash value first prevents
		 * expensive byte-by-byte strcmp().  Seems to improve performance. */
		if (chain->value == hashval &&
			!strcmp(chain->string, string_to_lookup))
		{
			r = chain;
			break;
		}

		chain = chain->next;
	}

	return r;
}

 void Adler32::Hash(const uint8 * buffer, uint32 size)
 {
     if (!buffer || !size) return;
     
     len += size;
     // Faster implementation
     const uint32 nMax = 5552;
     while (size >= nMax)
     {
         for (int i = 0; i < nMax / 16; i++)
         {
             DO16(buffer); buffer += 16;
         }
         MOD(a); MOD(b); size -= nMax;
     }
     
     while (size >= 16)
     {
         DO16(buffer); buffer += 16; size -= 16;
     }
     
     while (size > 0)
     {
         DO1(buffer, 0); buffer ++; size --;
     }
     MOD(a); MOD(b);
 }

void fasper(float x[], float y[], unsigned long n, float ofac, float hifac,
	float wk1[], float wk2[], unsigned long nwk, unsigned long *nout,
	unsigned long *jmax, float *prob)
{
	void avevar(float data[], unsigned long n, float *ave, float *var);
	void realft(float data[], unsigned long n, int isign);
	void spread(float y, float yy[], unsigned long n, float x, int m);
	unsigned long j,k,ndim,nfreq,nfreqt;
	float ave,ck,ckk,cterm,cwt,den,df,effm,expy,fac,fndim,hc2wt;
	float hs2wt,hypo,pmax,sterm,swt,var,xdif,xmax,xmin;

	*nout=0.5*ofac*hifac*n;
	nfreqt=ofac*hifac*n*MACC;
	nfreq=64;
	while (nfreq < nfreqt) nfreq <<= 1;
	ndim=nfreq << 1;
	if (ndim > nwk) nrerror("workspaces too small in fasper");
	avevar(y,n,&ave,&var);
	if (var == 0.0) nrerror("zero variance in fasper");
	xmin=x[1];
	xmax=xmin;
	for (j=2;j<=n;j++) {
		if (x[j] < xmin) xmin=x[j];
		if (x[j] > xmax) xmax=x[j];
	}
	xdif=xmax-xmin;
	for (j=1;j<=ndim;j++) wk1[j]=wk2[j]=0.0;
	fac=ndim/(xdif*ofac);
	fndim=ndim;
	for (j=1;j<=n;j++) {
		ck=(x[j]-xmin)*fac;
		MOD(ck,fndim)
		ckk=2.0*(ck++);
		MOD(ckk,fndim)
		++ckk;
		spread(y[j]-ave,wk1,ndim,ck,MACC);
		spread(1.0,wk2,ndim,ckk,MACC);
	}
	realft(wk1,ndim,1);
	realft(wk2,ndim,1);
	df=1.0/(xdif*ofac);
	pmax = -1.0;
	for (k=3,j=1;j<=(*nout);j++,k+=2) {
		hypo=sqrt(wk2[k]*wk2[k]+wk2[k+1]*wk2[k+1]);
		hc2wt=0.5*wk2[k]/hypo;
		hs2wt=0.5*wk2[k+1]/hypo;
		cwt=sqrt(0.5+hc2wt);
		swt=SIGN(sqrt(0.5-hc2wt),hs2wt);
		den=0.5*n+hc2wt*wk2[k]+hs2wt*wk2[k+1];
		cterm=SQR(cwt*wk1[k]+swt*wk1[k+1])/den;
		sterm=SQR(cwt*wk1[k+1]-swt*wk1[k])/(n-den);
		wk1[j]=j*df;
		wk2[j]=(cterm+sterm)/(2.0*var);
		if (wk2[j] > pmax) pmax=wk2[(*jmax=j)];
	}
	expy=exp(-pmax);
	effm=2.0*(*nout)/ofac;
	*prob=effm*expy;
	if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm);
}

void annealing(TSP *tsp)
{
    Path   p;
    int    i, j, pathchg;
    int    numOnPath, numNotOnPath;
    DTYPE    pathlen;
    int    n = tsp->n;
    double energyChange, T;

    pathlen = pathLength (tsp); 

    for (T = T_INIT; T > FINAL_T; T *= COOLING)  /* annealing schedule */
    {
        pathchg = 0;
        for (j = 0; j < TRIES_PER_T; j++)
        {
            do {
                p[0] = unifRand (n);
                p[1] = unifRand (n);
                /* non-empty path */
                if (p[0] == p[1]) p[1] = MOD(p[0]+1,n);
                numOnPath = MOD(p[1]-p[0],n) + 1;
                numNotOnPath = n - numOnPath;
            } while (numOnPath < 2 || numNotOnPath < 2); /* non-empty path */
            
            if (RANDOM() % 2) /*  threeWay */
            {
                do {
                    p[2] = MOD(unifRand (numNotOnPath)+p[1]+1,n);
                } while (p[0] == MOD(p[2]+1,n)); /* avoids a non-change */

                energyChange = getThreeWayCost (tsp, p);
                if (energyChange < 0 || RREAL < exp(-energyChange/T) )
                {
                    pathchg++;
                    pathlen += energyChange;
                    doThreeWay (tsp, p);
                }
            }
            else            /* path Reverse */
            {
                energyChange = getReverseCost (tsp, p);
                if (energyChange < 0 || RREAL < exp(-energyChange/T))
                {
                    pathchg++;
                    pathlen += energyChange;
                    doReverse(tsp, p); 
                }
            }
            // if the new length is better than best then save it as best
            if (pathlen < tsp->bestlen) {
                tsp->bestlen = pathlen;
                for (i=0; i<tsp->n; i++) tsp->border[i] = tsp->iorder[i];
            }
            if (pathchg > IMPROVED_PATH_PER_T) break; /* finish early */
        }   
        DBG("T:%f L:%f B:%f C:%d", T, pathlen, tsp->bestlen, pathchg);
        if (pathchg == 0) break;   /* if no change then quit */
    }
}

int _Dasm_ReadModRmAndSid (int nCur, BYTE * pInst)
{
  int nSize = 0 ;

  TRACE_INFO (TEXT("inst=0x%02X, mod=%d, opcode=%d, rm=%d\n"), 
	      pInst[0], OPCODE(pInst[nCur]), MOD(pInst[nCur]), RM(pInst[nCur])) ;
  
  switch( MOD(pInst[nCur]) )
    {
    case 0: // Mod==00 => No disp

      // /!\ not valid if R/M is 4
      switch( RM(pInst[nCur]) )
	{   
	case 4: // an SIB follows
	  nSize = nCur + 2 ;
	  break ;   
	case 5:  // has disp32
	  nSize = nCur + 5 ;
	  break ;
	default:
	  nSize = nCur + 1 ;
	}

      break ; 
      
    case 1: // Mod==01 => 8 bits disp

      switch( RM(pInst[nCur]) )
	{ 
	case 4: // an SIB follows
	  nSize = nCur + 3 ;
	  break ;

	default:
	  nSize = nCur + 2 ;
	}

      break ;
      
    case 2: // Mod==10 => 32 bits disp

      // /!\ not valid if R/M is 4
      nSize = nCur + 5 ;

      if( RM(pInst[nCur])==4 )
	TRACE_WARNING (TEXT("Not tested instruction decoding\n")) ;
      
      break ;
      
    case 3:

      nSize = nCur + 1 ;

      break ; 
    }

  return nSize ;
}

int is_abs_walkable(Entity** entities, int x, int y) {
    int walkable = is_walkable(entities[DIV(x)][DIV(y)]);
    if(MOD(x) != 0 && MOD(y) != 0)
        walkable = walkable && is_walkable(entities[DIV(x) + 1][DIV(y) + 1]);
    if(MOD(x) != 0)
        walkable = walkable && is_walkable(entities[DIV(x) + 1][DIV(y)]);
    if(MOD(y) != 0)
        walkable = walkable && is_walkable(entities[DIV(x)][DIV(y) + 1]);
    return walkable;
}

/**
 * @brief Mirrors a vector off another, stores results in vector.
 *
 *    @param r Resulting vector of the reflection.
 *    @param v Vector to reflect.
 *    @param n Normal to reflect off of.
 */
void vect_reflect( Vector2d* r, Vector2d* v, Vector2d* n )
{
   double dot;

   dot      = vect_dot( v, n );
   r->x     = v->x - ((2. * dot) * n->x);
   r->y     = v->y - ((2. * dot) * n->y);
   r->mod   = MOD(r->x,r->y);
   r->angle = MOD(r->x,r->y);
}

void
TSP::doThreeWay(Path p) {
    size_t count, m1, m2, m3, a, b, c, d, e, f;

    a = MOD(p[0]-1,n);
    b = p[0];
    c = p[1];
    d = MOD(p[1]+1,n);
    e = p[2];
    f = MOD(p[2]+1,n);    

    m1 = MOD(n + c - b, n) + 1;  /* num cities from b to c */
    m2 = MOD(n + a - f, n) + 1;  /* num cities from f to a */
    m3 = MOD(n + e - d, n) + 1;  /* num cities from d to e */

    count = 0;
    /* [b..c] */
    for (size_t i = 0; i < m1; i++)
        jorder[count++] = iorder[MOD(i + b, n)];

    /* [f..a] */
    for (size_t i = 0; i < m2; i++)
        jorder[count++] = iorder[MOD(i+f,n)];

    /* [d..e] */
    for (size_t i = 0; i < m3; i++)
        jorder[count++] = iorder[MOD(i+d,n)];

    /* copy segment back into iorder */
    for (size_t i = 0; i < n; i++) iorder[i] = jorder[i];
}

void
rehash_hashtable(struct hashtable *h)
{
	struct hashnode *l;
	struct hashnode *oldbucket, *oldtail;
	int newindex;

	++h->rehash_cnt;

	if (h->currentsize >= h->allocated)
		reallocate_buckets(h);

	oldbucket = h->buckets[h->p];
	oldtail = h->sentinels[h->p];
	l = oldbucket->next;

	oldbucket->next = oldtail;
	oldtail->prev = oldbucket;

	newindex = h->p + h->maxp;

	oldbucket->nodes_in_chain = 0;  /* may not be any lines in chain after rehash */
	++oldbucket->value;      /* number of splits */

	++h->p;

	if (h->p == h->maxp)
	{
		h->maxp *= 2;
		h->p = 0;
	}

	++h->currentsize;

	while (oldtail != l)
	{
		struct hashnode *t = l->next;
		int idx = MOD(l->value, h->maxp);

		if (idx < h->p)
			idx = MOD(l->value, (2*h->maxp));

		if (idx == newindex)
		{
			insert_node(h->buckets[newindex], l);
			++h->buckets[newindex]->nodes_in_chain;
		} else {
			insert_node(oldbucket, l);
			++oldbucket->nodes_in_chain;
		}

		l = t;
	}
}

/*
 *   c..b       c..b
 *    \/    =>  |  |
 *    /\        |  |
 *   a  d       a  d
 *
 *   a b 1  2   .. n-1 n c d
 *   a c n n-1  ..  2  1 c d
 */
double
TSP::getReverseCost(Path p) {

    auto a = iorder[MOD(p[0] - 1, n)];
    auto b = iorder[p[0]];
    auto c = iorder[p[1]];
    auto d = iorder[MOD(p[1] + 1, n)];

    return (D(d,b) + D(c,a) - D(a,b) - D(c,d));
    /* add cost between c and b if non symetric TSP */ 
}

/* Adjust time T by adding SECONDS.
   The absolute value of SECONDS cannot exceed 59 * INT_MAX,
   and also cannot exceed one month's worth of seconds;
   this is enough to handle any POSIX or real-life daylight-saving offset.
   Adjust only T's year, mon, mday, hour, min and sec members;
   plus adjust wday if it is defined.  */
void
adjzone (register struct tm *t,
         long seconds)
{
   int days = 0;

   /* This code can be off by a second if SECONDS is not a multiple of 60,
      if T is local time, and if a leap second happens during this minute.
      But this bug has never occurred, and most likely will not ever occur.
      Liberia, the last country for which SECONDS % 60 was nonzero,
      switched to UTC in May 1972; the first leap second was in June 1972.  */
   int leap_second = t->tm_sec == 60;
   long sec = seconds + (t->tm_sec - leap_second);
   if (sec < 0)
   {
      if ((t->tm_min -= (59 - sec) / 60) < 0
          && (t->tm_hour -= (59 - t->tm_min) / 60) < 0)
      {
         days = - ((23 - t->tm_hour) / 24);
         if ((t->tm_mday += days) <= 0)
         {
            if (--t->tm_mon < 0)
            {
               --t->tm_year;
               t->tm_mon = 11;
            }
            t->tm_mday += month_days (t);
         }
      }
   }
   else
   {
      if (60 <= (t->tm_min += sec / 60)
          && (24 <= (t->tm_hour += t->tm_min / 60)))
      {
         days = t->tm_hour / 24;
         if (month_days (t) < (t->tm_mday += days))
         {
            if (11 < ++t->tm_mon)
            {
               ++t->tm_year;
               t->tm_mon = 0;
            }
            t->tm_mday = 1;
         }
      }
   }
   if (TM_DEFINED (t->tm_wday))
      t->tm_wday = MOD (t->tm_wday + days, 7);
   t->tm_hour = MOD (t->tm_hour, 24);
   t->tm_min = MOD (t->tm_min, 60);
   t->tm_sec = (int) MOD (sec, 60) + leap_second;
}

int CALC(CF_ZCBondSG1D)(void *Opt,void *Mod,PricingMethod *Met)
{
  TYPEOPT* ptOpt=(TYPEOPT*)Opt;
  TYPEMOD* ptMod=(TYPEMOD*)Mod;

  return zcb_quad1d(    ptMod->flat_flag.Val.V_INT,
                        ptMod->a.Val.V_DOUBLE,
                        ptMod->Sigma.Val.V_PDOUBLE,
                        MOD(GetYield)(ptMod),
                        MOD(GetCurve)(ptMod),
                        ptOpt->BMaturity.Val.V_DATE-ptMod->T.Val.V_DATE,
                        &(Met->Res[0].Val.V_DOUBLE));
}

int CALC(CF_ZCCallBondEuroHW1DG)(void *Opt,void *Mod,PricingMethod *Met)
{
    TYPEOPT* ptOpt=(TYPEOPT*)Opt;
    TYPEMOD* ptMod=(TYPEMOD*)Mod;

    return cf_zbc1d(ptMod->flat_flag.Val.V_INT,
                    ptMod->a.Val.V_DOUBLE,
                    ptMod->CapletCurve.Val.V_ENUM.value,
                    MOD(GetYield)(ptMod),
                    MOD(GetCurve)(ptMod),
                    ptOpt->OMaturity.Val.V_DATE-ptMod->T.Val.V_DATE,
                    ptOpt->BMaturity.Val.V_DATE-ptMod->T.Val.V_DATE,
                    ptOpt->PayOff.Val.V_NUMFUNC_1,
                    &(Met->Res[0].Val.V_DOUBLE));
}

/**
 * @brief Set the vector value using cartesian coordinates
 *
 *    @param v Vector to set.
 *    @param x X value for vector.
 *    @param y Y value for vector.
 */
void vect_cset( Vector2d* v, const double x, const double y )
{
   v->x     = x;
   v->y     = y;
   v->mod   = MOD(x,y);
   v->angle = ANGLE(x,y);
}

int Shash_del(T *hash, char *key)
{
    int index, flag = 0;
    NODE_T **p, *q;
    BUCKET_T *bucket;

    assert(key);
    assert(hash);
    index  = MOD(str_hash(key), hash->len);
    bucket = &hash->buckets[index];

    if (!bucket->list)
        return -1;

    p = &bucket->list;
    do {
        q = (*p)->next;
        if (str_equal((*p)->key, key)) {
            FREE(*p);
            *p   = q;
            flag = 1;
            break;
        }
        p = &(*p)->next;
    } while (q);
    return flag? 0 : -1;
}

/**
 * @brief Adds a polar 2d vector to the current vector.
 *
 *    @param v Vector to add x and y to.
 *    @param m Module of vector to add.
 *    @param a Angle of vector to add.
 */
void vect_padd( Vector2d* v, const double m, const double a )
{
   v->x    += m*cos(a);
   v->y    += m*sin(a);
   v->mod   = MOD(v->x,v->y);
   v->angle = ANGLE(v->x,v->y);
}

/**
 * @brief Adds x and y to the current vector
 *
 *    @param v Vector to add x and y to.
 *    @param x X value to add to vector.
 *    @param y Y value to add to vector.
 */
void vect_cadd( Vector2d* v, const double x, const double y )
{
   v->x    += x;
   v->y    += y;
   v->mod   = MOD(v->x,v->y);
   v->angle = ANGLE(v->x,v->y);
}