How to divide two NUMERIC values and get a result that has the maximum number of digits the NUMERIC type supports?

The documentation of the NUMERIC type states:

Specifying: NUMERIC without any precision or scale creates an “unconstrained numeric” column in which numeric values of any length can be stored, up to the implementation limits. A column of this kind will not coerce input values to any particular scale, ...

However, when I execute the following query

SELECT (1::numeric / 3::numeric);

the result I get (using PostgreSQL version 17.2) is 0.33333333333333333333. The result of the division was coerced to a scale of 20. The documentation of the division operation does not mention any scale changes.

I assume that this is not a display issue as casting the result to a NUMERIC with bigger scale shows:

SELECT (1::numeric / 3::numeric)::numeric(26,25);
-----------------------------
 0.3333333333333333333300000

When I specifiy a scale > 20 for at least one the two operators, e.g. 1::numeric(1000,50), then the result also has more digits.

SELECT (1::numeric / 3::numeric(1000,50));
------------------------------------------------------
 0.33333333333333333333333333333333333333333333333333

According to the docs linked first, the maximum value that can be specified for the scale is 1000 which is far from the supported maximum of 16383.

My question is: How can I divide two NUMERIC values and get a result that has the maximum number of digits the NUMERIC type supports?

Unfortunately, you'll have to recompile to do that. Unconstrained numeric lets you store, process and operate on 16383 digits after decimal point just fine - here's a zero-point-16383-zeroes as a numeric constant in exponential notation:
_{demo at db<>fiddle}

create table test(x numeric);
insert into test values(0e-16383)returning scale(x),x;

Problem is, numeric/numeric division caps the result scale at a 1000 in a standard installation:

select scale(1e-16383/3e-16383)
      ,(1e-16383)/(3e-16383) as x;

You can look up the function behind numeric/numeric operator in pg_operator:

select oid::regoperator
      ,oprcode 
from pg_operator
where oprname='/'
and oprleft='numeric'::regtype
and oprright='numeric'::regtype;

Whenever documentation fails, you can always take a look at the source to see what it really does and how. The numeric_div in pg_operator.oprcode leads you to src/backend/utils/adt/numeric.c:3133:

/*
 * numeric_div() -
 *
 *  Divide one numeric into another
 */
Datum
numeric_div(PG_FUNCTION_ARGS)
{
    Numeric     num1 = PG_GETARG_NUMERIC(0);
    Numeric     num2 = PG_GETARG_NUMERIC(1);
    Numeric     res;

    res = numeric_div_opt_error(num1, num2, NULL);

    PG_RETURN_NUMERIC(res);
}

And in numeric_div_opt_error you can see it tries to estimate what scale is required for the division:

    /*
     * Select scale for division result
     */
    rscale = select_div_scale(&arg1, &arg2);

That function ends up capping the scale at the current value of NUMERIC_MAX_DISPLAY_SCALE:

    /* Select result scale */
    rscale = NUMERIC_MIN_SIG_DIGITS - qweight * DEC_DIGITS;
    rscale = Max(rscale, var1->dscale);
    rscale = Max(rscale, var2->dscale);
    rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
    rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);

To find how it came up with the 20-digit scale for your 1::numeric/3::numeric, you can follow the code that establishes the variables above, used by select_div_scale to cap the result - in your case it landed a bit above the minimum of 16. In src/include/utils /numeric.h:32 you can find NUMERIC_MAX_PRECISION, the maximum, is 1000:

#define NUMERIC_MAX_PRECISION       1000

And that's what you'd have to increase to get larger maximum scale in numeric division results. You're tweaking it at your own risk but once you do, you can just toss in a very precise zero to any calculation in order to move all the way up to the maximum scale:

SELECT (1+0e-16383)/3;