int SelectHeaviestSequence( CStudioHdr *pstudiohdr, int activity )
{
	if ( !pstudiohdr )
		return 0;

	VerifySequenceIndex( pstudiohdr );

	int maxweight = 0;
	int seq = ACTIVITY_NOT_AVAILABLE;
	int weight = 0;
	for (int i = 0; i < pstudiohdr->GetNumSeq(); i++)
	{
		int curActivity = GetSequenceActivity( pstudiohdr, i, &weight );
		if (curActivity == activity)
		{
			if ( iabs(weight) > maxweight )
			{
				maxweight = iabs(weight);
				seq = i;
			}
		}
	}

	return seq;
}

void CalculateSumOfSquares(void)
{

unsigned long GxSquare, GySquare, GzSquare;
long temp;

 // Calculate The Total Force
 // GForce = sqrt(x**2 + y**2 + z**2)
 // First  set them all positive. Sum of three 15 bits square won't overflow 32 bits unsigned
  
     
// iabs  ~9us
    GxSquare = (unsigned long) iabs( CurrentData.Gx);
    GySquare = (unsigned long) iabs(CurrentData.Gy);
    GzSquare = (unsigned long) iabs(CurrentData.Gz);



  GxSquare *= GxSquare;
  GySquare *= GySquare;
  GzSquare *= GzSquare;
  CurrentData.SumSquare =  GxSquare + GySquare + GzSquare;

if(CurrentData.SumSquare > PeakData.SumSquare)
   PeakData = CurrentData;
}

void w3j(mpf_t w, long j1, long j2, long j3, long m1, long m2, long m3)
{
  mpq_t delta_sq,r;
  mpz_t i;
  mpf_t h;

  mpq_init(delta_sq);
  mpq_init(r);
  mpz_init(i);
  mpf_init(h);
  mpq_set_si(r,0,1);

  if(m1+m2+m3!=0) return;
  if((iabs(m1)>j1) || (iabs(m2)>j2) || (iabs(m3)>j3)) return;
  if((j3<iabs(j1-j2)) || ((j1+j2)<j3)) return;

  w3j_Delta_sq(delta_sq, j1, j2, j3);
  w3j_intterm(i, j1, j2, j3, m1, m2, m3);
  if(iabs(j1-j2-m3)%2 == 1) mpz_neg(i,i);

  w3j_sqrt_sq(r, j1, j2, j3, m1, m2, m3);

  mpq_mul(r,r,delta_sq);
  mpf_set_q(w,r);
  mpf_sqrt(w,w);
  mpf_set_z(h,i);
  mpf_mul(w,w,h);

  mpf_clear(h);
  mpz_clear(i);
  mpq_clear(r);
  mpq_clear(delta_sq);
}

void ui::line(int x0, int y0, int x1, int y1, unsigned char color)
{
	if(y0==y1)
	{
		if(!correct(x0, y0, x1, y1, clipping))
			return;
		rmset(ptr(x0, y0), color, x1-x0+1);
	}
	else if(x0==x1)
	{
		if(!correct(x0, y0, x1, y1, clipping))
			return;
		linev(ptr(x0, y0), scanline, y1-y0+1, color);
	}
	else
	{
		int dx = iabs(x1-x0), sx = x0<x1 ? 1 : -1;
		int dy = -iabs(y1-y0), sy = y0<y1 ? 1 : -1;
		int err = dx+dy, e2;
		for(;;)
		{
			pixel(x0, y0, color);
			if(x0==x1 && y0==y1)
				break;
			e2 = 2*err;
			if(e2 >= dy) { err += dy; x0 += sx; }
			if(e2 <= dx) { err += dx; y0 += sy; }
		}
	}
}

char FindCrash(void)
{
   // if all absolute acceleration axis are bigger than 1.5 G this is a hit
  // 1.5 G is equivalent to 1.5 *  32768/16 => 3072
     const short OneAndHalfG = 3072;
  if((iabs(CurrentData.Gx) > OneAndHalfG)  || (iabs(CurrentData.Gy) > OneAndHalfG) || (iabs(CurrentData.Gz) > OneAndHalfG))
            return 1;
   return 0;

}

int isqrt(int n){
  int sign = signOf(n);
  n = iabs(n);
  int cur = n, next = (n+1)/2;
  while(iabs(next-cur) > 1){
    cur = next;
    next = (cur + n/cur)/2;
  }
  return sign*next;
}

