C++ Mult函数代码示例

BigInteger Karatsuba3Multiplication::Mult(const BigInteger&A, const BigInteger&B){

    int s1=A.GetSize();
    int s2=B.GetSize();

    int n,k;
    //for(n=1, m=0; n<s1 || n<s2; n*=3,++m); k=n/3;
    if(s1>s2) n=s1; else n=s2;
    k=n/3;


    if(n<100) return Multiplies::Mult(A,B);

    BigInteger A1 = A, A2 = A>>k, A3=A>>(2*k);
    A1.Resize(k);
    A1.ClearFirstZeros();
    A2.Resize(k);
    A2.ClearFirstZeros();

    BigInteger B1 = B, B2 = B>>k, B3=B>>(2*k);
    B1.Resize(k);
    B1.ClearFirstZeros();
    B2.Resize(k);
    B2.ClearFirstZeros();

    BigInteger A1B1 = Mult(A1,B1);
    BigInteger A2B2 = Mult(A2,B2);
    BigInteger A3B3 = Mult(A3,B3);

    BigInteger AB=A1B1+  (( Mult(A1+A2,B1+B2)-(A1B1+A2B2))<<k)+  (( Mult(A1+A3,B1+B3)-(A1B1+A3B3)+A2B2)<<2*k) + (( Mult(A2+A3,B2+B3)-(A2B2+A3B3))<<3*k) + (A3B3<<4*k);

    AB.ClearFirstZeros();
    return AB;
};

void Combine( int Positions[], FILE *InputFiles[], int NumberInputs )
{
    int OddNibble ;
    BOOLEAN Odd ;
    int i ;
    int j ;
    G13 C[15] ;
    G13 c ;
    G13 Y[15] ;
    BOOLEAN FillNibbles( G13 *Nibbles, FILE *InputFiles[], int NumberInputs );

    setmode( fileno( stdout ), O_BINARY ) ;

    /* Compute the coefficients by which the nibbles from the
     * various input files can be combined to produce the output
     * file.
     *     If X(i) is the "position" of the ith input file, and
     * Y(i) is the value of a particular nibble in the ith file,
     * we will find coefficients C(i) such that
     * 
     *     p = C(1) * Y(1) + C(2) * Y(2) + ...
     * 
     * where p is the appropriate nibble for the output file.
     * 
     *     The formula for the Cs is:
     * 
     *     C(i) = product of (   X(j) / ( X(i) - X(j) )   ) for all j != i .
     * 
     * 
     */

    for ( i = 0; i < NumberInputs; i++ )
    {
     c = 1 ;
     for ( j = 0; j < NumberInputs; j++ )
         if ( j != i )
          c = Mult( c,
               Div( Positions[j], Positions[i] ^ Positions[j] ) ) ;
     C[i] = c ;
    }

    /*
     *  Now, process the input files:
     */

    Odd = TRUE ;
    while ( FillNibbles( Y, InputFiles, NumberInputs ) )
    {
     c = 0 ;
     for ( i = 0; i < NumberInputs; i++ )
         c ^= Mult( C[i], Y[i] ) ;

     if ( Odd )
         OddNibble = c ;
     else
         putchar( ( OddNibble << 4 ) | c ) ;

     Odd = !Odd ;
    }
}

