cubicr.c

calc大数库

		X1 = INPUTI(&u);
		printf("Enter Y1: ");
		Y1 = INPUTI(&u);
		printf("Enter X2: ");
		X2 = INPUTI(&u);
		printf("Enter Y2: ");
		Y2 = INPUTI(&u);
		ADD_CUBICM(X1, Y1, X2, Y2, &X3, &Y3, A1, A2, A3, A4, A6, MODULUS);
		if (X3)
		{
			printf("(X3, Y3) = (");
			PRINTI(X3);
			printf(", ");
			PRINTI(Y3);
			printf(")\n");
			FREEMPI(X3);
			FREEMPI(Y3);
		}
		else
			printf("(X3, Y3) = IDENTITY\n");
		FREEMPI(X1);
		FREEMPI(Y1);
		FREEMPI(X2);
		FREEMPI(Y2);
	}
	FREEMPI(A1);
	FREEMPI(A2);
	FREEMPI(A3);
	FREEMPI(A4);
	FREEMPI(A6);
	return;
 
}

void POWER_CUBICM(MPI *X1, MPI *Y1, MPI **Xptr, MPI **Yptr, MPI *A1, MPI *A2, MPI *A3, MPI *A4, MPI *A6, unsigned int n, MPI *MODULUS)
/*
 * (*Xptr, *Yptr)= n(X1, Y1), where 0 <= n < R0 * R0 and (X1, Y1) is on the
 * elliptic curve y^2+A1*x*y+A3*y=X^3+A2*X^2+A4*x+A6 mod MODULUS.
 * See D. Husemoller, Elliptic curves, page 25.
 */
{
	MPI *BXtmp, *CXtmp, *BYtmp, *CYtmp;

	*Xptr = (MPI *)NULL;
	*Yptr = (MPI *)NULL;
	if (n == 0)
	{
		printf(" about to return\n");
		return;
	}
	BXtmp = COPYI(X1);
	BYtmp = COPYI(Y1);
	while (1)
	{
		if (n & 1)
		{
			CXtmp = *Xptr;
			CYtmp = *Yptr;
		ADD_CUBICM(BXtmp, BYtmp, CXtmp, CYtmp, Xptr, Yptr, A1, A2, A3, A4, A6, MODULUS);
			FREEMPI(CXtmp);
			FREEMPI(CYtmp);
			if (n == 1)
			{
				FREEMPI(BXtmp);
				FREEMPI(BYtmp);
				return;
			}
		}
		CXtmp = BXtmp;
		CYtmp = BYtmp;
		ADD_CUBICM(CXtmp, CYtmp, CXtmp, CYtmp, &BXtmp, &BYtmp, A1, A2, A3, A4, A6, MODULUS);
		FREEMPI(CXtmp);
		FREEMPI(CYtmp);
		n = n >> 1;
	}
}

