C#，数值计算——Globals的计算方法与源程序

 

1 文本格式

using System;
using System.Text;

namespace Legalsoft.Truffer
{
    public static partial class Globals
    {
        //const int FLT_RADIX = 2;
        //const int DBL_MANT_DIG = 53;
        //const int INT_DIGITS = 32;
        //const float FLT_EPSILON = 1.19209290E-07F;
        //const double DBL_EPSILON = 2.2204460492503131E-16;

        public static double SQR(double a)
        {
            return a * a;
        }

        /// <summary>
        /// 浮点数取余数的函数
        /// https://blog.csdn.net/Hanford/article/details/53633937
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        public static double fmod(double x, double y)
        {
            y = Math.Abs(y);
            if (x >= 0.0)
            {
                y = x - y * Math.Floor(x / y);
            }
            else
            {
                y = x - y * Math.Ceiling(x / y);
            }
            return y;
        }

        public static double SIGN(double a, double b)
        {
            return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
        }

        public static void SWAP(ref int a, ref int b)
        {
            (a, b) = (b, a);
        }

        public static void SWAP(ref double a, ref double b)
        {
            (a, b) = (b, a);
        }

        /// <summary>
        /// 数组复制(int)
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public static int[] CopyFrom(int[] b)
        {
            int[] v = new int[b.Length];
            for (int i = 0; i < v.Length; i++)
            {
                v[i] = b[i];
            }
            return v;
        }

        /// <summary>
        /// 数组复制(double)
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public static double[] CopyFrom(double[] b)
        {
            double[] v = new double[b.Length];
            for (int i = 0; i < v.Length; i++)
            {
                v[i] = b[i];
            }
            return v;
        }

        /// <summary>
        /// 数组转（方）矩阵(int)
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public static double[,] CopyFrom(int row, int col, double[] b)
        {
            double[,] v = new double[row, col];
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                {
                    v[i, j] = b[i * col + j];
                }
            }
            return v;
        }

        /// <summary>
        /// 矩阵复制(double)
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public static double[,] CopyFrom(double[,] b)
        {
            int nn = b.GetLength(0);
            int mm = b.GetLength(1);
            double[,] v = new double[nn, mm];
            for (int i = 0; i < nn; i++)
            {
                for (int j = 0; j < mm; j++)
                {
                    v[i, j] = b[i, j];
                }
            }
            return v;
        }

        /// <summary>
        /// 复制矩阵b的第k行
        /// </summary>
        /// <param name="k"></param>
        /// <param name="b"></param>
        /// <returns>数组</returns>
        public static double[] CopyFrom(int k, double[,] b)
        {
            int mm = b.GetLength(1);
            double[] v = new double[mm];
            for (int j = 0; j < mm; j++)
            {
                v[j] = b[k, j];
            }
            return v;
        }

        /// <summary>
        /// 提取三元矩阵的第k层的矩阵
        /// </summary>
        /// <param name="k"></param>
        /// <param name="b"></param>
        /// <returns>矩阵</returns>
        public static double[,] CopyFrom(int k, double[,,] b)
        {
            int nn = b.GetLength(1);
            int mm = b.GetLength(2);
            double[,] v = new double[nn, mm];
            for (int i = 0; i < nn; i++)
            {
                for (int j = 0; j < mm; j++)
                {
                    v[i, j] = b[k, i, j];
                }
            }
            return v;
        }

        public static double dist(double[] p1, double[] p2)
        {
            double sum = 0.0;
            for (int i = 0; i < p1.Length; i++)
            {
                sum += Globals.SQR(p1[i] - p2[i]);
            }
            if (sum <= float.Epsilon)
            {
                return 0.0;
            }
            return Math.Sqrt(sum);
        }

        public static double ldexp(double v, int exp)
        {
            return v * Math.Pow(2.0, exp);
        }

        public static double frexp(double x, out int exp)
        {
            if (Math.Abs(x) <= float.Epsilon)
            {
                exp = 0;
                return 0.0;
            }
            long bits = BitConverter.DoubleToInt64Bits(x);
            ulong u1 = 0x800fffffffffffffL;
            ulong u2 = (ulong)bits;
            long r = (long)(u1 & u2);
            double d = BitConverter.Int64BitsToDouble(r | 0x3fe0000000000000L);
            exp = (int)((0x7ff0000000000000L & bits) >> 52) - 1022;
            return d;
        }

