I Love Square Number
Source: aizu

I Love Square Number
Problem Statement

Consider a graph with $\frac{N(N+1)}{2}$ vertices and $\frac{3N(N-1}{2}$ edges, where $N$ is an integer greater than or equals to 2.

The set of vertices is $\{(i,j) | 1 \leq i \leq N, 1 \leq j \leq i\}$.
There is an edge with weight $a_{i,j}$ between $(i,j)$ and $(i+1,j)$ (for $1 \leq i \leq N-1$ and $1 \leq j \leq i$).
There is an edge with weight $b_{i,j}$ between $(i,j)$ and $(i+1,j+1)$ (for $1 \leq i \leq N-1$ and $1 \leq j \leq i$).
There is an edge with weight $c_{i,j}$ between $(i,j)$ and $(i,j+1)$ (for $2 \leq i \leq N$ and $1 \leq j \leq i-1$).

For a simple path in this graph, the weight of the path is defined as the product of the weights of the edges that the path traverses.

Determine the number of unordered pairs of distinct vertices $\{s, t\}$ such that any simple path from $s$ to $t$ has a weight that is a square number.

Input

The input is given in the following format:

$N$
$a_{1,1}$
$a_{2,1}$ $a_{2,2}$
$\vdots$
$a_{N-1,1}$ ... $a_{N-1,N-1}$
$b_{1,1}$
$b_{2,1}$ $b_{2,2}$
$\vdots$
$b_{N-1,1}$ ... $b_{N-1,N-1}$
$c_{2,1}$
$c_{3,1}$ $c_{3,2}$
$\vdots$
$c_{N,1}$ ... $c_{N,N-1}$

$2 \leq N \leq 1,000$
$1 \leq a_{i,j}, b_{i,j}, c_{i,j} \leq 10^6$
All input values are integers.
Output

Output the answer.

Sample Input 1
2
1
2
2

Sample Output 1
1

Sample Input 2
3
1
2 3
4
5 6
7
8 9

Sample Output 2
0


C++ solution:
// AOJ #3631 I Love Square Number
// 2025.3.31

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ull = unsigned long long;

#define gc() getchar_unlocked()
#define pc(c) putchar_unlocked(c)

int Cin() {
    int n = 0, c = gc();
    do n = 10*n + (c & 0xf), c = gc(); while (c >= '0');
    return n;
}

void Cout(ll n) {
    char b[30];
    if (!n) pc('0');
    else {
        if (n < 0) pc('-'), n = -n;
        int i = 0; while (n) b[i++] = n % 10 + '0', n /= 10;
        while (i--) pc(b[i]);
    }
    pc('\n');
}

const int MAXV = 1000000;
vector<int> spf(MAXV+5);

void buildSpf(){
    for (int i = 1; i <= MAXV; i++) spf[i] = i;
    for (int i = 2; i * i <= MAXV; i++)
        if (spf[i] == i)
            for (int j = i * i; j <= MAXV; j += i) if (spf[j] == j) spf[j] = i;
}

unordered_map<int, ull> H;
mt19937_64 rng((unsigned)chrono::steady_clock::now().time_since_epoch().count());

ull sf(int x) {
    ull r = 0;
    while(x > 1){
        int p = spf[x], c = 0;
        while(x % p == 0){ c++; x /= p; }
        if(c & 1) {
            if(H.find(p) == H.end()) H[p] = rng();
            r ^= H[p];
        }
    }
    return r;
}

inline int idx(int r, int p){ return r*(r+1)/2 + p; }

int main(){
    buildSpf();

    int n = Cin();
    vector<vector<int>> A(n-1), B(n-1);
    for (int r = 0; r < n-1; r++){
        A[r].resize(r+1);
        for (int j = 0; j <= r; j++) A[r][j] = Cin();
    }
    for (int r = 0; r < n-1; r++){
        B[r].resize(r+1);
        for (int j = 0; j <= r; j++) B[r][j] = Cin();
    }
    vector<vector<int>> C(n-1);
    for (int r = 1; r < n; r++){
        C[r-1].resize(r);
        for (int j = 0; j < r; j++) C[r-1][j] = Cin();
    }

    int tot = n*(n+1)/2;
    vector<vector<pair<int,ull>>> g(tot);

    for (int r = 0; r < n-1; r++){
        for (int j = 0; j <= r; j++){
            int u = idx(r, j), v = idx(r+1, j);
            ull w = sf(A[r][j]);
            g[u].push_back({v, w});
            g[v].push_back({u, w});
        }
    }
    for (int r = 0; r < n-1; r++){
        for (int j = 0; j <= r; j++){
            int u = idx(r, j), v = idx(r+1, j+1);
            ull w = sf(B[r][j]);
            g[u].push_back({v, w});
            g[v].push_back({u, w});
        }
    }
    for (int r = 1; r < n; r++){
        for (int j = 0; j < r; j++){
            int u = idx(r, j), v = idx(r, j+1);
            ull w = sf(C[r-1][j]);
            g[u].push_back({v, w});
            g[v].push_back({u, w});
        }
    }

    vector<int> vis(tot,0);
    vector<ull> f(tot, 0);
    ll ans = 0;
    for (int s = 0; s < tot; s++){
        if(!vis[s]){
            vector<int> comp;
            bool ok = true;
            deque<int> dq;
            dq.push_back(s);
            vis[s] = 1;
            f[s] = 0;
            comp.push_back(s);
            while(!dq.empty()){
                int u = dq.front(); dq.pop_front();
                for(auto &e : g[u]){
                    int v = e.first; ull w = e.second;
                    ull nv = f[u] ^ w;
                    if(!vis[v]){
                        vis[v] = 1;
                        f[v] = nv;
                        dq.push_back(v);
                        comp.push_back(v);
                    } else if(f[v] != nv) ok = false;
                }
            }
            if(ok){
                unordered_map<ull, ll> mp;
                for(auto v : comp) mp[f[v]]++;
                for(auto &p : mp){
                    ll m = p.second;
                    ans += m*(m-1LL)/2;
                }
            }
        }
    }
    Cout(ans);
    return 0;
}

Go solution:
package main

import (
	"bufio"
	"fmt"
	"math"
	"math/rand"
	"os"
	"strconv"
	"strings"
	"time"
)

const MAXV = 1000000

var spf [MAXV + 1]int

func buildSpf() {
	for i := 1; i <= MAXV; i++ {
		spf[i] = i
	}
	for i := 2; i*i <= MAXV; i++ {
		if spf[i] == i {
			for j := i * i; j <= MAXV; j += i {
				if spf[j] == j {
					spf[j] = i
				}
			}
		}
	}
}

var H = make(map[int]uint64)
var rng = rand.New(rand.NewSource(time.Now().UnixNano()))

func squareFree(x int) uint64 {
	res := uint64(0)
	for x > 1 {
		p := spf[x]
		cnt := 0
		for x%p == 0 {
			cnt++
			x /= p
		}
		if cnt%2 == 1 {
			if _, exists := H[p]; !exists {
				H[p] = rng.Uint64()
			}
			res ^= H[p]
		}
	}
	return res
}

func idx(r, p int) int {
	return r*(r+1)/2 + p
}

func main() {
	buildSpf()

	scanner := bufio.NewScanner(os.Stdin)
	scanner.Split(bufio.ScanWords)

	// Read N
	scanner.Scan()
	n, _ := strconv.Atoi(scanner.Text())

	// Read A
	A := make([][]int, n-1)
	for r := 0; r < n-1; r++ {
		A[r] = make([]int, r+1)
		for j := 0; j <= r; j++ {
			scanner.Scan()
			A[r][j], _ = strconv.Atoi(scanner.Text())
		}
	}

	// Read B
	B := make([][]int, n-1)
	for r := 0; r < n-1; r++ {
		B[r] = make([]int, r+1)
		for j := 0; j <= r; j++ {
			scanner.Scan()
			B[r][j], _ = strconv.Atoi(scanner.Text())
		}
	}

	// Read C
	C := make([][]int, n-1)
	for r := 1; r < n; r++ {
		C[r-1] = make([]int, r)
		for j := 0; j < r; j++ {
			scanner.Scan()
			C[r-1][j], _ = strconv.Atoi(scanner.Text())
		}
	}

	tot := n * (n + 1) / 2
	g := make([][]struct {
		v int
		w uint64
	}, tot)

	// Build graph edges for A
	for r := 0; r < n-1; r++ {
		for j := 0; j <= r; j++ {
			u := idx(r, j)
			v := idx(r+1, j)
			w := squareFree(A[r][j])
			g[u] = append(g[u], struct {
				v int
				w uint64
			}{v, w})
			g[v] = append(g[v], struct {
				v int
				w uint64
			}{u, w})
		}
	}

	// Build graph edges for B
	for r := 0; r < n-1; r++ {
		for j := 0; j <= r; j++ {
			u := idx(r, j)
			v := idx(r+1, j+1)
			w := squareFree(B[r][j])
			g[u] = append(g[u], struct {
				v int
				w uint64
			}{v, w})
			g[v] = append(g[v], struct {
				v int
				w uint64
			}{u, w})
		}
	}

	// Build graph edges for C
	for r := 1; r < n; r++ {
		for j := 0; j < r; j++ {
			u := idx(r, j)
			v := idx(r, j+1)
			w := squareFree(C[r-1][j])
			g[u] = append(g[u], struct {
				v int
				w uint64
			}{v, w})
			g[v] = append(g[v], struct {
				v int
				w uint64
			}{u, w})
		}
	}

	vis := make([]bool, tot)
	f := make([]uint64, tot)
	ans := int64(0)

	for s := 0; s < tot; s++ {
		if !vis[s] {
			var comp []int
			ok := true
			q := []int{s}
			vis[s] = true
			f[s] = 0
			comp = append(comp, s)

			for len(q) > 0 {
				u := q[0]
				q = q[1:]
				for _, e := range g[u] {
					v, w := e.v, e.w
					nv := f[u] ^ w
					if !vis[v] {
						vis[v] = true
						f[v] = nv
						q = append(q, v)
						comp = append(comp, v)
					} else if f[v] != nv {
						ok = false
					}
				}
			}

			if ok {
				freq := make(map[uint64]int64)
				for _, v := range comp {
					freq[f[v]]++
				}
				for _, m := range freq {
					ans += m * (m - 1) / 2
				}
			}
		}
	}

	fmt.Println(ans)
}

Java solution:
import java.io.*;
import java.util.*;

public class Main {
    static final int MAXV = 1000000;
    static int[] spf = new int[MAXV + 5];
    
    static void buildSpf() {
        for (int i = 1; i <= MAXV; i++) spf[i] = i;
        for (int i = 2; i * i <= MAXV; i++)
            if (spf[i] == i)
                for (int j = i * i; j <= MAXV; j += i) 
                    if (spf[j] == j) spf[j] = i;
    }
    
    static Map<Integer, Long> H = new HashMap<>();
    static Random random = new Random();
    
    static long sf(int x) {
        long r = 0;
        while (x > 1) {
            int p = spf[x], c = 0;
            while (x % p == 0) { c++; x /= p; }
            if ((c & 1) == 1) {
                if (!H.containsKey(p)) H.put(p, random.nextLong());
                r ^= H.get(p);
            }
        }
        return r;
    }
    
    static int idx(int r, int p) { return r * (r + 1) / 2 + p; }
    
    public static void main(String[] args) throws IOException {
        buildSpf();
        
        BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
        int n = Integer.parseInt(br.readLine());
        
        // Read A array
        int[][] A = new int[n-1][];
        for (int r = 0; r < n-1; r++) {
            String[] parts = br.readLine().split(" ");
            A[r] = new int[r+1];
            for (int j = 0; j <= r; j++) 
                A[r][j] = Integer.parseInt(parts[j]);
        }
        
        // Read B array
        int[][] B = new int[n-1][];
        for (int r = 0; r < n-1; r++) {
            String[] parts = br.readLine().split(" ");
            B[r] = new int[r+1];
            for (int j = 0; j <= r; j++)
                B[r][j] = Integer.parseInt(parts[j]);
        }
        
        // Read C array (note the different indexing)
        int[][] C = new int[n+1][]; // Using n+1 to match 1-based indexing
        for (int r = 2; r <= n; r++) {
            String[] parts = br.readLine().split(" ");
            C[r] = new int[r];
            for (int j = 1; j < r; j++)
                C[r][j] = Integer.parseInt(parts[j-1]);
        }
        
        int tot = n * (n + 1) / 2;
        List<List<Pair>> g = new ArrayList<>(tot);
        for (int i = 0; i < tot; i++) g.add(new ArrayList<>());
        
        // Build graph edges for A
        for (int r = 0; r < n-1; r++) {
            for (int j = 0; j <= r; j++) {
                int u = idx(r, j), v = idx(r+1, j);
                long w = sf(A[r][j]);
                g.get(u).add(new Pair(v, w));
                g.get(v).add(new Pair(u, w));
            }
        }
        
        // Build graph edges for B
        for (int r = 0; r < n-1; r++) {
            for (int j = 0; j <= r; j++) {
                int u = idx(r, j), v = idx(r+1, j+1);
                long w = sf(B[r][j]);
                g.get(u).add(new Pair(v, w));
                g.get(v).add(new Pair(u, w));
            }
        }
        
        // Build graph edges for C
        for (int r = 1; r <= n; r++) {
            if (C[r] == null) continue;
            for (int j = 1; j < r; j++) {
                int u = idx(r-1, j-1), v = idx(r-1, j);
                long w = sf(C[r][j]);
                g.get(u).add(new Pair(v, w));
                g.get(v).add(new Pair(u, w));
            }
        }
        
        int[] vis = new int[tot];
        long[] f = new long[tot];
        long ans = 0;
        
        for (int s = 0; s < tot; s++) {
            if (vis[s] == 0) {
                List<Integer> comp = new ArrayList<>();
                boolean ok = true;
                Queue<Integer> q = new LinkedList<>();
                q.add(s);
                vis[s] = 1;
                f[s] = 0;
                comp.add(s);
                
                while (!q.isEmpty()) {
                    int u = q.poll();
                    for (Pair e : g.get(u)) {
                        int v = e.first;
                        long w = e.second;
                        long nv = f[u] ^ w;
                        if (vis[v] == 0) {
                            vis[v] = 1;
                            f[v] = nv;
                            q.add(v);
                            comp.add(v);
                        } else if (f[v] != nv) {
                            ok = false;
                        }
                    }
                }
                
                if (ok) {
                    Map<Long, Long> mp = new HashMap<>();
                    for (int v : comp) {
                        mp.put(f[v], mp.getOrDefault(f[v], 0L) + 1);
                    }
                    for (long m : mp.values()) {
                        ans += m * (m - 1) / 2;
                    }
                }
            }
        }
        
        System.out.println(ans);
    }
    
    static class Pair {
        int first;
        long second;
        Pair(int first, long second) {
            this.first = first;
            this.second = second;
        }
    }
}

JavaScript solution:
const { Queue } = require('datastructures-js');
const crypto = require('crypto');
const fs = require('fs');

function buildSpf(maxV) {
    const spf = new Array(maxV + 1);
    for (let i = 1; i <= maxV; i++) spf[i] = i;
    for (let i = 2; i * i <= maxV; i++) {
        if (spf[i] === i) {
            for (let j = i * i; j <= maxV; j += i) {
                if (spf[j] === j) spf[j] = i;
            }
        }
    }
    return spf;
}

function createHashFunction() {
    const cache = new Map();
    return function(x) {
        if (cache.has(x)) return cache.get(x);
        const hash = crypto.randomBytes(8).toString('hex');
        const numHash = BigInt('0x' + hash);
        cache.set(x, numHash);
        return numHash;
    };
}

function squareFree(x, spf, getHash) {
    let result = 0n;
    while (x > 1) {
        const p = spf[x];
        let c = 0;
        while (x % p === 0) {
            c++;
            x = Math.floor(x / p);
        }
        if (c % 2 === 1) {
            result ^= getHash(p);
        }
    }
    return result;
}

function main() {
    // Read all input at once for efficiency
    const input = fs.readFileSync(0, 'utf-8').trim().split(/\s+/).map(Number);
    let ptr = 0;
    const n = input[ptr++];
    
    // Read input arrays
    const A = [];
    for (let r = 0; r < n - 1; r++) {
        const row = [];
        for (let j = 0; j <= r; j++) {
            row.push(input[ptr++]);
        }
        A.push(row);
    }
    
    const B = [];
    for (let r = 0; r < n - 1; r++) {
        const row = [];
        for (let j = 0; j <= r; j++) {
            row.push(input[ptr++]);
        }
        B.push(row);
    }
    
    const C = [];
    for (let r = 1; r < n; r++) {
        const row = [];
        for (let j = 0; j < r; j++) {
            row.push(input[ptr++]);
        }
        C.push(row);
    }
    
    // Precompute square-free hashes
    const MAXV = 1000000;
    const spf = buildSpf(MAXV);
    const getHash = createHashFunction();
    
    const idx = (r, p) => r * (r + 1) / 2 + p;
    const total = n * (n + 1) / 2;
    const graph = Array.from({ length: total }, () => []);
    
    // Build graph edges with weights
    for (let r = 0; r < n - 1; r++) {
        for (let j = 0; j <= r; j++) {
            const u = idx(r, j);
            const v = idx(r + 1, j);
            const w = squareFree(A[r][j], spf, getHash);
            graph[u].push([v, w]);
            graph[v].push([u, w]);
        }
    }
    
    for (let r = 0; r < n - 1; r++) {
        for (let j = 0; j <= r; j++) {
            const u = idx(r, j);
            const v = idx(r + 1, j + 1);
            const w = squareFree(B[r][j], spf, getHash);
            graph[u].push([v, w]);
            graph[v].push([u, w]);
        }
    }
    
    for (let r = 1; r < n; r++) {
        for (let j = 0; j < r; j++) {
            const u = idx(r, j);
            const v = idx(r, j + 1);
            const w = squareFree(C[r - 1][j], spf, getHash);
            graph[u].push([v, w]);
            graph[v].push([u, w]);
        }
    }
    
    // Process components
    const visited = new Array(total).fill(false);
    const f = new Array(total).fill(0n);
    let answer = 0n;
    
    for (let s = 0; s < total; s++) {
        if (!visited[s]) {
            const component = [];
            const queue = new Queue();
            queue.enqueue(s);
            visited[s] = true;
            f[s] = 0n;
            component.push(s);
            let isValidComponent = true;
            
            while (!queue.isEmpty()) {
                const u = queue.dequeue();
                for (const [v, w] of graph[u]) {
                    const nv = f[u] ^ w;
                    if (!visited[v]) {
                        visited[v] = true;
                        f[v] = nv;
                        queue.enqueue(v);
                        component.push(v);
                    } else if (f[v] !== nv) {
                        isValidComponent = false;
                    }
                }
            }
            
            if (isValidComponent) {
                const freq = new Map();
                for (const v of component) {
                    const key = f[v].toString();
                    freq.set(key, (freq.get(key) || 0n) + 1n);
                }
                
                for (const count of freq.values()) {
                    answer += count * (count - 1n) / 2n;
                }
            }
        }
    }
    
    console.log(answer.toString());
}

main();

Python solution:
import sys
import math
import random
from collections import deque

def main():
    input = sys.stdin.read().split()
    ptr = 0

    MAXV = 10**6
    spf = list(range(MAXV + 1))
    for i in range(2, int(math.isqrt(MAXV)) + 1):
        if spf[i] == i:
            for j in range(i * i, MAXV + 1, i):
                if spf[j] == j:
                    spf[j] = i

    H = {}
    rng = random.Random()

    def square_free(x):
        r = 0
        if x == 1:
            return r
        while x > 1:
            p = spf[x]
            c = 0
            while x % p == 0:
                c += 1
                x = x // p
            if c % 2 == 1:
                if p not in H:
                    H[p] = rng.getrandbits(64)
                r ^= H[p]
        return r

    N = int(input[ptr])
    ptr += 1

    A = []
    for r in range(1, N):
        row = list(map(int, input[ptr:ptr + r]))
        ptr += r
        A.append(row)

    B = []
    for r in range(1, N):
        row = list(map(int, input[ptr:ptr + r]))
        ptr += r
        B.append(row)

    C = []
    for r in range(2, N + 1):
        row = list(map(int, input[ptr:ptr + (r - 1)]))
        ptr += (r - 1)
        C.append(row)

    tot = N * (N + 1) // 2

    def idx(r, p):
        return r * (r + 1) // 2 + p

    g = [[] for _ in range(tot)]

    for r in range(N - 1):
        for j in range(r + 1):
            u = idx(r, j)
            v = idx(r + 1, j)
            w = square_free(A[r][j])
            g[u].append((v, w))
            g[v].append((u, w))

    for r in range(N - 1):
        for j in range(r + 1):
            u = idx(r, j)
            v = idx(r + 1, j + 1)
            w = square_free(B[r][j])
            g[u].append((v, w))
            g[v].append((u, w))

    for r in range(1, N):
        for j in range(r):
            u = idx(r, j)
            v = idx(r, j + 1)
            w = square_free(C[r - 1][j])
            g[u].append((v, w))
            g[v].append((u, w))

    vis = [False] * tot
    f = [0] * tot
    ans = 0

    for s in range(tot):
        if not vis[s]:
            comp = []
            ok = True
            q = deque()
            q.append(s)
            vis[s] = True
            f[s] = 0
            comp.append(s)
            while q:
                u = q.popleft()
                for (v, w) in g[u]:
                    nv = f[u] ^ w
                    if not vis[v]:
                        vis[v] = True
                        f[v] = nv
                        q.append(v)
                        comp.append(v)
                    elif f[v] != nv:
                        ok = False
            if ok:
                cnt = {}
                for v in comp:
                    cnt[f[v]] = cnt.get(f[v], 0) + 1
                for m in cnt.values():
                    ans += m * (m - 1) // 2

    print(ans)

if __name__ == '__main__':
    main()

Ruby solution:
require 'set'

MAXV = 1_000_000

# Build smallest prime factor table
$spf = (0..MAXV).to_a
(2..Math.sqrt(MAXV).to_i).each do |i|
  if $spf[i] == i
    (i*i..MAXV).step(i) do |j|
      $spf[j] = i if $spf[j] == j
    end
  end
end

# Generate random hashes for primes
$prime_hash = {}
$rng = Random.new

def square_free(x)
  res = 0
  while x > 1
    p = $spf[x]
    cnt = 0
    while x % p == 0
      cnt += 1
      x /= p
    end
    if cnt.odd?
      $prime_hash[p] ||= $rng.rand(2**64)
      res ^= $prime_hash[p]
    end
  end
  res
end

def idx(r, p)
  r * (r + 1) / 2 + p
end

# Read input
n = gets.to_i

a = []
(n-1).times do |r|
  a << gets.split.map(&:to_i)
end

b = []
(n-1).times do |r|
  b << gets.split.map(&:to_i)
end

c = []
(n-1).times do |r|
  c << gets.split.map(&:to_i)
end

tot = n * (n + 1) / 2
g = Array.new(tot) { [] }

# Build graph
(0..n-2).each do |r|
  (0..r).each do |j|
    u = idx(r, j)
    v = idx(r+1, j)
    w = square_free(a[r][j])
    g[u] << [v, w]
    g[v] << [u, w]
  end
end

(0..n-2).each do |r|
  (0..r).each do |j|
    u = idx(r, j)
    v = idx(r+1, j+1)
    w = square_free(b[r][j])
    g[u] << [v, w]
    g[v] << [u, w]
  end
end

(1..n-1).each do |r|
  (0..r-1).each do |j|
    u = idx(r, j)
    v = idx(r, j+1)
    w = square_free(c[r-1][j])
    g[u] << [v, w]
    g[v] << [u, w]
  end
end

vis = Array.new(tot, false)
f = Array.new(tot, 0)
ans = 0

(0...tot).each do |s|
  unless vis[s]
    comp = []
    ok = true
    q = [s]
    vis[s] = true
    f[s] = 0
    comp << s
    
    until q.empty?
      u = q.shift
      g[u].each do |(v, w)|
        nv = f[u] ^ w
        unless vis[v]
          vis[v] = true
          f[v] = nv
          q << v
          comp << v
        else
          ok = false if f[v] != nv
        end
      end
    end
    
    if ok
      freq = Hash.new(0)
      comp.each { |v| freq[f[v]] += 1 }
      freq.each_value do |m|
        ans += m * (m - 1) / 2
      end
    end
  end
end

puts ans