C++ Mul函数代码示例

int main()
{
    // printf("hello, world");

    char result[BUFFER_SIZE];
    char remainder[BUFFER_SIZE];

    // please make sure the bit length is enough before calculate
    // especially when you do a large power operation
    // you can change it by define: BIG_INT_BIT_LEN
    // the default bit length for BigInt is 1024

    // routine test
    puts(Add("2010", "4", result));
    puts(Sub("0", "2014", result));
    puts(Mul("2", "43", result));
    puts(Div("86", "10", result, remainder));
    puts(remainder);
    puts(Mod("-86", "10", result));
    puts(PowMod("7", "80", "86", result));

    // BigInt test
    puts(Sub("233333333333333333333333333333333333333333333333", "33", result));
    puts(Mul("2333333333333333333333333333333", "2333333333333333333", result));
    puts(Div("2333333333333333333333333333333", "2333333333333333332", result, remainder));
    puts(remainder);
    puts(Pow("8", "86", result));

    return 0;
}

bool RopeJoint::SolvePositionConstraints(const SolverData& data)
{
	Vec2 cA = data.positions[m_indexA].c;
	float32 aA = data.positions[m_indexA].a;
	Vec2 cB = data.positions[m_indexB].c;
	float32 aB = data.positions[m_indexB].a;

	Rot qA(aA), qB(aB);

	Vec2 rA = Mul(qA, m_localAnchorA - m_localCenterA);
	Vec2 rB = Mul(qB, m_localAnchorB - m_localCenterB);
	Vec2 u = cB + rB - cA - rA;

	float32 length = u.Normalize();
	float32 C = length - m_maxLength;

	C = Clamp(C, 0.0f, maxLinearCorrection);

	float32 impulse = -m_mass * C;
	Vec2 P = impulse * u;

	cA -= m_invMassA * P;
	aA -= m_invIA * Cross(rA, P);
	cB += m_invMassB * P;
	aB += m_invIB * Cross(rB, P);

	data.positions[m_indexA].c = cA;
	data.positions[m_indexA].a = aA;
	data.positions[m_indexB].c = cB;
	data.positions[m_indexB].a = aB;

	return length - m_maxLength < linearSlop;
}

int  ConjGradientMod(float3N &X, float3Nx3N &A, float3N &B,int3 hack)
{
// obsolete!!!
       // Solves for unknown X in equation AX=B
       conjgrad_loopcount=0;
       int n=B.count;
       float3N q(n),d(n),tmp(n),r(n);
       r = B - Mul(tmp,A,X);    // just use B if X known to be zero
       r[hack[0]] = r[hack[1]] = r[hack[2]] = float3(0,0,0);
       d = r;
       float s = dot(r,r);
       float starget = s * squared(conjgrad_epsilon);
       while( s>starget && conjgrad_loopcount++ < conjgrad_looplimit)
       {
               Mul(q,A,d); // q = A*d;
               q[hack[0]] = q[hack[1]] = q[hack[2]] = float3(0,0,0);
               float a = s/dot(d,q);
               X = X + d*a;
               if(conjgrad_loopcount%50==0)
               {
                       r = B - Mul(tmp,A,X);
                       r[hack[0]] = r[hack[1]] = r[hack[2]] = float3(0,0,0);
               }
               else
               {
                       r = r - q*a;
               }
               float s_prev = s;
               s = dot(r,r);
               d = r+d*(s/s_prev);
               d[hack[0]] = d[hack[1]] = d[hack[2]] = float3(0,0,0);
       }
       conjgrad_lasterror = s;
       return conjgrad_loopcount<conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
}

int  ConjGradientFiltered(float3N &X, const float3Nx3N &A, const float3N &B,const float3Nx3N &S)
{
       // Solves for unknown X in equation AX=B
       conjgrad_loopcount=0;
       int n=B.count;
       float3N q(n),d(n),tmp(n),r(n);
       r = B - Mul(tmp,A,X);    // just use B if X known to be zero
       filter(r,S);
       d = r;
       float s = dot(r,r);
       float starget = s * squared(conjgrad_epsilon);
       while( s>starget && conjgrad_loopcount++ < conjgrad_looplimit)
       {
               Mul(q,A,d); // q = A*d;
               filter(q,S);
               float a = s/dot(d,q);
               X = X + d*a;
               if(conjgrad_loopcount%50==0)
               {
                       r = B - Mul(tmp,A,X);
                       filter(r,S);
               }
               else
               {
                       r = r - q*a;
               }
               float s_prev = s;
               s = dot(r,r);
               d = r+d*(s/s_prev);
               filter(d,S);
       }
       conjgrad_lasterror = s;
       return conjgrad_loopcount<conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
}

