@Iyatsu

Are you sure you want to delete the question?

If your question is resolved, you may close it.

Leaving a resolved question undeleted may help others!

We hope you find it useful!

大規模データ実行のために処理を高速化させたい

Q&A

解決したいこと

以下の入力データ生成プログラムを用いて、
・与えられた始点と終点の1~5番目まで短い経路の距離と経由地点
・新たに与えられた点について、その地点と道路網に最小の長さの道を作成したときの交差地点の座標
・幹線道路の両端の地点番号
を出力するプログラムを作成したのですが、現状では実行時間が非常に長く処理が厳しいため、この入力プログラムに対応できるように実行プログラムの改善をしたいです。また、配列の動的確保を行っていますが、容量が充分に足りているかも確認したいと考えています。

入力データ生成プログラム

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define MAX_N 200000
#define MAX_M 100000
#define MAX_P 1000
#define MAX_Q 1000

// 座標範囲
#define COORD_MIN 0
#define COORD_MAX  1000000

int rand_in_range(int l, int r) {
    return l + rand() % (r - l + 1);
}

int main() {
    FILE *fp = fopen("testcase.in", "w");
    if (!fp) {
        perror("Cannot open file");
        return 1;
    }

    srand(time(NULL));

    int N = rand_in_range(2, 200000);           
    int M = rand_in_range(1, 100000);           
    int P = rand_in_range(0, 1000);           
    int Q = rand_in_range(0, 1000);            

    fprintf(fp, "%d %d %d %d\n", N, M, P, Q);

    // 頂点座標
    for (int i = 0; i < N; i++) {
        int x = rand_in_range(COORD_MIN, COORD_MAX);
        int y = rand_in_range(COORD_MIN, COORD_MAX);
        fprintf(fp, "%d %d\n", x, y);
    }

    // 線分接続
    for (int i = 0; i < M; i++) {
        int u = rand_in_range(1, N);
        int v = rand_in_range(1, N);
        while (v == u) {
            v = rand_in_range(1, N); // 自己ループ回避
        }
        fprintf(fp, "%d %d\n", u, v);
    }

    // クエリ
    int cross_point_count = rand_in_range(N / 2, N * 2); // クエリ内でC+番号を入れる用

    for (int i = 0; i < Q; i++) {
        int use_cross_s = rand_in_range(0, 1);
        int use_cross_t = rand_in_range(0, 1);

        if (use_cross_s && cross_point_count > 0) {
            fprintf(fp, "C%d ", rand_in_range(1, cross_point_count));
        } else {
            fprintf(fp, "%d ", rand_in_range(1, N));
        }

        if (use_cross_t && cross_point_count > 0) {
            fprintf(fp, "C%d ", rand_in_range(1, cross_point_count));
        } else {
            fprintf(fp, "%d ", rand_in_range(1, N));
        }

        int k = rand_in_range(1, 5);
        fprintf(fp, "%d\n", k);
    }

    // 新たなP個の頂点の座標
    for (int i = 0; i < P; i++) {
        int x = rand_in_range(COORD_MIN, COORD_MAX);
        int y = rand_in_range(COORD_MIN, COORD_MAX);
        fprintf(fp, "%d %d\n", x, y);
    }

    fclose(fp);
    return 0;
}

実行プログラム

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define MAXN 200000  // 最大ノード数
#define MAXM 100000  // 最大クエリ数
#define MAXQ 1000   // 最大クエリ数
#define MAX_K 5    // 最多経路数
#define INF 1e18   // 無限大
#define EPS 1e-8   // 浮動小数点比較の誤差許容範囲

// 座標点構造体
typedef struct {
    double a, b;
} Point;

// 線分
typedef struct {
    int a, b;
} Segment;

// グラフの隣接リスト用エッジ構造体
typedef struct Edge {
    int to;
    double cost;
    struct Edge* next;
} Edge;

typedef struct {
    double cost;
    int path[MAXN];
    int length;
} Path;

