How to include Vref in thermistor temperature calculation?
The formula to get the analog voltage from the ADC reading is:

&nbsp;&nbsp;&nbsp;&nbsp;V = Vref × reading ÷ 1024

(yes, it's 1024, not 1023). The formula to get the thermistor resistance
from the measured voltage is:

&nbsp;&nbsp;&nbsp;&nbsp;R = Rref × V ÷ (Vref − V)

where Rref is your 10&nbsp;kΩ pull-up. You can combine these
two formulas in order to get the resistance directly from the ADC
reading:

&nbsp;&nbsp;&nbsp;&nbsp;R = Rref × reading ÷ (1024 − reading)

You may notice that, while doing this simplification, Vref
cancels out. It's not an error. The readings you get with this setup are
indeed independent of your reference voltage.



Update 1: After the OP's edit, we now have a completely different
question. The question is now: how come the readings do depend on the
reference voltage that is used, whereas in theory they should not.

The answer lies most likely in the ADC calibration. If high accuracy is
a requirement, the correct formula for the measured voltage is

&nbsp;&nbsp;&nbsp;&nbsp;V = Vref ×
(reading + Eoff) ÷ (1024 + Egain)

where Eoff and Egain are the offset error and the
gain error respectively. These errors are not knows a priori, and they
are usually neglected when the accuracy requirement is low. The
datasheet only gives constraints on how big these errors can be. The
actual values depend on the specific microcontroller and can somewhat
depend on the reference voltage. The only solution to get rid of these
errors is to calibrate your own ADC. See Response of the Arduino
ADC for an example on how this could be done.



Update 2: If I had to program this in an Arduino, I would rather
simplify the calculations as much as possible before coding:


do not compute the measured analog voltage, nor the resistance, as
they are not needed
compute instead the ratio R/Rref directly from the analog
reading
rewrite the calibration polynomial as a function of
log(R/Rref) instead of log(R) (which actually means
log(R/(1&nbsp;Ω)))
optimize the evaluation of the polynomial by using
Horner's method.


With these optimizations, we get:

// Coefficients of polynomial fit of log(R/Rref) -&gt; 1/T.
const float c0 = 3.3540165e-3;
const float c1 = 2.5416075e-4;
const float c2 = -2.5959644e-6;
const float c3 = -9.3951087e-8;

// Convert ADC reading to temperature in deg. C.
static float temperature(int reading)
{
    float R_ratio = reading / (1024.0 - reading);  // = R / Rref
    float x = log(R_ratio);
    float inverse_T = c0 + x*(c1 + x*(c2 + x*c3));
    return 1/inverse_T - 273.15;
}

void loop()
{
    int reading = analogRead(A3);
    Serial.println(temperature(reading));
    deay(500);
}