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();
        }
    }
}