import numpy as npfrom numba import cuda# 定义CUDA Kernel函数@cuda.jitdef vector_add_kernel(a, b, c):    """    向量加法CUDA Kernel    参数:        a: 输入数组A        b: 输入数组B          c: 输出数组C = A + B    """    # 计算当前线程的全局索引    idx = cuda.grid(1)    # 边界检查    if idx < c.size:        c[idx] = a[idx] + b[idx]# 主机端调用代码def vector_add(a, b):    """    向量加法主机端函数    """    # 确保输入数组是连续的    a = np.ascontiguousarray(a)    b = np.ascontiguousarray(b)    # 创建输出数组    c = np.zeros_like(a)    # 将数组传输到GPU内存    d_a = cuda.to_device(a)    d_b = cuda.to_device(b)    d_c = cuda.device_array_like(c)    # 计算线程块和网格配置    threads_per_block = 256    blocks_per_grid = (c.size + (threads_per_block - 1)) // threads_per_block    # 启动CUDA Kernel    vector_add_kernel[blocks_per_grid, threads_per_block](d_a, d_b, d_c)    # 将结果传输回主机    c = d_c.copy_to_host()    return c![Python CUDA数据流](https://s.coze.cn/image/y0bvAy4DI7Y/)# 性能测试if __name__ == "__main__":    # 创建测试数据    size = 10_000_000    a = np.random.randn(size).astype(np.float32)    b = np.random.randn(size).astype(np.float32)    # GPU计算    c_gpu = vector_add(a, b)    # NumPy计算(对比)    c_numpy = a + b    # 验证结果    print(f"最大误差: {np.max(np.abs(c_gpu - c_numpy))}")    print(f"结果匹配: {np.allclose(c_gpu, c_numpy)}")

@cuda.jitdef demonstrate_thread_hierarchy(output):    """    演示CUDA线程层次结构    """    # 获取线程和块的信息    tx = cuda.threadIdx.x    # 线程在线程块中的索引    ty = cuda.threadIdx.y    bx = cuda.blockIdx.x     # 线程块在网格中的索引    by = cuda.blockIdx.y    bw = cuda.blockDim.x     # 线程块的维度    bh = cuda.blockDim.y    # 计算全局索引    global_x = bx * bw + tx    global_y = by * bh + ty    # 计算扁平化索引    flat_idx = global_y * (cuda.gridDim.x * bw) + global_x    if flat_idx < output.size:        output[flat_idx] = flat_idx# 使用二维网格配置def test_2d_grid(size_x, size_y):    output = np.zeros(size_x * size_y, dtype=np.int32)    # 配置二维线程块和网格    threads_per_block = (16, 16)    blocks_per_grid_x = (size_x + 15) // 16    blocks_per_grid_y = (size_y + 15) // 16    blocks_per_grid = (blocks_per_grid_x, blocks_per_grid_y)    d_output = cuda.device_array_like(output)    demonstrate_thread_hierarchy[blocks_per_grid, threads_per_block](d_output)    return d_output.copy_to_host()