void Intermediate ( Quaternion* q0, Quaternion* q1, Quaternion* q2,
   Quaternion* a, Quaternion* b)
{
   /* assert:  q0, q1, q2 are unit quaternions */
   /*Quaternion q0inv = q0.UnitInverse();
   Quaternion q1inv = q1.UnitInverse();
   Quaternion p0 = q0inv*q1;
   Quaternion p1 = q1inv*q2;
   Quaternion arg = 0.25*(p0.Log()-p1.Log());
   Quaternion marg = -arg;
 
   a = q1*arg.Exp();
   b = q1*marg.Exp();
   */

   Quaternion p0, p1;
   UnitInverse(&p0, q0);
   UnitInverse(&p1, q1);
   MulSelf(&p0, q1);
   MulSelf(&p1, q2);
   LogSelf(&p0);
   LogSelf(&p1);
   Sub(a,&p0, &p1);
   MulScalSelf(a, 0.25);
   Neg(b, a);

   ExpSelf(a);
   ExpSelf(b);
   Mul(a, q1, a);
   Mul(b, q1, b);
}

u64 SumOfDivisorsOfBinomialCoefficient(u64 n, u64 k) {
    if (n == 2) {
        return k == 1 ? 3 : 1;
    } else if (n < 2) {
        return 1;
    }

    Vector primes = GetPrimes(n);
    Vector powers;
    for (auto x: primes) {
        powers.push_back(
            PowerOfDivisor(x, n) -
            PowerOfDivisor(x, k) -
            PowerOfDivisor(x, n - k) 
        );
    }

    u64 result = 1;
    for (u64 i = 0; i < primes.size(); ++i) {
        u64 sum_of_powers = 0;
        u64 x = 1;
        for (u64 p = 0; p <= powers[i]; ++p) {
            sum_of_powers = Sum(sum_of_powers, x);
            x = Mul(x, primes[i]);
        }
        result = Mul(result, sum_of_powers);
    }
    return result;
}

void mul(Ring_Element& ans,const Ring_Element& a,const Ring_Element& b)
{
  if (a.rep!=b.rep)   { throw rep_mismatch(); }
  if (a.FFTD!=b.FFTD) { throw pr_mismatch();  }
  ans.partial_assign(a);
  if (ans.rep==evaluation)
    { // In evaluation representation, so we can just multiply componentwise
      for (int i=0; i<(*ans.FFTD).phi_m(); i++)
        { Mul(ans.element[i],a.element[i],b.element[i],(*a.FFTD).get_prD()); }
    }
  else if ((*ans.FFTD).get_twop()!=0)
    { // This is the case where m is not a power of two

      // Here we have to do a poly mult followed by a reduction
      // We could be clever (e.g. use Karatsuba etc), but instead
      // we do the school book method followed by term re-writing

      // School book mult
      vector<modp> aa(2*(*a.FFTD).phi_m());
      for (int i=0; i<2*(*a.FFTD).phi_m(); i++)
        { assignZero(aa[i],(*a.FFTD).get_prD()); }
      modp temp;
      for (int i=0; i<(*a.FFTD).phi_m(); i++)
        { for (int j=0; j<(*a.FFTD).phi_m(); j++)
	    { Mul(temp,a.element[i],b.element[j],(*a.FFTD).get_prD()); 
              int k=i+j;
              Add(aa[k],aa[k],temp,(*a.FFTD).get_prD());
            }
        }
      // Now apply reduction, assumes Ring.poly is monic
      for (int i=2*(*a.FFTD).phi_m()-1; i>=(*a.FFTD).phi_m(); i--)
        { reduce_step(aa,i,*a.FFTD);
          assignZero(aa[i],(*a.FFTD).get_prD()); 
        }
     // Now stick into answer
     for (int i=0; i<(*ans.FFTD).phi_m(); i++)
       { ans.element[i]=aa[i]; }
    }
  else if ((*ans.FFTD).get_twop()==0)
    { // m a power of two case
      Ring_Element aa(*ans.FFTD,ans.rep);
      modp temp;
      for (int i=0; i<(*ans.FFTD).phi_m(); i++)
        { for (int j=0; j<(*ans.FFTD).phi_m(); j++)
            { Mul(temp,a.element[i],b.element[j],(*a.FFTD).get_prD());
              int k=i+j;
              if (k>=(*ans.FFTD).phi_m())
                 { k-=(*ans.FFTD).phi_m();
                   Negate(temp,temp,(*a.FFTD).get_prD());
                 }
              Add(aa.element[k],aa.element[k],temp,(*a.FFTD).get_prD());
            }
        }
      ans=aa;
    }
  else
    { throw not_implemented(); }
}

