E - Evaluate Matrix Sum

本文介绍了一种高效的矩阵子区域求和算法,通过预处理矩阵元素并利用累积和技巧来快速计算任意指定子区域内的元素平方和。适用于处理大量查询的情况。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

Description

Given a matrix, the elements of which are all integer number from 0 to 50, you are required to evaluate the square sum of its specified sub-matrix.


Input

The first line of the input contains a single integer T (1 <= T <= 5), the number of test cases.

For each test case, the first line contains two integers m and n (1 <= m, n <= 500), which are the row and column sizes of the matrix, respectively. The next m lines with n numbers each gives the elements of the matrix.

The next line contains a single integer N (1 <= N <= 100,000), the number of queries. The next N lines give one query on each line, with four integers r1, c1, r2, c2 (1 <= r1 <= r2 <= m, 1 <= c1 <= c2 <= n), which are the indices of the upper-left corner and lower-right corner of the sub-matrix in question.


Output

For each test case, first print the number of the test case, then N lines with one number on each line, the required square sum. Refer to the sample output for details.


Sample Input

2
2 3
1 2 3
4 5 6
2
1 1 2 2
1 3 2 3
3 3
4 2 3
2 5 1
7 9 2
1
1 1 3 3


Sample Output

Case 1:
46
45
Case 2:
193

题意:
给一个n*m矩阵 下表(1,1)~~~(n,m)
给两个点做对角线 求中间矩阵的值
朴素算法也可以过 但是这个在不断输入的同时也将值存在数组中 查询更快
  #include <iostream>
 #include <string.h>
 #include <stdio.h>

 using namespace std;

 int main()
 {
     int t;
     int n,m;
     int k;
     int tt,x1,x2,y1,y2;
     int num;
     int map[505][505];
     scanf("%d",&t);
     for(num=1;num<=t;num++)
     {
         scanf("%d%d",&n,&m);
         memset(map,0,sizeof(map));
         for(int i=1;i<=n;i++)
         {
             for(int j=1;j<=m;j++)
             {
                 scanf("%d",&k);
                 map[i][j]=k*k+map[i-1][j]+map[i][j-1]-map[i-1][j-1];
             }
         }
         printf("Case %d:\n",num);
         scanf("%d",&tt);
         while(tt--)
         {
             scanf("%d%d%d%d",&x1,&y1,&x2,&y2);
             int ans=map[x1-1][y1-1]+map[x2][y2]-map[x1-1][y2]-map[x2][y1-1];
             printf("%d\n",ans);

         }
     }
     return 0;
 }
 
 
 
 
 
 
 
 #include<iostream>
#include <string.h>
#include <stdio.h>

using namespace std;

int a[ 510 ][ 510 ];
int sum[ 510 ][ 510 ];

int main()
{
    int t,n,m,tt;
    int ss;
    int k;
    int r1,c1,r2,c2;
    scanf("%d",&t);
    for(k=1; k<=t; k++)
    {
        scanf("%d%d",&n,&m);
        memset(sum,0,sizeof(sum));
        for(int i=1; i<=n; i++)
        {
            ss = 0;
            for(int j=1; j<=m; j++)
            {
                scanf("%d",&a[i][j]);
                ss += a[i][j] * a[i][j];
                sum[i][j] = sum[ i - 1 ][ j ] + ss;
            }
        }
        printf("Case %d:\n",k);
        scanf("%d",&tt);
        while(tt--)
        {
            scanf("%d%d%d%d",&r1,&c1,&r2,&c2);
            printf("%d\n",sum[ r2 ][ c2 ] - sum[ r1-1 ][ c2 ] -  sum[ r2 ][ c1 - 1 ] + sum[ r1 - 1 ][ c1 - 1 ]  );
        }
    }
    return 0;

}


转载于:https://www.cnblogs.com/zhangying/p/3924337.html