void POWER_CUBICMX()
/*
 * Front end for POWER_CUBICM().
 */ 
{
	MPI  *X1, *X2, *Y1, *Y2, *A1, *A2, *A3, *A4, *A6;
	MPI *Tmp, *Tmp1,*Tmp2,*Tmp3,*Tmp4,*Tmp5, *MODULUS;
	MPI *TWO, *THREE, *FOUR, *NINE, *DELTA, *B2, *B4, *B6, *B8;
	unsigned int u, n;
	
	printf("Calculating n(X1,Y1) for the nonsingular elliptic curve \n");
	printf("y^2+A1*x*y+A3*y=X^3+A2*X^2+A4*x+A6 mod MODULUS.\n");
	printf("Enter A1: ");
	A1 = INPUTI(&u);
	printf("Enter A2: ");
	A2 = INPUTI(&u);
	printf("Enter A3: ");
	A3 = INPUTI(&u);
	printf("Enter A4: ");
	A4 = INPUTI(&u);
	printf("Enter A6: ");
	A6 = INPUTI(&u);
	printf("Enter MODULUS: ");
	MODULUS = INPUTI(&u);
 
	TWO = CHANGE(2);
	THREE = CHANGE(3);
	FOUR = MULTM(TWO, TWO, MODULUS);
	NINE = MULTM(THREE, THREE, MODULUS);

	Tmp = MULTM(FOUR, A2, MODULUS);
	B2 = MULTABCM(Tmp, A1, A1, MODULUS);
	FREEMPI(Tmp);

	Tmp = MULTM(TWO, A4, MODULUS);
	B4 = MULTABCM(Tmp, A1, A3, MODULUS);
	FREEMPI(Tmp);

	Tmp = MULTM(FOUR, A6, MODULUS);
	B6 = MULTABCM(Tmp, A3, A3, MODULUS);
	FREEMPI(Tmp);

	Tmp1 = MULTM(B2, B6, MODULUS);
	Tmp2 = MULTM(B4, B4, MODULUS);
	Tmp3 = SUBM(Tmp1, Tmp2, MODULUS);
	B8 = DIVM(Tmp3, FOUR, MODULUS);
	FREEMPI(Tmp1);
	FREEMPI(Tmp2);
	FREEMPI(Tmp3);

	/* Finally DELTA = -B2^2*B8 -8*B4^3 - 27*B6^2 + 9*B2*B4*B6.  */

	Tmp1 = MULTM3(B2, B2, B8, MODULUS);
	Tmp2 = MULTM(TWO, B4, MODULUS);
	Tmp3 = MULTM3(Tmp2, Tmp2, Tmp2, MODULUS);
	Tmp = Tmp3;
	Tmp3 = ADDM(Tmp1, Tmp3, MODULUS);
	FREEMPI(Tmp);
	FREEMPI(Tmp1);
	FREEMPI(Tmp2);

	Tmp2 = MULTM(B2, B4, MODULUS);
	Tmp4 = MULTM(THREE, B6, MODULUS);
	Tmp5 = SUBM(Tmp2, Tmp4, MODULUS);
	FREEMPI(Tmp2);
	FREEMPI(Tmp4);
	Tmp2 = MULTM3(NINE, B6, Tmp5, MODULUS);
	DELTA = SUBM(Tmp2, Tmp3, MODULUS);
	printf("DELTA = "); PRINTI(DELTA); printf("\n");
	FREEMPI(Tmp2);
	FREEMPI(Tmp3);
	FREEMPI(Tmp5);
	FREEMPI(TWO);
	FREEMPI(THREE);
	FREEMPI(FOUR);
	FREEMPI(NINE);
	FREEMPI(B2);
	FREEMPI(B4);
	FREEMPI(B6);
	FREEMPI(B8);

	if (EQZEROI(DELTA))
	{
		fprintf(stderr, "Discriminant is zero.\n");
		exit (1);
	}
	FREEMPI(DELTA);

	while (GetYN())
	{
		printf("Enter X1: ");
		X1 = INPUTI(&u);
		printf("Enter Y1: ");
		Y1 = INPUTI(&u);
		printf("Enter n: ");
		scanf("%u", &n);
		printf("n = %u\n", n);
		POWER_CUBICM(X1, Y1, &X2, &Y2, A1, A2, A3, A4, A6, n, MODULUS);
		if (X2)
		{
			printf("(X2, Y2) = (");
			PRINTI(X2);
			printf(", ");
			PRINTI(Y2);
			printf(")\n");
			FREEMPI(X2);
			FREEMPI(Y2);
		}
		else
			printf("(X2, Y2) = IDENTITY\n");
		FREEMPI(X1);
		FREEMPI(Y1);
	}
	FREEMPI(A1);
	FREEMPI(A2);
	FREEMPI(A3);
	FREEMPI(A4);
	FREEMPI(A6);
	FREEMPI(MODULUS);
	return;
}

MPI *MULTABCM(MPI *A, MPI *B, MPI *C, MPI *M)
/* Returns A + B * C. */
{
	MPI *Temp1, *Temp2;

	Temp1 = MULTM(B,C, M);
	Temp2 = ADDM(A, Temp1, M);
	FREEMPI(Temp1);
	return (Temp2);
}

MPI *MULTM3(MPI *A, MPI *B, MPI *C, MPI *M)
/* Returns A * B * C mod M. */
{
	MPI *Tmp1, *Tmp2;

	Tmp1 = MULTM(A, B, M);
	Tmp2 = MULTM(Tmp1, C, M);
	FREEMPI(Tmp1);
	return (Tmp2);
}


unsigned int ORDER_CUBICM(MPI *X1, MPI *Y1, MPI *A1, MPI *A2, MPI *A3, MPI *A4, MPI *A6, MPI *MODULUS)
/* Returns order of (X1,Y1) on the elliptic curve 
 * y^2+A1*x*y+A3*y=X^3+A2*X^2+A4*x+A6 mod MODULUS.
 */
{
	unsigned int k;
	MPI *X2, *Y2, *X3, *Y3;
	k = 1;
	if (X1 == NULL)
		return 1;
	else
		X2= COPYI(X1);
		Y2= COPYI(Y1);
	while (1){
		ADD_CUBICM(X1, Y1, X2, Y2, &X3, &Y3, A1, A2, A3, A4, A6, MODULUS);
		FREEMPI(X2);
		FREEMPI(Y2);
			k++;
		if (X3 == NULL)
			break;
		else
		{
			X2 = X3;
			Y2 = Y3;
		}
	}
	return (k);
}