        public static double gammln(double xx)
        {
            double[] cof = { 57.1562356658629235, -59.5979603554754912, 14.1360979747417471, -0.491913816097620199, .339946499848118887e-4, .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3, -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3, .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5 };
            if (xx <= 0)
            {
                throw new Exception("bad arg in gammln");
            }
            double x = xx;
            double y = xx;
            double tmp = x + 5.24218750000000000;
            tmp = (x + 0.5) * Math.Log(tmp) - tmp;
            double ser = 0.999999999999997092;
            for (int j = 0; j < 14; j++)
            {
                ser += cof[j] / ++y;
            }
            return tmp + Math.Log(2.5066282746310005 * ser / x);
        }

        private static double[] factrl_a;
        private static bool factrl_init = true;

        public static double factrl(int n)
        {
            if (factrl_init)
            {
                factrl_init = false;
                factrl_a = new double[171];
                factrl_a[0] = 1.0;
                for (int i = 1; i < 171; i++)
                {
                    factrl_a[i] = i * factrl_a[i - 1];
                }
            }
            if (n < 0 || n > 170)
            {
                throw new Exception("factrl out of range");
            }
            return factrl_a[n];
        }

        private static double[] factln_a = new double[2000];
        private static bool factln_init = true;

        public static double factln(int n)
        {
            const int NTOP = 2000;
            if (factln_init)
            {
                factln_init = false;
                for (int i = 0; i < NTOP; i++)
                {
                    factln_a[i] = gammln(i + 1.0);
                }
            }
            if (n < 0)
            {
                throw new Exception("negative arg in factln");
            }
            if (n < NTOP)
            {
                return factln_a[n];
            }
            return gammln(n + 1.0);
        }

        public static double bico(int n, int k)
        {
            if (n < 0 || k < 0 || k > n)
            {
                throw new Exception("bad args in bico");
            }
            if (n < 171)
            {
                return Math.Floor(0.5 + factrl(n) / (factrl(k) * factrl(n - k)));
            }
            return Math.Floor(0.5 + Math.Exp(factln(n) - factln(k) - factln(n - k)));
        }

        public static double beta(double z, double w)
        {
            return Math.Exp(gammln(z) + gammln(w) - gammln(z + w));
        }

        public static string ToString(double[] x)
        {
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < x.Length; i++)
            {
                sb.AppendFormat("{0:F12},", x[i]);
            }
            sb.AppendLine("<br>");
            return sb.ToString();
        }

        public static string ToString(double[,] x)
        {
            StringBuilder sb = new StringBuilder();
            sb.AppendLine("<style>td { padding:5px; } </style>");            
            sb.AppendLine("<table border=1 bordercolor='#888888' style='border-collapse:collapse;'>");
            for (int i = 0; i < x.GetLength(0); i++)
            {
                sb.AppendLine("<tr>");
                for (int j = 0; j < x.GetLength(1); j++)
                {
                    sb.AppendFormat("<td>{0:F12}</td>", x[i, j]);
                }
                sb.AppendLine("</tr>");
            }
            sb.AppendLine("</table>");
            sb.AppendLine("<br>");
            return sb.ToString();
        }
    }
}
 

2 代码格式

