int makeMSP_v3(MSP* msp, attr_tree_node* node){
   int** matrix = msp -> matrix;
   char** label = msp -> label;
   int cols = msp -> cols;
   int rows = msp -> rows;

   attr_tree_node* queue[100] = {0};
   int* vectors [100] = {0};
   int write_i = -1;
   int read_i = -1;
   int* v = calloc(cols, sizeof(int));
   v[0] = 1;

   queue[++write_i] = node;
   vectors[write_i] = v;

   int count = 1;
   int index = 0;

   while( read_i < write_i ){

      attr_tree_node* node = queue[++read_i];
      int* v = vectors[read_i];

      if( node -> node_type == ATTR_NODE_OR ){
	 int i = 0;
         for( i; i < 2; i++){
	    queue[++write_i] = (node -> subnode)[i];
	    vectors[write_i] = v;
	 }
      }
      else if( node -> node_type == ATTR_NODE_AND ){
         int* l_v = calloc(cols, sizeof(int));
         int* r_v = calloc(cols, sizeof(int));
	 memcpy(l_v, v, cols*sizeof(int));
	 l_v[count] = 1;
	 {
	    int i = 0;
	    for(i; i < cols; i++)
	       r_v[i] = v[i] - l_v[i];
	 }
	 count++;

	 queue[++write_i] = (node -> subnode)[0];
	 vectors[write_i] = l_v;

	 queue[++write_i] = (node -> subnode)[1];
	 vectors[write_i] = r_v;
      }
      else if( node -> node_type == ATTR_NODE_THRESHOLD ){
	 int i = 0;
	 for( i; i < node->num_subnodes; i++ ){
            int* l_v = calloc(cols, sizeof(int));
	    memcpy(l_v, v, cols*sizeof(int));
            int j = 0;
	    for( j; j < node->threshold_k - 1; j++){
	       l_v[count+j] = power(i+1, j+1);
	       //printf("power:%d\n", power(i+1, j+1));
	    }

	    queue[++write_i] = (node -> subnode)[i];
	    vectors[write_i] = l_v;
            //free(l_v);
	 }
	 count = count + (node->threshold_k) - 1;

      }
      else if( node -> node_type == ATTR_NODE_LEAF){
         memcpy( matrix[index], v, cols*sizeof(int) );
         label[index] = node->attribute;
         index++;
      }

   }
}

inline
time_duration
str_from_delimited_time_duration(const std::basic_string<char_type>& s)
{
    unsigned short min=0, sec =0;
    int hour =0;
    bool is_neg = (s.at(0) == '-');
    lslboost::int64_t fs=0;
    int pos = 0;

    typedef typename std::basic_string<char_type>::traits_type traits_type;
    typedef lslboost::char_separator<char_type, traits_type> char_separator_type;
    typedef lslboost::tokenizer<char_separator_type,
            typename std::basic_string<char_type>::const_iterator,
            std::basic_string<char_type> > tokenizer;
    typedef typename lslboost::tokenizer<char_separator_type,
            typename std::basic_string<char_type>::const_iterator,
            typename std::basic_string<char_type> >::iterator tokenizer_iterator;

    char_type sep_chars[5] = {'-',':',',','.'};
    char_separator_type sep(sep_chars);
    tokenizer tok(s,sep);
    for(tokenizer_iterator beg=tok.begin(); beg!=tok.end(); ++beg) {
        switch(pos) {
        case 0: {
            hour = lslboost::lexical_cast<int>(*beg);
            break;
        }
        case 1: {
            min = lslboost::lexical_cast<unsigned short>(*beg);
            break;
        }
        case 2: {
            sec = lslboost::lexical_cast<unsigned short>(*beg);
            break;
        };
        case 3: {
            int digits = static_cast<int>(beg->length());
            //Works around a bug in MSVC 6 library that does not support
            //operator>> thus meaning lexical_cast will fail to compile.
#if (defined(BOOST_MSVC) && (_MSC_VER < 1300))
            // msvc wouldn't compile 'time_duration::num_fractional_digits()'
            // (required template argument list) as a workaround a temp
            // time_duration object was used
            time_duration td(hour,min,sec,fs);
            int precision = td.num_fractional_digits();
            // _atoi64 is an MS specific function
            if(digits >= precision) {
                // drop excess digits
                fs = _atoi64(beg->substr(0, precision).c_str());
            }
            else {
                fs = _atoi64(beg->c_str());
            }
#else
            int precision = time_duration::num_fractional_digits();
            if(digits >= precision) {
                // drop excess digits
                fs = lslboost::lexical_cast<lslboost::int64_t>(beg->substr(0, precision));
            }
            else {
                fs = lslboost::lexical_cast<lslboost::int64_t>(*beg);
            }
#endif
            if(digits < precision) {
                // trailing zeros get dropped from the string,
                // "1:01:01.1" would yield .000001 instead of .100000
                // the power() compensates for the missing decimal places
                fs *= power(10, precision - digits);
            }

            break;
        }
        default:
            break;
        }//switch
        pos++;
    }
    if(is_neg) {
        return -time_duration(hour, min, sec, fs);
    }
    else {
        return time_duration(hour, min, sec, fs);
    }
}