static void LogPyramid(int PyramidIndex, const TanPile &Pile, FamData *Fam)
	{
	Log("\n");
	Log("Pyramid %d\n", PyramidIndex);

	Log(
"     Label   QStart1     QEnd1   TStart1     TEnd1   QStart1     QEnd1   TStart1     TEnd1    Diag   Start     End\n");
	Log(
"----------  --------  --------  --------  --------  --------  --------  --------  --------  ------  ------  ------\n");
	for (PtrFamData q1 = Fam->begin(); q1 != Fam->end(); ++q1)
		{
		FamMemberData &FamMember1 = *q1;
		int HitIndex1 = FamMember1.PileIndex;
		if (HitIndex1 < 0 || HitIndex1 >= Pile.HitCount)
			Quit("Hit index out of range");
		const HitData &Hit1 = Pile.Hits[HitIndex1];
		int Diag1Start = Hit1.TargetFrom - Hit1.QueryFrom;
		int Diag1End = Hit1.TargetTo - Hit1.QueryTo;
		int Diag1 = (Diag1Start + Diag1End)/2;

		PtrFamData q2 = q1;
		for (++q2; q2 != Fam->end(); ++q2)
			{
			FamMemberData &FamMember2 = *q2;
			int HitIndex2 = FamMember2.PileIndex;
			if (HitIndex2 < 0 || HitIndex2 >= Pile.HitCount)
				Quit("Hit index out of range");
			const HitData &Hit2 = Pile.Hits[HitIndex2];

			int Diag2Start = Hit2.TargetFrom - Hit2.QueryFrom;
			int Diag2End = Hit2.TargetTo - Hit2.QueryTo;
			int Diag2 = (Diag2Start + Diag2End)/2;

			int DiagDist = iabs(Diag1 - Diag2);

			int StartDist = iabs(Hit1.TargetFrom - Hit2.TargetFrom);
			int EndDist = iabs(Hit1.QueryTo - Hit2.QueryTo);
			Log("%10.10s  %8d  %8d  %8d  %8d  %8d  %8d  %8d  %8d  %6d",
				Pile.Label,
				Hit1.QueryFrom,
				Hit1.QueryTo,
				Hit1.TargetFrom,
				Hit1.TargetTo,
				Hit2.QueryFrom,
				Hit2.QueryTo,
				Hit2.TargetFrom,
				Hit2.TargetTo,
				DiagDist);
			Log("  %6d", StartDist);
			Log("  %6d", EndDist);
			Log("\n");
			}
		}
	}

/*************************************************************************
* This function performs a coarse decomposition and determines the 
* min-cover.
* REF: Pothen ACMTrans. on Amth Software
**************************************************************************/
void MinCover_Decompose(idx_t *xadj, idx_t *adjncy, idx_t asize, idx_t bsize, idx_t *mate, idx_t *cover, idx_t *csize) {
    idx_t i, k;
    idx_t *where;
    idx_t card[10];

    where = imalloc(bsize, "MinCover_Decompose: where");
    for (i = 0; i < 10; i++) {
        card[i] = 0;
    }

    for (i = 0; i < asize; i++) {
        where[i] = SC;
    }
    for (; i < bsize; i++) {
        where[i] = SR;
    }

    for (i = 0; i < asize; i++) {
        if (mate[i] == -1) {
            MinCover_ColDFS(xadj, adjncy, i, mate, where, INCOL);
        }
    }
    for (; i < bsize; i++) {
        if (mate[i] == -1) {
            MinCover_RowDFS(xadj, adjncy, i, mate, where, INROW);
        }
    }

    for (i = 0; i < bsize; i++) {
        card[where[i]]++;
    }

    k = 0;
    if (iabs(card[VC] + card[SC] - card[HR]) < iabs(card[VC] - card[SR] - card[HR])) {  /* S = VC+SC+HR */
        /* printf("%"PRIDX" %"PRIDX" ",vc+sc, hr); */
        for (i = 0; i < bsize; i++) {
            if (where[i] == VC || where[i] == SC || where[i] == HR) {
                cover[k++] = i;
            }
        }
    } else {  /* S = VC+SR+HR */
        /* printf("%"PRIDX" %"PRIDX" ",vc, hr+sr); */
        for (i = 0; i < bsize; i++) {
            if (where[i] == VC || where[i] == SR || where[i] == HR) {
                cover[k++] = i;
            }
        }
    }

    *csize = k;
    gk_free((void **) &where, LTERM);

}