void display(void)
{
	// clear the screen
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	GLfloat modelView[16], camMatrix[16], cam_Matrix_skybox[16], total[16];
	
	printError("pre display");
	
	check_keys();

        t += 1;

	glUseProgram(program);

	// Build matrix
	
	lookAt(&cam_pos, &obj_pos, up.x, up.y, up.z, camMatrix);
        CopyMatrix(camMatrix, cam_Matrix_skybox);
        cam_Matrix_skybox[3] = 0;
        cam_Matrix_skybox[7] = 0;
        cam_Matrix_skybox[11] = 0;
        cam_Matrix_skybox[15] = 1;
        // disable z-buffer for skybox
        glDisable(GL_DEPTH_TEST);

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        // Draw skybox
        // Set the skybox variable 
        glUniform1i(glGetUniformLocation(program, "skybox"), 1);
        T(0,-0.5,0,modelView);
        Mult(cam_Matrix_skybox, modelView, total);
        glUniformMatrix4fv(glGetUniformLocation(program, "mdlMatrix"), 1, GL_TRUE, total);
        glBindTexture(GL_TEXTURE_2D, skytex);
        DrawModel(skybox, program, "inPosition", "inNormal", "inTexCoord");

        glUniform1i(glGetUniformLocation(program, "skybox"), 0);
        glEnable(GL_DEPTH_TEST);
        // Bind terrain
	IdentityMatrix(modelView);
	Mult(camMatrix, modelView, total);
	glUniformMatrix4fv(glGetUniformLocation(program, "mdlMatrix"), 1, GL_TRUE, total);
	// Bind Our Texture tex1
	glBindTexture(GL_TEXTURE_2D, tex1);		
	DrawModel(tm, program, "inPosition", "inNormal", "inTexCoord");

        // Draw sphere
        sphere_pos = update_sphere(t, total);
        sphere_pos.y = calculate_height(sphere_pos.x, sphere_pos.z, ttex.width, vertexArray);
        T(sphere_pos.x, sphere_pos.y, sphere_pos.z, trans);
        Mult(camMatrix, trans, total);
        //Ry(0.01*t, roty);
	glUniformMatrix4fv(glGetUniformLocation(program, "mdlMatrix"), 1, GL_TRUE, total);
	glBindTexture(GL_TEXTURE_2D, spheretex);		
	DrawModel(sphere, program, "inPosition", "inNormal", "inTexCoord");

	printError("display 2");
	
	glutSwapBuffers();
}

	void Work(void) {
		int64 P[5] = {1, 2};
		while (!Zero(a) && !Zero(b)) {
			bool Odda = a[1] & 1, Oddb = b[1] & 1;
			if (Odda && Oddb) {
				Decreas(a, b);
			}else
			if (!Odda && !Oddb) {
				Div(a);
				Div(b);
				Mult(delta, P);
			}else
			if (Odda && !Oddb) {
				Div(b);
			}else
			if (!Odda && Oddb) {
				Div(a);
			}
		}
		if (Zero(a)) {
			Mult(b, delta);
			Print(b);
		}
		else {
			Mult(a, delta);
			Print(a);
		}
	}

/*
 *	Takes care of orbit and rotation
 */
void UpdatePlanetMovement(GLint t)
{
	GLuint i;
	mat4 ModelToWorld;
	for(i = 0; i < numberOfPlanets; i++)
	{
		t = glutGet(GLUT_ELAPSED_TIME) -  planetsList[i].timeOfCreation;

		vec3 arbAxis = {0,1,0};
		ModelToWorld = ArbRotate(planetsList[i].rotationalAxis, planetsList[i].rotationalSpeed*(GLfloat)t); //IdentityMatrix(); //Rotation around own axis
		ModelToWorld = Mult(T(planetsList[i].startingPosition.x, planetsList[i].startingPosition.y, planetsList[i].startingPosition.z), ModelToWorld); //Offset

		if(!(planetsList[i].startingPosition.x == planetsList[0].startingPosition.x 
		  && planetsList[i].startingPosition.y == planetsList[0].startingPosition.y 
		  && planetsList[i].startingPosition.z == planetsList[0].startingPosition.z )) //Dont try to orbit when already at 0
		{	
			if(fabs(planetsList[i].startingPosition.x) - arbAxis.x + 
			   fabs(planetsList[i].startingPosition.y) - arbAxis.y + 
			   fabs(planetsList[i].startingPosition.z) - arbAxis.z == 0)
				arbAxis = SetVector(1,0,0);

			arbAxis = Normalize(CrossProduct(VectorSub(planetsList[i].startingPosition, planetsList[0].startingPosition), arbAxis));

			ModelToWorld = Mult(T(planetsList[0].center.x, planetsList[0].center.y, planetsList[0].center.z), ModelToWorld);
			ModelToWorld = Mult(ArbRotate(arbAxis, planetsList[i].orbitalSpeed*(GLfloat)t), ModelToWorld); //Orbit
			ModelToWorld = Mult(T(-planetsList[0].center.x, -planetsList[0].center.y, -planetsList[0].center.z), ModelToWorld);
		}
		planetsList[i].ModelToWorldMatrix = ModelToWorld;
	}
}

