Re: creating (non)zero vector

[Karim Belabas <Karim.Belabas@u-bordeaux.fr>]

* Max Alekseyev [2026-08-09 17:17]:
> I thought that creating nonzero vector (say, of all 1s) like
> 
> allocatemem(2^31);
> V = vector(10^8,i,1);
> 
> is more efficient than first creating a zero vector and then filling it
> with the required values:
> 
> allocatemem(2^31);
> V = vector(10^8);
> for(i=1,#V,V[i]=1);
> 
> In reality, the latter code works fine, but the former one results in "the
> PARI stack overflows" error, which is quite counterintuitive.
> Is such behavior expected?

Yes.

1) The integer 1 requires more space than the integer 0 (50% more),
there's not enough room in your PARI stack of size 2^31 to store the former:

? V = vector(10^6,i,0); sizebyte(V)
%8 = 24000008
? V = vector(10^6,i,1); sizebyte(V)
%1 = 32000008

(using lower values to avoid sending my old laptop to swapping hell)

So I wouldn't expect the first command to fit in your constrained stack
   10^8 * (8 bytes in a word) * (3 words for a small nonzero t_INT)
is > 2^31. With 0 it would work because 10^8 * 8 * 2 < 2^31

(The above is a white lie: I neglected 10^8 + 1 words for t_VEC
overhead, but it doesn't change the inequalities.)

2) Why does the loop version work ?

The PARI stack is a scratch space where computations occur, not the
total memory footprint of the application. It's empty when the prompt is
shown. Contents of variables, history entries, internal shenanigans
related to memory optimizations (and constraints) do not happen in the
PARI stack ... or only temporarily.

PARI has an additional "heap" of memory blocks (actually a linked list
and an AVL tree), beyond user control, to store such objets. You can
visualize its size with the debugging command 'getheap' or using \s

This is quite technical but you can read about it in the libpari user's
manual or directly

  ?? "Introduction: initializations, universal objects"@4

(found using ??? heap@ = look for all sections in GP + libpari manual
where the pattern 'heap' occurs)

Here's what it looks like in practice:

? getheap()   \\ almost empty when starting: one object, 9 words.
%1 = [1, 9]

? V = vector(10^6,i,0); sizebyte(V) / 8 \\ convert bytesize to number of words
%2 = 3000001

? getheap()   \\ increased because both V is now stored 
%3 = [4, 1000046]  

N.B. %3[2] is smaller than %2 because V is somewhat compressed (entries
equal to 0 are not duplicated).

? for(i=1,#V,V[i]=1); sizebyte(V)/8 \\ V is filled with 1s, hence bigger !
%4 = 4000001

? getheap() \\ ... but the V[i] now reside on the heap :-)
%5 = [1000006, 12000074]

Long story short: in the loop case, part of the computation does not
occur in the stack because the V[i] are computed there then moved to the
heap. In the 'vector' case, V must reside entirely in the stack at one momet
and it fails.

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/