mp_integer bv_arithmetict::pack() const
{
  if(value>=0) return value;
  return value+power(2, spec.width);
}

void main()
{
	int x=5,y=2;
	printf("%d\n",power(x,y));
	printf("%d\n",poweropti(5,4));
}

// returns a pointer to the loaded bitmap
// only supports 24-bit (and 32-bit?) images
wBitmap* wLoadBitmap(char *filename)
{
	FILE *bmpFile;
	wBitmap *bmp;
	int i;
	int fileLineSize, memLineSize, fileScanLines, memScanLines, dataSize;
	int pad, generic;
	BYTE *line;

	bmp = (wBitmap*)malloc(sizeof(wBitmap));
	memset(bmp, 0, sizeof(wBitmap));
	bmp->name = filename;

	bmpFile = fopen(filename, "rb");
	if (bmpFile)
	{
		printf("Loading bitmap %s\n", filename);

		fread(&bmp->Header, sizeof(bmp->Header), 1, bmpFile);
		fread(&bmp->Info, sizeof(bmp->Info), 1, bmpFile);

		// only supports 24-bit+ color bitmaps
		if (bmp->Info.biBitCount < 24)
		{
			printf("cannot load %s:  bitmap is less than 24-bit color\n", filename);
			free(bmp);
			return NULL;
		}
		if (bmp->Info.biCompression)
		{
			printf("cannot load %s:  bitmap is compressed\n", filename);
			free(bmp);
			return NULL;
		}

		fileLineSize = bmp->Info.biWidth*bmp->Info.biBitCount / 8;
		pad = 4 - fileLineSize % 4;
		if (pad == 4)
			pad = 0;
		fileScanLines = bmp->Info.biHeight;

		// check to make sure that image is 2^x * 2^y
		// image in memory can be 2^x * 2^y even if disk file isn't, but there will be blank pixels
		for (i = 1; i <= MAX_EXP; i++)
		{
			generic = (int)power(2, i);
			if (generic >= bmp->Info.biWidth)
			{
				bmp->Info.biWidth = generic;
				break;
			}
		}
		if (i > MAX_EXP)
		{
			printf("cannot load %s:  bitmap is too large\n", filename);
			free(bmp);
			return NULL;
		}
		for (i = 1; i <= MAX_EXP; i++)
		{
			generic = power(2, i);
			if (generic >= bmp->Info.biHeight)
			{
				bmp->Info.biHeight = generic;
				break;
			}
		}
		if (i > MAX_EXP)
		{
			printf("cannot load %s:  bitmap is too large\n", filename);
			free(bmp);
			return NULL;
		}

		memLineSize = bmp->Info.biWidth*bmp->Info.biBitCount / 8;
		memScanLines = bmp->Info.biHeight;
		dataSize = memLineSize*memScanLines;
		bmp->pixels = (BYTE*)malloc(dataSize);
		memset(bmp->pixels, 0, dataSize);

		// end 2^n check
		/*		printf("image is %i by %i", fileLineSize*8/bmp->Info.biBitCount, fileScanLines);
		if ( fileLineSize != memLineSize || fileScanLines != memScanLines )
		printf(", expanded to %i by %i", memLineSize*8/bmp->Info.biBitCount, memScanLines);
		printf("\n");
		*/

		// bmps are stored last line first
		fseek(bmpFile, bmp->Header.bfOffBits, 0);
		line = bmp->pixels + (fileScanLines - 1)*memLineSize;
		for (i = 0; i < fileScanLines; i++)
		{
			fread(line, fileLineSize, 1, bmpFile);
			// lines are padded to be 32-bit divisible
			if (pad)
				fread(&generic, pad, 1, bmpFile);
			line -= memLineSize;
		}
//.........这里部分代码省略.........

BaseIF* makeGeometry(Box&      a_domain,
                     RealVect& a_origin,
                     Real&     a_dx)
{
    RealVect center;
    RealVect radii;
    bool     insideRegular;

    // parse input file
    ParmParse pp;

    Vector<int> n_cell(SpaceDim);
    pp.getarr("n_cell",n_cell,0,SpaceDim);

    CH_assert(n_cell.size() == SpaceDim);

    IntVect lo = IntVect::Zero;
    IntVect hi;

    for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
        if (n_cell[ivec] <= 0)
        {
            pout() << "Bogus number of cells input = " << n_cell[ivec];
            exit(1);
        }

        hi[ivec] = n_cell[ivec] - 1;
    }

    a_domain.setSmall(lo);
    a_domain.setBig(hi);

    Vector<Real> prob_lo(SpaceDim,1.0);
    Real prob_hi;

    pp.getarr("prob_lo",prob_lo,0,SpaceDim);
    pp.get("prob_hi",prob_hi);

    a_dx = (prob_hi-prob_lo[0])/n_cell[0];

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        a_origin[idir] = prob_lo[idir];
    }

    // ParmParse doesn't get RealVects, so work-around with Vector<Real>
    Vector<Real> vectorCenter;
    pp.getarr("center",vectorCenter,0,SpaceDim);

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        center[idir] = vectorCenter[idir];
    }

    Vector<Real> vectorRadii;
    pp.getarr("radii",vectorRadii,0,SpaceDim);

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        radii[idir] = vectorRadii[idir];
    }

    int halfTurns;
    pp.get("halfturns",halfTurns);

    // Parm Parse doesn't get bools, so work-around with int
    int intInsideRegular;
    pp.get("insideRegular",intInsideRegular);

    if (intInsideRegular != 0) insideRegular = true;
    if (intInsideRegular == 0) insideRegular = false;

    IntVect zero=IntVect::Zero;
    bool inside = true;
    EllipsoidIF ellipsoid(radii,center,inside);

    Real coef = halfTurns / 2.0;
    IntVect power(D_DECL(1,0,0));

    Vector<PolyTerm> term;
    term.resize(1);
    term[0].coef   = coef;
    term[0].powers = power;

    PolynomialIF rotate(term,inside);

    LatheIF implicit(ellipsoid,rotate,center,insideRegular);

    RealVect vectDx = RealVect::Unit;
    vectDx *= a_dx;

    GeometryShop workshop(implicit,0,vectDx);

    // This generates the new EBIS
    EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
    ebisPtr->define(a_domain,a_origin,a_dx,workshop);

    return implicit.newImplicitFunction();
}

