文档章节

用盛金公式解三次方程(ansi c版)

wangxuwei
 wangxuwei
发布于 09/22 00:37
字数 1121
阅读 12
收藏 0

ansi风格: 

/*
   cc cubic.c -lm
   gcc cubic.c -lm

Shengjin's Formulas Univariate cubic equation aX ^ 3 + bX ^ 2 + cX + d = 0, (a, b, c, d < R, and a!= 0).
       Multiple root discriminant: delta1 = b^2-3*a*c; delta2 = b*c-9*a*d; delta3 = c^2-3*b*d,
           The total discriminant is delta=delta2^2-4*delta1*delta3.
           When delta1 = delta2 = 0, Shengjin Formula (1): X1=X2=X3=-b/(3*a)=-c/b=-3d/c.
           When delta=B^2-4*A*C>0, Shengjin Formula II:
                 Y1= delta1*b + 3*a *((-B + (delta)^1/2))/ 2.
                 Y2= delta1*b + 3*a *((-B - (delta)^1/2))/ 2.
                 x1=(-b-Y1^(1/3) - Y1^(1/3)/(3*a);
                 X2=(-2*b+Y1^(1/3)+Y2^(1/3)/(6*a)+3^(1/2)* (Y1^(1/3)-Y2^(1/3)/(6a)i,
                 X3=(-2*b+Y1^(1/3)+Y2^(1/3)/(6*a)-3^(1/2)* (Y1^(1/3)-Y2^(1/3)/(6a)i,

           When delta=B^2-4AC=0, Shengjin Formula 3:
                 X1=-b/a+K; X2=X3=-K/2,            K = delta2/delta1, (A<>0).
           When delta=B^2-4AC<0, Shengjin Formula 4:
                 X1= (-b-2*sqrt(delta1)*cos(theta/3))/(3*a);
                 X2= (-b+sqrt(delta1)*(cos(theta/3)+sqrt(3)*sin(theta/3)))/(3*a);
                 X3= (-b+sqrt(delta1)*(cos(theta/3)-sqrt(3)*sin(theta/3)))/(3*a)
                 theta=arccosT,T=(2Ab-3aB)/(2A^(3/2))
    Shengjin's Distinguishing Means
        (1)A = B = 0, the equation has a triple real root.
        (2) When delta =B^2-4AC>0, the equation has a real root and a pair of conjugate complex roots.
        (3) When delta=B^2-4AC=0, the equation has three real roots, one of which has two double roots.
        (4) When delta=B^2-4AC<0, the equation has three unequal real roots.
*/

#include <stdio.h>
#include <math.h>
enum roottype{UniReal,OneRPairComplex,TwoReal,UnequalReal,NotCubic};

