C#，数值计算——Sobol拟随机序列的计算方法与源程序

1 文本格式

using System;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Sobol quasi-random sequence
    /// </summary>
    public class Sobol
    {

        public Sobol() { }

        public static void sobseq(int n, double[] x)
        {
            const int MAXBIT = 30;
            const int MAXDIM = 6;

            uint xin = 0;
            int[] mdeg = new int[] { 1, 2, 3, 3, 4, 4 };
            uint[] ix = new uint[MAXDIM];
            List<uint[]> iu = new List<uint[]>(MAXBIT);
            uint[] ip = new uint[] {
                0, 1, 1, 2, 1, 4
            };
            uint[] iv = new uint[] {
                1, 1, 1, 1, 1, 1,
                3, 1, 3, 3, 1, 1,
                5, 7, 7, 3, 3, 5,
                15, 11, 5, 15, 13, 9
            };
            double fac = 0.0;

            if (n < 0)
            {
                for (int k = 0; k < MAXDIM; k++)
                {
                    ix[k] = 0;
                }
                xin = 0;
                if (iv[0] != 1)
                {
                    return;
                }
                fac = 1.0 / (1 << MAXBIT);
                for (int j = 0, k = 0; j < MAXBIT; j++, k += MAXDIM)
                {
                    iu[j] = new uint[iv[k]];
                }
                for (int k = 0; k < MAXDIM; k++)
                {
                    for (int j = 0; j < mdeg[k]; j++)
                    {
                        iu[j][k] <<= (MAXBIT - 1 - j);
                    }
                    for (int j = mdeg[k]; j < MAXBIT; j++)
                    {
                        uint ipp = ip[k];
                        uint i = iu[j - mdeg[k]][k];
                        i ^= (i >> mdeg[k]);
                        for (int l = mdeg[k] - 1; l >= 1; l--)
                        {
                            if ((ipp & 1) != 0)
                            {
                                i ^= iu[j - l][k];
                            }
                            ipp >>= 1;
                        }
                        iu[j][k] = i;
                    }
                }
            }
            else
            {
                uint im = xin++;
                int j = 0;
                for (; j < MAXBIT; j++)
                {
                    if ((im & 1) == 0)
                    {
                        break;
                    }
                    im >>= 1;
                }
                if (j >= MAXBIT)
                {
                    throw new Exception("MAXBIT too small in sobseq");
                }
                im = (uint)(j * MAXDIM);
                for (int k = 0; k < Math.Min(n, MAXDIM); k++)
                {
                    ix[k] ^= iv[im + k];
                    x[k] = ix[k] * fac;
                }
            }
        }

    }
}
 

2 代码格式

using System;
using System.Collections.Generic;namespace Legalsoft.Truffer
{/// <summary>/// Sobol quasi-random sequence/// </summary>public class Sobol{public Sobol() { }public static void sobseq(int n, double[] x){const int MAXBIT = 30;const int MAXDIM = 6;uint xin = 0;int[] mdeg = new int[] { 1, 2, 3, 3, 4, 4 };uint[] ix = new uint[MAXDIM];List<uint[]> iu = new List<uint[]>(MAXBIT);uint[] ip = new uint[] {0, 1, 1, 2, 1, 4};uint[] iv = new uint[] {1, 1, 1, 1, 1, 1,3, 1, 3, 3, 1, 1,5, 7, 7, 3, 3, 5,15, 11, 5, 15, 13, 9};double fac = 0.0;if (n < 0){for (int k = 0; k < MAXDIM; k++){ix[k] = 0;}xin = 0;if (iv[0] != 1){return;}fac = 1.0 / (1 << MAXBIT);for (int j = 0, k = 0; j < MAXBIT; j++, k += MAXDIM){iu[j] = new uint[iv[k]];}for (int k = 0; k < MAXDIM; k++){for (int j = 0; j < mdeg[k]; j++){iu[j][k] <<= (MAXBIT - 1 - j);}for (int j = mdeg[k]; j < MAXBIT; j++){uint ipp = ip[k];uint i = iu[j - mdeg[k]][k];i ^= (i >> mdeg[k]);for (int l = mdeg[k] - 1; l >= 1; l--){if ((ipp & 1) != 0){i ^= iu[j - l][k];}ipp >>= 1;}iu[j][k] = i;}}}else{uint im = xin++;int j = 0;for (; j < MAXBIT; j++){if ((im & 1) == 0){break;}im >>= 1;}if (j >= MAXBIT){throw new Exception("MAXBIT too small in sobseq");}im = (uint)(j * MAXDIM);for (int k = 0; k < Math.Min(n, MAXDIM); k++){ix[k] ^= iv[im + k];x[k] = ix[k] * fac;}}}}
}