/*----------------------------------------------------------------------------------------------
	Core routine to map between local and Utc. This assumes that DST doesn't kick in
	near a year boundary so that when the switch happens, the DST and STD times are in the
	same year.
----------------------------------------------------------------------------------------------*/
int64 TimeMapper::MapMsec(int64 msecSrc, bool fToUtc)
{
	int64 msecT;
	int msecDay;
	int day;
	int dayMin;
	int dyear;
	int yday;
	int yt;
	int ymin;

	// First determines the year type and minute within the year.

	// Mod to 400 year period.
	msecT = ModPos(msecSrc, kmsecPerPeriod);

	// Find the day within the 400 year period.
	day = (int)(msecT / kmsecPerDay);
	Assert(0 <= day && day < kdayPerPeriod);

	// Get the time within the day.
	msecDay = (int)(msecT - day * (int64)kmsecPerDay);
	Assert(0 <= msecDay && msecDay < kmsecPerDay);

	// Find the year within the 400 year period (dyear) and the day within that year (yday).
	Assert(0 <= day && day < 146097);
	dyear = YearFromDayInPeriod(day, &yday);
	Assert(0 <= dyear && dyear < 400);
	Assert(0 <= yday && (yday < 365 || yday == 365 && SilTime::IsLeapYear(dyear + kyearBase)));

	// Calculate the day number (within the 400 year period) of the first day of
	// the year (dayMin).
	dayMin = day - yday;
	Assert(dayMin == DayFromYearInPeriod(dyear));

	// Get the year type. kdayMonday is used because T0 is a Monday.
	yt = ModPos(dayMin + kwdayMonday, kdayPerWeek);
	if (SilTime::IsLeapYear(dyear + kyearBase))
		yt += kdayPerWeek;

	// Calculate the minute within the year.
	ymin = yday * kminPerDay + msecDay / kmsecPerMin;

	if (!fToUtc)
	{
		// Mapping from utc to local.
		if (m_rgyminMin[yt] <= ymin && ymin < m_rgyminLim[yt])
			return msecSrc + Mul(m_dminTz2, kmsecPerMin);
		return msecSrc + Mul(m_dminTz1, kmsecPerMin);
	}

	// The overlap is always ambiguous and there's nothing we can do about it.
	ymin -= m_dminTz2;
	if (m_rgyminMin[yt] <= ymin && ymin < m_rgyminLim[yt])
		return msecSrc - Mul(m_dminTz2, kmsecPerMin);
	return msecSrc - Mul(m_dminTz1, kmsecPerMin);
}

//If you ever call something with 1 octave, call this instead
 void plainSIMD3d(SIMD* out, Settings*  S,ISIMDNoise3d noise )
{

	SIMD vfx = Mul(S->x.m, S->frequency);
	SIMD vfy = Mul(S->y.m, S->frequency);
	SIMD vfz = Mul(S->z.m, S->frequency);

	*out = noise(&vfx, &vfy, &vfz);
}

	Rect Rect::operator*(mat4 matrix) const
	{
		matrix = Transpose(matrix);
		vec4 returnVec1 = Mul(vec4(m_LeftBottom.x, m_LeftBottom.y, 0, 1), matrix);
		vec4 returnVec2 = Mul(vec4(m_RightBottom.x, m_RightBottom.y, 0, 1), matrix);
		vec4 returnVec3 = Mul(vec4(m_LeftTop.x, m_LeftTop.y, 0, 1), matrix);
		vec4 returnVec4 = Mul(vec4(m_RightTop.x, m_RightTop.y, 0, 1), matrix);

		return Rect(vec2(returnVec1.x , returnVec1.y),
					vec2(returnVec2.x , returnVec2.y),
					vec2(returnVec3.x , returnVec3.y),
					vec2(returnVec4.x , returnVec4.y));
	}