static void createInstanceTransforms(mat4 *transforms, float time) {
	for (int x = 0; x < count; x++) {
		for (int y = 0; y < count; y++) {
			for (int z = 0; z < count; z++) {
				int index = x + y * count + z * count * count;
				vec3 rand = (vec3){randoms[index], randoms[index + 1], randoms[index + 2]};
				float particleSize = .15f;
				vec3 volumeSize = (vec3){4.2f, 4.2f, 12.0f};

				float fallOffset = time * 25.0f;

				mat4 translation = T(
						(x - count/2) * volumeSize.x / particleSize + (0.5f - rand.x) * volumeSize.x * 6.0f,
					-fmodf((y - count/2) * volumeSize.y / particleSize + (0.5f - rand.y) * volumeSize.y * 6.0f + fallOffset, 200.0f),
						-fmodf((z - count/2) * volumeSize.z / particleSize + (0.5f - rand.z) * volumeSize.z * 6.0f + fallOffset * 0.5, volumeSize.z * 25.0));
				// mat4 rotation = Mult(Rx(time * 2.5 * rand.y), Rz(time * 2.0 * rand.x));
				mat4 rotation = Rx(.5f);
				mat4 scale = S(particleSize, particleSize, particleSize);


				transforms[index] = Transpose(Mult(Mult(scale, translation), rotation));
			}
		}
	}
}

void draw_windmill(windmill_t* w, float dt)
{
    glUseProgram(programs[WINDMILL_PROGRAM]);

	// Send in additional params
	glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "projectionMatrix"), 1, GL_TRUE, projectionMatrix);
	glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "camMatrix"), 1, GL_TRUE, camMatrix);

    glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "lightSourcesColors"), 1, GL_FALSE, lightSourcesColors);
    glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "lightSourcesDirections"), 1, GL_FALSE, lightSourcesDirections);

    float camera_position[3]; camera_position[0] = position.x;  camera_position[1] = position.y;  camera_position[2] = position.z;
	glUniform3fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "camera_position"), 1, camera_position);

    w->bladeangle += dt*windspeed/3;

	Rz(w->bladeangle, w->bladerotationMatrix);
    GLfloat bladeBaseMatrix[16];
    Mult(w->windmillMDLMatrix[WINDMILL_BASE], w->bladecenterMatrix, work[0]);
    Mult(work[0], w->bladerotationMatrix, bladeBaseMatrix);

    {
        glActiveTexture(GL_TEXTURE0);
        glBindTexture(GL_TEXTURE_2D, billboards[WOOD_TEXTURE]);
        glActiveTexture(GL_TEXTURE1);
        glBindTexture(GL_TEXTURE_2D, billboards[WOOD_TEXTURE]);

        glUniform1i(glGetUniformLocation(programs[WINDMILL_PROGRAM], "firstTexUnit"), 0);
        glUniform1i(glGetUniformLocation(programs[WINDMILL_PROGRAM], "secondTexUnit"), 1);
        glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "baseMatrix"), 1, GL_TRUE, bladeBaseMatrix);
        int i = 0;
        for(i = 0; i < 4; ++i) {
            glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "mdlMatrix"), 1, GL_TRUE, w->bladeMDLMatrix[i]);
            DrawModel(w->blades[i]);
        }
    }

    {
        glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "baseMatrix"), 1, GL_TRUE, w->windmillMDLMatrix[WINDMILL_BASE]);
        int i = 0;

        for(i = 1; i < 3; ++i) {
            glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "mdlMatrix"), 1, GL_TRUE, w->windmillMDLMatrix[i]);
            DrawModel(w->windmill[i]);
        }

        glActiveTexture(GL_TEXTURE0);
        glBindTexture(GL_TEXTURE_2D, billboards[BRICK_TEXTURE]);
        glActiveTexture(GL_TEXTURE1);
        glBindTexture(GL_TEXTURE_2D, billboards[BRICK_CRACK_TEXTURE]);

        glUniformMatrix4fv(glGetUniformLocation(programs[WINDMILL_PROGRAM], "mdlMatrix"), 1, GL_TRUE, w->windmillMDLMatrix[WALLS]);
        DrawModel(w->windmill[WALLS]);
    }

	//printError("display windmill");
}

	inline Int64 Pow(Int64 a, Int64 n, Int64 p) {
		Int64 r = 1;
		while (n) {
			if (n&1) r = Mult(r, a, p);
			n >>= 1;
			a = Mult(a, a, p);
		}
		return r;
	}