// グローバル変数
int N,M,P;   // ノード数、クエリ数、新たに加えられる点の数
Point* points;   // 座標点一覧
Segment* segs;   // 線分一覧
Edge** graph;   // グラフの隣接リスト
Path* paths;
int node_count = 0;   // ノードの総数
int* cur_path;
int cur_len = 0;
int* vis;
int path_count_total = 0;
Point current_a, current_b;

// 2点が同じ座標かを比較(誤差考慮)
int equal(Point a, Point b) {
    return fabs(a.a - b.a) < EPS && fabs(a.b - b.b) < EPS;
}

// 2点間のユークリッド距離
double distance(Point a, Point b) {
    double dx = a.a - b.a, dy = a.b - b.b;
    return sqrt(dx * dx + dy * dy);
}

// 点pが線分ab上にあるか判定
int point_on_segment(Point p, Point a, Point b) {
    
    double dx1 = p.a - a.a;
    double dy1 = p.b - a.b;
    double dx2 = b.a - a.a;
    double dy2 = b.b - a.b;
  
    // 外積で一直線上にあるか判定し、内積で範囲内か確認
    if (fabs(dx1 * dy2 - dy1 * dx2) < EPS) {
        double dot = dx1 * dx2 + dy1 * dy2;
        double len_sq = dx2 * dx2 + dy2 * dy2;
        return (0 <= dot && dot <= len_sq) || fabs(dot) < EPS || fabs(dot - len_sq) < EPS;
    }
    return 0;
}

// 線分同士の交差判定と交点計算
int intersect(Point a, Point q1, Point b, Point q2, double* ix, double* iy) {
    
    // 端点同士が重なっているか確認
    if (point_on_segment(a, b, q2)) {
        *ix = a.a;
        *iy = a.b;
        return 1;
    }

    if (point_on_segment(q1, b, q2)) {
        *ix = q1.a;
        *iy = q1.b;
        return 1;
    }

    if (point_on_segment(b, a, q1)) {
        *ix = b.a;
        *iy = b.b;
        return 1;
    }

    if (point_on_segment(q2, a, q1)) {
        *ix = q2.a;
        *iy = q2.b;
        return 1;
    }

    // 行列式で交点計算
    double a_val = q1.a - a.a;
    double b_val = -(q2.a - b.a);
    double c = q1.b - a.b;
    double d = -(q2.b - b.b);
    double e = b.a - a.a;
    double f = b.b - a.b;
    double det = a_val * d - b_val * c;

    if (fabs(det) < EPS) return 0;   // 平行または一致

    double s = (d * e - b_val * f) / det;
    double t = (-c * e + a_val * f) / det;

    // 交差範囲のチェック
    if (s < EPS || s > 1 - EPS || t < EPS || t > 1 - EPS) return 0;

    *ix = a.a + s * (q1.a - a.a);
    *iy = a.b + s * (q1.b - a.b);

    // 既存点との重複回避
    for (int i = 0; i < node_count; i++) {
        if (equal((Point){*ix, *iy}, points[i])) {
            *ix = points[i].a;
            *iy = points[i].b;
            break;
        }
    }

    return 1;
}

// グラフに無向エッジを追加
void add_edge(int u, int v, double cost) {
    Edge* e = malloc(sizeof(Edge));
    e->to = v;
    e->cost = cost;
    e->next = graph[u];
    graph[u] = e;
}

Segment current_seg;

int cmp_point_along_segment(const void *a, const void *b) {
    const Point *pa = (const Point *)a;
    const Point *pb = (const Point *)b;

    double dx = current_b.a - current_a.a;
    double dy = current_b.b - current_a.b;

    double da = (pa->a - current_a.a) * dx + (pa->b - current_a.b) * dy;
    double db = (pb->a - current_a.a) * dx + (pb->b - current_a.b) * dy;

    if (da < db) return -1;
    if (da > db) return 1;
    return 0;
}

// 既存ノードのインデックス取得(なければ追加)
int get_node_index(Point p) {
    for (int i = 0; i < node_count; i++) {
        if (equal(points[i], p)) return i;
    }
    if (node_count >= (N + P + M * 20 + 1000)) {
        fprintf(stderr, "Error: node_count exceeded MAX_TOTAL_NODE (%d)\n", (N + P + M * 20 + 1000));
        exit(1);
    }
    points[node_count++] = p;
    return node_count - 1;
}

