Длинная арифметика: найти сумму 1^2+2^2+3^2+...+n^2 при n > 20000

#include <iostream>
#include <iomanip>
 
const int BASIS = 1000;
const int N = 1000;
 
using namespace std;
 
void  Transform(int n, int a[], int &n_a);
void Print(int a[], int n_a);
int Max(int a, int b);
void Sum(int a[], int n_a, int b[], int n_b, int c[], int &n_c);
void Product(int a[], int n_a, int b[], int n_b, int c[], int &n_c);
 
int main(){
    
    int n;
    do{
        cout << "n = ";
        cin >> n;
    }while(n < 1 || 20000 < n);
    
    int sum[N] = {};
    int n_sum = 1;
    for (int i = 1; i <= n; ++i){
        
        int a[N] = {};
        int n_a;
        Transform(i, a, n_a);
        
        int c[N] = {};
        int n_c = 1;
        Product(a, n_a, a, n_a, c, n_c);
        
        for (int j = 0; j < N; ++j)
            if (j < n_sum)
                a[j] = sum[j];
            else
                a[j] = 0;
        n_a = n_sum;
        
        Sum(a, n_a, c, n_c, sum, n_sum);
    }
    cout << "sum = ";
    Print(sum, n_sum);
    
    long long s = 0;
    for (int i = 1; i <= n; ++i)
        s += i*i;
    cout << "s   = " << s << endl;
}
 
void  Transform(int n, int a[], int &n_a){
    
    n_a = 0;
    while (n > 0){
        
        a[n_a] = n%BASIS;
        ++n_a;
        n /= BASIS;
    }
}
 
void Print(int a[], int n_a){
    
    int h = 0;
    int buf = BASIS - 1;
    while(buf > 0){
        buf /= 10;
        ++h;
    }
    
    for (int i = n_a - 1; i >= 0; --i){
        cout << setfill('0') << setw(h) << a[i];
    }
    cout << endl;
}
 
int Max(int a, int b){
    return a > b ? a : b;
}
 
void Sum(int a[], int n_a, int b[], int n_b, int c[], int &n_c){
    
    int p = 0;
    for (int i = 0; i < Max(n_a, n_b); ++i){
        
        c[i] = (a[i] + b[i] + p)%BASIS;
        p = (a[i] + b[i] + p)/BASIS;
    }
    if (p == 0)
        n_c = Max(n_a, n_b);
    else{
        n_c = Max(n_a, n_b);
        c[n_c] = p;
        ++n_c;
    }
}
 
void Product(int a[], int n_a, int b[], int n_b, int c[], int &n_c){
    
    c[0] = 0;
    n_c = 1;
    
    for (int i = 0; i < n_a; ++i){
        
        int buf[N] = {};
        int p = 0;
        for (int j = 0; j < n_b; ++j){
            buf[j] = (a[i]*b[j] + p)%BASIS;
            p = (a[i]*b[j] + p)/BASIS;
        }
        int n_buf = n_b;
        if (p > 0){
            buf[n_buf] = p;
            ++n_buf;
        }
        
        int buf2[N] = {};
        for (int j = 0; j < n_c; ++j)
            buf2[j] = c[j]; 
        int n_buf2 = n_c;
        
        int buf1[N] = {};
        for (int j = i; j < n_buf + i; ++j)
            buf1[j] = buf[j - i];
        int n_buf1 = n_buf + i;
        
        Sum(buf1, n_buf1, buf2, n_buf2, c, n_c);
    }
}