@cuda.jitdef matrix_multiply_shared(a, b, c):    """    使用共享内存优化的矩阵乘法    """    # 定义块大小    TILE_SIZE = 16    # 声明共享内存    tile_a = cuda.shared.array((TILE_SIZE, TILE_SIZE), dtype=np.float32)    tile_b = cuda.shared.array((TILE_SIZE, TILE_SIZE), dtype=np.float32)    # 计算全局行列索引    row = cuda.blockIdx.y * TILE_SIZE + cuda.threadIdx.y    col = cuda.blockIdx.x * TILE_SIZE + cuda.threadIdx.x    # 累积结果    sum_val = 0.0    # 遍历所有块    for i in range((a.shape[1] + TILE_SIZE - 1) // TILE_SIZE):        # 加载数据到共享内存        if row < a.shape[0] and i * TILE_SIZE + cuda.threadIdx.x < a.shape[1]:            tile_a[cuda.threadIdx.y, cuda.threadIdx.x] = a[row, i * TILE_SIZE + cuda.threadIdx.x]        else:            tile_a[cuda.threadIdx.y, cuda.threadIdx.x] = 0.0        if i * TILE_SIZE + cuda.threadIdx.y < b.shape[0] and col < b.shape[1]:            tile_b[cuda.threadIdx.y, cuda.threadIdx.x] = b[i * TILE_SIZE + cuda.threadIdx.y, col]        else:            tile_b[cuda.threadIdx.y, cuda.threadIdx.x] = 0.0        # 同步所有线程        cuda.syncthreads()        # 计算部分乘积        for j in range(TILE_SIZE):            sum_val += tile_a[cuda.threadIdx.y, j] * tile_b[j, cuda.threadIdx.x]        # 再次同步        cuda.syncthreads()    # 写入结果    if row < c.shape[0] and col < c.shape[1]:        c[row, col] = sum_val# 性能测试def benchmark_matrix_multiplication(size):    """    矩阵乘法性能基准测试    """    import time    # 创建随机矩阵    a = np.random.randn(size, size).astype(np.float32)    b = np.random.randn(size, size).astype(np.float32)    c = np.zeros((size, size), dtype=np.float32)    # 准备GPU内存    d_a = cuda.to_device(a)    d_b = cuda.to_device(b)    d_c = cuda.device_array_like(c)    # 配置执行参数    TILE_SIZE = 16    threads_per_block = (TILE_SIZE, TILE_SIZE)    blocks_per_grid = ((size + TILE_SIZE - 1) // TILE_SIZE,                        (size + TILE_SIZE - 1) // TILE_SIZE)    # 预热    matrix_multiply_shared[blocks_per_grid, threads_per_block](d_a, d_b, d_c)    cuda.synchronize()    # 测试GPU性能    start_time = time.time()    for _ in range(10):        matrix_multiply_shared[blocks_per_grid, threads_per_block](d_a, d_b, d_c)    cuda.synchronize()    gpu_time = (time.time() - start_time) / 10    # 测试NumPy性能    start_time = time.time()    for _ in range(10):        c_numpy = np.dot(a, b)    numpy_time = (time.time() - start_time) / 10    print(f"矩阵大小: {size}x{size}")    print(f"GPU时间: {gpu_time*1000:.2f} ms")    print(f"NumPy时间: {numpy_time*1000:.2f} ms")    print(f"加速比: {numpy_time/gpu_time:.2f}x")

def demonstrate_streams(data_chunks):    """    演示使用CUDA Stream进行异步操作    """    # 创建多个Stream    num_streams = 4    streams = [cuda.stream() for _ in range(num_streams)]    # 准备输出    results = []    # 定义处理函数    @cuda.jit    def process_data(data, output):        idx = cuda.grid(1)        if idx < output.size:            output[idx] = data[idx] * 2.0 + 1.0    # 在不同流中并发处理数据块    for i, chunk in enumerate(data_chunks):        stream = streams[i % num_streams]        # 在流中分配GPU内存        d_chunk = cuda.to_device(chunk, stream=stream)        d_output = cuda.device_array_like(chunk, stream=stream)        # 配置Kernel        threads_per_block = 256        blocks_per_grid = (chunk.size + threads_per_block - 1) // threads_per_block        # 在流中启动Kernel        process_data[blocks_per_grid, threads_per_block, stream](d_chunk, d_output)        # 异步复制结果        result = d_output.copy_to_host(stream=stream)        results.append(result)    # 同步所有流    for stream in streams:        stream.synchronize()    return results# 性能对比测试def benchmark_streams():    """    Stream性能基准测试    """    import time    # 创建测试数据    chunk_size = 10_000_000    num_chunks = 8    data_chunks = [np.random.randn(chunk_size).astype(np.float32)                    for _ in range(num_chunks)]    # 无Stream版本    start_time = time.time()    for chunk in data_chunks:        result = chunk * 2.0 + 1.0    serial_time = time.time() - start_time    # 使用Stream版本    start_time = time.time()    stream_results = demonstrate_streams(data_chunks)    stream_time = time.time() - start_time    print(f"串行处理时间: {serial_time:.4f} 秒")    print(f"Stream处理时间: {stream_time:.4f} 秒")    print(f"性能提升: {serial_time/stream_time:.2f}x")