// 線分の交差を考慮しながらグラフ構築
void build_graph(int M) {

    // グラフ初期化
    for (int i = 0; i < (N + P + M * 20 + 1000); i++) {
        graph[i] = NULL;
    }
    
    for (int i = 0; i < M; i++) {
        Point a = points[segs[i].a];
        Point b = points[segs[i].b];
        Point* seg_pts = malloc(sizeof(Point) * (N + M * 10));   // 分割された点列
        int count = 0;
        seg_pts[count++] = a;
        seg_pts[count++] = b;

        // 他の線分との交点を探す
        for (int j = 0; j < M; j++) {
            if (i == j) continue;
            double ix, iy;
            if (intersect(a, b, points[segs[j].a], points[segs[j].b], &ix, &iy)) {
                Point p = {ix, iy};
                int found = 0;
                for (int k = 0; k < count; k++) {
                    if (equal(p, seg_pts[k])) {
                        found = 1;
                        break;
                    }
                }  
                if (!found) {
                    seg_pts[count++] = p;
                }
            }
        }

        current_a = a;
        current_b = b;
        // 点をソートして順番に辺を張る
        qsort(seg_pts, count, sizeof(Point), cmp_point_along_segment);

        // 重複を除去
        int unique_count = 1;

        for (int j = 1; j < count; j++) {
            if (!equal(seg_pts[j], seg_pts[unique_count-1])) {
                seg_pts[unique_count++] = seg_pts[j];
            }
        }

        count = unique_count;

        for (int k = 0; k < count - 1; k++) {
            int u = get_node_index(seg_pts[k]);
            int v = get_node_index(seg_pts[k + 1]);
            if (u == v) continue;  
            double d = distance(seg_pts[k], seg_pts[k + 1]);
            // 重複エッジチェック
            int edge_exists = 0;
            for (Edge* e = graph[u]; e; e = e->next) {
                if (e->to == v) {
                    edge_exists = 1;
                    break;
                }
            }
            if (!edge_exists) {
                add_edge(u, v, d);
                add_edge(v, u, d);   // 無向グラフ
            }
        }

        free(seg_pts);
        
    }

}

// 幹線道路検出用変数
int* tin; 
int* low; 
int* visited;
int timer = 0;

// 橋検出用DFS
void bridge_dfs(int v, int parent) {
    visited[v] = 1;
    tin[v] = low[v] = ++timer;

    for (Edge* e = graph[v]; e != NULL; e = e->next) {
        int to = e->to;
        if (to == parent) continue;

        if (visited[to]) {
            if (low[v] > tin[to]) low[v] = tin[to];
        } else {
            bridge_dfs(to, v);
            if (low[v] > low[to]) low[v] = low[to];

            if (low[to] > tin[v]) {
                printf("%d %d\n", v + 1, to + 1);  // 両端のノード番号を出力
            }
        }
    }
}


// ノード名(数値 or C+数値)からノードインデックスを取得
int parse_node(char* s, int N) {
    if (s[0] == 'C') {
        int id = atoi(s + 1);
        if (id < 1 || id > node_count - N) return -1;
        return N + id - 1;
    } else {
        int id = atoi(s);
        if (id < 1 || id > N) return -1;
        return id - 1;
    }
}

// 昇順にdoubleを比較するための関数(誤差考慮)
int comp_double(const void* A, const void* B) {
    double a = *(double*)A, b = *(double*)B;
    if (fabs(a - b) > EPS) return a < b ? -1 : 1;
    return 0;
}

// 経路を距離の昇順にソート
int compare_path(const void* a, const void* b) {
    double d1 = ((Path*)a)->cost;
    double d2 = ((Path*)b)->cost;
    if (fabs(d1 - d2) > EPS) return d1 < d2 ? -1 : 1;
    return 0;
}

