3_5.cpp

用VC编写的求解多元方程的算法

/*第5题	求解方程--源代码及关键源代码注解如下：*/
/*Basic theory:
	1. # of roots = highest power
	2. All rational roots will be factors of k / factors of a in the general equation of ax^n + bx^(n-1) + ... + cx + k
	3. Quadratic formula can solve for irrational / imaginary roots (i dont know cubic / quartic formulas)
	4. Irrational roots and imaginaries come in pairs (quadratic formula guarantees this)
	5. Given polynomial P(x) = ax^n + bx^(n-1) + ... + cx + k and given that q is a root, P(q) = 0*/

#include "stdio.h"
#include "conio.h"
#include "iostream.h"
#include "string.h"
#include "stdlib.h"
#include "ctype.h"
#include "math.h"
#include "process.h"
void quadraticequation(int,int,int,char*,char*);
void getterm(char*,char*,int&);
void getfactorsof(int*,int&,int);
void reduceequbyfactor(char*,char*,float);
void gettermforpower(char*,char*,int);
int getcoeff(char*);
int getdegreeofterm(char*);
int remaindertheorem(char*,float);
int exp(int,int);
float forallfactors(int*,int*,int,int,char*,char*);
float exp(float,int);
float calcequval(char*,float);
float getcoeff_f(char*);
char terms[6][40];