char FindDrop(void)
{
   // if all absolute acceleration axis are less than 0.5 G this is a drop
  // 0.5 G is equivalent to 0.5 *  32768/16 => 1024
 const short HalfG = 1024;

  if(iabs(CurrentData.Gx) < HalfG)
    if(iabs(CurrentData.Gy) < HalfG)
      if(iabs(CurrentData.Gz) < HalfG)
            return 1;
return 0;
} 

void DiagramItem::shift(QPoint & p, QPoint other, bool contains_other) const
{
    QRect r = rect();

    p = r.center();

    int x = p.x();
    int y = p.y();
    int ox = other.x();
    int oy = other.y();

    if (contains_other) {
        // move p outside
        int dx = x - ox;
        int dy = y - oy;

        if ((dx == 0) && (dy == 0))
            dx = 16;
        else if ((iabs(dx) < 3) && (iabs(dy) < 3)) {
            dx *= 16;
            dy *= 16;
        }

        do {
            x -= dx;
            y -= dy;
        }
        while (contains(x, y));
    }

    // move p to border
    for (;;) {
        int mx = (x + ox) / 2;
        int my = (y + oy) / 2;

        if (((mx == x) || (mx == ox)) &&
            ((my == y) || (my == oy))) {
            p.setX(mx);
            p.setY(my);
            return;
        }

        if (contains(mx, my) == contains_other) {
            ox = mx;
            oy = my;
        }
        else {
            x = mx;
            y = my;
        }
    }
}

void FrameBuffer::bresenham(int x1, int y1, int x2, int y2, int red, int green, int blue, int line) {
  int Fx[] = { 1,  0, -1,  0};
  int Fy[] = { 0,  1,  0, -1};

  int Gx[] = { 1, -1, -1,  1};
  int Gy[] = { 1,  1, -1, -1};

  int dx = x2 - x1;
  int dy = y2 - y1;

  int X = dx > 0;
  int Y = dy > 0;
  int Z = (iabs(dx) - iabs(dy)) > 0;

  int f = F(X, Y, Z);
  int g = G(X, Y);

  int m1x = Fx[f];
  int m1y = Fy[f];
  int m2x = Gx[g];
  int m2y = Gy[g];

  int da, db;

  if (Z) {
    da = iabs(dx); db = iabs(dy);
  } else {
    da = iabs(dy); db = iabs(dx);
  }

  int D = 2 * db - da;

  plot(x1, y1, red, green, blue);

  int x = x1;
  int y = y1;
  
  int** oosElement = &oosMap[line];
  *oosElement = (int*) malloc((iabs(y2-y1)+2)*sizeof(int));
  (*oosElement)[0] = x1;

  while ((Z && (x != x2)) || (!Z && (y != y2))) {

    D = D + 2*db;

    if (D >= 0) {
      x += m2x;
      y += m2y;
      D = D - 2*da;
    } else {
      x += m1x;
      y += m1y;
    }
  
    (*oosElement)[iabs(y - y1)] = x;
    plot(x, y, red, green, blue); 

  }
}

void w3j_intterm(mpz_t sum, long j1, long j2, long j3, long m1, long m2, long m3)
{
  mpz_t term,h;
  mpz_init(term);
  mpz_init(h);

  long I1 = 0;
  if(j1-j3+m2>I1) I1=j1-j3+m2;
  if(j2-j3-m1>I1) I1=j2-j3-m1;
  long I2 = j1+j2-j3;
  if(j1-m1<I2) I2=j1-m1;
  if(j2+m2<I2) I2=j2+m2;

  mpz_set_ui(sum,0);

  long sgn=1;
  if(iabs(I1)%2==1) sgn=-1;

  for(long k=I1; k<=I2; k++)
  {
//    sum+=sgn*binomialCoefficient(j1+j2-j3,k)*binomialCoefficient(j1-j2+j3,j1-m1-k)
//            *binomialCoefficient(-j1+j2+j3,j2+m2-k);
    binomialCoefficient(term, j1+j2-j3, k);
    binomialCoefficient(h,j1-j2+j3,j1-m1-k);
    mpz_mul(term,term,h);
    binomialCoefficient(h,-j1+j2+j3,j2+m2-k);
    mpz_mul(term,term,h);
    mpz_mul_si(term,term,sgn);
    mpz_add(sum,sum,term);
    sgn*= -1;
  }
  mpz_clear(h);
  mpz_clear(term);
}