void Power( Matrix A, int P )
{
  Matrix R;
  memset(R, 0, sizeof(R));
  for (int i = 0; i < N; i++)
    R[i][i] = 1;
  for (; P; P >>= 1)
  {
    if (P & 1)
      Mul(R, A);
    Mul(A, A);
  }
  A = R;
}

void f ()
{
   Matrix<2, 3> q;
   Matrix<2, 4> a;
   Matrix<4, 3> b;
   q = Mul (q, a, b);
}

/****************************************************
	FFT()

	参数：

		TD为时域值
		FD为频域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数实现快速傅立叶变换
****************************************************/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;

	/*计算傅立叶变换点数*/
	count=1<<power;
	
	/*分配运算所需存储器*/
	W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	
	/*计算加权系数*/
	for(i=0;i<count/2;i++)
	{
		angle=-i*PI*2/count;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}
	
	/*将时域点写入存储器*/
	memcpy(X1,TD,sizeof(COMPLEX)*count);
	
	/*蝶形运算*/
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for(i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}
	
	/*重新排序*/
	for(j=0;j<count;j++)
	{
		p=0;
		for(i=0;i<power;i++)
		{
			if (j&(1<<i)) p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	
	/*释放存储器*/
	free(W);
	free(X1);
	free(X2);
}

void FastMul(const BigInt &A,const BigInt &B, BigInt &C) {
	ulong Len = 2;
	while ( Len <  A.Size + B.Size ) Len *=2;

	if ( Len < 40 ) {
		BigInt Atmp(A), Btmp(B);
		Mul(Atmp,Btmp,C);
		return;
	}

	ulong x;
	const ushort *a=A.Coef, *b=B.Coef;

	for (x = 0; x < A.Size; x++)    LongNum1[x] = a[x];
	for (; x < Len; x++) 		LongNum1[x] = 0.0;

	FHT_F(LongNum1,Len);

	if (a == b) {
		FHTConvolution(LongNum1,LongNum1,Len);
	} else {
		for (x = 0; x < B.Size; x++)    LongNum2[x] = b[x];
		for (; x < Len; x++) 		LongNum2[x] = 0.0;
		FHT_F(LongNum2,Len);
		FHTConvolution(LongNum1,LongNum2,Len);
	 }

	FHT_T(LongNum1,Len);

	CarryNormalize(LongNum1, Len, C);
}

char *Ricerca_e_deriva(char *funzione_1, char *funzione_2, char *operatore, char *output)
	{
		output = (char *)malloc(sizeof(char)*10000);
		if ( strcmp(operatore,"x") == 0 || ( strcmp(operatore,"plus") != 0 && strcmp(operatore,"mul") != 0 && strcmp(operatore,"sot") != 0 && strcmp(operatore,"pow") != 0 ))
			output = D_Fondamentali(operatore); // Se il primo operando è una x oppure non è nessuna delle funzioni previste allora chiama D_Fondamentali 	
	
		if(strcmp(operatore,"pow") == 0)	//	Potenza
			output = Pow(funzione_1,funzione_2, output);
	
		if(strcmp(operatore,"mul") ==  0)	//	Prodotto
			output = Mul(funzione_1 ,funzione_2, output);	
		
		if(strcmp(operatore,"div") == 0)	//	Rapporto
			output = Div(funzione_1 ,funzione_2, output);	
		
		if(strcmp(operatore,"plus") == 0)	//	Somma
			output = Sum(funzione_1,funzione_2, output);
		
		if(strcmp(operatore,"sot") == 0)	//	Differenza
			output = Sot(funzione_1,funzione_2, output);

		if(strcmp(operatore,"sin") == 0)	//	Seno
			output = Sin(funzione_1, "", output);
		
		if(strcmp(operatore,"cos") == 0)	//	Coseno
			output = Cos(funzione_1, "", output);
		
		
		return output;		// Ritorno la funzione già derivata	
}