void main()
{
	char *equ,*equ2,roots[5][40];
	int loc=0,current_rootnum=0,a=0,b=0,c=0,numterms;
	int deg,f_a[1000],f_k[1000],num_a,num_k;
	float result;
	equ=(char*)malloc(5000);
	equ2=(char*)malloc(5000);

printf("Quadratic, cubic, most quartic & some quintic solver\nBy Albert Chaulk\n\nDos & Don'ts\n\tUse ^n for exponents\n\tDon't do x ^ n, put it like x^n\n\tDo not use x^1 or x^0. These really screw up the searches\n\tDon't use decimal or fractional coefficients\n\tCubics must have at least 1 rational root. Quartics, 2, and quintics, 3\n\tTerms must be listed in decending order (highest power to lowest)\n\tTerms may be skipped by not including them or using a coefficient of 0\n\n");
	printf("Equation: ");
	gets(equ);
	if(strlen(equ)==0){
		printf("No equation!\n");
		return;
	}

	for(int i=0;i<6&&equ[loc]!='='&&equ[loc]!='\0'&&equ[loc]!='\n';i++)
		getterm(equ,terms[i],loc);
	numterms=i;
	if((deg=getdegreeofterm(terms[0]))>5||deg<2){
		printf("Invalid term \"%s\" Exponent is outside range (2 to 5)\n",terms[0],deg);
		return;
	}
	switch(deg){//there are no breaks after 3-5 because after factoring, its 1 power smaller
	case 5:
		getfactorsof(f_a,num_a,getdegreeofterm(terms[0]));
		c=0;
		for(i=1;i<numterms;i++){
			if(getdegreeofterm(terms[i])==0){
				c=getcoeff(terms[i]);
			}
		}
		getfactorsof(f_k,num_k,c);
		result=forallfactors(f_a,f_k,num_a,num_k,equ,roots[current_rootnum]);
		if(result<0.0f&&result>-0.05f){
			printf("Error! Couldn't divide any factors\n");
			return;
		}
		current_rootnum++;
		reduceequbyfactor(equ,equ2,result);
		strcpy(equ,equ2);
		loc=0;
		for(i=0;i<6&&equ[loc]!='='&&equ[loc]!='\0'&&equ[loc]!='\n';i++)
			getterm(equ,terms[i],loc);
		numterms=i;
	case 4:
		getfactorsof(f_a,num_a,getdegreeofterm(terms[0]));
		c=0;
		for(i=1;i<numterms;i++){
			if(getdegreeofterm(terms[i])==0){
				c=getcoeff(terms[i]);
			}
		}
		getfactorsof(f_k,num_k,c);
		result=forallfactors(f_a,f_k,num_a,num_k,equ,roots[current_rootnum]);
		if(result<0.0f&&result>-0.05f){
			printf("Error! Couldn't divide any factors\n");
			return;
		}
		current_rootnum++;
		reduceequbyfactor(equ,equ2,result);
		strcpy(equ,equ2);
		loc=0;
		for(i=0;i<6&&equ[loc]!='='&&equ[loc]!='\0'&&equ[loc]!='\n';i++)
			getterm(equ,terms[i],loc);
		numterms=i;
	case 3:
		getfactorsof(f_a,num_a,getdegreeofterm(terms[0]));
		c=0;
		for(i=1;i<numterms;i++){
			if(getdegreeofterm(terms[i])==0){
				c=getcoeff(terms[i]);
			}
		}
		getfactorsof(f_k,num_k,c);
		if((result=forallfactors(f_a,f_k,num_a,num_k,equ,roots[current_rootnum]))==0.0f);
		if(result<0.0f&&result>-0.05f){
			printf("Error! Couldn't divide any factors\n");
			return;
		}
		current_rootnum++;
		reduceequbyfactor(equ,equ2,result);
		strcpy(equ,equ2);
		loc=0;
		for(i=0;i<6&&equ[loc]!='='&&equ[loc]!='\0'&&equ[loc]!='\n';i++)
			getterm(equ,terms[i],loc);
		numterms=i;
	case 2:
		b=c=0;
		for(i=1;i<numterms;i++)
			if(getdegreeofterm(terms[i])==1)
				b=getcoeff(terms[i]);
		for(i=1;i<numterms;i++)
			if(getdegreeofterm(terms[i])==0)
				c=getcoeff(terms[i]);
	     quadraticequation(getcoeff(terms[0]),b,c,roots[current_rootnum],roots[current_rootnum+1]);
		current_rootnum+=2;
		break;
	}
	printf("Roots: ");
	for(i=0;i<current_rootnum;i++){
		printf("%s",roots[i]);
		if(i+1<current_rootnum)
			printf(", ");
	}
	printf("\n");
	system("pause");
}
void quadraticequation(int a,int b,int c,char* r1,char* r2)
{
	// (-b +- sqrt( b*b - 4*a*c )) / 2*a
	if(b*b-4*a*c<0){
		sprintf(r1,"(%d+i*sqrt(%d))/%d",b*-1,(b*b-4*a*c)*-1,2*a);
		sprintf(r2,"(%d-i*sqrt(%d))/%d",b*-1,(b*b-4*a*c)*-1,2*a);
	}
	else{
		if(sqrt(b*b-4*a*c)-(int)sqrt(b*b-4*a*c)==0){
			sprintf(r1,"%.3f",(float)((float)-b+(float)sqrt(b*b-4*a*c))/(float)(2*a));
			sprintf(r2,"%.3f",(float)((float)-b-(float)sqrt(b*b-4*a*c))/(float)(2*a));
		}
		else{
			sprintf(r1,"(%d+sqrt(%d))/%d",b*-1,b*b-4*a*c,2*a);
			sprintf(r2,"(%d-sqrt(%d))/%d",b*-1,b*b-4*a*c,2*a);
		}
	}
}
void getterm(char *src,char *dest,int &loc)
{
	int loc2=0;
	if(src[loc]=='-'){
		dest[loc2]='-';
		loc2++;
	}
for(;src[loc+loc2]!='+'&&src[loc+loc2]!='-'&&src[loc+loc2]!='\n'&&src[loc+loc2]!='\0'&&src[loc+loc2]!='=';loc2++)
	dest[loc2]=src[loc+loc2];
	dest[loc2]='\0';
	loc+=loc2;
if(src[loc]!='-'&&src[loc]!='\0')
		loc++;//+1 passes terminating character, but a - sign is part of a term
}
void getfactorsof(int *factor_list,int &num,int val)
{
	int fcount=0;
	num=0;
	val=abs(val);
	for(int i=1;i<=val;i++)
		if((float)((float)((float)val/(float)i)-(int)((float)val/(float)i))==0.0f)
			num++;
//	factor_list=(int*)malloc(1000);
	for(i=1;i<=val;i++)
		if((float)((float)((float)val/(float)i)-(int)((float)val/(float)i))==0.0f){
			factor_list[fcount]=i;
			fcount++;
		}
}
int getcoeff(char *term)
{
	char *temp;
	int mult=1;
	temp=(char*)malloc(strlen(term));
	for(int i=0;term[i]==' ';i++);//bypass whitespace
	if(term[i]=='-'){
		mult=-1;
		i++;
	}
	for(;!isdigit(term[i])&&!isalpha(term[i]);i++);
	if(isalpha(term[i]))
		return mult;
	for(int j=0;isdigit(term[i+j]);j++)
		temp[j]=term[i+j];
	return atoi(temp)*mult;
}
float getcoeff_f(char *term)
{
	char *temp;
	int mult=1;
	if(term[0]=='\0')
		return 0.0f;
	temp=(char*)malloc(strlen(term));
	for(int i=0;term[i]==' ';i++);//bypass whitespace
	if(term[i]=='-'){
		mult=-1;
		i++;
	}
	for(;!isdigit(term[i])&&!isalpha(term[i]);i++);
	if(isalpha(term[i]))
		return (float)mult;
	for(int j=0;isdigit(term[i+j]);j++)
		temp[j]=term[i+j];
	return (float)atof(temp)*(float)mult;
}
int getdegreeofterm(char *term)
{
	if(strchr(term,'x')!=NULL){//variable term
		if(strchr(term,'^')!=NULL){//exponent term
			for(int i=0;term[i]!='^';i++);
//no error checking needed because ^ must be in the string
			i++;
			switch(term[i]){
			case '2':	return 2;
			case '3':	return 3;
			case '4':	return 4;
			case '5':	return 5;
			}
		}
		else
			return 1;
	}
	else
		return 0;
	return -1;
}
int remaindertheorem(char *equ,float root)
{//P(x) = ax^n + bx^(n-1) + ... + cx + k. If q is a root, P(q) = 0
	char terms[5][40];
	float val[5],totalval=0.0f;
	int loc=0;
	for(int i=0;i<6&&equ[loc]!='='&&equ[loc]!='\0'&&equ[loc]!='\n';i++)
		getterm(equ,terms[i],loc);
	int numterms=i;
	for(i=0;i<numterms;i++){
		val[i]=calcequval(terms[i],root);
		totalval+=val[i];
	}
	if(totalval*0.99f>-0.002f&&totalval*0.99f<0.002f)	return 0;
	else  return 1;
}
int exp(int val,int e)
{
	int v2=val;
	if(e==0)	return 1;
	for(int i=1;i<e;i++)	v2*=val;
	return v2;
}
float exp(float val,int e)
{
	float v2=val;
	if(e==0)	return 1.0f;
	for(int i=1;i<e;i++)	v2*=val;	
return v2;
}
float calcequval(char *equ,float val)
{
	int coeff=getcoeff(equ);
	int power=getdegreeofterm(equ);
	return coeff*exp(val,power);
}
void reduceequbyfactor(char *base,char *newequ,float factor)
{
	float synthterms[6]={0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};//synthetic division vars
	int newterms[5]={0,0,0,0,0};
	int j=0,olddegree;
char temp[50]="",exp[10]="",tempterms[6][40];
	if(remaindertheorem(base,factor))		return;
	olddegree=getdegreeofterm(terms[0]);
	for(int i=/*olddegree*/0;i<=olddegree;i++){
		gettermforpower(base,temp,i);
		synthterms[i]=getcoeff_f(temp);
	}
	newterms[4]=(int)synthterms[5];
	for(i=4;i>=1;i--)	newterms[i-1]=(int)((float)(newterms[i]*factor+synthterms[i]));
	for(i=4;newterms[i]==0;i--);
	strcpy(newequ,"");
	for(;i>=0;i--){
		if(newterms[i]!=0){
		switch(i){
		case 0:
			strcpy(exp,"");
			break;
		case 1:
			strcpy(exp,"x");
			break;
		case 2:
			strcpy(exp,"x^2");
			break;
		case 3:
			strcpy(exp,"x^3");
			break;
		case 4:
			strcpy(exp,"x^4");
			break;
		}
		if(newterms[i]!=1){
			if(newterms[i]<0)
				sprintf(tempterms[i],"%d%s",newterms[i],exp);
			else{
				if(i<olddegree-1){
					sprintf(tempterms[i],"+%d%s",newterms[i],exp);
				}
				else
					sprintf(tempterms[i],"%d%s",newterms[i],exp);
			}
		}
		else{
			if(newterms[i]<0)
				sprintf(tempterms[i],"%s",exp);
			else{
				if(i<olddegree-1){
					sprintf(tempterms[i],"+%s",exp);
				}
				else
					sprintf(tempterms[i],"%s",exp);
			}
		}
		if(newterms[i]==1&&exp[0]=='\0'){
			if(newterms[i]>0)
				sprintf(tempterms[i],"+%d",newterms[i]);
			else
				sprintf(tempterms[i],"%d",newterms[i]);
		}
		strcat(newequ,tempterms[i]);
		}
	}
}

void gettermforpower(char *equ,char *dest,int power)
{
	for(int i=0;i<6;i++)
		if(getdegreeofterm(terms[i])==power){
			strcpy(dest,terms[i]);
			return;
		}
//	if(i==6)
	strcpy(dest,"0");
}

float forallfactors(int *a,int *k,int n_a,int n_k,char *equ,char *root)
{
	for(int i=0;i<n_a;i++)
		for(int j=0;j<n_k;j++){
			if(remaindertheorem(equ,(float)((float)k[j]/(float)a[i]))==0){
				int r=k[j]/a[i];
				sprintf(root,"%d",r);
				return (float)r;
			}
			if(remaindertheorem(equ,(float)((float)k[j]/(float)a[i])*-1.0f)==0){
				int r=-k[j]/a[i];
				sprintf(root,"%d",r);
				return (float)r;
			}
		}

	if(remaindertheorem(equ,0.0f)==0){
		sprintf(root,"0");
		return 0.0f;
	}

	return -0.0001f;//value too small to be returned, used for errors
}