@cuda.jitdef parallel_sum_with_atomic(data, result):    """    使用原子操作实现并行求和    """    idx = cuda.grid(1)    if idx < data.size:        # 原子加法        cuda.atomic.add(result, 0, data[idx])@cuda.jitdef parallel_histogram(data, histogram, bins, min_val, max_val):    """    并行直方图计算    """    idx = cuda.grid(1)    if idx < data.size:        # 计算bin索引        value = data[idx]        if min_val <= value < max_val:            bin_idx = int((value - min_val) / (max_val - min_val) * bins)            bin_idx = min(bin_idx, bins - 1)            # 原子递增            cuda.atomic.add(histogram, bin_idx, 1)def test_atomic_operations():    """    测试原子操作    """    # 并行求和测试    data = np.random.randn(1_000_000).astype(np.float32)    d_data = cuda.to_device(data)    result = np.array([0.0], dtype=np.float32)    d_result = cuda.to_device(result)    threads_per_block = 256    blocks_per_grid = (data.size + threads_per_block - 1) // threads_per_block    parallel_sum_with_atomic[blocks_per_grid, threads_per_block](d_data, d_result)    cuda.synchronize()    sum_result = d_result.copy_to_host()[0]    expected_sum = np.sum(data)    print(f"原子求和结果: {sum_result:.6f}")    print(f"期望结果: {expected_sum:.6f}")    print(f"误差: {abs(sum_result - expected_sum):.6e}")    # 直方图测试    data = np.random.randn(100_000) * 10 + 50  # 正态分布,均值50,标准差10    histogram = np.zeros(50, dtype=np.int32)    d_data = cuda.to_device(data)    d_histogram = cuda.to_device(histogram)    threads_per_block = 256    blocks_per_grid = (data.size + threads_per_block - 1) // threads_per_block    parallel_histogram[blocks_per_grid, threads_per_block](        d_data, d_histogram, 50, 0.0, 100.0)    cuda.synchronize()    histogram = d_histogram.copy_to_host()    print(f"直方图统计: {histogram[:10]}...")  # 打印前10个bin

@cuda.jit(device=True)def merge_device(arr, temp, left, mid, right):    """    设备端归并函数    """    i = left    j = mid + 1    k = left    while i <= mid and j <= right:        if arr[i] <= arr[j]:            temp[k] = arr[i]            i += 1        else:            temp[k] = arr[j]            j += 1        k += 1    while i <= mid:        temp[k] = arr[i]        i += 1        k += 1    while j <= right:        temp[k] = arr[j]        j += 1        k += 1    for i in range(left, right + 1):        arr[i] = temp[i]@cuda.jitdef merge_sort_device(arr, temp, left, right, depth, max_depth):    """    设备端归并排序    """    idx = cuda.grid(1)    # 使用线程ID来决定处理的子数组    # 这里简化实现,实际应用中需要更复杂的调度策略    if idx == 0 and depth < max_depth:        if left < right:            mid = (left + right) // 2            # 递归调用(注意:这不是真正的动态并行)            # 在实际应用中,需要使用其他策略            merge_sort_device(arr, temp, left, mid, depth + 1, max_depth)            merge_sort_device(arr, temp, mid + 1, right, depth + 1, max_depth)            merge_device(arr, temp, left, mid, right)def hybrid_merge_sort(arr):    """    混合归并排序:主机端递归 + 设备端归并    """    import math    if len(arr) <= 1:        return arr    # 将数组复制到设备    d_arr = cuda.to_device(arr)    d_temp = cuda.device_array_like(arr)    # 在主机端进行递归分解,在设备端执行归并    def sort_range(left, right):        if left >= right:            return        mid = (left + right) // 2        sort_range(left, mid)        sort_range(mid + 1, right)        # 在设备端执行归并        threads_per_block = min(256, right - left + 1)        blocks_per_grid = 1        # 定义设备端归并Kernel        @cuda.jit        def merge_kernel(arr, temp, left, mid, right):            idx = cuda.grid(1)            global_idx = left + idx            if global_idx > right:                return            i = left            j = mid + 1            k = left            # 找到global_idx位置应该放置的元素            count = 0            while i <= mid and j <= right:                if arr[i] <= arr[j]:                    if count == idx:                        temp[global_idx] = arr[i]                        break                    i += 1                else:                    if count == idx:                        temp[global_idx] = arr[j]                        break                    j += 1                count += 1            if count == idx:                # 已经找到位置                pass            elif i <= mid:                temp[global_idx] = arr[i + idx - count]            else:                temp[global_idx] = arr[j + idx - count]        merge_kernel[blocks_per_grid, threads_per_block](d_arr, d_temp, left, mid, right)        cuda.synchronize()        # 复制回原数组        d_arr[left:right+1] = d_temp[left:right+1]    sort_range(0, len(arr) - 1)    return d_arr.copy_to_host()

