How to solve the error ' [not a vector ]'

Question

I ran this code to find the norm of some fundamnetal units of a biqaudratic number field, but I faced following problem

for (q=5, 200,  for(p=q+1, 200,  if (isprime(p)==1 && isprime(q)==1 ,k1=bnfinit(y^2-2*p,1); k2=bnfinit(y^2-q,1);  k3=bnfinit(y^2-2*p*q,1); ep1=k1[8][5][1]; ep2=k2[8][5][1]; ep3=k3[8][5][1]; normep1=nfeltnorm(k1,ep1); normep2=nfeltnorm(k2,ep2); normep3=nfeltnorm(k3,ep3); li=[[q,p], [normep1, normep2, normep3]]; lis4=concat(lis4,[li]))))

and it works for small p and q. However, when I ran that for p and q greater than 150, it gives the following error:

First, I didn't use the flag=1 for bnf, but after adding that, still I get the same error.

Answer 1

Please, do not use indexing like ...[8][5][1] to get the fundamental units (FU). It seems that bnfinit omits FU matrix for some p and q. Instead, use the member function fu to receive FU. Please, find the example below:

> [q, p] = [23, 109];
> k = bnfinit(y^2 - 2*p*q, 1);
> k[8][5]
[;]
> k[8][5][1]  \\ you will get the error here trying to index the empty matrix.
...
incorrect type in _[_] OCcompo1 [not a vector] (t_MAT).

> k.fu
[Mod(-355285121749346859670064114879166870*y - 25157598731408198132266996072608016699, y^2 - 5014)]
> norm(k.fu[1])
1

Answer 2

1. The reason you get the error is that by default bnfinit will not try very hard to compute fundamental units exactly, since the result can be huge (its size may be exponential in the size of the discriminant), hence very expensive to compute and even overflow the possibilities of the implementation of algebraic numbers. When units get large, the simple-minded default method will fail and the bnfinit result will not contain the units but a sentinel value (an empty matrix). Thus you get an exception.

2. There are a few naive ways out, none of which I would recommend, to fix your code. First, you can check for the sentinel value, or even trap the exception with iferr, and do something different in that case. Or, as in Piotr's (more helpful) answer you can add the flag 1 to bnfinit. This uses more complicated internal code, in particular exact representations of quantities instead of faster but less reliable floating point approximations. This will still not compute large fundamental units (so using k[5][8][1] would still fail), but will allow the .fu method to do so (in which case the result is cached and will not be recomputed).

3. Why is this still not good ? Well, as mentioned above, the result might be so large that this fails again: because PARI can not represent integers large enough to write down the result you requested. And also, it's completely useless: you do not need that unit, only its norm ! PARI always computes units with bnfinit(,1) as exact algebraic numbers, but in factored form, not as expanded values. Think of the difference between a nice and compact answer like 2^2^100 and its expanded decimal representation (by the way, if you try it, it will overflow; PARI can't represent integers whose length has 2^100 bits).

4. Units in factored form are available from bnfunits, from which we can easily compute their norm since the latter is multiplicative. The function nfeltnorm(bnf, x) handles all that transparently for you and is better than norm for that reason. You can give x in any representation (as a rational number, a column vector in terms of the integer basis, a polynomial, or a factorization matrix representing a product of such elements raised to some powers). Here's the best solution:

    [q, p] = [23, 109];
    k = bnfinit(y^2 - 2*p*q, 1);
    u = bnfunits(k)[1]; \\ trailing elements are technical, for bnfisunit()
    nfeltnorm(k, u[1])  \\ u[2] is the torsion unit, so -1 
   

And the answer is 1 as in Piotr's answer.

P.S. Note that your double loop is inefficient. The data corresponding to q should be computed outside the p loop. You should use a list as well and PARI provides loops over primes natively. Here's a suggestion

  Normu(D) = my(k = bnfinit(y^2-D, 1), u = bnfunits(k)[1][1]); nfeltnorm(k, u); 
  {
    L = List();
    forprime (q = 5, 200,
      N2 = Normu(q);
      forprime(p = q+1, 200,
        N1 = Normu(2*p);
        N3 = Normu(2*p*q);
        listput(L, [[q,p], [N1, N2, N3]])));
  }