int main(){
    const double PI=3.14159265359;
    double a,b,c,d;    /*Coefficient of cubic Equation*/
    double delta1,delta2,delta3,delta;
    double Y1,Y2,expY1,expY2;
    double K,theta,T;
    enum roottype rt;

    printf("Input Coefficient of cubic Equation:a b c d.\n");
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);

    delta1 = b*b-3*a*c;
    delta2 = b*c-9*a*d;
    delta3 = c*c-3*b*d;
    delta  = delta2*delta2-4*delta1*delta3;

    if(delta1==0&&delta2==0) rt=UniReal;
    else if(delta >0 )  rt=OneRPairComplex;
    else if(delta ==0)  rt=TwoReal;
    else if(delta <0 )  rt=UnequalReal;
    if (a==0) rt=NotCubic;
    switch(rt){
      case UniReal:
         printf("the equation has a triple real root.\n");
         printf("x1=x2=x3=%f\n",-b/3/a);
         break;
      case OneRPairComplex:
         printf("the equation has a real root and a pair of conjugate complex roots.\n");
         Y1 = delta1*b + 3*a *((-delta2 + sqrt(delta)))/2;
         Y2 = delta1*b + 3*a *((-delta2 - sqrt(delta)))/2;

         if (Y1>0)  expY1 = pow(Y1,1.0/3.0); else expY1 = (-1)*pow(fabs(Y1),1.0/3.0);
         if (Y2>0)  expY2 = pow(Y2,1.0/3.0); else expY2 = (-1)*pow(fabs(Y2),1.0/3.0);
         printf("x1=%f\n",(-b-expY1-expY2)/(3*a));
         printf("x2=%f+%fi\n",(-2*b+expY1+expY2)/(6*a),sqrt(3.0)* (expY1-expY2)/(6*a));
         printf("x3=%f-%fi\n",(-2*b+expY1+expY2)/(6*a),sqrt(3.0)* (expY1-expY2)/(6*a));
         break;
      case TwoReal:
         printf("the equation has three real roots, one of which has two double roots.\n");
         K =delta2/delta1;
         printf("X1=X2=%f\n",-K/2);
         printf("X3=%f\n",-b/a+K);
         break;
      case UnequalReal:
         printf("the equation has three unequal real roots.\n");
         T=(2.0*delta1*b-3.0*a*delta2)/(2.0*sqrt(delta1*delta1*delta1));
         theta=acos(T);
         printf("X1=%f\n",(-b-2.0*sqrt(delta1)*cos(theta/3.0))/(3.0*a));
         printf("X2=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)+sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         printf("X3=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)-sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         break;
      case NotCubic:
         printf("the equation is  not a cubic.\n");
         break;
    } /* end of switch*/

   return(0);
}

K&R风格

#include <stdio.h>
#include <math.h>
#define PI 3.14159265359

enum roottype{UniReal,OneRPairComplex,TwoReal,UnequalReal,NotCubic};

main(){
    double a,b,c,d;    /*Coefficient of cubic Equation*/
    double delta1,delta2,delta3,delta;
    double Y1,Y2,expY1,expY2;
    double K,theta,T;
    enum roottype rt;

    printf("Input Coefficient of cubic Equation:a b c d.\n");
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);

    delta1 = b*b-3*a*c;
    delta2 = b*c-9*a*d;
    delta3 = c*c-3*b*d;
    delta  = delta2*delta2-4*delta1*delta3;

    if(delta1==0&&delta2==0) rt=UniReal;
    else if(delta >0 )  rt=OneRPairComplex;
    else if(delta ==0)  rt=TwoReal;
    else if(delta <0 )  rt=UnequalReal;
    if (a==0) rt=NotCubic;

    switch(rt){
      case UniReal:
         printf("the equation has a triple real root.\n");
         printf("x1=x2=x3=%f\n",-b/3/a);
         break;
      case OneRPairComplex:
         printf("the equation has a real root and a pair of conjugate complex roots.\n");
         Y1 = delta1*b + 3*a *((-delta2 + sqrt(delta)))/2;
         Y2 = delta1*b + 3*a *((-delta2 - sqrt(delta)))/2;

         if (Y1>0)  expY1 = pow(Y1,1.0/3.0); else expY1 = (-1)*pow(fabs(Y1),1.0/3.0);
         if (Y2>0)  expY2 = pow(Y2,1.0/3.0); else expY2 = (-1)*pow(fabs(Y2),1.0/3.0);
         /*printf("expY1=%f,expY2=%f,x3=%f\n",expY1,expY2,a);*/		 
         printf("x1=%f\n",(-b-expY1-expY2)/(3*a));
         printf("x2=%f+%fi\n",(-2*b+expY1+expY2)/(6.0*a),sqrt(3.0)* (expY1-expY2)/(6.0*a));
         printf("x3=%f-%fi\n",(-2*b+expY1+expY2)/(6.0*a),sqrt(3.0)* (expY1-expY2)/(6.0*a));
         break;
      case TwoReal:
         printf("the equation has three real roots, one of which has two double roots.\n");
         K =delta2/delta1;
         printf("X1=X2=%f\n",-K/2);
         printf("X3=%f\n",-b/a+K);
         break;
      case UnequalReal:
         printf("the equation has three unequal real roots.\n");
         T=(2.0*delta1*b-3.0*a*delta2)/(2.0*sqrt(delta1*delta1*delta1));
         theta=acos(T);
         printf("X1=%f\n",(-b-2.0*sqrt(delta1)*cos(theta/3.0))/(3.0*a));
         printf("X2=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)+sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         printf("X3=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)-sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         break;
      case NotCubic:
         printf("the equation is  not a cubic.\n");
         break;
    } /* end of switch*/

}