import cupy as cpdef demonstrate_cupy_numba_integration():    """    演示CuPy与Numba的集成    """    # 使用CuPy进行数据传输和基础操作    size = 10_000_000    a_cpu = np.random.randn(size).astype(np.float32)    b_cpu = np.random.randn(size).astype(np.float32)    # 使用CuPy创建GPU数组    a_gpu = cp.asarray(a_cpu)    b_gpu = cp.asarray(b_cpu)    # CuPy基础操作(高度优化)    c_gpu_cupy = a_gpu + b_gpu    # Numba CUDA Kernel    @cuda.jit    def custom_kernel(a, b, c):        idx = cuda.grid(1)        if idx < c.size:            # 自定义计算逻辑            c[idx] = a[idx] * b[idx] + (a[idx] + b[idx]) * 0.5    # 将CuPy数组转换为Numba可用的格式    d_a = cuda.as_cuda_array(a_gpu)    d_b = cuda.as_cuda_array(b_gpu)    d_c = cuda.device_array_like(d_a)    threads_per_block = 256    blocks_per_grid = (size + threads_per_block - 1) // threads_per_block    custom_kernel[blocks_per_grid, threads_per_block](d_a, d_b, d_c)    cuda.synchronize()    # 将结果转换回CuPy数组    c_gpu_numba = cp.asarray(d_c)    # 性能对比    import time    # CuPy性能    start = time.time()    for _ in range(100):        c_temp = a_gpu + b_gpu    cp.cuda.Stream.null.synchronize()    cupy_time = time.time() - start    # Numba性能    start = time.time()    for _ in range(100):        custom_kernel[blocks_per_grid, threads_per_block](d_a, d_b, d_c)    cuda.synchronize()    numba_time = time.time() - start    print(f"CuPy计算时间: {cupy_time*1000:.2f} ms")    print(f"Numba计算时间: {numba_time*1000:.2f} ms")    return c_gpu_cupy, c_gpu_numbadef advanced_cupy_operations():    """    高级CuPy操作示例    """    # 创建GPU数组    x = cp.random.randn(1000, 1000, dtype=np.float32)    # 矩阵运算    result = cp.dot(x, x.T)    # FFT    fft_result = cp.fft.fft2(x)    # 线性代数    eigenvalues, eigenvectors = cp.linalg.eigh(result)    # 广播和通用函数    y = cp.sin(x) * cp.exp(-x**2)    # 归约操作    mean_val = cp.mean(x)    std_val = cp.std(x)    print(f"矩阵形状: {result.shape}")    print(f"特征值数量: {len(eigenvalues)}")    print(f"均值: {mean_val:.4f}, 标准差: {std_val:.4f}")