// 深さ優先探索で全経路を探索し、目的地までの距離と経路を記録
void dfs(int u, int goal, double cost, int k) {
    if (path_count_total >= k && paths[k - 1].cost < cost) {
        return;
    }

    cur_path[cur_len++] = u;

    if (u == goal) {
        if (path_count_total >= k) {
            path_count_total--;
        }
        paths[path_count_total].cost = cost;
        paths[path_count_total].length = cur_len;
        for (int i = 0; i < cur_len; i++) {
            paths[path_count_total].path[i] = cur_path[i];
        }
        path_count_total++;
        qsort(paths, path_count_total, sizeof(Path), compare_path);
        cur_len--;
        return;
    }

    for (Edge* e = graph[u]; e; e = e->next) {
        int v = e->to;
        if (vis[v]) continue;
        vis[v] = 1;
        dfs(v, goal, cost + e->cost, k);
        vis[v] = 0;
    }

    cur_len--;
}


// DFSを呼び出し、結果をソートして準備
void solve(int start, int goal, int k, int output_path) {
    memset(vis, 0, sizeof(vis));
    path_count_total = 0;
    cur_len = 0;

    vis[start] = 1;
    dfs(start, goal, 0, k);
    vis[start] = 0;

    qsort(paths, path_count_total, sizeof(Path), compare_path);

    for (int i = 0; i < k && i < path_count_total; i++) {
        printf("%.5f\n", paths[i].cost);
        if (output_path) {
            for (int j = 0; j < paths[i].length; j++) {
                int node = paths[i].path[j];
                // 数字ノードとCラベルノードの区別
                if (node < N) {
                    printf("%d ", node + 1);   // ノード1〜Nはそのまま出力
                } else {
                    printf("C%d ", node - N + 1);   // ノードN以降は "C" + 番号で出力
                }
            }
            printf("\n");
        }
    }

    // 経路が存在しない場合
    if (path_count_total == 0) {
        printf("NA\n");
    }
}

void connect_new_points(int P) {
    for (int i = 0; i < P; i++) {
        Point p;
        scanf("%lf %lf", &p.a, &p.b);

        double min_dist = INF;
        Point best_point;

        for (int j = 0; j < M; j++) {
            Point a = points[segs[j].a];
            Point b = points[segs[j].b];

            double dx = b.a - a.a;
            double dy = b.b - a.b;
            double t = ((p.a - a.a) * dx + (p.b - a.b) * dy) / (dx * dx + dy * dy);
            if (t < 0) t = 0;
            if (t > 1) t = 1;

            Point proj = { a.a + t * dx, a.b + t * dy };
            double d = distance(p, proj);

            if (d + EPS < min_dist) {
                min_dist = d;
                best_point = proj;
            }
        }

        // 出力
        printf("%.5f %.5f\n", best_point.a, best_point.b);

        // グラフに接続
        int pi = get_node_index(p);
        int ci = get_node_index(best_point);
        double d = distance(p, best_point);
        add_edge(pi, ci, d);
        add_edge(ci, pi, d);

        // 線分にも記録(順番どおりにMを増やす)
        segs[M].a = pi;
        segs[M].b = ci;
        M++;
    }
}