static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
{
	const ex x_pt = x.subs(rel, subs_options::no_pattern);
	if (x_pt.info(info_flags::numeric)) {
		// First special case: x==0 (derivatives have poles)
		if (x_pt.is_zero()) {
			// method:
			// The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot 
			// simply substitute x==0.  The limit, however, exists: it is 1.
			// We also know all higher derivatives' limits:
			// (d/dx)^n Li2(x) == n!/n^2.
			// So the primitive series expansion is
			// Li2(x==0) == x + x^2/4 + x^3/9 + ...
			// and so on.
			// We first construct such a primitive series expansion manually in
			// a dummy symbol s and then insert the argument's series expansion
			// for s.  Reexpanding the resulting series returns the desired
			// result.
			const symbol s;
			ex ser;
			// manually construct the primitive expansion
			for (int i=1; i<order; ++i)
				ser += pow(s,i) / pow(numeric(i), *_num2_p);
			// substitute the argument's series expansion
			ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
			// maybe that was terminating, so add a proper order term
			epvector nseq { expair(Order(_ex1), order) };
			ser += pseries(rel, std::move(nseq));
			// reexpanding it will collapse the series again
			return ser.series(rel, order);
			// NB: Of course, this still does not allow us to compute anything
			// like sin(Li2(x)).series(x==0,2), since then this code here is
			// not reached and the derivative of sin(Li2(x)) doesn't allow the
			// substitution x==0.  Probably limits *are* needed for the general
			// cases.  In case L'Hospital's rule is implemented for limits and
			// basic::series() takes care of this, this whole block is probably
			// obsolete!
		}
		// second special case: x==1 (branch point)
		if (x_pt.is_equal(_ex1)) {
			// method:
			// construct series manually in a dummy symbol s
			const symbol s;
			ex ser = zeta(_ex2);
			// manually construct the primitive expansion
			for (int i=1; i<order; ++i)
				ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i));
			// substitute the argument's series expansion
			ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
			// maybe that was terminating, so add a proper order term
			epvector nseq { expair(Order(_ex1), order) };
			ser += pseries(rel, std::move(nseq));
			// reexpanding it will collapse the series again
			return ser.series(rel, order);
		}
		// third special case: x real, >=1 (branch cut)
		if (!(options & series_options::suppress_branchcut) &&
			ex_to<numeric>(x_pt).is_real() && ex_to<numeric>(x_pt)>1) {
			// method:
			// This is the branch cut: assemble the primitive series manually
			// and then add the corresponding complex step function.
			const symbol &s = ex_to<symbol>(rel.lhs());
			const ex point = rel.rhs();
			const symbol foo;
			epvector seq;
			// zeroth order term:
			seq.push_back(expair(Li2(x_pt), _ex0));
			// compute the intermediate terms:
			ex replarg = series(Li2(x), s==foo, order);
			for (size_t i=1; i<replarg.nops()-1; ++i)
				seq.push_back(expair((replarg.op(i)/power(s-foo,i)).series(foo==point,1,options).op(0).subs(foo==s, subs_options::no_pattern),i));
			// append an order term:
			seq.push_back(expair(Order(_ex1), replarg.nops()-1));
			return pseries(rel, std::move(seq));
		}
	}
	// all other cases should be safe, by now:
	throw do_taylor();  // caught by function::series()
}