/*unsigned int ORDER_CUBICM(MPI *X1, MPI *Y1, MPI *A1, MPI *A2, MPI *A3, MPI *A4, MPI *A6, MPI *MODULUS)*/
/* Returns order of (X1,Y1) on the elliptic curve 
 * y^2+A1*x*y+A3*y=X^3+A2*X^2+A4*x+A6 mod MODULUS.
 */
/*{
	unsigned int k;
	MPI *X2, *Y2;
	k = 1;
	if (X1 == NULL)
		return 1;
	while (1){
		POWER_CUBICM(X1, Y1, &X2, &Y2, A1, A2, A3, A4, A6, k, MODULUS);
		if (X2 == NULL)
			break;
		else
		{
			FREEMPI(X2);
			FREEMPI(Y2);
			k++;
		}
	}
	return (k);
}
*/
void ORDER_CUBICMX()
/*
 * Front end for ORDER_CUBICM().
 */ 
{
	MPI  *X1, *Y1, *A1, *A2, *A3, *A4, *A6;
	MPI *Tmp, *Tmp1,*Tmp2,*Tmp3,*Tmp4,*Tmp5, *MODULUS;
	MPI *TWO, *THREE, *FOUR, *NINE, *DELTA, *B2, *B4, *B6, *B8;
	unsigned int k, u;
	
	printf("Calculating o(X1,Y1) for the nonsingular elliptic curve \n");
	printf("y^2+A1*xy+A3*y=X^3+A2*X^2+A4*x+A6 modulo a prime.\n");
	printf("Enter A1: ");
	A1 = INPUTI(&u);
	printf("Enter A2: ");
	A2 = INPUTI(&u);
	printf("Enter A3: ");
	A3 = INPUTI(&u);
	printf("Enter A4: ");
	A4 = INPUTI(&u);
	printf("Enter A6: ");
	A6 = INPUTI(&u);
	printf("Enter MODULUS: ");
	MODULUS = INPUTI(&u);
 
	TWO = CHANGE(2);
	THREE = CHANGE(3);
	FOUR = MULTM(TWO, TWO, MODULUS);
	NINE = MULTM(THREE, THREE, MODULUS);

	Tmp = MULTM(FOUR, A2, MODULUS);
	B2 = MULTABCM(Tmp, A1, A1, MODULUS);
	FREEMPI(Tmp);

	Tmp = MULTM(TWO, A4, MODULUS);
	B4 = MULTABCM(Tmp, A1, A3, MODULUS);
	FREEMPI(Tmp);

	Tmp = MULTM(FOUR, A6, MODULUS);
	B6 = MULTABCM(Tmp, A3, A3, MODULUS);
	FREEMPI(Tmp);

	Tmp1 = MULTM(B2, B6, MODULUS);
	Tmp2 = MULTM(B4, B4, MODULUS);
	Tmp3 = SUBM(Tmp1, Tmp2, MODULUS);
	B8 = DIVM(Tmp3, FOUR, MODULUS);
	FREEMPI(Tmp1);
	FREEMPI(Tmp2);
	FREEMPI(Tmp3);

	/* Finally DELTA = -B2^2*B8 -8*B4^3 - 27*B6^2 + 9*B2*B4*B6.  */

	Tmp1 = MULTM3(B2, B2, B8, MODULUS);
	Tmp2 = MULTM(TWO, B4, MODULUS);
	Tmp3 = MULTM3(Tmp2, Tmp2, Tmp2, MODULUS);
	Tmp = Tmp3;
	Tmp3 = ADDM(Tmp1, Tmp3, MODULUS);
	FREEMPI(Tmp);
	FREEMPI(Tmp1);
	FREEMPI(Tmp2);

	Tmp2 = MULTM(B2, B4, MODULUS);
	Tmp4 = MULTM(THREE, B6, MODULUS);
	Tmp5 = SUBM(Tmp2, Tmp4, MODULUS);
	FREEMPI(Tmp2);
	FREEMPI(Tmp4);
	Tmp2 = MULTM3(NINE, B6, Tmp5, MODULUS);
	DELTA = SUBM(Tmp2, Tmp3, MODULUS);
	printf("DELTA = "); PRINTI(DELTA); printf("\n");
	FREEMPI(Tmp2);
	FREEMPI(Tmp3);
	FREEMPI(Tmp5);
	FREEMPI(TWO);
	FREEMPI(THREE);
	FREEMPI(FOUR);
	FREEMPI(NINE);
	FREEMPI(B2);
	FREEMPI(B4);
	FREEMPI(B6);
	FREEMPI(B8);

	if (EQZEROI(DELTA))
	{
		fprintf(stderr, "Discriminant is zero.\n");
		exit (1);
	}
	FREEMPI(DELTA);

	while (GetYN())
	{
		printf("Enter X1: ");
		X1 = INPUTI(&u);
		printf("Enter Y1: ");
		Y1 = INPUTI(&u);
		k=ORDER_CUBICM(X1, Y1, A1, A2, A3, A4, A6, MODULUS);
		printf("Order of (");
		PRINTI(X1); 
		printf(",");
		PRINTI(Y1); 
		printf(") on the curve\n");
		printf("y^2");
		if (A1->S)
		{
			printf("+");
			if (!EQONEI(A1))
				PRINTI(A1);
			printf("xy");
		}
		if (A3->S)
		{
			printf("+");
			if (!EQONEI(A3))
				PRINTI(A3);
			printf("y");
		}
		printf("=x^3");
		if (A2->S > 0)
		{
			printf("+");
			if (!EQONEI(A2))
				PRINTI(A2);
			printf("x^2");
		}
		if (A4->S > 0)
		{
			printf("+");
			if (!EQONEI(A4))
				PRINTI(A4);
			printf("x");
		}
		if (A6->S > 0)
		{
			printf("+");
			PRINTI(A6);
		}
		printf(" mod ");
		PRINTI(MODULUS);
		printf("\n");
		printf(" is %u\n ", k);
		FREEMPI(X1);
		FREEMPI(Y1);
	}
	FREEMPI(A1);
	FREEMPI(A2);
	FREEMPI(A3);
	FREEMPI(A4);
	FREEMPI(A6);
	FREEMPI(MODULUS);
	return;
}

