Let us get an explicit field $L$, i will ignore any restriction on $a,b$.
a, b = 7, 3 # chosen so that we have a nice discriminant below
L.<u> = NumberField(sqrt(a + b*sqrt(5)).minpoly())
print(f"L is <{L}> with discriminant {L.disc()}")
So our test situation is:
L is <Number Field in u with defining polynomial x^4 - 14*x^2 + 4> with discriminant 1600
Numerically i tried then with a precision equal to $300$ (bits):
sage: LZ = L.zeta_function(prec=300)
sage: LZ(3)
1.01897545804910243658865624113355956226224991920025923733956382609692854856490756716927418
sage: LZ(5)
1.00101170969216360222562623040290551783199837820002159868621703090382967126817104154134429
sage: LZ(7)
1.00006145739122202145223364281491843791554787649331378421635463305298859721049129811714303
Note that the code goes to pari per default. So let us see what happens in there exactly:
? lfun(x^4 - 14*x^2 + 4, 3)
%42 = 1.0189754580491024365886562411335595623
? \p 100
realprecision = 115 significant digits (100 digits displayed)
? lfun(x^4 - 14*x^2 + 4, 3)
%43 = 1.018975458049102436588656241133559562262249919200259237339563826096928548564907567169274179589141566
So an alternative is
sage: f = (sqrt(7 + 3*sqrt(5))).minpoly()
sage: pari.lfun(f, 3)
1.01897545804910
(And we can go also with similar code to have a better control of the precision.)