// 入力読み込みと各クエリに対する出力
int main() {
    int P, Q;

    scanf("%d %d %d %d", &N, &M, &P, &Q);   // 点数、線分数、新たに加えられる点の数、クエリ数

    points = malloc(sizeof(Point) * (N + P + M * 20 + 1000)); 
    segs = malloc(sizeof(Segment) * (M + P + 100));     
    graph = malloc(sizeof(Edge*) * (N + P + M * 20 + 1000));
    paths = malloc(sizeof(Path) * (Q * MAX_K + 100));
    cur_path = malloc(sizeof(int) * (N + P + M * 20 + 1000));
    vis = malloc(sizeof(int) * (N + P + M * 20 + 1000));
    tin = malloc(sizeof(int) * (N + P + M * 20 + 1000));
    low = malloc(sizeof(int) * (N + P + M * 20 + 1000));
    visited = malloc(sizeof(int) * (N + P + M * 20 + 1000));

    if (!graph || !points || !segs) {
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
    }
    
    for (int i = 0; i < N; i++) {
        scanf("%lf %lf", &points[i].a, &points[i].b);
    }

    node_count = N;

    for (int i = 0; i < M; i++) {
        scanf("%d %d", &segs[i].a, &segs[i].b);
        segs[i].a--; segs[i].b--;   // 0-indexed に変換
    }

    build_graph(M);   // グラフ構築

    memset(visited, 0, sizeof(visited));   // 幹線道路検出
    timer = 0;

    for (int i = 0; i < node_count; i++) {
        if (!visited[i]) {
        bridge_dfs(i, -1);
        }
    }

    for (int q = 0; q < Q; q++) {
        char s1[10], s2[10];
        int k;
        scanf("%s %s %d", s1, s2, &k);
        int start = parse_node(s1, N);
        int goal = parse_node(s2, N);
        if (start == -1 || goal == -1) {
            printf("NA\n");
            continue;
        }

        solve(start, goal,k,1);   
        
    }

    connect_new_points(P);   

    free(points);
    free(segs);
    free(graph);
    free(paths);
    free(cur_path);
    free(vis);
    free(tin);
    free(low);
    free(visited);

    return 0;
}

自分で試したこと

どちらのプログラムもコンパイル済みです。
入力ファイルについては問題なく作成できました。
実行プログラムについてはまだ出力結果が出せていない状況です。

1 likes

2Answer

アルゴリズムやプログラムに詳しくないので、Windows copilotに聞いてみました。
以下、得られた回答を貼ってみます。


1.セグメント交差判定前の AABB チェック
build_graph 内の他セグメントとの二重ループで、まず長方形(AABB)の重なりを確認して早期スキップします。

// build_graph の for(j) の直前に追加
double min_ax = fmin(a.a, b.a), max_ax = fmax(a.a, b.a);
double min_ay = fmin(a.b, b.b), max_ay = fmax(a.b, b.b);

for(int j = 0; j < M; j++){
    if(i == j) continue;
    Point c = points[segs[j].a], d = points[segs[j].b];
    double min_cx = fmin(c.a, d.a), max_cx = fmax(c.a, d.a);
    double min_cy = fmin(c.b, d.b), max_cy = fmax(c.b, d.b);

    // AABB が重ならなければ交差判定をスキップ
    if(max_ax + EPS < min_cx  ||  max_cx + EPS < min_ax ||
       max_ay + EPS < min_cy  ||  max_cy + EPS < min_ay) {
        continue;
    }

    double ix, iy;
    if(intersect(a, b, c, d, &ix, &iy)){
        // …既存の処理…
    }
}

これだけで「ほとんどぶつからない組み合わせ」が定数時間でバイパスでき、交差関数呼び出しを大幅削減できます。


2.seg_pts の malloc/free → スタティックバッファ化
各セグメントごとに毎回 malloc/free すると非常に重いので、最大点数を予測してスタティック配列を再利用します。

// グローバル or build_graph 外部で一度だけ
#define MAX_SEG_PTS 4000   // M=100k なら交点平均 40 点でも安心余裕
static Point seg_pts_buf[MAX_SEG_PTS];

// build_graph 内
// Point* seg_pts = malloc(…);  を
Point* seg_pts = seg_pts_buf;
int count = 0;
seg_pts[count++] = a;
seg_pts[count++] = b;

// …最後
// free(seg_pts); を削除

→ malloc/free オーバーヘッドが消えます。


3.add_edge の malloc → メモリプール化
エッジを都度 malloc せず、あらかじめ大きな配列を作ってそこから振り分けます。

// グローバル
#define MAX_EDGES ( (N + P + M*20 + 1000) * 4 )
static Edge edge_pool[MAX_EDGES];
static int edge_pool_ptr = 0;

// add_edge を以下に置き換え
void add_edge(int u, int v, double cost) {
    Edge* e = &edge_pool[edge_pool_ptr++];
    e->to   = v;
    e->cost = cost;
    e->next = graph[u];
    graph[u] = e;
}

