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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/124955.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

基于SpringBoot的小区物业管理系统

基于SpringBoot的小区物业管理系统的设计与实现 开发语言&#xff1a;Java数据库&#xff1a;MySQL技术&#xff1a;SpringBootMyBatis工具&#xff1a;IDEA/Ecilpse、Navicat、Maven 系统展示 首页 管理员界面 摘要 基于SpringBoot的小区物业管理系统是一款为小区物业管理提…

进程调度算法之先来先服务(FCFS),短作业优先(SJF)以及高响应比优先(HRRN)

1.先来先服务&#xff08;FCFS&#xff09; first come first service 1.算法思想 主要从“公平”的角度考虑(类似于我们生活中排队买东西的例子) 2.算法规则 按照作业/进程到达的先后顺序进行服务。 3.用于作业/进程调度 用于作业调度时&#xff0c;考虑的是哪个作业先…

plt 画图不显示label

没写 plt.legend() 这个 ! # 效果模拟-------------- import matplotlib.pyplot as plt import matplotlib as mpl # matplotlib其实是不支持显示中文的 显示中文需要一行代码设置字体 mpl.rcParams[font.family] = STKAITI # STKAITI——字体 plt.rcParams[axes.unicode_m…

微信小程序点单左右联动的效果实现

微信小程序点单左右联动的效果实现 原理解析&#xff1a;   点击左边标签会跳到右边相应位置&#xff1a;点击改变rightCur值&#xff0c;转跳相应位置滑动右边&#xff0c;左边标签会跳到相应的位置&#xff1a;监听并且设置每个右边元素的top和bottom&#xff0c;再判断当…

Visopsys 0.92 发布

Visopsys 是一个 PC 机的操作系统&#xff0c;系统小型、快速而且开源。有着丰富的图形界面、抢先式多任务机制以及支持虚拟内存。Visopsys 视图兼容很多操作系统&#xff0c;但并不是他们的克隆版本。Visopsys 0.92 现已发布&#xff0c;此维护版本引入了多任务处理程序、文件…

MySQL 多表关联查询优化实践和原理解析

目录 一、前言二、表数据准备三、表关联查询原理和两种算法3.1、研究关联查询算法必备知识点3.2、嵌套循环连接 Nested-Loop Join(NLJ) 算法3.3、基于块的嵌套循环连接 Block Nested-Loop Join(BNL)算法3.4、被驱动表的关联字段没索引为什么要选择使用 BNL 算法而不使用 Nested…

Spring Cloud Loadbalancer 实现客户端负载均衡

针对 ribbon 负载均衡组件&#xff0c; 官方提出的替换解决方案是 Spring Cloud Loadbalancer。本次主要通过学习示例介绍了 Spring Cloud Loadbalancer 的基础使用。 1&#xff0c;引入pom <dependency><groupId>org.springframework.cloud</groupId><…

Pikachu靶场——文件包含漏洞(File Inclusion)

文章目录 1. File Inclusion1.2 File Inclusion(local)1.2.1 源代码分析1.2.2 漏洞防御 1.3 File Inclusion(remote)1.3.1 源代码分析1.3.2 漏洞防御 1.4 文件包含漏洞防御 1. File Inclusion 还可以参考我的另一篇文章&#xff1a;文件包含漏洞及漏洞复现。 File Inclusion(…

Cortex-A9 架构

一、Cortex-A 处理器运行模式 Cortex-A9处理器有 9中处理模式&#xff0c;如下表所示&#xff1a; 九种运行模式 在上表中&#xff0c;除了User(USR)用户模式以外&#xff0c;其它8种运行模式都是特权模式&#xff0c;在特权模式下&#xff0c;程序可以访问所有的系统资源。这…

2023年中国半导体IP行业发展概况及趋势分析:半导体IP的市场空间广阔[图]

半导体指IP指芯片设计中预先没计、验证好的功能模块&#xff0c;处于半导体产业链最上游&#xff0c;为芯片设计厂商提供设计模块。半导体IP按交付方式可分为软核、硬核和固核&#xff1b;按产品类型可分为处理器IP、接口IP、其他物理IP及其他数字IP。 半导体IP分类 资料来源&…

CompletableFuture异步回调

CompletableFuture异步回调 CompletableFutureFuture模式CompletableFuture详解1.CompletableFuture的UML类关系2.CompletionStage接口3.使用runAsync和supplyAcync创建子任务4.设置子任务回调钩子5.调用handle()方法统一处理异常和结果6.线程池的使用 异步任务的串行执行thenA…

二叉树题目:路径总和 II

文章目录 题目标题和出处难度题目描述要求示例数据范围 前言解法一思路和算法代码复杂度分析 解法二思路和算法代码复杂度分析 题目 标题和出处 标题&#xff1a;路径总和 II 出处&#xff1a;113. 路径总和 II 难度 4 级 题目描述 要求 给你二叉树的根结点 root \tex…