[心得] C# 代數處理的類

[Lucemia (生の直感、死の予感)]

今天臨時需要用到解方程式，牛頓法
找不到C#的類別。
python的 sympy又跑太慢了,
只好自己硬幹了一個能跑牛頓法的。

請問C#有這類的lib嗎?

// 用法
// 定義變數
Furmula x = new Furmula("x");
Furmula y = new Furmula("y");
Furmula z = new Furmula("z");

// 定義方程式
Furmula f = x*x + y*y+z*z + x*y*z;

// 變數代入
double v =f.subs(x,3).subs(y,2).subs(z,1).Value;

// 微分
Furmula df = f.Diff(x);



//  -- source code --

 class Furmula
    {
        // Furmula
        public String mName = "";
        public bool mIsValue = false;
        public double mValue = 0;

        //
        public OPERATORS mOperator = OPERATORS.NULL;

        public Furmula mLeftChild = null;
        public Furmula mRightChild = null;

        public enum OPERATORS { ADD, SUB, MUL, DIV, SIN, COS, TAN, POW, NULL };

        public Furmula(String name)
        {
            mName = name;
        }

        public Furmula(double value)
        {
            mValue = value;
            mIsValue = true;
        }

        public Furmula(Furmula b)
        {
            mName = b.mName;
            mIsValue = b.mIsValue;
            mValue = b.mValue;

            mOperator = b.mOperator;
            mLeftChild = b.mLeftChild;
            mRightChild = b.mRightChild;
        }

        public Furmula(OPERATORS oper, Furmula LeftChild, Furmula RightChild)
        {
            mOperator = oper;
            mLeftChild = LeftChild;
            mRightChild = RightChild;
        }

        public static Furmula operator +(Furmula a, Furmula b)
        {
            if (a.mIsValue && b.mIsValue)
                return a.mValue + b.mValue;
            else if (a == 0)
                return new Furmula(b);
            else if (b == 0)
                return new Furmula(a);


            return new Furmula(OPERATORS.ADD, a, b);
        }

        public static Furmula operator -(Furmula a, Furmula b)
        {
            if (a.mIsValue && b.mIsValue)
                return a.mValue - b.mValue;
            else if (a == 0)
                return -1 * new Furmula(b);
            else if (b == 0)
                return new Furmula(a);

            return new Furmula(OPERATORS.SUB, a, b);
        }

        public static Furmula operator *(Furmula a, Furmula b)
        {
            if (a.mIsValue && b.mIsValue)
                return a.mValue * b.mValue;
            else if (a == 0)
                return 0;
            else if (b == 0)
                return 0;
            else if (a == 1)
                return new Furmula(b);
            else if (b == 1)
                return new Furmula(a);

            if (a.mIsValue && b.mIsValue)
                return a.mValue * b.mValue;

            return new Furmula(OPERATORS.MUL, a, b);
        }

        public static Furmula operator /(Furmula a, Furmula b)
        {
            if (a.mIsValue && b.mIsValue)
                return a.mValue / b.mValue;
            else if (b == 0)
                throw new Exception();
            else if (a == 0)
                return 0;
            else if (b == 1)
                return new Furmula(a);

            return new Furmula(OPERATORS.DIV, a, b);
        }
        public static Furmula Pow(Furmula a, Furmula b)
        {
            if (a.mIsValue && b.mIsValue)
                return Math.Pow(a.mValue, b.mValue);
            else if (b == 0)
                return 1;
            else if (a == 0)
                return 0;
            else if (a == 1)
                return 1;

            return new Furmula(OPERATORS.POW, a, b);
        }

        public static Furmula Sin(Furmula a)
        {
            if (a.mIsValue)
                return Math.Sin(a.mValue);

            return new Furmula(OPERATORS.SIN, a, null);
        }
        public static Furmula Cos(Furmula a)
        {
            if (a.mIsValue)
                return Math.Cos(a.mValue);
            return new Furmula(OPERATORS.COS, a, null);
        }
        public static Furmula Tan(Furmula a)
        {
            if (a.mIsValue)
                return Math.Tan(a.mValue);
            return new Furmula(OPERATORS.TAN, a, null);
        }
        public static bool operator ==(Furmula a, double b)
        {
            if (a.mIsValue && a.mValue == b)
                return true;
            return false;
        }
        public static bool operator !=(Furmula a, double b)
        {
            return !(a == b);
        }

        public static implicit operator Furmula(double v)
        {
            return new Furmula(v);
        }

        public double Value
        {
            get
            {
                if (mIsValue)
                {
                    return mValue;
                }
                throw new Exception();
            }

        }

        public Furmula Subs(Dictionary<Furmula, double> values)
        {
            if (mOperator == OPERATORS.NULL)
            {
                foreach (var v in values)
                {
                    if (mName == v.Key.mName)
                    {
                        return new Furmula(v.Value);
                    }
                }
                return new Furmula(this);
            }
            else
            {
                Furmula NewLeft = null, NewRight = null;
                if (mLeftChild != null)
                {
                    NewLeft = mLeftChild.Subs(values);
                }
                if (mRightChild != null)
                {
                   NewRight= mRightChild.Subs(values);
                }
                switch (mOperator)
                {
                    case OPERATORS.ADD:
                        return NewLeft + NewRight;
                    case OPERATORS.SUB:
                        return NewLeft - NewRight;
                    case OPERATORS.MUL:
                        return NewLeft * NewRight;
                    case OPERATORS.DIV:
                        return NewLeft / NewRight;
                    case OPERATORS.SIN:
                        return Sin(NewLeft);
                    case OPERATORS.COS:
                        return Cos(NewLeft);
                    case OPERATORS.TAN:
                        return Tan(NewLeft);
                    case OPERATORS.POW:
                        return Pow(NewLeft, NewRight);
                }
                throw new Exception();
            }
        }

        public Furmula Subs(String name, double value)
        {
            if (mOperator == OPERATORS.NULL)
            {
                if (mName == name)
                {
                    return new Furmula(value);
                }
                else
                {
                    return new Furmula(this);
                }
            }
            else
            {
                Furmula NewLeft = null, NewRight = null;
                if (mLeftChild != null)
                {
                   NewLeft= mLeftChild.Subs(name, value);
                }
                if (mRightChild != null)
                {
                    NewRight= mRightChild.Subs(name, value);
                }
                switch (mOperator)
                {
                    case OPERATORS.ADD:
                        return NewLeft + NewRight;
                    case OPERATORS.SUB:
                        return NewLeft - NewRight;
                    case OPERATORS.MUL:
                        return NewLeft * NewRight;
                    case OPERATORS.DIV:
                        return NewLeft / NewRight;
                    case OPERATORS.SIN:
                        return Sin(NewLeft);
                    case OPERATORS.COS:
                        return Cos(NewLeft);
                    case OPERATORS.TAN:
                        return Tan(NewLeft);
                    case OPERATORS.POW:
                        return Pow(NewLeft, NewRight);
                }
                throw new Exception();
            }
        }
        public Furmula Diff(Furmula furmula)
        {
            return Diff(furmula.mName);
        }

        public Furmula Diff(String name)
        {
            if (mOperator == OPERATORS.NULL)
            {
                if (mName == name)
                {
                    return new Furmula(1);
                }
                else
                {
                    return new Furmula(0);
                }
            }
            else
            {
                switch (mOperator)
                {
                    case OPERATORS.ADD:
                        return mLeftChild.Diff(name) + mRightChild.Diff(name);
                    case OPERATORS.SUB:
                        return mLeftChild.Diff(name) - mRightChild.Diff(name);
                    case OPERATORS.MUL:
                        return mLeftChild.Diff(name) * mRightChild + mLeftChild * mRightChild.Diff(name);
                    case OPERATORS.DIV:
                        return ((mLeftChild.Diff(name) * mRightChild) - (mLeftChild * mRightChild.Diff(name))) / Pow(mRightChild, 2);
                    case OPERATORS.SIN:
                        return Cos(mLeftChild) * mLeftChild.Diff(name);
                    case OPERATORS.COS:
                        return -1 * Sin(mLeftChild) * mLeftChild.Diff(name);
                    case OPERATORS.TAN:
                        return mLeftChild.Diff(name) / Pow(Cos(mLeftChild), 2);
                    case OPERATORS.POW:
                        return mRightChild * Pow(mLeftChild, mRightChild - 1) * mLeftChild.Diff(name);
                }
            }
            throw new Exception();
        }
        public override string ToString()
        {
            switch(mOperator){
                case OPERATORS.NULL:
                    if (mIsValue)
                        return mValue.ToString();
                    else
                        return mName;
                case OPERATORS.ADD:
                    return String.Format("({0}+{1})", mLeftChild, mRightChild);
                case OPERATORS.SUB:
                    return String.Format("({0}-{1})", mLeftChild, mRightChild);
                case OPERATORS.MUL:
                    return String.Format("({0}*{1})", mLeftChild, mRightChild);
                case OPERATORS.DIV:
                    return String.Format("({0}/{1})", mLeftChild, mRightChild);
                case OPERATORS.POW:
                    return String.Format("({0}**{1})", mLeftChild, mRightChild);
                case OPERATORS.SIN:
                    return String.Format("Sin({0})", mLeftChild);
                case OPERATORS.COS:
                    return String.Format("Cos({0})", mLeftChild);
                case OPERATORS.TAN:
                    return String.Format("Tan({0})", mLeftChild);

            }
            throw new Exception();
        }

        public static Dictionary<Furmula, double> Newton(Furmula f, Dictionary<Furmula, double> Guass)
        {
            List<Furmula> vars = new List<Furmula>(Guass.Keys);
            int K = vars.Count;

            Furmula[] Gradient = new Furmula[K];
            Furmula[,] Hessian = new Furmula[K, K];

            for (int i = 0; i < K; i++)
                Gradient[i] = f.Diff(vars[i]);

            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < K; j++)
                {
                    Hessian[i, j] = f.Diff(vars[i]).Diff(vars[j]);
                }
            }

            CSML.Matrix Xn = new CSML.Matrix(K, 1);
            CSML.Matrix Fn = new CSML.Matrix(K, 1);
            CSML.Matrix Hn = new CSML.Matrix(K, K);
            System.IO.StreamWriter writer = new System.IO.StreamWriter("iter.log");
            for (int iter = 0; iter < 100; iter++)
            {
                writer.WriteLine("{0} {1}", iter, f.Subs(Guass));
                writer.Flush();

                for (int i = 0; i < K; i++)
                {
                    Xn[i + 1, 1] = new Complex(Guass[vars[i]]);
                }

                for (int i = 0; i < K; i++)
                {
                    Fn[i + 1, 1] = new Complex(Gradient[i].Subs(Guass).Value);
                }

                for (int i = 0; i < K; i++)
                {
                    for (int j = 0; j < K; j++)
                    {
                        Hn[i + 1, j + 1] = new Complex(Hessian[i, j].Subs(Guass).Value);
                    }
                }
                Xn = Xn - Hn.Inverse() * Fn;

                for (int i = 0; i < K; i++)
                {
                    Guass[vars[i]] = Xn[i + 1, 1].Re;
                }
            }
            return Guass;
        }
}

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 154.20.172.64
推 Steven0422:IMSL for .Net
推 pyrochlore:Open Source的有Math.Net
→ Lucemia:感謝