自学教程：C++ subtract函数代码示例

Vector2D computeVelocity(Vector2D start, Vector2D end, COIN *cToBeHit,COIN *cHitBy,Vector2D velRequired){Vector2D fVel;double fVelMag = 0;Vector2D dir;dir = subtract(&end,&start);double mag = magnitude(dir.X,dir.Y);dir.X = dir.X / mag;dir.Y = dir.Y / mag;Vector2D vReqDir;if(cToBeHit == NULL){fVelMag = sqrt(2 * friction * mag * 1.6f );fVel.X = dir.X * fVelMag;fVel.Y = dir.Y * fVelMag;}else{double reqVel = magnitude(velRequired.X,velRequired.Y);vReqDir.X = velRequired.X / reqVel;vReqDir.Y = velRequired.Y / reqVel;double xpVel = (magnitude(velRequired.X,velRequired.Y) * (cToBeHit->mass + cHitBy->mass))/(cHitBy->mass * (1.0f+eCS));double cosT = dotProduct(&dir,&vReqDir);if(cosT < 0) cosT = -cosT;double iniVelReq = sqrt((xpVel / cosT) * (xpVel / cosT) + 2 * friction * mag * 1.6f);fVel.X = iniVelReq * dir.X; fVel.Y = iniVelReq * dir.Y;if(iniVelReq > maxPower){fVel.X = maxPower * dir.X; fVel.Y = maxPower * dir.Y;}}return fVel;}

/*GLuint load_texture(const char* file){   SDL_Surface* surface = IMG_Load(file);   if(!surface)   {       printf("Error while loading texture %s/n", file);       exit(1);   }      GLuint texture;   glPixelStorei(GL_UNPACK_ALIGNMENT,4);   glGenTextures(1, &texture);   glBindTexture(GL_TEXTURE_2D, texture);   SDL_PixelFormat *format = surface->format;   if (format->Amask)   {      gluBuild2DMipmaps(GL_TEXTURE_2D, 4,         surface->w, surface->h, GL_RGBA,GL_UNSIGNED_BYTE, surface->pixels);   }   else   {      gluBuild2DMipmaps(GL_TEXTURE_2D, 3,         surface->w, surface->h, GL_RGB, GL_UNSIGNED_BYTE, surface->pixels);   }   SDL_FreeSurface(surface);   return texture;}*/void add_section(    tree_type_t *t, tree_section_type_t *type, GLfloat w,    GLfloat p1x, GLfloat p1y, GLfloat p1z){    tree_section_t *s = &t->section[t->section_count];    s->type = type;    s->width = w;    s->pos[0] = p1x;    s->pos[1] = p1y;    s->pos[2] = p1z;    if(t->section_count > 0 && s->type)    {	tree_section_t *p = &t->section[t->section_count-1];	subtract(	    p->pos,	    s->pos,	    s->normal,	    3);		s->length = sqrt(dot_prod(s->normal, s->normal, 3));		normalize(	    s->normal,	    s->normal,	    3);	for(int i=0; i<3; i++)	{	    assert(fabs(s->pos[i]+s->normal[i]*s->length - p->pos[i])< 0.01);	    	}	    }    t->section_count++;}