clc; n = 40; % number of elements finalT = 0.2; % final time boundary = 1; % switch between Dirichlet conditions (1) and Absorbing (0) k = 3; % polynomial degree gamma = 1.4; % specific heat ratio CFL = 0.1; % Courant number for stability flux = 1; % flux type, 0 = LF, 1 = HDG nc = k+1; % number of quadrature points left = -5; % left end of the domain right = 5; % right end of the domain plot_accurate = 1; % plot resolution %% Quadrature setup [xg, wg] = get_gauss_quadrature(nc); xunit = get_gauss_lobatto_quadrature(k+1); %% Mesh creation kp1 = length(xunit); y = zeros(2*n, 1); h = zeros(n, 1); for e = 1:n y(2*e-1) = left + (right-left)*(e-1)/n; y(2*e) = left + (right-left)*e/n; h(e) = y(2*e) - y(2*e-1); end dx = (right-left)/n; % Node points x = zeros(kp1, n); for e = 1:n for i = 1:kp1 x(i,e) = y(2*e-1) + (y(2*e)-y(2*e-1))*(0.5 + 0.5*xunit(i)); end end % Time step calculation a_max = sqrt(gamma*1.0/1.0); % max wave speed estimate dt = CFL*dx/a_max; if dt <= 0 dt = finalT/100; end num_steps = ceil(finalT/dt); dt = finalT/num_steps; %% Initial conditions - Sod shock tube rho = zeros(kp1, n); u = zeros(kp1, n); p = zeros(kp1, n); for e = 1:n for j = 1:kp1 x_pos = x(j,e); if x_pos < 0.0 rho(j,e) = 1.0; u(j,e) = 0.0; p(j,e) = 1.0; else rho(j,e) = 0.125; u(j,e) = 0.0; p(j,e) = 0.1; end end end % Conservation variables E = p./(gamma-1) + 0.5*rho.*u.^2; U = zeros(3, kp1, n); U(1,:,:) = rho; U(2,:,:) = rho.*u; U(3,:,:) = E; %% Basis functions and mass matrix [psi, dpsi] = evaluate_lagrange_basis(xunit, xg); Me = psi * diag(wg) * psi'; M1inv = zeros(kp1, kp1, n); for e = 1:n M1inv(:,:,e) = inv(0.5*dx*Me); end %% Time stepping t = 0; for step = 1:num_steps % Runge-Kutta 3 k1 = rhs_corrected(U, dx, gamma, xunit, psi, dpsi, wg, M1inv); U1 = U + dt*k1; k2 = rhs_corrected(U1, dx, gamma, xunit, psi, dpsi, wg, M1inv); U2 = 0.75*U + 0.25*U1 + 0.25*dt*k2; k3 = rhs_corrected(U2, dx, gamma, xunit, psi, dpsi, wg, M1inv); U = (1/3)*U + (2/3)*U2 + (2/3)*dt*k3; t = t + dt; end %% 可视化部分修正 figure; % 统一使用[]自动计算维度 x_plot = reshape(x, [], 1); rho_plot = reshape(U(1,:,:), [], 1); % 修正速度计算,先重塑再除法 rho_reshaped = reshape(U(1,:,:), [], 1); u_reshaped = reshape(U(2,:,:), [], 1); u_plot = u_reshaped ./ rho_reshaped; % 修正压力计算 E_reshaped = reshape(U(3,:,:), [], 1); p_plot = reshape((gamma-1)*(E_reshaped - 0.5*u_reshaped.^2./rho_reshaped), [], 1); subplot(3,1,1); plot(x_plot, rho_plot, 'b.', 'MarkerSize', 10); title('Density Distribution'); xlabel('x'); ylabel('\rho'); grid on; subplot(3,1,2); plot(x_plot, u_plot, 'r.', 'MarkerSize', 10); title('Velocity Distribution'); xlabel('x'); ylabel('u'); grid on; subplot(3,1,3); plot(x_plot, p_plot, 'g.', 'MarkerSize', 10); title('Pressure Distribution'); xlabel('x'); ylabel('p'); grid on; %% Corrected right-hand side function function rhs = rhs_corrected(U, dx, gamma, xunit, psi, dpsi, weights, M1inv) [~, kp1, n] = size(U); rhs = zeros(size(U)); J = dx/2; % Jacobian % Compute Euler fluxes at solution points F = zeros(3, kp1, n); for e = 1:n for i = 1:kp1 rho = U(1,i,e); u_val = U(2,i,e)/rho; E_val = U(3,i,e); p_val = (gamma-1)*(E_val - 0.5*rho*u_val^2); F(1,i,e) = rho*u_val; F(2,i,e) = rho*u_val^2 + p_val; F(3,i,e) = u_val*(E_val + p_val); end end % Element integral term: -∫dφ/dx·F dx for e = 1:n % Interpolate flux to quadrature points F_q = zeros(3, length(weights)); for q = 1:length(weights) F_interp = [0; 0; 0]; for i = 1:kp1 F_interp = F_interp + psi(i,q)*F(:,i,e); end F_q(:,q) = F_interp; % Corrected assignment end % Compute integral for each component for comp = 1:3 integral = zeros(kp1,1); for i = 1:kp1 for q = 1:length(weights) integral(i) = integral(i) + weights(q)*dpsi(i,q)*F_q(comp,q); end end rhs(comp,:,e) = -J * integral; end % Apply inverse mass matrix for comp = 1:3 rhs(comp,:,e) = (M1inv(:,:,e) * rhs(comp,:,e)')'; end end % Boundary flux terms for e = 1:n % Left boundary if e == 1 % Left boundary condition rho_L = 1.0; u_L = 0.0; p_L = 1.0; E_L = p_L/(gamma-1) + 0.5*rho_L*u_L^2; U_L = [rho_L; rho_L*u_L; E_L]; % Current element left state rho_R = U(1,1,e); u_R = U(2,1,e)/rho_R; E_R = U(3,1,e); p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2); U_R = [rho_R; rho_R*u_R; E_R]; % Numerical flux F_num = roe_flux(U_L, U_R, gamma); % Apply flux (using temporary variables) Minv = M1inv(:,:,e); boundary_vector = Minv(:,1); contribution = F_num * boundary_vector'; rhs(:,:,e) = rhs(:,:,e) - contribution; else % Internal left boundary U_L = [U(1,end,e-1); U(2,end,e-1); U(3,end,e-1)]; U_R = [U(1,1,e); U(2,1,e); U(3,1,e)]; F_num = roe_flux(U_L, U_R, gamma); % Apply to left element Minv_left = M1inv(:,:,e-1); boundary_vector_left = Minv_left(:,end); contribution_left = F_num * boundary_vector_left'; rhs(:,:,e-1) = rhs(:,:,e-1) + contribution_left; % Apply to current element Minv_current = M1inv(:,:,e); boundary_vector_current = Minv_current(:,1); contribution_current = F_num * boundary_vector_current'; rhs(:,:,e) = rhs(:,:,e) - contribution_current; end % Right boundary if e == n % Right boundary condition rho_R = 0.125; u_R = 0.0; p_R = 0.1; E_R = p_R/(gamma-1) + 0.5*rho_R*u_R^2; U_R = [rho_R; rho_R*u_R; E_R]; % Current element right state rho_L = U(1,end,e); u_L = U(2,end,e)/rho_L; E_L = U(3,end,e); p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2); U_L = [rho_L; rho_L*u_L; E_L]; F_num = roe_flux(U_L, U_R, gamma); % Apply flux Minv = M1inv(:,:,e); boundary_vector = Minv(:,end); contribution = F_num * boundary_vector'; rhs(:,:,e) = rhs(:,:,e) + contribution; else % Internal right boundary U_L = [U(1,end,e); U(2,end,e); U(3,end,e)]; U_R = [U(1,1,e+1); U(2,1,e+1); U(3,1,e+1)]; F_num = roe_flux(U_L, U_R, gamma); % Apply to current element Minv_current = M1inv(:,:,e); boundary_vector_current = Minv_current(:,end); contribution_current = F_num * boundary_vector_current'; rhs(:,:,e) = rhs(:,:,e) + contribution_current; % Apply to right element Minv_right = M1inv(:,:,e+1); boundary_vector_right = Minv_right(:,1); contribution_right = F_num * boundary_vector_right'; rhs(:,:,e+1) = rhs(:,:,e+1) - contribution_right; end end end %% Roe flux function function F = roe_flux(U_L, U_R, gamma) % Unpack variables rho_L = U_L(1); u_L = U_L(2)/rho_L; E_L = U_L(3); rho_R = U_R(1); u_R = U_R(2)/rho_R; E_R = U_R(3); % Compute pressures p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2); p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2); % Left and right fluxes F_L = [rho_L*u_L; rho_L*u_L^2 + p_L; u_L*(E_L + p_L)]; F_R = [rho_R*u_R; rho_R*u_R^2 + p_R; u_R*(E_R + p_R)]; % Roe averages sqrt_rho_L = sqrt(rho_L); sqrt_rho_R = sqrt(rho_R); sum_sqrt = sqrt_rho_L + sqrt_rho_R; u_roe = (sqrt_rho_L*u_L + sqrt_rho_R*u_R)/sum_sqrt; H_L = (E_L + p_L)/rho_L; H_R = (E_R + p_R)/rho_R; H_roe = (sqrt_rho_L*H_L + sqrt_rho_R*H_R)/sum_sqrt; a_roe = sqrt((gamma-1)*(H_roe - 0.5*u_roe^2)); % Wave strengths drho = rho_R - rho_L; du = u_R - u_L; dp = p_R - p_L; dU = [drho; du; dp]; lambda = [u_roe-a_roe, u_roe, u_roe+a_roe]; % Entropy fix epsilon = 0.1*a_roe; alpha = abs(lambda); for i = 1:3 if abs(lambda(i)) < epsilon alpha(i) = (lambda(i)^2 + epsilon^2)/(2*epsilon); end end % Numerical flux F = 0.5*(F_L + F_R) - 0.5*max(alpha)*dU; end %% Gaussian quadrature function [x, w] = get_gauss_quadrature(n) beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2)); T = diag(beta,1) + diag(beta,-1); [V, D] = eig(T); x = diag(D); w = 2*V(1,:).^2; [x, idx] = sort(x); w = w(idx)'; end %% Gauss-Lobatto quadrature function x = get_gauss_lobatto_quadrature(n) if n == 2 x = [-1; 1]; else M = zeros(n-2); for i = 1:n-3 M(i,i+1) = 0.5*sqrt(i*(i+2)/((i+0.5)*(i+1.5))); M(i+1,i) = M(i,i+1); end [V, D] = eig(M); x = [-1; diag(D); 1]; end end %% Lagrange basis evaluation function [psi, dpsi] = evaluate_lagrange_basis(points_lagrange, points_evaluation) n = length(points_lagrange); m = length(points_evaluation); psi = zeros(n, m); dpsi = zeros(n, m); % Compute weights w = ones(1, n); for i = 1:n for j = 1:n if i ~= j w(i) = w(i)/(points_lagrange(i) - points_lagrange(j)); end end end % Evaluate basis for k = 1:m x = points_evaluation(k); for i = 1:n % Compute product product = 1; for j = 1:n if j ~= i product = product * (x - points_lagrange(j)); end end psi(i,k) = w(i)*product; % Compute derivative sum_deriv = 0; for j = 1:n if j ~= i term = 1; for l = 1:n if l ~= i && l ~= j term = term * (x - points_lagrange(l)); end end sum_deriv = sum_deriv + term; end end dpsi(i,k) = w(i)*sum_deriv; end end end 分析上述代码并进行修正,给出修正后的代码
06-28
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值