参考

matlab实现一元三次方程求解

© 著作权归作者所有

wangxuwei
粉丝 27
博文 343
码字总数 137316
作品 0
杭州
其他
私信 提问
10071 - Back to High School Physics

解题思路: 1. 两个基本物理公式 V = V0 + at (V0为初速度, a 为加速度, t 为时间, V 为当前速度) S = (V + V0) * t / 2 (V0为初速度, V 为末速度, t 为时间) 2. 由题目可知, 输入的是 t 和 ...

sailtseng
2013/06/04
216
0
Eran/BezierMathUtils

#BezierMathUtils 介绍 BezierMathUtils是用来提供计算贝塞尔曲线时候需要用到的数学公式。 比如存在如下需求: 有两点坐标A和B,希望某一个物品从A点移动到B点. 移动方式希望采用贝塞尔曲线...

Eran
2015/02/18
0
0
合取范式的可满足性判定算法和谓词逻辑不可判定性

作为本系列的最后一篇文章,我们来看被广为研究的SAT问题。 SAT问题是第一个被证明为NP问题的判定问题。更多信息可以去百度或者维基一下。 前面我们看到了Horn公式可满足性的判定算法,现在把...

davelet
2014/05/11
0
0
学Objective-C后如何进行IOS软件开发呢

这两天也开始学习了objective-c, 因为自己有C基础(应该说还算扎实)同时有一点java经验. objective-c虽然是C的衍生但感觉还是差别很大. 通过这两天对xcode4.4.1的初步了解后发现其界面开发能力...

fcsong000833
2013/01/24
375
1
error C2143: syntax error : missing ';' before ...

今天偶然写了下面的程序(原来我写的程序不一样,下面的只是为了把问题简化) void foo() { int p = 0; if ( p == 0 ) { int i = 0; } int a; } int main() { foo(); } 不幸的是偶然将这个文件...

dalu~
2012/12/06
297
0

没有更多内容

加载失败,请刷新页面

加载更多

iOS Xcode升级包地址(感谢大神)

下载地址:DeviceSupport

_____1____
10分钟前
4
0
Qt编写自定义控件71-圆弧进度条

一、前言 现在web形式的图表框架非常流行,国产代表就是echart,本人用过几次,三个字屌爆了来形容,非常强大,而且易用性也非常棒,还是开源免费的,使用起来不要太爽,内置的各种图表和仪表...

飞扬青云
11分钟前
3
0
润乾报表与 ActiveReport JS 功能对比

简介 润乾报表是用于报表制作的大型企业级报表软件,核心特点在于开创性地提出了非线性报表数学模型,采用了革命性的多源关联分片、不规则分组、自由格间运算、行列对称等技术,使得复杂报表...

泡泡糖儿
22分钟前
4
0
【1015】LNMP架构二

【1015】LNMP架构二 三、PHP安装 PHP安装和LAMP安装PHP方法有差别,需要开启php-fpm服务 1、下载PHP7至/usr/local/src/ 切换目录:cd /usr/local/src 2、解压缩 tar -jxvf php-7.3.0.tar.gz...

飞翔的竹蜻蜓
56分钟前
5
0
浅谈Visitor访问者模式

一、前言 什么叫访问,如果大家学过数据结构,对于这点就很清晰了,遍历就是访问的一般形式,单独读取一个元素进行相应的处理也叫作访问,读取到想要查看的内容+对其进行处理就叫作访问,那么...

青衣霓裳
今天
6
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部