unsigned int ORDER_CUBICR(MPR *X1, MPR *Y1, MPR *A1, MPR *A2, MPR *A3, MPR *A4, MPR *A6)
/* Returns order of (X1,Y1) on the elliptic curve 
 * y^2+A1*x*y+A3*y=X^3+A2*X^2+A4*x+A6. Uses fact that |tor(E)|<= 12.
 */
{
	unsigned int k;
	MPR *X2, *Y2;
	k = 1;
	if (X1 == NULL)
		return 1;
	while (k <= 12){
		POWER_CUBICR(X1, Y1, &X2, &Y2, A1, A2, A3, A4, A6, k);
		if (X2 == NULL)
			break;
		else
		{
			FREEMPR(X2);
			FREEMPR(Y2);
			k++;
		}
	}
	if (k == 13)
		return (0);
	else
		return (k);
}

void ORDER_CUBICRX()
/*
 * Front end for ORDER_CUBICR().
 */ 
{
	MPR  *X1, *Y1, *A1, *A2, *A3, *A4, *A6;
	MPR *Tmp, *Tmp1,*Tmp2,*Tmp3,*Tmp4,*Tmp5;
	MPR *TWO, *THREE, *FOUR, *NINE, *DELTA, *B2, *B4, *B6, *B8;
	unsigned int k, u;
	
	printf("Calculating o(X1,Y1) for the nonsingular elliptic curve \n");
	printf("y^2+A1*xy+A3y=X^3+A2*X^2+A4*x+A6.\n");
	printf("Enter A1: ");
	A1 = INPUTR(&u);
	printf("Enter A2: ");
	A2 = INPUTR(&u);
	printf("Enter A3: ");
	A3 = INPUTR(&u);
	printf("Enter A4: ");
	A4 = INPUTR(&u);
	printf("Enter A6: ");
	A6 = INPUTR(&u);
	TWO = BUILDMPR();
	TWO->N = CHANGE(2);
	TWO->D = ONEI();
	THREE = BUILDMPR();
	THREE->N = CHANGE(3);
	THREE->D = ONEI();
	FOUR = MULTR(TWO, TWO);
	NINE = MULTR(THREE, THREE);

	Tmp = MULTR(FOUR, A2);
	B2 = MULTABCR(Tmp, A1, A1);
	FREEMPR(Tmp);

	Tmp = MULTR(TWO, A4);
	B4 = MULTABCR(Tmp, A1, A3);
	FREEMPR(Tmp);

	Tmp = MULTR(FOUR, A6);
	B6 = MULTABCR(Tmp, A3, A3);
	FREEMPR(Tmp);

	Tmp1 = MULTR(B2, B6);
	Tmp2 = MULTR(B4, B4);
	Tmp3 = SUBR(Tmp1, Tmp2);
	B8 = RATIOR(Tmp3, FOUR);
	FREEMPR(Tmp1);
	FREEMPR(Tmp2);
	FREEMPR(Tmp3);

	/* Finally DELTA = -B2^2*B8 -8*B4^3 - 27*B6^2 + 9*B2*B4*B6.  */

	Tmp1 = MULTR3(B2, B2, B8);
	Tmp2 = MULTR(TWO, B4);
	Tmp3 = MULTR3(Tmp2, Tmp2, Tmp2);
	Tmp = Tmp3;
	Tmp3 = ADDR(Tmp1, Tmp3);
	FREEMPR(Tmp);
	FREEMPR(Tmp1);
	FREEMPR(Tmp2);

	Tmp2 = MULTR(B2, B4);
	Tmp4 = MULTR(THREE, B6);
	Tmp5 = SUBR(Tmp2, Tmp4);
	FREEMPR(Tmp2);
	FREEMPR(Tmp4);
	Tmp2 = MULTR3(NINE, B6, Tmp5);
	DELTA = SUBR(Tmp2, Tmp3);
	printf("DELTA = "); PRINTR(DELTA); printf("\n");
	FREEMPR(Tmp2);
	FREEMPR(Tmp3);
	FREEMPR(Tmp5);
	FREEMPR(TWO);
	FREEMPR(THREE);
	FREEMPR(FOUR);
	FREEMPR(NINE);
	FREEMPR(B2);
	FREEMPR(B4);
	FREEMPR(B6);
	FREEMPR(B8);

	if (EQZEROR(DELTA))
	{
		fprintf(stderr, "Discriminant is zero.\n");
		exit (1);
	}
	FREEMPR(DELTA);

	while (GetYN())
	{
		printf("Enter X1: ");
		X1 = INPUTR(&u);
		printf("Enter Y1: ");
		Y1 = INPUTR(&u);
		k=ORDER_CUBICR(X1, Y1, A1, A2, A3, A4, A6);
		printf("Order of (");
		PRINTR(X1); 
		printf(",");
		PRINTR(Y1); 
		printf(") on the curve\n");
		printf("y^2");
		if ((A1->N)->S > 0)
		{
			printf("+");
			if (!EQONER(A1))
				PRINTR(A1);
			printf("xy");
		}
		else if ((A1->N)->S < 0)
		{
			if (!EQMINUSONER(A1))
				PRINTR(A1);
			else
				printf("-");
			printf("xy");
		}
		if ((A3->N)->S > 0)
		{
			printf("+");
			if (!EQONER(A3))
				PRINTR(A3);
			printf("y");
		}
		else if ((A3->N)->S < 0)
		{
			if (!EQMINUSONER(A3))
				PRINTR(A3);
			else
				printf("-");
			printf("y");
		}
		printf("=x^3");
		if ((A2->N)->S > 0)
		{
			printf("+");
			if (!EQONER(A2))
				PRINTR(A2);
			printf("x^2");
		}
		else if ((A2->N)->S < 0)
		{
			if (!EQMINUSONER(A2))
				PRINTR(A2);
			else
				printf("-");
			printf("x^2");
		}
		if ((A4->N)->S > 0)
		{
			printf("+");
			if (!EQONER(A4))
				PRINTR(A4);
			printf("x");
		}
		else if ((A4->N)->S < 0)
		{
			if (!EQMINUSONER(A4))
				PRINTR(A4);
			else
				printf("-");
			printf("x");
		}
		if ((A6->N)->S > 0)
		{
			printf("+");
			PRINTR(A6);
		}
		else if ((A6->N)->S < 0)
			PRINTR(A6);
		printf("\n");
		if (k)
			printf(" is %u\n ", k);
		else 
			printf(" is infinite\n ");
		FREEMPR(X1);
		FREEMPR(Y1);
	}
	FREEMPR(A1);
	FREEMPR(A2);
	FREEMPR(A3);
	FREEMPR(A4);
	FREEMPR(A6);
	return;
}