inlinetypename ICR::EnsembleLearning::Discrete<T>::data_tICR::EnsembleLearning::Discrete<T>::CalcLogNorm(vector_data_parameter unLogProbs) {  /* Unnormalised   *If all the probs are very small then can easily get servere numerical errors,   * eg. norm = 0.   * To solve this we subtract most significant log before exponetating    *   (and add it again after).   */  const data_t LogMax = *PARALLEL_MAX(unLogProbs.begin(),unLogProbs.end());  std::vector<T> unLogProbsTmp(unLogProbs.size());    PARALLEL_TRANSFORM( unLogProbs.begin(),unLogProbs.end(), unLogProbsTmp.begin(), 		      subtract(LogMax));  //exponentiate log (prob/max)  std::vector<T> unProbs(unLogProbsTmp.size());  PARALLEL_TRANSFORM( unLogProbsTmp.begin(),unLogProbsTmp.end(), unProbs.begin(), 		      exponentiate());  const data_t norm = PARALLEL_ACCUMULATE(unProbs.begin(), unProbs.end(),0.0) ;    return -std::log(norm)- LogMax;}

void AdaptiveManifoldFilterN::gatherResult(InputArray src_, OutputArray dst_){    int dDepth = src_.depth();    vector<Mat> dstCn(srcCnNum);    if (!adjust_outliers_)    {        for (int i = 0; i < srcCnNum; i++)            divide(sum_w_ki_Psi_blur_[i], sum_w_ki_Psi_blur_0_, dstCn[i], 1.0, dDepth);        merge(dstCn, dst_);    }    else    {        Mat1f& alpha = minDistToManifoldSquared;        double sigmaMember = -0.5 / (sigma_r_*sigma_r_);        multiply(minDistToManifoldSquared, sigmaMember, minDistToManifoldSquared);        cv::exp(minDistToManifoldSquared, alpha);        for (int i = 0; i < srcCnNum; i++)        {            Mat& f = srcCn[i];            Mat& g = dstCn[i];            divide(sum_w_ki_Psi_blur_[i], sum_w_ki_Psi_blur_0_, g);            subtract(g, f, g);            multiply(alpha, g, g);            add(g, f, g);            g.convertTo(g, dDepth);        }        merge(dstCn, dst_);    }}

/* Sigma (stdev) probabilities:1σ 	68.2689492% 	31.7310508% 			1 / 3.15148722σ 	95.4499736% 	4.5500264% 				1 / 21.9778953σ 	99.7300204% 	0.2699796% 				1 / 370.3984σ 	99.993666% 		0.006334% 				1 / 157875σ 	99.9999426697% 	0.0000573303% 			1 / 17442786σ 	99.9999998027% 	0.0000001973% 			1 / 5067973467σ 	99.9999999997440% 	0.000000000256% 	1 / 390682215445"if GYRO way different from ACCEL"			Compares the Gyro reading to the computed Accelerometer DeltasIn other words:  at the angular rate level(Not the angles level)*/BOOL significant_gyro_accel_discrepancy(){	// if GYRO way different from ACCEL (ie. +- 3 sigma)	struct fXYZ g_sigma3 ;	struct fXYZ ad_sigma3;	struct fXYZ sigma3	 ;	struct fXYZ Difference;	// Compute Discrepancy:	subtract( cast_f &GyroAdjusted, cast_f &AccelAngularVelocity, cast_f &Difference );	// Compute Combined 3 Sigma:			scale( cast_f &GyroVariance, 	    3.0, cast_f &g_sigma3  );	scale( cast_f &AccelDeltaVariance,  3.0, cast_f &ad_sigma3 );	add  ( cast_f &g_sigma3, 	cast_f &ad_sigma3,  cast_f &sigma3 );			// if any 1 axis is > 3sigma, then significant	if ( (Difference.x > sigma3.x) || 		 (Difference.y > sigma3.y) || 		 (Difference.z > sigma3.z) ) 		 return TRUE;		else 		 return FALSE;	}

DateTime& DateTime::operator-=(const Duration &duration) {    subtract(duration);    return *this;}

bool emulator::runProgram(){	int accumulator = 0;	// Start at location 100	int loc = 100;	while (true) {		// The value at the memory location		int content = m_memory[loc];		// The operation code		int opCode = content / 10000;		// The address		int address = content % 10000;		if (opCode == 13) {			// Halt.			break;		}		if (opCode == 8) {			// Write.			write(address);		}		if (opCode == 7) {			// Read.			read(address);		}		if (opCode == 6) {			// Store.			store(address, accumulator);		}		if (opCode == 5) {			// Load.			load(address, accumulator);		}		if (opCode == 4) {			// Divide.			divide(address, accumulator);		}		if (opCode == 3) {			// Multiply.			multiply(address, accumulator);		}		if (opCode == 2) {			// Subtract.			subtract(address, accumulator);		}		if (opCode == 1) {			// Add.			add(address, accumulator);		}		// Increment location.		loc++;		if (opCode == 12) {			// Branch if accumulator is positive.			branchPositive(address, accumulator, loc);		}		if (opCode == 11) {			// Branch if accumulator is zero.			branchZero(address, accumulator, loc);		}		if (opCode == 10) {			// Branch if accumulator is negative			branchMinus(address, accumulator, loc);		}		if (opCode == 9) {			// Branch to location			branch(address, loc);		}	}	return true;} /* bool emulator::runProgram(); */

Vector3 Vector3::operator - (Vector3 a){    return subtract(a);}

int InverseElementTransformation::NewtonSolve(const Vector &pt,                                              IntegrationPoint &ip){   MFEM_ASSERT(pt.Size() == T->GetSpaceDim(), "invalid point");   const double phys_tol = phys_rtol*pt.Normlinf();   const int geom = T->GetGeometryType();   const int dim = T->GetDimension();   const int sdim = T->GetSpaceDim();   IntegrationPoint xip, prev_xip;   double xd[3], yd[3], dxd[3], dx_norm = -1.0, err_phys, real_dx_norm = -1.0;   Vector x(xd, dim), y(yd, sdim), dx(dxd, dim);   bool hit_bdr = false, prev_hit_bdr = false;   // Use ip0 as initial guess:   xip = *ip0;   xip.Get(xd, dim); // xip -> x   if (print_level >= 3)   {      NewtonPrint(1, 0.); // iter 0      NewtonPrintPoint(",    ref_pt", x, "/n");   }   for (int it = 0; true; )   {      // Remarks:      // If f(x) := 1/2 |pt-F(x)|^2, then grad(f)(x) = -J^t(x) [pt-F(x)].      // Linearize F(y) at y=x: F(y) ~ L[x](y) := F(x) + J(x) [y-x].      // Newton iteration for F(y)=b is given by L[x_old](x_new) = b, i.e.      // F(x_old) + J(x_old) [x_new-x_old] = b.      //      // To minimize: 1/2 |F(y)-b|^2, subject to: l(y) >= 0, we may consider the      // iteration: minimize: |L[x_old](x_new)-b|^2, subject to l(x_new) >= 0,      // i.e. minimize: |F(x_old) + J(x_old) [x_new-x_old] - b|^2.      // This method uses:      // Newton iteration:    x := x + J(x)^{-1} [pt-F(x)]      // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)]      // Compute the physical coordinates of the current point:      T->Transform(xip, y);      if (print_level >= 3)      {         NewtonPrint(11, 0.); // continuation line         NewtonPrintPoint("approx_pt", y, ", ");         NewtonPrintPoint("exact_pt", pt, "/n");      }      subtract(pt, y, y); // y = pt-y      // Check for convergence in physical coordinates:      err_phys = y.Normlinf();      if (err_phys < phys_tol)      {         if (print_level >= 1)         {            NewtonPrint(1, (double)it);            NewtonPrint(3, dx_norm);            NewtonPrint(30, err_phys);         }         ip = xip;         if (solver_type != Newton) { return Inside; }         return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside;      }      if (print_level >= 1)      {         if (it == 0 || print_level >= 2)         {            NewtonPrint(1, (double)it);            NewtonPrint(3, dx_norm);            NewtonPrint(18, err_phys);         }      }      if (hit_bdr)      {         xip.Get(xd, dim); // xip -> x         if (prev_hit_bdr || it == max_iter || print_level >= 2)         {            prev_xip.Get(dxd, dim); // prev_xip -> dx            subtract(x, dx, dx);    // dx = xip - prev_xip            real_dx_norm = dx.Normlinf();            if (print_level >= 2)            {               NewtonPrint(41, real_dx_norm);            }            if (prev_hit_bdr && real_dx_norm < ref_tol)            {               if (print_level >= 0)               {                  if (print_level <= 1)                  {                     NewtonPrint(1, (double)it);                     NewtonPrint(3, dx_norm);                     NewtonPrint(18, err_phys);                     NewtonPrint(41, real_dx_norm);                  }                  mfem::out << "Newton: *** stuck on boundary!/n";               }               return Outside;//.........这里部分代码省略.........

Decimal & Decimal::modulus(const Decimal & other){    Decimal left(*this);    left.divide(other).truncate(0).multiply(other);    return subtract(left);}

//Overload operator '-' for subtractionVector3 Vector3::operator - (const Vector3 &v){	return subtract(v);}

QBitmap basicSegmentation::segmentImage(const QString &method){  if(currentIndex < 0) {  	QMessageBox::warning(0, "No images open.", "Please open an image prior to segmentation.");  	return(QBitmap());  }  QImage *im = images.at(currentIndex);  int upper;  int lower;  if(method == "Intensity threshold") {  	   	int ath;   	ath = autoThreshold(im, currentChannel);    thresholdRange = new rangeDialog();    thresholdRange->setAutoTh(ath);    thresholdRange->exec();    int r = thresholdRange->result();    if (r==0) {    	QBitmap a = QBitmap();    	return(a);   	}    upper = thresholdRange->getRange()[1];    lower = thresholdRange->getRange()[0];    return(segmentByThreshold(im, upper, lower));  } else if(method == "Halo") {   	int ath;   	ath = autoThreshold(im, currentChannel);    haloSettings = new haloDialog();    haloSettings->setAutoTh(ath);    haloSettings->exec();    int r = haloSettings->result();    if (r==0) {    	return(QBitmap());   	}    upper = haloSettings->getIntRange()[1];    lower = haloSettings->getIntRange()[0];    int offset = haloSettings->getOffRange()[0];    int width = haloSettings->getOffRange()[1];    QBitmap bm = segmentByThreshold(im, upper, lower);  	QBitmap dbm = nonMergeDilate(bm, offset+1);  	QBitmap dbm2 = nonMergeDilate(bm, width);    bm = subtract(dbm, dbm2);    return(bm);  } else if(method == "Threshold and dilate (non merging)") {  	   	int ath;   	ath = autoThreshold(im, currentChannel);    haloSettings = new haloDialog();    haloSettings->setAutoTh(ath);    haloSettings->exec();    int r = haloSettings->result();    if (r==0) {    	QBitmap a = QBitmap();    	return(a);   	}    upper = haloSettings->getIntRange()[1];    lower = haloSettings->getIntRange()[0];    int width = haloSettings->getOffRange()[1];    QBitmap bm = segmentByThreshold(im, upper, lower);  	QBitmap dbm = nonMergeDilate(bm, width);    //QBitmap dbm2 = erode(dbm, 1);    return(dbm);      } else if(method == "Two channel threshold (non merging)") {  	   	int ath1, ath2;   	ath1 = autoThreshold(im, currentChannel);   	//ath2 = autoThreshold(im, 1);    ath2 = 45;    QBitmap nuc = segmentByThreshold(im, 255, ath1, 2);    QBitmap nucg = nonMergeDilate(nuc, 100);    QBitmap cyt = segmentByThreshold(im, 255, ath2, 1);    QBitmap cytr = xand(nucg, cyt);    //QBitmap dbm2 = erode(dbm, 1);    return(cytr);  } else {  	return(QBitmap()); }}

vector2d vector2d::subtract(float x, float y) const{	vector2d subtract(afData[0] - x, afData[1] - y);	return subtract;}

		vec3& vec3::operator-=(const vec3& other)		{			return subtract(other);		}

 currency &currency::operator--() {     return subtract(1); }

int main(){    int i;    FILE *fp;    big K,rid,id,w,a,b,n,q1;    miracl *mip=mirsys(200,256);    for (i=0;i<NPRIMES;i++)    {        pp[i]=mirvar(0);        rem[i]=mirvar(0);    }    w=mirvar(0);    n=mirvar(0);    a=mirvar(0);    b=mirvar(0);    p=mirvar(0);    p1=mirvar(0);         q1=mirvar(0);    K=mirvar(0);    lim1=mirvar(0);    lim2=mirvar(0);    id=mirvar(0);    rid=mirvar(0);    order=mirvar(0);    printf("Enter ID= ");    innum(rid,stdin);    getprime("trap1.dat");    copy(p,n);    getprime("trap2.dat");       multiply(n,p,n);    printf("/ncomposite =/n");    cotnum(n,stdout);    premult(rid,256,id);       while (jack(id,n)!=1)    { /* bad identity - id=256*rid+i */        printf("No Discrete Log. for this ID -- incrementing/n");        incr(id,1,id);    }    getprime("trap1.dat");    copy(p1,q1);    pollard(id,b);    getprime("trap2.dat");    pollard(id,a);    xgcd(p1,q1,K,K,K);     subtract(b,a,w);    mad(w,K,w,q1,q1,w);    if(size(w)<0) add_r(w,q1,w);    subdiv(w,2,w);    multiply(w,p1,w);    add_r(w,a,w);    fp=fopen("secret.dat","w");    otnum(rid,fp);    cotnum(w,fp);    cotnum(n,fp);    fclose(fp);    printf("/nDiscrete log (secret key) /n");    cotnum(w,stdout);    powltr(PROOT,w,n,id);    subdiv(id,256,id);    otstr(id,mip->IOBUFF);    printf("Check Identity= %s/n",mip->IOBUFF);    return 0;}

void Matrix::subtract(const Matrix& mat){    subtract(*this, mat, this);}

// Ensure subtract() subtracts two integers correctly.void test_subtract(void **state) {	assert_int_equal(subtract(3, 3), 0);	assert_int_equal(subtract(3, -3), 6);}

Vec3 operator-(const Vec3 & lhs, const Vec3 & rhs) {	return subtract(lhs, rhs);}

Result VideoHandle::getDirectionPoints(){    Mat frame;    Mat prev;    vector<Point> List;	for (int i=30;i--;) getFrame();    prev = getFrame();    while(true)    {        int prev_clock = clock();        frame = getFrame();        if(frame.empty()) break;        Mat temp;        subtract(prev, frame, temp);        const Size size = Size(160, 120);        resize(temp, temp, size, 0, 0, CV_INTER_LINEAR);        cvtColor(temp, temp, CV_BGR2GRAY);        threshold(temp, temp, 20, 255, CV_THRESH_BINARY);        morphologyEx(temp.clone(), temp, MORPH_OPEN, Mat::ones(3, 3, CV_8U));        Moments m = ::moments(temp);        Point p = Point(m.m10/m.m00, m.m01/m.m00);		// Mat_<Point2f> points(1,1), dst(1,1);		// points(0) = Point2f(p.x,p.y);		// undistortPoints(points, dst, distortmtx, distortdist);		// p.x = - dst(0).y * size.width;		// p.y = - dst(0).x * size.height;        cout << "Point : " << p.x << " " << p.y << endl;        List.push_back(p);        if (p.x < 5 || p.y < 5 || p.x > size.width - 6 || p.y > size.height - 6) {            cout << "455555555555555" << endl;            List.clear();        }        const int TIMES = 1;        Point sum = Point(0, 0);        for(int i=1;i<=TIMES;i++)        {            if(List.size() < TIMES + 1) continue;            Point a = List[List.size() - i];            Point b = List[List.size() - i - 1];            Point sub = Point(a.x-b.x, a.y-b.y);            sum.x += sub.x;            sum.y += sub.y;        }        sum.x /= TIMES;        sum.y /= TIMES;        cout << "vector : " << sum.x << " " << sum.y << endl;        if(abs(sum.x) >= 2 || abs(sum.y) >= 2) {            Result ret = generateOutput(p, Point(p.x+sum.x, p.y+sum.y));            ret.angle *= -1;            return ret;        }        int now_clock = clock();        double speed = double(now_clock - prev_clock) / CLOCKS_PER_SEC;        cout << "speed : " << speed << " " << (1.0/speed) << endl;    }}

int main(){   polynomial a, b, c, sum, difference, product, quotient, normal;  double coeff;  int order;  //these will be converted into dynamic arrays later on.  //they just store the coeffs to be used for the polynomials.  double data1[3];  data1[2] = 5;  double data2[3] = {5.2, -2, 4};  double data3[7]= {5.0, 2.0, 20, 0, 5, 2, 40.2};  //clear the screen to focus on this single session's tests  system("clear");  //////////////////////////////////////////////////////////////////  //creating base polynomials we'll be using for the entirity of  //these test cases.  //  //NOTE: second parameter of initalisePolynomial(poly, size)   //      MUST match the test case arrays above.  //////////////////////////////////////////////////////////////////    createPolynomial(&a, 3, data1);  createPolynomial(&b, 3, data2);  createPolynomial(&c, 7, data3);    printf("--------------------------------------------------------------------------------/n");  printf("These are the starting polynomials which will be used for the later operations/n");  printf("--------------------------------------------------------------------------------/n");  printf("Polynomial (A): ");  verify(printPolynomial(&a));  printf("Polynomial (B): ");  verify(printPolynomial(&b));  printf("Polynomial (C): ");  verify(printPolynomial(&c));  //////////////////////////////////////////////////////////////////  //testing add/subtract polynomials of same length  //////////////////////////////////////////////////////////////////  printf("--------------------------------------------------------------------------------/n");  printf("testing add/subtract polynomials of same length/n");  printf("--------------------------------------------------------------------------------/n");  printf("Polynomial (A): ");  verify(printPolynomial(&a));  printf("Polynomial (B): ");  verify(printPolynomial(&b));  //add  printf("A + B =  ");  verify(add(&a, &b, &sum));  verify(printPolynomial(&sum));  //subtract  printf("A - B =  ");  verify(subtract(&a, &b, &difference));  verify(printPolynomial(&difference));  deletePolynomial(&sum);  deletePolynomial(&difference);  //////////////////////////////////////////////////////////////////  //testing add/subtract polynomials of different length    //////////////////////////////////////////////////////////////////  printf("--------------------------------------------------------------------------------/n");  printf("testing add/subtract polynomials of different length/n");  printf("--------------------------------------------------------------------------------/n");  printf("Polynomial (A): ");  verify(printPolynomial(&a));  printf("Polynomial (C): ");  verify(printPolynomial(&c));  printf("A + C =  ");  verify(add(&a, &c, &sum));  verify(printPolynomial(&sum));  printf("A - C =  ");  verify(subtract(&a, &c, &difference));  verify(printPolynomial(&difference));  deletePolynomial(&sum);  deletePolynomial(&difference);  //////////////////////////////////////////////////////////////////  //testing multiplcation/division by a scalar    //////////////////////////////////////////////////////////////////  printf("--------------------------------------------------------------------------------/n");  printf("testing multiplcation/division by a scalar/n");  printf("--------------------------------------------------------------------------------/n");  printf("Polynomial (A): ");  printPolynomial(&a);  printf("Polynomial (C): ");  printPolynomial(&c);  coeff = 5;  //multiply  printf("A * %.2lf =  ", coeff);  multiply(&a, coeff, &product);  printPolynomial(&product);//.........这里部分代码省略.........

Result VideoHandle::getDirection(){    vector<Point> List;    Mat frame;    Mat prev;    for(int i=10;i--;)        getUndistortFrame();    prev = getUndistortFrame();    while(true)    {        int prev_clock = clock();        frame = getUndistortFrame();        if(frame.empty()) exit(1);        Mat temp;        subtract(prev, frame, temp);        const Size size = Size(160, 120);        resize(temp, temp, size, 0, 0, CV_INTER_LINEAR);// FIXME        cvtColor(temp, temp, CV_BGR2GRAY);        threshold(temp, temp, 20, 255, CV_THRESH_BINARY);        Moments m = ::moments(temp);        Point p = Point(m.m10/m.m00, m.m01/m.m00);        cout << "Point : " << p.x << " " << p.y << endl;        List.push_back(p);        if (p.x < 5 || p.y < 5 || p.x > size.width-5 || p.y > size.height-5) { // FIXME            cout << "233333333333333" << endl;            List.clear();        }        const int TIMES = 1;// FIXME        Point sum = Point(0, 0);        for(int i=1;i<=TIMES;i++)        {            if(List.size() < TIMES + 1) continue;            Point a = List[List.size() - i];            Point b = List[List.size() - i - 1];            Point sub = Point(a.x-b.x, a.y-b.y);            sum.x += sub.x;            sum.y += sub.y;        }        sum.x /= TIMES;        sum.y /= TIMES;        cout << "vector : " << sum.x << " " << sum.y << endl;        if(sum.x*sum.x + sum.y*sum.y >= 30) {    //FIXME 建议修改为||关系            Result ret = generateOutput(p, Point(p.x+sum.x, p.y+sum.y));            ret.angle *= -1;            return ret;        }        int now_clock = clock();        double speed = double(now_clock - prev_clock) / CLOCKS_PER_SEC;        cout << "speed : " << speed << " " << (1.0/speed) << endl;    }}

Vector4 Vector4::operator - (Vector4 a){    return subtract(a);}

//.........这里部分代码省略.........                {                    u=muldvm(mr_mip->w1->w[n-1],mr_mip->w1->w[n-2],m,&sr);                    v=muldvm(mr_mip->w2->w[n-1],mr_mip->w2->w[n-2],m,&sr);                }#endif#endif#ifndef MR_SIMPLE_BASE            }            else            {#ifdef mr_dltype                if (n>2)                { /* squeeze out as much significance as possible */                    u=muldiv(mr_mip->w1->w[n-1],mr_mip->base,mr_mip->w1->w[n-2],m,&sr);                    u=u*mr_mip->base+muldiv(sr,mr_mip->base,mr_mip->w1->w[n-3],m,&sr);                    v=muldiv(mr_mip->w2->w[n-1],mr_mip->base,mr_mip->w2->w[n-2],m,&sr);                    v=v*mr_mip->base+muldiv(sr,mr_mip->base,mr_mip->w2->w[n-3],m,&sr);                }                else                {                    u=(mr_large)mr_mip->base*mr_mip->w1->w[n-1]+mr_mip->w1->w[n-2];                    v=(mr_large)mr_mip->base*mr_mip->w2->w[n-1]+mr_mip->w2->w[n-2];                    last=TRUE;                }#else                u=muldiv(mr_mip->w1->w[n-1],mr_mip->base,mr_mip->w1->w[n-2],m,&sr);                v=muldiv(mr_mip->w2->w[n-1],mr_mip->base,mr_mip->w2->w[n-2],m,&sr);#endif            }#endif        }        dplus=TRUE;        a=1; b=0; c=0; d=1;        forever        { /* work only with most significant piece */            if (last)            {                if (v==0) break;                q=qdiv(u,v);                if (q==0) break;            }            else            {                if (dplus)                {                     if ((mr_small)(v-c)==0 || (mr_small)(v+d)==0) break;                    q=qdiv(u+a,v-c);                    if (q==0) break;                    if (q!=qdiv(u-b,v+d)) break;                }                else                 {                    if ((mr_small)(v+c)==0 || (mr_small)(v-d)==0) break;                    q=qdiv(u-a,v+c);                    if (q==0) break;                    if (q!=qdiv(u+b,v-d)) break;                }            }            if (q==1)            {                if ((mr_small)(b+d) >= MAXBASE) break;                 r=a+c;  a=c; c=r;                r=b+d;  b=d; d=r;                lr=u-v; u=v; v=lr;                  }            else            {                 if (q>=MR_DIV(MAXBASE-b,d)) break;                r=a+q*c;  a=c; c=r;                r=b+q*d;  b=d; d=r;                lr=u-q*v; u=v; v=lr;            }            iter++;            dplus=!dplus;        }        iter%=2;    }    if (s==MINUS) iter++;    if (iter%2==1) subtract(_MIPP_ y,mr_mip->w3,mr_mip->w3);    if (xd!=yd)    {        negify(x,mr_mip->w2);        mad(_MIPP_ mr_mip->w2,mr_mip->w3,mr_mip->w1,y,mr_mip->w4,mr_mip->w4);        copy(mr_mip->w4,yd);    }    copy(mr_mip->w3,xd);    if (z!=xd && z!=yd) copy(mr_mip->w1,z);    MR_OUT    return (size(mr_mip->w1));}

int main(){  /*  factoring program using Brents method */    long k,r,i,m,iter;    big x,y,z,n,q,ys,c3;    miracl *mip=mirsys(50,0);    x=mirvar(mip,0);    y=mirvar(mip,0);    ys=mirvar(mip,0);    z=mirvar(mip,0);    n=mirvar(mip,0);    q=mirvar(mip,0);    c3=mirvar(mip,3);    printf("input number to be factored/n");    cinnum(mip,n,stdin);    if (isprime(mip,n))    {        printf("this number is prime!/n");        return 0;    }    m=10L;    r=1L;    iter=0L;    do    {        printf("iterations=%5ld",iter);        convert(mip,1,q);        do        {            copy(y,x);            for (i=1L;i<=r;i++)                mad(mip,y,y,c3,n,n,y);            k=0;            do            {                iter++;                if (iter%10==0) printf("/b/b/b/b/b%5ld",iter);                fflush(stdout);                  copy(y,ys);                for (i=1L;i<=mr_min(m,r-k);i++)                {                    mad(mip,y,y,c3,n,n,y);                    subtract(mip,y,x,z);                    mad(mip,z,q,q,n,n,q);                }                egcd(mip,q,n,z);                k+=m;            } while (k<r && size(z)==1);            r*=2;        } while (size(z)==1);        if (compare(z,n)==0) do         { /* back-track */            mad(mip,ys,ys,c3,n,n,ys);            subtract(mip,ys,x,z);        } while (egcd(mip,z,n,z)==1);        if (!isprime(mip,z))             printf("/ncomposite factor ");        else printf("/nprime factor     ");        cotnum(mip,z,stdout);        if (compare(z,n)==0) return 0;        divide(mip,n,z,n);        divide(mip,y,n,n);    } while (!isprime(mip,n));    printf("prime factor     ");    cotnum(mip,n,stdout);    return 0;}

ATF_TC_BODY(subtractrdataset, tc) {	UNUSED(tc);	subtract(NULL);}

booloperator-(const mat<N,S> &a, const mat<N,S> &b) { 	return subtract(a, b); }