def performance_analysis_demo():    """    性能分析演示    """    from numba import cuda    # 创建测试数据    size = 10_000_000    data = np.random.randn(size).astype(np.float32)    # 定义待分析的Kernel    @cuda.jit    def analyze_kernel(data, output):        idx = cuda.grid(1)        if idx < output.size:            # 复杂计算逻辑            temp = data[idx]            for i in range(10):                temp = temp * 1.1 + 0.5                if temp > 100:                    temp = temp / 2.0            output[idx] = temp    d_data = cuda.to_device(data)    d_output = cuda.device_array_like(data)    threads_per_block = 256    blocks_per_grid = (size + threads_per_block - 1) // threads_per_block    # 使用CUDA事件进行精确计时    start = cuda.event()    end = cuda.event()    start.record()    analyze_kernel[blocks_per_grid, threads_per_block](d_data, d_output)    end.record()    end.synchronize()    elapsed_time = cuda.event_elapsed_time(start, end)    print(f"Kernel执行时间: {elapsed_time:.3f} ms")    # 计算性能指标    num_operations = size * 10  # 每个元素10次操作    gflops = (num_operations / elapsed_time / 1e6)    print(f"性能: {gflops:.2f} GFLOPS")    # 内存带宽分析    bytes_transferred = data.nbytes * 3  # 读输入,写输出,读输入    bandwidth = (bytes_transferred / elapsed_time / 1e6)    print(f"内存带宽: {bandwidth:.2f} GB/s")def optimization_checklist():    """    CUDA优化检查清单    """    print("CUDA性能优化检查清单:")    print("=" * 50)    checklist = {        "内存合并访问": "确保连续线程访问连续内存",        "共享内存使用": "利用共享内存减少全局内存访问",        "避免分支分歧": "最小化Warp内的条件分支",        "合适的块大小": "选择32的倍数作为块大小",        "足够的并行度": "确保有足够的块来隐藏延迟",        "内存传输最小化": "减少主机与设备间的数据传输",        "使用Stream": "实现计算和传输的重叠",        "原子操作优化": "减少原子操作的使用频率",        "数据类型选择": "使用适当的数据类型(如float32)",        "编译器优化": "使用适当的编译器优化标志"    }    for i, (item, description) in enumerate(checklist.items(), 1):        print(f"{i}. {item}: {description}")

import cv2from numba import cudaimport numpy as npclass GPUImageProcessor:    """    GPU图像处理器    """    def __init__(self):        self.stream = cuda.stream()    @staticmethod    @cuda.jit    def rgb_to_grayscale_kernel(rgb, gray):        """        RGB转灰度Kernel        """        x, y = cuda.grid(2)        if x < gray.shape[1] and y < gray.shape[0]:            # 使用标准灰度转换公式            gray[y, x] = (0.299 * rgb[y, x, 0] +                          0.587 * rgb[y, x, 1] +                          0.114 * rgb[y, x, 2])    @staticmethod    @cuda.jit    def gaussian_blur_kernel(input_img, output_img, kernel, radius):        """        高斯模糊Kernel        """        x, y = cuda.grid(2)        height, width = input_img.shape        if x < width and y < height:            sum_val = 0.0            kernel_sum = 0.0            # 应用卷积核            for ky in range(-radius, radius + 1):                for kx in range(-radius, radius + 1):                    px = min(max(x + kx, 0), width - 1)                    py = min(max(y + ky, 0), height - 1)                    k_idx_x = kx + radius                    k_idx_y = ky + radius                    k_val = kernel[k_idx_y, k_idx_x]                    sum_val += input_img[py, px] * k_val                    kernel_sum += k_val            output_img[y, x] = sum_val / kernel_sum    @staticmethod    @cuda.jit    def sobel_edge_detection_kernel(input_img, output_img):        """        Sobel边缘检测Kernel        """        x, y = cuda.grid(2)        height, width = input_img.shape        if x > 0 and x < width - 1 and y > 0 and y < height - 1:            # Sobel算子            gx = (-1 * input_img[y-1, x-1] + 1 * input_img[y-1, x+1] +                  -2 * input_img[y, x-1]   + 2 * input_img[y, x+1] +                  -1 * input_img[y+1, x-1] + 1 * input_img[y+1, x+1])            gy = (-1 * input_img[y-1, x-1] + -2 * input_img[y-1, x] + -1 * input_img[y-1, x+1] +                   1 * input_img[y+1, x-1] +  2 * input_img[y+1, x] +  1 * input_img[y+1, x+1])            output_img[y, x] = np.sqrt(gx**2 + gy**2)    def process_image_pipeline(self, image_path):        """        完整的图像处理流水线        """        # 读取图像        image = cv2.imread(image_path)        if image is None:            raise ValueError("无法读取图像")        # 转换为float32        image = image.astype(np.float32)        # 配置线程块和网格        threads_per_block = (16, 16)        blocks_per_grid = ((image.shape[1] + 15) // 16,                           (image.shape[0] + 15) // 16)        # 步骤1: RGB转灰度        gray = np.zeros((image.shape[0], image.shape[1]), dtype=np.float32)        d_image = cuda.to_device(image, stream=self.stream)        d_gray = cuda.device_array_like(gray, stream=self.stream)        self.rgb_to_grayscale_kernel[blocks_per_grid, threads_per_block, self.stream](            d_image, d_gray)        # 步骤2: 高斯模糊        # 创建高斯核        radius = 1        kernel_size = 2 * radius + 1        kernel = np.zeros((kernel_size, kernel_size), dtype=np.float32)        sigma = 1.0        for y in range(kernel_size):            for x in range(kernel_size):                dx = x - radius                dy = y - radius                kernel[y, x] = np.exp(-(dx**2 + dy**2) / (2 * sigma**2))        kernel = kernel / np.sum(kernel)        d_blurred = cuda.device_array_like(d_gray, stream=self.stream)        d_kernel = cuda.to_device(kernel, stream=self.stream)        self.gaussian_blur_kernel[blocks_per_grid, threads_per_block, self.stream](            d_gray, d_blurred, d_kernel, radius)        # 步骤3: 边缘检测        d_edges = cuda.device_array_like(d_blurred, stream=self.stream)        self.sobel_edge_detection_kernel[blocks_per_grid, threads_per_block, self.stream](            d_blurred, d_edges)        # 获取结果        edges = d_edges.copy_to_host(stream=self.stream)        self.stream.synchronize()        # 归一化到0-255范围        edges = np.clip(edges, 0, 255).astype(np.uint8)        return edges# 性能测试def benchmark_image_processing(image_path):    """    图像处理性能基准测试    """    import time    processor = GPUImageProcessor()    # GPU版本    start_time = time.time()    gpu_result = processor.process_image_pipeline(image_path)    gpu_time = time.time() - start_time    # CPU版本(使用OpenCV)    image = cv2.imread(image_path)    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)    blurred = cv2.GaussianBlur(gray, (3, 3), 1.0)    cpu_edges = cv2.Sobel(blurred, cv2.CV_64F, 1, 1, ksize=3)    cpu_edges = np.abs(cpu_edges).astype(np.uint8)    cpu_time = time.time() - start_time - gpu_time  # 减去GPU时间    print(f"GPU处理时间: {gpu_time*1000:.2f} ms")    print(f"CPU处理时间: {cpu_time*1000:.2f} ms")    print(f"加速比: {cpu_time/gpu_time:.2f}x")    return gpu_result, cpu_edges