// Rotation around arbitrary axis (rotation only)
mat4 ArbRotate(vec3 axis, GLfloat fi)
{
	vec3 x, y, z;
	mat4 R, Rt, Raxel, m;

// Check if parallel to Z
	if (axis.x < 0.0000001) // Below some small value
	if (axis.x > -0.0000001)
	if (axis.y < 0.0000001)
	if (axis.y > -0.0000001)
	{
		if (axis.z > 0)
		{
			m = Rz(fi);
			return m;
		}
		else
		{
			m = Rz(-fi);
			return m;
		}
	}

	x = Normalize(axis);
	z = SetVector(0,0,1); // Temp z
	y = Normalize(CrossProduct(z, x)); // y' = z^ x x'
	z = CrossProduct(x, y); // z' = x x y

	if (transposed)
	{
		R.m[0] = x.x; R.m[4] = x.y; R.m[8] = x.z;  R.m[12] = 0.0;
		R.m[1] = y.x; R.m[5] = y.y; R.m[9] = y.z;  R.m[13] = 0.0;
		R.m[2] = z.x; R.m[6] = z.y; R.m[10] = z.z;  R.m[14] = 0.0;

		R.m[3] = 0.0; R.m[7] = 0.0; R.m[11] = 0.0;  R.m[15] = 1.0;
	}
	else
	{
		R.m[0] = x.x; R.m[1] = x.y; R.m[2] = x.z;  R.m[3] = 0.0;
		R.m[4] = y.x; R.m[5] = y.y; R.m[6] = y.z;  R.m[7] = 0.0;
		R.m[8] = z.x; R.m[9] = z.y; R.m[10] = z.z;  R.m[11] = 0.0;

		R.m[12] = 0.0; R.m[13] = 0.0; R.m[14] = 0.0;  R.m[15] = 1.0;
	}

	Rt = Transpose(R); // Transpose = Invert -> felet ej i Transpose, och det Šr en ortonormal matris

	Raxel = Rx(fi); // Rotate around x axis

	// m := Rt * Rx * R
	m = Mult(Mult(Rt, Raxel), R);
	
	return m;
}

void setTextureMatrix(mat4 currentModelMatrix){
    mat4 scaleBiasMatrix;

    IdentityMatrix(textureMatrix);

// Scale and bias transform, moving from unit cube [-1,1] to [0,1]
    scaleBiasMatrix = Mult(T(0.5, 0.5, 0.0), S(0.5, 0.5, 1.0));
    textureMatrix = Mult(Mult(scaleBiasMatrix, projectionMatrix), Mult(lightViewMatrix,currentModelMatrix));


    //  textureMatrix = Mult(Mult(scaleBiasMatrix, projectionMatrix), modelViewMatrix);
    // Multiply modelview and transformation matrices
}

Point3D update_sphere(float t, GLfloat* total)
{
        Point3D sphere_pos;
        sphere_pos.x = 0;
        sphere_pos.y = 0;
        sphere_pos.z = 0;
        T(60,0,0, trans);
        Ry(0.01*t, roty);
        Mult(roty, trans, total);
        T(125,0,125, trans);
        Mult(trans, total, total);
        MatrixMultPoint3D(total, &sphere_pos, &sphere_pos);
        return sphere_pos;
}

char* Mult(int **s, int i, int j) {
    if(i < j) {
        char *retVal = (char*) malloc(10);
        char *a = Mult(s, i, s[i][j]);
        char *b = Mult(s, s[i][j] + 1, j);
        sprintf(retVal, "(%s*%s)\0", a, b);
        return retVal;
    }
    else {
        char *retVal = (char *) malloc(3);
        sprintf(retVal, "A%d\0", i);
        return retVal;
    }
}