/*!
 ************************************************************************
 * \brief
 *    2x2 transform of chroma DC
 * \par Input:
 *    level and run for coefficients
 * \par Output:
 *    length and info
 * \note
 *    see ITU document for bit assignment
 ************************************************************************
 */
void levrun_linfo_c2x2(int level,int run,int *len,int *info)
{
  static const byte NTAB[2][2]=
  {
    {1,5},
    {3,0}
  };
  static const byte LEVRUN[4]=
  {
    2,1,0,0
  };

  if (level == 0) //  check if the coefficient sign EOB (level=0)
  {
    *len=1;
    return;
  }
  else
  {
    int levabs = iabs(level);
    int sign = level <= 0 ? 1 : 0;
    int n = (levabs <= LEVRUN[run]) ? NTAB[levabs - 1][run] + 1 : (levabs - LEVRUN[run]) * 8 + run * 2;
    int nn = n >> 1;
    int i;

    for (i=0; i < 16 && nn != 0; i++)
    {
      nn >>= 1;
    }
    *len  = (i << 1) + 1;
    *info = n - (1 << i) + sign;
  }
}

bool TandemPair(const HitData &Hit1, const HitData &Hit2)
	{
	int Length1 = (Hit1.QueryTo - Hit1.QueryFrom + Hit1.TargetTo- Hit1.TargetFrom)/2;
	int Length2 = (Hit2.QueryTo - Hit2.QueryFrom + Hit2.TargetTo - Hit2.TargetFrom)/2;
	int ShorterLength = min(Length1, Length2);
	int LongerLength = max(Length1, Length2);
	if ((double) ShorterLength / (double) LongerLength < MIN_RATIO)
		return false;

	int StartDist = iabs(Hit1.TargetFrom - Hit2.TargetFrom);
	int EndDist = iabs(Hit1.QueryTo - Hit2.QueryTo);

	double StartMargin = (double) StartDist / (double) ShorterLength;
	double EndMargin = (double) EndDist / (double) ShorterLength;
	return StartMargin <= MAX_FRACT_MARGIN && EndMargin <= MAX_FRACT_MARGIN;
	}

int SelectWeightedSequence( CStudioHdr *pstudiohdr, int activity, int curSequence )
{
	VPROF( "SelectWeightedSequence" );

	if (! pstudiohdr)
		return 0;

	if (!pstudiohdr->SequencesAvailable())
		return 0;

	VerifySequenceIndex( pstudiohdr );

#if STUDIO_SEQUENCE_ACTIVITY_LOOKUPS_ARE_SLOW
	int weighttotal = 0;
	int seq = ACTIVITY_NOT_AVAILABLE;
	int weight = 0;
	for (int i = 0; i < pstudiohdr->GetNumSeq(); i++)
	{
		int curActivity = GetSequenceActivity( pstudiohdr, i, &weight );
		if (curActivity == activity)
		{
			if ( curSequence == i && weight < 0 )
			{
				seq = i;
				break;
			}
			weighttotal += iabs(weight);
			
			int randomValue;

			if ( IsInPrediction() )
				randomValue = SharedRandomInt( "SelectWeightedSequence", 0, weighttotal - 1, i );
			else
				randomValue = RandomInt( 0, weighttotal - 1 );
			
			if (!weighttotal || randomValue < iabs(weight))
				seq = i;
		}
	}

	return seq;
#else
	return pstudiohdr->SelectWeightedSequence( activity, curSequence );
#endif
}

/*!
 ************************************************************************
 * \brief
 *    Quantization process for All coefficients for a 4x4 block
 *
 * \par Input:
 *
 * \par Output:
 *
 ************************************************************************
 */