@cuda.jitdef monte_carlo_pi_kernel(samples, results):    """    蒙特卡洛方法计算π的CUDA Kernel    """    idx = cuda.grid(1)    if idx < results.size:        count = 0        # 每个线程处理多个样本        samples_per_thread = samples // results.size        # 使用简单的伪随机数生成器        seed = idx * 12345 + 67890        for i in range(samples_per_thread):            # 线性同余生成器            seed = (1664525 * seed + 1013904223) & 0xFFFFFFFF            x = (seed & 0xFFFF) / 65535.0            seed = (1664525 * seed + 1013904223) & 0xFFFFFFFF            y = (seed & 0xFFFF) / 65535.0            # 检查是否在单位圆内            if x*x + y*y <= 1.0:                count += 1        results[idx] = countdef monte_carlo_pi_estimation(total_samples, num_threads=1024):    """    使用蒙特卡洛方法估算π值    """    import time    # 分配结果数组    results = np.zeros(num_threads, dtype=np.int32)    d_results = cuda.to_device(results)    # 配置Kernel    threads_per_block = 256    blocks_per_grid = (num_threads + threads_per_block - 1) // threads_per_block    # GPU计算    start_time = time.time()    monte_carlo_pi_kernel[blocks_per_grid, threads_per_block](        total_samples, d_results)    cuda.synchronize()    gpu_time = time.time() - start_time    # 获取结果    results = d_results.copy_to_host()    total_inside = np.sum(results)    # 计算π值    pi_estimate = 4.0 * total_inside / total_samples    error = abs(pi_estimate - np.pi)    print(f"总样本数: {total_samples:,}")    print(f"GPU计算时间: {gpu_time*1000:.2f} ms")    print(f"π估算值: {pi_estimate:.10f}")    print(f"π真实值: {np.pi:.10f}")    print(f"误差: {error:.2e}")    print(f"相对误差: {error/np.pi*100:.4f}%")    return pi_estimatedef option_pricing_monte_carlo(spot_price, strike_price, risk_free_rate,                                volatility, time_to_maturity, num_simulations):    """    使用蒙特卡洛方法进行欧式期权定价    """    @cuda.jit    def option_pricing_kernel(spot, strike, rate, vol, t, num_paths, results):        idx = cuda.grid(1)        if idx < results.size:            paths_per_thread = num_paths // results.size            # 预计算常数            drift = (rate - 0.5 * vol * vol) * t            diffusion = vol * np.sqrt(t)            total_payoff = 0.0            for i in range(paths_per_thread):                # 生成标准正态随机数(Box-Muller变换)                seed = idx * paths_per_thread + i + 12345                u1 = ((seed * 1664525 + 1013904223) & 0xFFFFFFFF) / 4294967295.0                u2 = ((seed * 1103515245 + 12345) & 0xFFFFFFFF) / 4294967295.0                z0 = np.sqrt(-2.0 * np.log(u1)) * np.cos(2.0 * np.pi * u2)                # 几何布朗运动                terminal_price = spot * np.exp(drift + diffusion * z0)                # 计算收益                payoff = max(terminal_price - strike, 0.0)                total_payoff += payoff            results[idx] = total_payoff    num_threads = 1024    results = np.zeros(num_threads, dtype=np.float64)    d_results = cuda.to_device(results)    threads_per_block = 256    blocks_per_grid = (num_threads + threads_per_block - 1) // threads_per_block    option_pricing_kernel[blocks_per_grid, threads_per_block](        spot_price, strike_price, risk_free_rate,         volatility, time_to_maturity, num_simulations, d_results)    cuda.synchronize()    results = d_results.copy_to_host()    # 计算期权价格    avg_payoff = np.sum(results) / num_simulations    option_price = np.exp(-risk_free_rate * time_to_maturity) * avg_payoff    return option_price# 测试期权定价def test_option_pricing():    """    测试欧式看涨期权定价    """    spot_price = 100.0      # 标的资产当前价格    strike_price = 105.0    # 执行价格    risk_free_rate = 0.05   # 无风险利率    volatility = 0.2        # 波动率    time_to_maturity = 1.0  # 到期时间(年)    num_simulations = 10_000_000  # 模拟路径数    import time    start_time = time.time()    option_price = option_pricing_monte_carlo(        spot_price, strike_price, risk_free_rate,        volatility, time_to_maturity, num_simulations    )    gpu_time = time.time() - start_time    print(f"欧式看涨期权定价结果:")    print(f"标的资产价格: ${spot_price:.2f}")    print(f"执行价格: ${strike_price:.2f}")    print(f"无风险利率: {risk_free_rate*100:.1f}%")    print(f"波动率: {volatility*100:.1f}%")    print(f"到期时间: {time_to_maturity:.1f}年")    print(f"模拟路径数: {num_simulations:,}")    print(f"期权价格: ${option_price:.4f}")    print(f"计算时间: {gpu_time*1000:.2f} ms")