unsigned long long int power(unsigned long long int a,unsigned long long int b) { if(b==0) return 1; unsigned long long int x=power(a,b/2);
return (((x*x)%mod)*((b&1)?a:1))%mod;}

//TODO remove pModel
CEvaluationNode* CDerive::deriveBranch(const CEvaluationNode* node, const CCopasiObject * pObject,
                                       std::vector<const CEvaluationNode*>& env,
                                       //std::vector<const CCopasiObject*>& objenv,
                                       const CEvaluationTree* pTree,
                                       bool simplify)
{
  CEvaluationNode * newNode = NULL;

  const CEvaluationNodeOperator * pENO = dynamic_cast<const CEvaluationNodeOperator*>(node);

  if (pENO)
    {
      if (!pENO->getLeft() || !pENO->getRight()) return NULL;

      CEvaluationNode * pLeftDeriv = deriveBranch(pENO->getLeft(), pObject, env, pTree, simplify);

      if (!pLeftDeriv) return NULL;

      CEvaluationNode * pRightDeriv = deriveBranch(pENO->getRight(), pObject, env, pTree, simplify);

      if (!pRightDeriv) {delete pLeftDeriv; return NULL;}

      // we now know that derivations of the left and right branch exist

      switch ((CEvaluationNodeOperator::SubType) CEvaluationNode::subType(pENO->getType()))
        {
          case CEvaluationNodeOperator::MULTIPLY:
          {
            CEvaluationNode * pLeftCopy = copyBranch_var2obj(pENO->getLeft(), env);
            CEvaluationNode * pRightCopy = copyBranch_var2obj(pENO->getRight(), env);

            CEvaluationNode * tmpNode1 = multiply(pRightCopy, pLeftDeriv, simplify);
            CEvaluationNode * tmpNode2 = multiply(pRightDeriv, pLeftCopy, simplify);

            return add(tmpNode1, tmpNode2, simplify);
          }
          break;

          case CEvaluationNodeOperator::DIVIDE:
          {
            CEvaluationNode * pLeftCopy = copyBranch_var2obj(pENO->getLeft(), env);
            CEvaluationNode * pRightCopy = copyBranch_var2obj(pENO->getRight(), env);

            //numerator
            CEvaluationNode * tmpNode1 = multiply(pRightCopy, pLeftDeriv, simplify);
            CEvaluationNode * tmpNode2 = multiply(pRightDeriv, pLeftCopy, simplify);

            CEvaluationNode * minusNode = subtract(tmpNode1, tmpNode2, simplify);

            minusNode->compile(NULL);

            //denominator
            CEvaluationNode * powerNode = power(copyBranch_var2obj(pENO->getRight(), env),
                                                new CEvaluationNodeNumber(CEvaluationNodeNumber::INTEGER, "2"),
                                                simplify);

            return divide(minusNode, powerNode, simplify);
          }
          break;

          case CEvaluationNodeOperator::PLUS:

            return add(pLeftDeriv, pRightDeriv, simplify);
            break;

          case CEvaluationNodeOperator::MINUS:

            return subtract(pLeftDeriv, pRightDeriv, simplify);
            break;

          case CEvaluationNodeOperator::POWER:
          {
            // b-1
            CEvaluationNode * tmpNode = subtract(copyBranch_var2obj(pENO->getRight(), env),
                                                 new CEvaluationNodeNumber(CEvaluationNodeNumber::INTEGER, "1"),
                                                 simplify);

            // a^(b-1)
            CEvaluationNode * powerNode = power(copyBranch_var2obj(pENO->getLeft(), env), tmpNode, simplify);

            // b*a'
            tmpNode = multiply(copyBranch_var2obj(pENO->getRight(), env),
                               pLeftDeriv, simplify);

            // ln a
            CEvaluationNodeFunction * funcNode = new CEvaluationNodeFunction(CEvaluationNodeFunction::LOG, "ln");
            funcNode->addChild(copyBranch_var2obj(pENO->getLeft(), env)); // add a

            // a * b' * ln a
            CEvaluationNode * tmpNode2 = multiply(copyBranch_var2obj(pENO->getLeft(), env),
                                                  multiply(pRightDeriv, funcNode, simplify),
                                                  simplify);

            // b*a + a*b * ln a
            CEvaluationNode * plusNode = add(tmpNode, tmpNode2, simplify);

            // a^(b-1)*(b*a + a*b * ln a)
            return multiply(powerNode, plusNode, simplify);
          }
//.........这里部分代码省略.........

int load_bmp(char *file, struct image **res_image)
{
    FILE *fp;
    struct bitmap_file_header hdr;
    struct bitmap_info_header infohdr;
    struct image *img;
    int num_pixels;
    int j, pad;


    if (!res_image) {
        printf("**res_image pointer to pointer is not initialized!!!\n");
        errno = -EAGAIN;
        return -EAGAIN;
    }
    
    img = (struct image *) malloc(sizeof(*img));
    if (!img)
        return errno;

    *res_image = img;

    /* open the file */
    if ((fp = fopen(file,"r")) == NULL)
        return errno;

    /* Read bitmap header */
    if (fread(&hdr, sizeof(hdr), 1, fp) != 1)
        goto error;

    /* Check format */
    if (hdr.id[0] != 'B' || hdr.id[1] != 'M') {
        printf("%s is not a bitmap file in windows-format.\n", file);
        errno = -EAGAIN;
        goto error;
    }

    /* Read bitmap info header */
    if (fread(&infohdr, sizeof(infohdr), 1, fp) != 1)
        goto error;

    if (infohdr.color_planes != 1 || infohdr.bits_per_pixel != 24 ||
        infohdr.compression != 0 || infohdr.num_colors != 0 ||
        infohdr.important_colors != 0) {
        printf("The bmp image is not in a supported format!\n");
        printf("The supported format requires: colorplanes == 0, 24 bits per "
               "pixel, no compression, num-colors == 0 and important-colors == 0\n");
        printf("But we got: colorplanes %u bits per pixel %u compression %u "
               "num-colors %u important-colors %u\n", infohdr.color_planes,
               infohdr.bits_per_pixel, infohdr.compression, infohdr.num_colors,
               infohdr.important_colors);
        errno = -EAGAIN;
        goto error;
    }

    if (infohdr.num_colors == 0)
        /* This means, the number of colors in the color-pallette is 2**bits_per_pixel */
        infohdr.num_colors = power(2, infohdr.bits_per_pixel);

    img->width = infohdr.width;
    img->height = infohdr.height;
    img->hor_res = infohdr.hor_res;
    img->ver_res = infohdr.ver_res;

    /* Now, move the pointer to the pixel-array */
    if (fseek(fp, infohdr.hdrsize - sizeof(infohdr), SEEK_CUR))
        goto error;

    num_pixels = img->width * img->height;
    img->pixels = (struct pixel *) malloc(sizeof(struct pixel)*num_pixels);

    if (!img->pixels)
        goto error;

    pad = (4 - (img->width * 3) % 4) % 4;
    for (j = 0; j < img->height; j++) {
        if (fread(&img->pixels[j*img->width], 3, img->width, fp) != img->width)
            goto error;

        if (fseek(fp, pad, SEEK_CUR))
            goto error;
    }

    fclose(fp);
    return 0;

error:
    fclose(fp);
    return errno;
}

 void power_mpz(elem &result, elem a, mpz_ptr n) const
 {
   int n1 = static_cast<int>(mpz_fdiv_ui(n, p1));
   power(result,a,n1);
 }

void Tester_68k::testSbcd() {
    u16 opcode;

    //MOVE.L    #$12345689, D1
    //MOVE.L    #$f745ff78, D2
    //MOVE #$4, CCR
    //SBCD D2, D1
    setUp();
    opcode = 0x8100 | (1 << 9) | 2;
    addWord(opcode);
    addWord(0xffff);
    power();
    setCCR(0x4);
    setRegD(1, 0x12345689);
    setRegD(2, 0xf745ff78);
    process();
    check("SBCD D2, D1");

    //MOVE.L    #$12345605, D1
    //MOVE.L    #$f745ff05, D2
    //MOVE #$4, CCR
    //SBCD D2, D1
    setUp();
    opcode = 0x8100 | (1 << 9) | 2;
    addWord(opcode);
    addWord(0xffff);
    power();
    setCCR(0x4);
    setRegD(1, 0x123456ff);
    setRegD(2, 0xf745ffff);
    process();
    check("SBCD D2, D1 zero");

    //MOVE.L    #$12345634, D1
    //MOVE.L    #$f745ff45, D2
    //MOVE #$10, CCR
    //SBCD D2, D1
    setUp();
    opcode = 0x8100 | (1 << 9) | 2;
    addWord(opcode);
    addWord(0xffff);
    power();
    setCCR(0x10);
    setRegD(1, 0x12345634);
    setRegD(2, 0xf745ff45);
    process();
    check("SBCD D2, D1 neg");

    //MOVE.L    #$12345634, D1
    //MOVE.L    #$f745ff45, D2
    //MOVE #$10, CCR
    //SBCD D2, D1
    setUp();
    opcode = 0x8100 | (1 << 9) | 2;
    addWord(opcode);
    addWord(0xffff);
    power();
    setCCR(0x10);
    setRegD(1, 0x123456a9);
    setRegD(2, 0xf745ffff);
    process();
    check("SBCD D2, D1 overflow");

    //MOVE.L    #$3001, A1
    //MOVE.L    #$4001, A2
    //MOVE.B    #$a2, $3000
    //MOVE.B    #$19, $4000
    //MOVE #$10, CCR
    //SBCD -(A2), -(A1)
    //MOVE.B    $3000, D1
    //MOVE.B    $4000, D2
    setUp();
    opcode = 0x8100 | (1 << 9) | (1 << 3) | 2;
    addWord(opcode);
    addWord(0xffff);
    power();
    memoryblock.write(0x3000, 0xa2);
    memoryblock.write(0x4000, 0x19);
    setCCR(0x10);
    setRegA(1, 0x3001);
    setRegA(2, 0x4001);
    process();
    check("SBCD -(A2), -(A1)");
}

ui get_head_num(ull n, ui d)
{
    return n / power(10, d - 2);
}

void
absval(void)
{
	int h;
	save();
	p1 = pop();

	if (istensor(p1)) {
		absval_tensor();
		restore();
		return;
	}

	if (isnum(p1)) {
		push(p1);
		if (isnegativenumber(p1))
			negate();
		restore();
		return;
	}

	if (iscomplexnumber(p1)) {
		push(p1);
		push(p1);
		conjugate();
		multiply();
		push_rational(1, 2);
		power();
		restore();
		return;
	}

	// abs(1/a) evaluates to 1/abs(a)

	if (car(p1) == symbol(POWER) && isnegativeterm(caddr(p1))) {
		push(p1);
		reciprocate();
		absval();
		reciprocate();
		restore();
		return;
	}

	// abs(a*b) evaluates to abs(a)*abs(b)

	if (car(p1) == symbol(MULTIPLY)) {
		h = tos;
		p1 = cdr(p1);
		while (iscons(p1)) {
			push(car(p1));
			absval();
			p1 = cdr(p1);
		}
		multiply_all(tos - h);
		restore();
		return;
	}

	if (isnegativeterm(p1) || (car(p1) == symbol(ADD) && isnegativeterm(cadr(p1)))) {
		push(p1);
		negate();
		p1 = pop();
	}

	push_symbol(ABS);
	push(p1);
	list(2);

	restore();
}

double calculate(int numInputTokens, char **inputString){
	int i;
	double result = 0.0;
	char *s;
	struct DynArr *stack;
	double num;
	//set up the stack
	stack = createDynArr(20);

	// start at 1 to skip the name of the calculator calc
	for(i=1;i < numInputTokens;i++)
	{
		s = inputString[i];

		// Hint: General algorithm:
		// (1) Check if the string s is in the list of operators.
		//   (1a) If it is, perform corresponding operations.
		//   (1b) Otherwise, check if s is a number.
		//     (1b - I) If s is not a number, produce an error.
		//     (1b - II) If s is a number, push it onto the stack

		if(strcmp(s, "+") == 0){
			add(stack);
			printf("Adding\n");
		}
		else if(strcmp(s,"-") == 0){
			subtract(stack);
			printf("Subtracting\n");
		}
		else if(strcmp(s, "/") == 0){
			divide(stack);
			printf("Dividing\n");
		}
		else if(strcmp(s, "x") == 0){
			multiply(stack);
			printf("Multiplying\n");
		}
        else if(strcmp(s,"^") == 0){
			power(stack);
			printf("Power\n");
        }
		else if(strcmp(s, "^2") == 0){
			squared(stack);
			printf("Squaring\n");
		}
		else if(strcmp(s, "^3") == 0){
			cubed(stack);
			printf("Cubing\n");
		}
		else if(strcmp(s, "abs") == 0){
			absoluteValue(stack);
			printf("Absolute value\n");
		}
		else if(strcmp(s, "sqrt") == 0){
			squareRoot(stack);
			printf("Square root\n");
		}
		else if(strcmp(s, "exp") == 0){
			exponential(stack);
			printf("Exponential\n");
		}
		else if(strcmp(s, "ln") == 0){
			naturalLog(stack);
			printf("Natural Log\n");
		}
		else if(strcmp(s, "log") == 0){
			logBase10(stack);
			printf("Log\n");
		}

        else{
    // FIXME: You need to develop the code here (when s is not an operator)
    // Remember to deal with special values ("pi" and "e")
    //check if not a number
            if (isNumber(s, &num) == 0){
                if (strcmp(s, "pi") == 0){
                    num = 3.14159265;
                }
                else if (strcmp(s, "e") == 0){
                    num = 2.7182818;
                }
                else{	//wrong
                    printf("%s is not valid (number or operator) \n", s);
                    break;
                }
            }
            pushDynArr(stack, num);
        }
    }	//end for
/* FIXME: You will write this part of the function (2 steps below)
* (1) Check if everything looks OK and produce an error if needed.
* (2) Store the final value in result and print it out.
*/
    if (sizeDynArr(stack) != 1) {
        printf("Incorrect count of numbers is detected! Calculations CANNOT be preformed. ");
        return 0;
    }
    else {
        result = topDynArr(stack);
    }
//.........这里部分代码省略.........

void
yybesselj(void)
{
	double d;
	int n;

	N = pop();
	X = pop();

	push(N);
	n = pop_integer();

	// numerical result

	if (isdouble(X) && n != (int) 0x80000000) {
		//d = jn(n, X->u.d);
		push_double(d);
		return;
	}

	// bessej(0,0) = 1

	if (iszero(X) && iszero(N)) {
		push_integer(1);
		return;
	}

	// besselj(0,n) = 0

	if (iszero(X) && n != (int) 0x80000000) {
		push_integer(0);
		return;
	}

	// half arguments

	if (N->k == NUM && MEQUAL(N->u.q.b, 2)) {

		// n = 1/2

		if (MEQUAL(N->u.q.a, 1)) {
			push_integer(2);
			push_symbol(PI);
			divide();
			push(X);
			divide();
			push_rational(1, 2);
			power();
			push(X);
			sine();
			multiply();
			return;
		}

		// n = -1/2

		if (MEQUAL(N->u.q.a, -1)) {
			push_integer(2);
			push_symbol(PI);
			divide();
			push(X);
			divide();
			push_rational(1, 2);
			power();
			push(X);
			cosine();
			multiply();
			return;
		}

		// besselj(x,n) = (2/x) (n-sgn(n)) besselj(x,n-sgn(n)) - besselj(x,n-2*sgn(n))

		push_integer(MSIGN(N->u.q.a));
		SGN = pop();

		push_integer(2);
		push(X);
		divide();
		push(N);
		push(SGN);
		subtract();
		multiply();
		push(X);
		push(N);
		push(SGN);
		subtract();
		besselj();
		multiply();
		push(X);
		push(N);
		push_integer(2);
		push(SGN);
		multiply();
		subtract();
		besselj();
		subtract();

		return;
	}

//.........这里部分代码省略.........

 bool get_num_branches(contains_app& x, expr* fml, rational& nb) override {
     unsigned sz = m_bv.get_bv_size(x.x());
     nb = power(rational(2), sz);
     return true;
 }

/* Am folosit Fibonacci generalizat */
int main(int args, char** argv) {
	
	int n, m, k, i, j;
	long long **a, *initial;
	
	FILE *fr, *fw;

	fr = fopen("date.in", "r");	
	fw = fopen("date.out", "w");
	
	fscanf(fr, "%d", &n);
	fscanf(fr, "%d", &m);
	fscanf(fr, "%d", &k);
	
	a = (long long**)calloc((k + 1), sizeof(long long*));
	
	for(i = 0; i <= k; i++) {
		a[i] = (long long *)calloc((k + 1), sizeof(long long));
	}
	
	initial = (long long *)calloc((k + 1), sizeof(long long));
	
	/* Pe vectorul termenilor initiali se pun puteri ale lui 2 */
	for(i = 0; i <= k; i++) {
		initial[i] = pow(2, i);
	}

	/* In matrice se pune 1 pe diagonala de deasupra celei principale (pentru a
	se avansa cu cate un termen la fiecare inmultire si a se lua in considerare
	doar ultimi k) si 1 pe ultima linie (al k + 1 termen se calculeaza ca suma 
	celor k + 1 termeni de la pasul anterior) */
	for(i = 0; i < k + 1; i++) {
		for(j = 0; j < k + 1; j++) {
			if(j == i + 1) {
				a[i][j] = 1;
			}
			if(i == k) {
				a[i][j] = 1;
			}
		}
	}

	/* Se ridica a la n - k in vederea calculului celui de-al n-lea termen */
	a = logPowMatrix(a, n - k, k);

	
	/* Se inmulteste vectorul intial cu matricea a */
	long long result = 0;
	for(i = 0; i <= k; i++) {
		result += ((a[k][i] * initial[i]) % MOD);
	}
	
	/* Rezultatul se ridica la puterea numarul de coloane pt numarul total de
	cadre */
	result = power(result, m);
	

	fprintf(fw, "%lli", result);
	
	/* Se elibereaza memoria */
	for(i = 0; i <= k; i++) {
		free(a[i]);
	}
	free(a);
	free(initial);
	fclose(fr);
	fclose(fw);	

	return 0;
}

int main(void)
{
    printf("2 ^ 4 = %.2f\n", power(2, 4));
    return 0;
}

double power(double x, long n)
{
   if(n == 0) return 1;
   if(n < 0) return power ( 1 / x, -n);
   return x * power(x, n - 1);
}