int quant_4x4_normal(int (*tblock)[16], int block_y, int block_x, int  qp,
                     int*  ACLevel, int*  ACRun, 
                     int **fadjust4x4, int **levelscale, int **invlevelscale, int **leveloffset,
                     int *coeff_cost, const byte (*pos_scan)[2], const byte *c_cost, int is_cavlc)
{
  static int i,j, coeff_ctr;

  static int *m7;
  static int scaled_coeff;

  int   level, run = 0;
  int   nonzero = FALSE;
  int   qp_per = qp_per_matrix[qp];
  int   q_bits = Q_BITS + qp_per;
  const byte *p_scan = &pos_scan[0][0];
  int*  ACL = &ACLevel[0];
  int*  ACR = &ACRun[0];

  // Quantization
  for (coeff_ctr = 0; coeff_ctr < 16; coeff_ctr++)
  {
    i = *p_scan++;  // horizontal position
    j = *p_scan++;  // vertical position

    m7 = &tblock[j][block_x + i];

    if (*m7 != 0)
    {
      scaled_coeff = iabs (*m7) * levelscale[j][i];
      level = (scaled_coeff + leveloffset[j][i]) >> q_bits;

      if (level != 0)
      {
        if (is_cavlc)
          level = imin(level, CAVLC_LEVEL_LIMIT);

        *coeff_cost += (level > 1) ? MAX_VALUE : c_cost[run];

        level  = isignab(level, *m7);
        *m7    = rshift_rnd_sf(((level * invlevelscale[j][i]) << qp_per), 4);
        // inverse scale can be alternative performed as follows to ensure 16bit
        // arithmetic is satisfied.
        // *m7 = (qp_per<4) ? rshift_rnd_sf((level*invlevelscale[j][i]),4-qp_per) : (level*invlevelscale[j][i])<<(qp_per-4);
        *ACL++ = level;
        *ACR++ = run; 
        // reset zero level counter
        run    = 0;
        nonzero = TRUE;        
      }
      else
      {
        run++;
        *m7 = 0;
      }
    }
    else
    {

long 
labs(long j)
{
#if	defined(__TCS__)
	return iabs((int)j);		/* N.B. assumes int==long */
#else	/* defined(__TCS__) */
	return (j < 0) ? -j : j;
#endif	/* defined(__TCS__) */
}

/*!
 ************************************************************************
 * \brief
 *    Quantization process for All coefficients for a 2x2 DC block
 *
 * \par Input:
 *
 * \par Output:
 *
 ************************************************************************
 */
int quant_dc2x2_around(Macroblock *currMB, int **tblock, int qp, int* DCLevel, int* DCRun, 
                       LevelQuantParams *q_params_4x4, int **fadjust, const byte (*pos_scan)[2])
{
  QuantParameters *p_Quant = currMB->p_Vid->p_Quant;
  Boolean is_cavlc = (Boolean) (currMB->p_Slice->symbol_mode == CAVLC);
  int coeff_ctr;

  int *m7;
  int scaled_coeff;

  int   level, run = 0;
  int   nonzero = FALSE;  
  int   qp_per = p_Quant->qp_per_matrix[qp];
  int   q_bits = Q_BITS + qp_per + 1;
  //const byte *p_scan = &pos_scan[0][0];
  int*  DCL = &DCLevel[0];
  int*  DCR = &DCRun[0];

  m7 = *tblock;

  // Quantization
  for (coeff_ctr=0; coeff_ctr < 4; ++coeff_ctr)
  {
    // we need to update q_params_4x4->OffsetComp to a 4x1 array that would contain offset info for 
    // every 2x2 DC position
    if (*m7)
    {
      scaled_coeff = iabs (*m7) * q_params_4x4->ScaleComp;
      level = (scaled_coeff + (q_params_4x4->OffsetComp << 1) ) >> q_bits;

      if (level  != 0)
      {
        if (is_cavlc)
          level = imin(level, CAVLC_LEVEL_LIMIT);

        level = isignab(level, *m7);

        *m7++ = ((level * q_params_4x4->InvScaleComp) << qp_per);

        *DCL++  = level;
        *DCR++  = run;
        // reset zero level counter
        run     = 0;
        nonzero = TRUE;
      }
      else
      {
        ++run;
        *m7++ = 0;
      }
    }
    else
    {
      ++run;
      ++m7;
    }
  }

/*!
***********************************************************************
* \brief
*    Calculate SAD
***********************************************************************
*/
distblk distortion4x4SAD(short* diff, distblk min_dist)
{
	int distortion = 0, k;
	for (k = 0; k < 16; k++)
	{
		distortion += iabs(*diff++);
	}
	return (dist_scale((distblk) distortion));
}

int findtorq(int btyp1, int btyp2) {
    int i;
    for (i=0; i < BENDATASIZE; i++) {
		
	if  (iabs(bendata[i].b1typ) == iabs(btyp1) &&
	     iabs(bendata[i].b2typ) == iabs(btyp2)) return i;
	if  (iabs(bendata[i].b1typ) == iabs(btyp2) &&
	     iabs(bendata[i].b2typ) == iabs(btyp1)) return i;
    }
    // fprintf(stderr, "Bend type %d-%d not found\n",btyp1,btyp2);
    return 0; // the default bend type
}