@cuda.jitdef monte_carlo_pi_kernel(samples, results):    """    蒙特卡洛方法计算π的CUDA Kernel    """    idx = cuda.grid(1)    if idx < results.size:        count = 0        # 每个线程处理多个样本        samples_per_thread = samples // results.size        # 使用简单的伪随机数生成器        seed = idx * 12345 + 67890        for i in range(samples_per_thread):            # 线性同余生成器            seed = (1664525 * seed + 1013904223) & 0xFFFFFFFF            x = (seed & 0xFFFF) / 65535.0            seed = (1664525 * seed + 1013904223) & 0xFFFFFFFF            y = (seed & 0xFFFF) / 65535.0            # 检查是否在单位圆内            if x*x + y*y <= 1.0:                count += 1        results[idx] = countdef monte_carlo_pi_estimation(total_samples, num_threads=1024):    """    使用蒙特卡洛方法估算π值    """    import time    # 分配结果数组    results = np.zeros(num_threads, dtype=np.int32)    d_results = cuda.to_device(results)    # 配置Kernel    threads_per_block = 256    blocks_per_grid = (num_threads + threads_per_block - 1) // threads_per_block    # GPU计算    start_time = time.time()    monte_carlo_pi_kernel[blocks_per_grid, threads_per_block](        total_samples, d_results)    cuda.synchronize()    gpu_time = time.time() - start_time    # 获取结果    results = d_results.copy_to_host()    total_inside = np.sum(results)    # 计算π值    pi_estimate = 4.0 * total_inside / total_samples    error = abs(pi_estimate - np.pi)    print(f"总样本数: {total_samples:,}")    print(f"GPU计算时间: {gpu_time*1000:.2f} ms")    print(f"π估算值: {pi_estimate:.10f}")    print(f"π真实值: {np.pi:.10f}")    print(f"误差: {error:.2e}")    print(f"相对误差: {error/np.pi*100:.4f}%")    return pi_estimatedef option_pricing_monte_carlo(spot_price, strike_price, risk_free_rate,                                volatility, time_to_maturity, num_simulations):    """    使用蒙特卡洛方法进行欧式期权定价    """    @cuda.jit    def option_pricing_kernel(spot, strike, rate, vol, t, num_paths, results):        idx = cuda.grid(1)        if idx < results.size:            paths_per_thread = num_paths // results.size            # 预计算常数            drift = (rate - 0.5 * vol * vol) * t            diffusion = vol * np.sqrt(t)            total_payoff = 0.0            for i in range(paths_per_thread):                # 生成标准正态随机数(Box-Muller变换)                seed = idx * paths_per_thread + i + 12345                u1 = ((seed * 1664525 + 1013904223) & 0xFFFFFFFF) / 4294967295.0                u2 = ((seed * 1103515245 + 12345) & 0xFFFFFFFF) / 4294967295.0                z0 = np.sqrt(-2.0 * np.log(u1)) * np.cos(2.0 * np.pi * u2)                # 几何布朗运动                terminal_price = spot * np.exp(drift + diffusion * z0)                # 计算收益                payoff = max(terminal_price - strike, 0.0)                total_payoff += payoff            results[idx] = total_payoff    num_threads = 1024    results = np.zeros(num_threads, dtype=np.float64)    d_results = cuda.to_device(results)    threads_per_block = 256    blocks_per_grid = (num_threads + threads_per_block - 1) // threads_per_block    option_pricing_kernel[blocks_per_grid, threads_per_block](        spot_price, strike_price, risk_free_rate,         volatility, time_to_maturity, num_simulations, d_results)    cuda.synchronize()    results = d_results.copy_to_host()    # 计算期权价格    avg_payoff = np.sum(results) / num_simulations    option_price = np.exp(-risk_free_rate * time_to_maturity) * avg_payoff    return option_price# 测试期权定价def test_option_pricing():    """    测试欧式看涨期权定价    """    spot_price = 100.0      # 标的资产当前价格    strike_price = 105.0    # 执行价格    risk_free_rate = 0.05   # 无风险利率    volatility = 0.2        # 波动率    time_to_maturity = 1.0  # 到期时间(年)    num_simulations = 10_000_000  # 模拟路径数    import time    start_time = time.time()    option_price = option_pricing_monte_carlo(        spot_price, strike_price, risk_free_rate,        volatility, time_to_maturity, num_simulations    )    gpu_time = time.time() - start_time    print(f"欧式看涨期权定价结果:")    print(f"标的资产价格: ${spot_price:.2f}")    print(f"执行价格: ${strike_price:.2f}")    print(f"无风险利率: {risk_free_rate*100:.1f}%")    print(f"波动率: {volatility*100:.1f}%")    print(f"到期时间: {time_to_maturity:.1f}年")    print(f"模拟路径数: {num_simulations:,}")    print(f"期权价格: ${option_price:.4f}")    print(f"计算时间: {gpu_time*1000:.2f} ms")