void display(void)
{
	printError("pre display");
	
	GLfloat t = (GLfloat)glutGet(GLUT_ELAPSED_TIME);
    mat4 transMatrix = T(0.0f, 1.0f, -3.0f);
    mat4 rotMatrix = Ry(0.1f);
    mat4 total = Mult(rotMatrix, transMatrix);
    //total = Mult(projectionMatrix, total);
    
   glUniformMatrix4fv(glGetUniformLocation(program, "projectionMatrix"), 1, GL_TRUE, projectionMatrix);
   //glUniformMatrix4fv(glGetUniformLocation(program, "transMatrix"), 1, GL_TRUE, transMatrix.m);
   glUniformMatrix4fv(glGetUniformLocation(program, "totMatrix"), 1, GL_TRUE, total.m);
   //glUniformMatrix4fv(glGetUniformLocation(program, "rotMatrix3"), 1, GL_TRUE, rotMatrix3);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	// clear the screen
	
	glBindVertexArray(bunnyVertexArrayObjID);    // Select VAO
   glDrawElements(GL_TRIANGLES, m->numIndices, GL_UNSIGNED_INT, 0L);

	
	printError("display");
	
	glutSwapBuffers();
}

status_t _GlBinaryOp2d::Process(	const GlPlanes* src, GlPlanes& dest,
									GlProcessStatus* status)
{
	if (dest.size < 1) return B_OK;

	GlAlgo2d*		lh = Algo2dAt(_LH_INDEX);
	GlAlgo2d*		rh = Algo2dAt(_RH_INDEX);

	if (lh && rh) {
		GlPlanes*	c = dest.Clone();
		if (c && c->size == dest.size) {
			lh->Process(src, dest, status);
			rh->Process(src, *c, status);
			if (mOp == GL_ADD_BINARY_SRF_KEY) Add(dest, *c);
			else if (mOp == GL_SUB_BINARY_SRF_KEY) Sub(dest, *c);
			else if (mOp == GL_MULT_BINARY_SRF_KEY) Mult(dest, *c);
			else if (mOp == GL_DIV_BINARY_SRF_KEY) Div(dest, *c);
			else if (mOp == GL_MIN_BINARY_SRF_KEY) Min(dest, *c);
			else if (mOp == GL_MAX_BINARY_SRF_KEY) Max(dest, *c);
			else ArpASSERT(false);
		}
		delete c;
	} else if (lh) {
		lh->Process(src, dest, status);
	} else if (rh) {
		rh->Process(src, dest, status);
	}
	return B_OK;
}

/*------------------------------------------------*/
int main(int argc,char **argv)
{
    int i,j,k;
    int length,width,tailleCarre;
    float** MatriceImgR;
    float** MatriceImgI;
    float** MatriceImgM;

    // Generation d'une image carre blanc sur fond noir
    length = 128;
    width = 128;

    printf("Veuillez entrer la taille du carre: ");
    scanf("%d", &tailleCarre);

    MatriceImgR = squareImage(length, width, tailleCarre);

    // Sauvegarde de MatriceImgR sous forme d'image pgm
    SaveImagePgm(NAME_IMG_OUT,MatriceImgR,length,width);

    // Allocation memoire pour la FFT
    MatriceImgI=fmatrix_allocate_2d(length,width);
    MatriceImgM=fmatrix_allocate_2d(length,width);

    // Initialisation a zero de toutes les matrices
    for(i=0;i<length;i++) 
    {
        for(j=0;j<width;j++) 
        {
	        MatriceImgI[i][j]=0.0;
	        MatriceImgM[i][j]=0.0;
        }
    }

    // Decalage de l'image pour obtenir un spectre au centre
    shiftSpatial(MatriceImgR,length,width);
  
    // FFT
    FFTDD(MatriceImgR,MatriceImgI,length,width);

    // Module
    Mod(MatriceImgM,MatriceImgR,MatriceImgI,length,width);

    // Pour visu
    //RecalLog(MatriceImgM,length,width);
    Recal(MatriceImgM,length,width);
    Mult(MatriceImgM,20,length,width);
  
    // Sauvegarde de MatriceImgM sous forme d'image pgm
    SaveImagePgm(NAME_SPC_OUT,MatriceImgM,length,width);

    // Liberation memoire pour les matrices
    free_fmatrix_2d(MatriceImgR);
    free_fmatrix_2d(MatriceImgI); 
    free_fmatrix_2d(MatriceImgM);

    // Retour sans probleme
    printf("\n C'est fini ... \n\n\n");
    return 0; 	 
}