using System;
using System.Text;namespace Legalsoft.Truffer
{public static partial class Globals{//const int FLT_RADIX = 2;//const int DBL_MANT_DIG = 53;//const int INT_DIGITS = 32;//const float FLT_EPSILON = 1.19209290E-07F;//const double DBL_EPSILON = 2.2204460492503131E-16;public static double SQR(double a){return a * a;}/// <summary>/// 浮点数取余数的函数/// https://blog.csdn.net/Hanford/article/details/53633937/// </summary>/// <param name="x"></param>/// <param name="y"></param>/// <returns></returns>public static double fmod(double x, double y){y = Math.Abs(y);if (x >= 0.0){y = x - y * Math.Floor(x / y);}else{y = x - y * Math.Ceiling(x / y);}return y;}public static double SIGN(double a, double b){return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}public static void SWAP(ref int a, ref int b){(a, b) = (b, a);}public static void SWAP(ref double a, ref double b){(a, b) = (b, a);}/// <summary>/// 数组复制(int)/// </summary>/// <param name="b"></param>/// <returns></returns>public static int[] CopyFrom(int[] b){int[] v = new int[b.Length];for (int i = 0; i < v.Length; i++){v[i] = b[i];}return v;}/// <summary>/// 数组复制(double)/// </summary>/// <param name="b"></param>/// <returns></returns>public static double[] CopyFrom(double[] b){double[] v = new double[b.Length];for (int i = 0; i < v.Length; i++){v[i] = b[i];}return v;}/// <summary>/// 数组转（方）矩阵(int)/// </summary>/// <param name="b"></param>/// <returns></returns>public static double[,] CopyFrom(int row, int col, double[] b){double[,] v = new double[row, col];for (int i = 0; i < row; i++){for (int j = 0; j < col; j++){v[i, j] = b[i * col + j];}}return v;}/// <summary>/// 矩阵复制(double)/// </summary>/// <param name="b"></param>/// <returns></returns>public static double[,] CopyFrom(double[,] b){int nn = b.GetLength(0);int mm = b.GetLength(1);double[,] v = new double[nn, mm];for (int i = 0; i < nn; i++){for (int j = 0; j < mm; j++){v[i, j] = b[i, j];}}return v;}/// <summary>/// 复制矩阵b的第k行/// </summary>/// <param name="k"></param>/// <param name="b"></param>/// <returns>数组</returns>public static double[] CopyFrom(int k, double[,] b){int mm = b.GetLength(1);double[] v = new double[mm];for (int j = 0; j < mm; j++){v[j] = b[k, j];}return v;}/// <summary>/// 提取三元矩阵的第k层的矩阵/// </summary>/// <param name="k"></param>/// <param name="b"></param>/// <returns>矩阵</returns>public static double[,] CopyFrom(int k, double[,,] b){int nn = b.GetLength(1);int mm = b.GetLength(2);double[,] v = new double[nn, mm];for (int i = 0; i < nn; i++){for (int j = 0; j < mm; j++){v[i, j] = b[k, i, j];}}return v;}public static double dist(double[] p1, double[] p2){double sum = 0.0;for (int i = 0; i < p1.Length; i++){sum += Globals.SQR(p1[i] - p2[i]);}if (sum <= float.Epsilon){return 0.0;}return Math.Sqrt(sum);}public static double ldexp(double v, int exp){return v * Math.Pow(2.0, exp);}public static double frexp(double x, out int exp){if (Math.Abs(x) <= float.Epsilon){exp = 0;return 0.0;}long bits = BitConverter.DoubleToInt64Bits(x);ulong u1 = 0x800fffffffffffffL;ulong u2 = (ulong)bits;long r = (long)(u1 & u2);double d = BitConverter.Int64BitsToDouble(r | 0x3fe0000000000000L);exp = (int)((0x7ff0000000000000L & bits) >> 52) - 1022;return d;}public static double gammln(double xx){double[] cof = { 57.1562356658629235, -59.5979603554754912, 14.1360979747417471, -0.491913816097620199, .339946499848118887e-4, .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3, -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3, .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5 };if (xx <= 0){throw new Exception("bad arg in gammln");}double x = xx;double y = xx;double tmp = x + 5.24218750000000000;tmp = (x + 0.5) * Math.Log(tmp) - tmp;double ser = 0.999999999999997092;for (int j = 0; j < 14; j++){ser += cof[j] / ++y;}return tmp + Math.Log(2.5066282746310005 * ser / x);}private static double[] factrl_a;private static bool factrl_init = true;public static double factrl(int n){if (factrl_init){factrl_init = false;factrl_a = new double[171];factrl_a[0] = 1.0;for (int i = 1; i < 171; i++){factrl_a[i] = i * factrl_a[i - 1];}}if (n < 0 || n > 170){throw new Exception("factrl out of range");}return factrl_a[n];}private static double[] factln_a = new double[2000];private static bool factln_init = true;public static double factln(int n){const int NTOP = 2000;if (factln_init){factln_init = false;for (int i = 0; i < NTOP; i++){factln_a[i] = gammln(i + 1.0);}}if (n < 0){throw new Exception("negative arg in factln");}if (n < NTOP){return factln_a[n];}return gammln(n + 1.0);}public static double bico(int n, int k){if (n < 0 || k < 0 || k > n){throw new Exception("bad args in bico");}if (n < 171){return Math.Floor(0.5 + factrl(n) / (factrl(k) * factrl(n - k)));}return Math.Floor(0.5 + Math.Exp(factln(n) - factln(k) - factln(n - k)));}public static double beta(double z, double w){return Math.Exp(gammln(z) + gammln(w) - gammln(z + w));}public static string ToString(double[] x){StringBuilder sb = new StringBuilder();for (int i = 0; i < x.Length; i++){sb.AppendFormat("{0:F12},", x[i]);}sb.AppendLine("<br>");return sb.ToString();}public static string ToString(double[,] x){StringBuilder sb = new StringBuilder();sb.AppendLine("<style>td { padding:5px; } </style>");            sb.AppendLine("<table border=1 bordercolor='#888888' style='border-collapse:collapse;'>");for (int i = 0; i < x.GetLength(0); i++){sb.AppendLine("<tr>");for (int j = 0; j < x.GetLength(1); j++){sb.AppendFormat("<td>{0:F12}</td>", x[i, j]);}sb.AppendLine("</tr>");}sb.AppendLine("</table>");sb.AppendLine("<br>");return sb.ToString();}}
}