malloc/free が消えて、メモリアクセスが連続化されキャッシュヒット率も改善します。


4.get_node_index の線形探索 → ハッシュテーブル化
毎回 O(node_count) で比較している部分を、開番地法のシンプルなハッシュテーブルに置き換えます。

// グローバル
#define HASH_SIZE (1<<19)  // 524,288 項目
static int   ht_used[HASH_SIZE];
static int   ht_idx[HASH_SIZE];

// 大まかな double→int ハッシュ
static inline int hash_point(Point p) {
    uint64_t x = (uint64_t)llround(p.a * 1e6);
    uint64_t y = (uint64_t)llround(p.b * 1e6);
    return (int)((x * 1315423911u + y) & (HASH_SIZE - 1));
}

int get_node_index(Point p) {
    int h = hash_point(p);
    while (ht_used[h]) {
        int idx = ht_idx[h];
        if (equal(points[idx], p)) return idx;
        h = (h + 1) & (HASH_SIZE - 1);
    }
    // 新規登録
    ht_used[h] = 1;
    ht_idx[h]  = node_count;
    points[node_count] = p;
    return node_count++;
}

これで「座標→ノード番号」がほぼ O(1) に


5.k 最短経路探索を DFS → 修正版 Dijkstra ベースへ
指数関数的に全単純パスを列挙する DFS は、ノード数・枝数が増えると全く追いつきません。
以下のように「各ノードにつき最大 k 回だけ到達を許す」Dijkstra 形式のプライオリティキュー探索に切り替えると、O(k·E·logV) 程度で収束します。

// 必要なヘッダー
#include <stdbool.h>
#include <float.h>

typedef struct {
    double cost;
    int    v;
} State;

// 比較関数
static int cmp_state(const void* pa, const void* pb) {
    double da = ((State*)pa)->cost;
    double db = ((State*)pb)->cost;
    return da < db ? -1 : da > db ? 1 : 0;
}

// solve() 内を以下に差し替え
void solve(int start, int goal, int K, int output_path) {
    // 各ノード何回まで取り出したか
    int *cnt = calloc(node_count, sizeof(int));

    // ミニヒープとしての配列バッファ (消費上限は K * node_count)
    State *heap = malloc(sizeof(State) * K * node_count);
    int    heap_sz = 0;

    // 解を格納
    double *answers = malloc(sizeof(double) * K);
    int    found   = 0;

    // 初期状態プッシュ
    heap[heap_sz++] = (State){0.0, start};
    qsort(heap, heap_sz, sizeof(State), cmp_state);

    while (heap_sz > 0 && found < K) {
        // pop_min
        State cur = heap[0];
        heap[0] = heap[--heap_sz];
        qsort(heap, heap_sz, sizeof(State), cmp_state);

        if (++cnt[cur.v] > K) continue;
        if (cur.v == goal) {
            answers[found++] = cur.cost;
            continue;
        }

        // 隣接ノード展開
        for (Edge* e = graph[cur.v]; e; e = e->next) {
            heap[heap_sz++] = (State){ cur.cost + e->cost, e->to };
            // 小さい方に維持するため逐次ソート
            qsort(heap, heap_sz, sizeof(State), cmp_state);
        }
    }

    // 出力
    if (found == 0) {
        printf("NA\n");
    } else {
        for(int i = 0; i < found; i++) {
            printf("%.5f\n", answers[i]);
        }
    }

    free(cnt);
    free(heap);
    free(answers);
}

最後に

  • まずは①~④を入れて「グラフ構築がどれだけ速くなるか」をプロファイラで見てみる
  • 次に⑤を入れて k-最短経路が実用範囲内か計測
  • 必要なら「GPU 並列化」「空間分割データ構造 (Quadtree/K-D tree)」など、さらなる手を検討
    これらを段階的に適用すれば、「どこを直せば何が変わるか」が可視化され、安心して大幅高速化が狙えます。まずは AABB チェック+メモリプールあたりから挑戦してみてください。

以上になります。

0Like

Comments

  1. @Iyatsu

    Questioner

    ありがとうございます!

Your answer might help someone💌