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