sub 0~+~tau ) |~+~{0.8~<<~epsilon sub e sup 2 ( tau )~>> } over sqrt {
<<~epsilon sub i sup 2 ( tau )~>> }> ,

where <$Eepsilon sub i ( tau )> and <$Eepsilon sub e ( tau )> are the
accumulated error of the ith clock and entire clock ensemble,
respectively. The accumulated error of the entire ensemble is

<$E<<~epsilon sub e sup 2 ( tau )~>>~=~left [ sum from i=1 to n~1 over {
<<~epsilon sub i sup 2 ( tau )~>> } right ] sup {~-1}>.

Finally, the weight factor for the ith clock is calculated as

<$Ew sub i ( tau )~=~ { <<~epsilon sub e sup 2 ( tau )~>> } over {
<<~epsilon sub i sup 2 ( tau )~>> }> .

When all estimators and weight factors have been updated, the origin of
the estimation interval is shifted and the new value of t0 becomes the
old value of <$Et sub 0 ~+~ tau>.

While not entering into the above calculations, it is useful to estimate
the frequency error, since the ensemble clocks can be located some
distance from each other and become isolated for some time due to
network failures. The frequency-offset error in Ri is equivalent to the
fractional frequency yi,

<$Ey sub i~=~{ nu sub i~-~nu sub I } over nu sub I>

measured between the ith timescale and the standard timescale I.
Temporarily dropping the subscript i for clarity, consider a sequence of
N independent frequency-offset samples <$Ey(j)~ (j~=~1,~2,~... ,~N)>
where the interval between samples is uniform and equal to T. Let
<$Etau> be the nominal interval over which these samples are averaged.
The Allan variance <$Esigma sub y sup 2 ( N,~T,~tau )> [ALL74a] is
defined as
<$E<< sigma sub y sup 2 ( N,~T,~tau )~>>~=~<< ~ 1 over { N~-~1 }~ left [
sum from j=1 to N~y (j) sup 2~-~1 over N~left ( sum from j=1 to N~y(j)
right ) sup 2 right ]~>>> ,

A particularly useful formulation is <$EN~=~2> and <$ET~=~tau>:

<$E<< sigma sub y sup 2 (N~=~2,~T~=~tau ,~tau )>>~==~sigma sub y sup 2 (
tau )~=~<< {[y(j~+~1)~-~y(j)] sup 2 } over 2 >>> ,

so that

<$Esigma sub y sup 2 ( tau )~=~1 over {2(N~-~1)}sum from { j = 1 } to
{n-1 }~[y(j~+~1)~-~y(j)] sup 2> .

While the Allan variance has found application when estimating errors in
ensembles of cesium clocks, its application to NTP is limited due to the

computation and storage burden. As described in the next section, it is
possible to estimate errors with some degree of confidence using normal
byproducts of NTP processing algorithms.

Application to NTP

The NTP clock model is somewhat less complex than the general model
described above. For instance, at the present level of development it is
not necessary to separately estimate the time and frequency of all peer
clocks, only the time and frequency of the local clock. If the
timekeeping reference is the local clock itself, then the offsets
available in the peer.offset peer variables can be used directly for the
<$ET sub ij> quantities above. In addition, the NTP local-clock model
incorporates a type-II phase-locked loop, which itself reliably
estimates frequency errors and corrects accordingly. Thus, the
requirement for estimating frequency is entirely eliminated.

There remains the problem of how to determine a robust and easily
computable error estimate <$Eepsilon sub i>. The method described above,
although analytically justified, is most difficult to implement.
Happily, as a byproduct of the NTP clock-filter algorithm, a useful
error estimate is available in the form of the dispersion. As described
in the NTP specification, the dispersion includes the absolute value of
the weighted average of the offsets between the chosen offset sample and
the <$En~-~1> other samples retained for selection. The effectiveness of
this estimator was compared with the above estimator by simulation using
observed timekeeping data and found to give quite acceptable results.

The NTP clock-combining algorithm can be implemented with only minor
modifications to the algorithms as described in the NTP specification.
Although elsewhere in the NTP specification the use of general-purpose
multiply/divide routines has been successfully avoided, there seems to
be no way to avoid them in the clock-combining algorithm. However, for
best performance the local-clock algorithm described elsewhere in this
document should be implemented as well, since the combining algorithms
result in a modest increase in phase noise which the revised local-clock
algorithm is designed to suppress.

Clock-Combining Procedure

The result of the NTP clock-selection procedure is a set of survivors
(there must be at least one) that represent truechimers, or correct
clocks. When clock combining is not implemented, one of these peers,
chosen as the most likely candidate, becomes the synchronization source
and its computed offset becomes the final clock correction.
Subsequently, the system variables are adjusted as described in the NTP
clock-update procedure. When clock combining is implemented, these
actions are unchanged, except that the final clock correction is
computed by the clock-combining procedure.

The clock-combining procedure is called from the clock-select procedure.
It constructs from the variables of all surviving peers the final clock
correction <$ETHETA>. The estimated error required by the algorithms
previously described is based on the synchronization distance <$ELAMBDA>
computed by the distance procedure, as defined in the NTP specification.
The reciprocal of <$ELAMBDA> is the weight of each clock-offset
contribution to the final clock correction. The following pseudo-code
describes the procedure.

begin clock-combining procedure
        <$Etemp1~<<-~0>;
        <$Etemp2~<<-~0>;
        for (each peer remaining on the candidate list)         /* scan
all survivors */
                <$ELAMBDA~<<-~roman distance (peer)>;
                <$Etemp~<<-~1 over roman
{peer.stratum~times~NTP.MAXDISPERSE~+~LAMBDA }>;
                <$Etemp1~<<-~temp1~+~temp>;             /* update weight
and offset */
                <$Etemp2~<<-~temp2~+~temp~times~roman peer.offset>;
                endif;
        <$ETHETA~<<-~temp2 over temp1>;                                 
/* compute final correction */
        end clock-combining procedure;

The value <$ETHETA> is the final clock correction used by the local-
clock procedure to adjust the clock.

Appendix G. Computer Clock Modelling and Analysis

A computer clock includes some kind of reference oscillator, which is
stabilized by a quartz crystal or some other means, such as the power
grid. Usually, the clock includes a prescaler, which divides the
oscillator frequency to a standard value, such as 1 MHz or 100 Hz, and a
counter, implemented in hardware, software or some combination of the
two, which can be read by the processor. For systems intended to be
synchronized to an external source of standard time, there must be some
means to correct the phase and frequency by occasional vernier
adjustments produced by the timekeeping protocol. Special care is
necessary in all timekeeping system designs to insure that the clock
indications are always monotonically increasing; that is, system time
never <169>runs backwards.<170>

Computer Clock Models

The simplest computer clock consists of a hardware latch which is set by
overflow of a hardware counter or prescaler, and causes a processor
interrupt or tick. The latch is reset when acknowledged by the
processor, which then increments the value of a software clock counter.
The phase of the clock is adjusted by adding periodic corrections to the

counter as necessary. The frequency of the clock can be adjusted by
changing the value of the increment itself, in order to make the clock
run faster or slower. The precision of this simple clock model is
limited to the tick interval, usually in the order of 10 ms; although in
some systems the tick interval can be changed using a kernel variable.

This software clock model requires a processor interrupt on every tick,
which can cause significant overhead if the tick interval is small, say
in the order less 1 ms with the newer RISC processors. Thus, in order to
achieve timekeeping precisions less than 1 ms, some kind of hardware
assist is required. A straightforward design consists of a voltage-
controlled oscillator (VCO), in which the frequency is controlled by a
buffered, digital/analog converter (DAC). Under the assumption that the
VCO tolerance is 10-4 or 100 parts-per-million (ppm) (a reasonable value
for inexpensive crystals) and the precision required is 100 <$Emu roman
s> (a reasonable goal for a RISC processor), the DAC must include at
least ten bits.

A design sketch of a computer clock constructed entirely of hardware
logic components is shown in Figure 10a<$&fig10>. The clock is read by
first pulsing the read signal, which latches the current value of the
clock counter, then adding the contents of the clock-counter latch and a
64-bit clock-offset variable, which is maintained in processor memory.
The clock phase is adjusted by adding a correction to the clock-offset
variable, while the clock frequency is adjusted by loading a correction
to the DAC latch. In principle, this clock model can be adapted to any
precision by changing the number of bits of the prescaler or clock
counter or changing the VCO frequency. However, it does not seem useful
to reduce precision much below the minimum interrupt latency, which is
in the low microseconds for a modern RISC processor.

If it is not possible to vary the oscillator frequency, which might be
the case if the oscillator is an external frequency standard, a design
such as shown in Figure 10b may be used. It includes a fixed-frequency
oscillator and prescaler which includes a dual-modulus swallow counter
that can be operated in either divide-by-10 or divide-by-11 modes as
controlled by a pulse produced by a programmable divider (PD). The PD is
loaded with a value representing the frequency offset. Each time the
divider overflows a pulse is produced which switches the swallow counter
from the divide-by-10 mode to the divide-by-11 mode and then back again,
which in effect <169>swallows<170> or deletes a single pulse of the
prescaler pulse train.

The pulse train produced by the prescaler is controlled precisely over a
small range by the contents of the PD. If programmed to emit pulses at a
low rate, relatively few pulses are swallowed per second and the
frequency counted is near the upper limit of its range; while, if
programmed to emit pulses at a high rate, relatively many pulses are
swallowed and the frequency counted is near the lower limit. Assuming
some degree of freedom in the choice of oscillator frequency and

prescaler ratios, this design can compensate for a wide range of
oscillator frequency tolerances.

In all of the above designs it is necessary to limit the amount of
adjustment incorporated in any step to insure that the system clock
indications are always monotonically increasing. With the software clock
model this is assured as long as the increment is never negative. When
the magnitude of a phase adjustment exceeds the tick interval (as
corrected for the frequency adjustment), it is necessary to spread the
adjustments over mulitple tick intervals. This strategy amounts to a
deliberate frequency offset sustained for an interval equal to the total
number of ticks required and, in fact, is a feature of the Unix clock
model discussed below.

In the hardware clock models the same considerations apply; however, in
these designs the tick interval amounts to a single pulse at the
prescaler output, which may be in the order of 1 ms. In order to avoid
decreasing the indicated time when a negative phase correction occurs,
it is necessary to avoid modifying the clock-offset variable in
processor memory and to confine all adjustments to the VCO or prescaler.
Thus, all phase adjustments must be performed by means of programmed
frequency adjustments in much the same way as with the software clock
model described previously.

It is interesting to conjecture on the design of a processor assist that
could provide all of the above functions in a compact, general-purpose
hardware interface. The interface might consist of a multifunction timer
chip such as the AMD 9513A, which includes five 16-bit counters, each
with programmable load and hold registers, plus an onboard crystal
oscillator, prescaler and control circuitry. A 48-bit hardware clock
counter would utilize three of the 16-bit counters, while the fourth
would be used as the swallow counter and the fifth as the programmable
divider. With the addition of a programmable-array logic device and
architecture-specific host interface, this compact design could provide
all the functions necessary for a comprehensive timekeeping system.

The Fuzzball Clock Model

The Fuzzball clock model uses a combination of hardware and software to
provide precision timing with a minimum of software and processor
overhead. The model includes an oscillator, prescaler and hardware
counter; however, the oscillator frequency remains constant and the
hardware counter produces only a fraction of the total number of bits
required by the clock counter. A typical design uses a 64-bit software
clock counter and a 16-bit hardware counter which counts the prescaler
output. A hardware-counter overflow causes the processor to increment
the software counter at the bit corresponding to the frequency <$E2 sup
N f sub p>, where N is the number of bits of the hardware counter and fp
is the counted frequency at the prescaler output. The processor reads
the clock counter by first generating a read pulse, which latches the

hardware counter, and then adding its contents, suitably aligned, to the
software counter.

The Fuzzball clock can be corrected in phase by adding a (signed)
adjustment to the software clock counter. In practice, this is done only
when the local time is substantially different from the time indicated
by the clock and may violate the monotonicity requirement. Vernier phase
adjustments determined in normal system operation must be limited to no
more than the period of the counted frequency, which is 1 kHz for LSI-11
Fuzzballs. In the Fuzzball model these adjustments are performed at
intervals of 4 s, called the adjustment interval, which provides a
maximum frequency adjustment range of 250 ppm. The adjustment
opportunities are created using the interval-timer facility, which is a
feature of most operating systems and independent of the time-of-day
clock. However,  if the counted frequency is increased from 1 kHz to 1
MHz for enhanced precision, the adjustment frequency must be increased
to 250 Hz, which substantially increases processor overhead. A modified
design suitable for high precision clocks is presented in the next
section.

In some applications involving the Fuzzball model, an external pulse-
per-second (pps) signal is available from a reference source such as a
cesium clock or GPS receiver. Such a signal generally provides much
higher accuracy than the serial character string produced by a radio
timecode receiver, typically in the low nanoseconds. In the Fuzzball
model this signal is processed by an interface which produces a hardware
interrupt coincident with the arrival of the pps pulse. The processor
then reads the clock counter and computes the residual modulo 1 s of the
clock counter. This represents the local-clock error relative to the pps
signal.

Assuming the seconds numbering of the clock counter has been determined
by a reliable source, such as a timecode receiver, the offset within the
second is determined by the residual computed above. In the NTP local-
clock model the timecode receiver or NTP establishes the time to within
<F128M>?<F255D>128 ms, called the aperture, which guarantees the seconds
numbering to within the second. Then, the pps residual can be used
directly to correct the oscillator, since the offset must be less than
the aperture for a correctly operating timecode receiver and pps signal.

The above technique has an inherent error equal to the latency of the
interrupt system, which in modern RISC processors is in the low tens of
microseconds. It is possible to improve accuracy by latching the
hardware time-of-day counter directly by the pps pulse and then reading
the counter in the same way as usual. This requires additional circuitry
to prioritize the pps signal relative to the pulse generated by the
program to latch the counter.
The Unix Clock Model

The Unix 4.3bsd clock model is based on two system calls, settimeofday

and adjtime, together with two kernel variables tick and tickadj. The
settimeofday call unceremoniously resets the kernel clock to the value
given, while the adjtime call slews the kernel clock to a new value
numerically equal to the sum of the present time of day and the (signed)
argument given in the adjtime call. In order to understand the behavior
of the Unix clock as controlled by the Fuzzball clock model described
above, it is helpful to explore the operations of adjtime in more
detail.

The Unix clock model assumes an interrupt produced by an onboard
frequency source, such as the clock counter and prescaler described
previously, to deliver a pulse train in the 100-Hz range. In priniciple,
the power grid frequency can be used, although it is much less stable
than a crystal oscillator. Each interrupt causes an increment called
tick to be added to the clock counter. The value of the increment is
chosen so that the clock counter, plus an initial offset established by
the settimeofday call, is equal to the time of day in microseconds.

The Unix clock can actually run at three different rates, one
corresponding to tick, which is related to the intrinsic frequency of
the particular oscillator used as the clock source, one to
<$Etick~+~tickadj> and the third to <$Etick~-~tickadj>. Normally the
rate corresponding to tick is used; but, if adjtime is called, the
argument <$Edelta> given is used to calculate an interval <$EDELTA
t~=~delta~tick over tickadj> during which one or the other of the two
rates are used, depending on the sign of <$Edelta>. The effect is to
slew the clock to a new value at a small, constant rate, rather than
incorporate the adjustment all at once, which could cause the clock to
be set backward. With common values of <$Etick~=~10> ms and
<$Etickadj~=~5~mu roman s>, the maximum frequency adjustment range is
<$E+- tickadj over tick~=~+- {5~roman x~10 sup -6} over {10 sup -2}> or
<F128M>æ<F255D>500 ppm. Even larger ranges may be required in the case
of some workstations (e.g., SPARCstations) with extremely poor component
tolerances.

When precisions not less than about 1 ms are required, the Fuzzball
clock model can be adapted to the Unix model by software simulation, as
described in Section 5 of the NTP specification, and calling adjtime at
each adjustment interval. When precisions substantially better than this
are required, the hardware microsecond clock provided in some
workstations can be used together with certain refinements of the
Fuzzball and Unix clock models. The particular design described below is
appropriate for a maximum oscillator frequency tolerance of 100 ppm
(.01%), which can be obtained using a relatively inexpensive quartz
crystal oscillator, but is readily scalable for other assumed
tolerances.

The clock model requires the capability to slew the clock frequency over
the range <F128M>æ<F255D>100 ppm with an intrinsic oscillator frequency
error as great as <F128M>æ<F255D>100 ppm. Figure 11<$&fig11> shows the

timing relationships at the extremes of the requirements envelope.
Starting from an assumed offset of nominal zero and an assumed error of
+100 ppm at time 0 s, the line AC shows how the uncorrected offset grows
with time. Let <$Esigma> represent the adjustment interval and a the
interval AB, in seconds, and let r be the slew, or rate at which
corrections are introduced, in ppm. For an accuracy specification of 100
<$Emu roman s>, then

<$Esigma~<<=~{100~mu roman s} over {100~roman ppm}~+~{100~mu roman s}
over {(r~-~100)~roman ppm}~=~r over {r~-~100}> .
The line AE represents the extreme case where the clock is to be steered
<F128M>-<F255D>100 ppm. Since the slew must be complete at the end of
the adjustment interval,

<$Ea~<<=~{(r~-~200)~sigma} over r>.

These relationships are satisfied only if <$Er~>>~200~roman ppm> and
<$Esigma~<<~2~roman s>. Using <$Er~=~300~roman ppm> for convenience,
<$Esigma~=~1.5~roman s> and <$Ea~<<=~0.5~roman s>. For the Unix clock
model with <$Etick~=~10~roman ms>, this results in the value of
<$Etickadj~=~3~mu roman s>.

One of the assumptions made in the Unix clock model is that the period
of adjustment computed in the adjtime call must be completed before the
next call is made. If not, this results in an error message to the
system log. However, in order to correct for the intrinsic frequency
offset of the clock oscillator, the NTP clock model requires adjtime to
be called at regular adjustment intervals of <$Esigma> s. Using the
algorithms described here and the architecture constants in the NTP
specification, these adjustments will always complete.

Mathematical Model of the NTP Logical Clock

The NTP logical clock can be represented by the feedback-control model
shown in Figure 12<$&fig12>. The model consists of an adaptive-
parameter, phase-lock loop (PLL), which continuously adjusts the phase
and frequency of an oscillator to compensate for its intrinsic jitter,
wander and drift. A mathematical analysis of this model developed along
the lines of [SMI86] is presented in following sections, along with a
design example useful for implementation guidance in operating-systems
environments such as Unix and Fuzzball. Table 9<$&tab9> summarizes the
quantities ordinarily treated as variables in the model. By convention,
<$Ev> is used for internal loop variables, <$Etheta> for phase,
<$Eomega> for frequency and <$Etau> for time. Table 10<$&tab10>
summarizes those quantities ordinarily fixed as constants in the model.
Note that these are all expressed as a power of two in order to simplify
the implementation.

In Figure 12 the variable <$Etheta sub r> represents the phase of the
reference signal and <$Etheta sub o> the phase of the voltage-controlled

oscillator (VCO). The phase detector (PD) produces a voltage <$Ev sub d>
representing the phase difference <$Etheta sub r~-~theta sub o> . The
clock filter functions as a tapped delay line, with the output <$Ev sub
s> taken at the tap selected by the clock-filter algorithm described in
the NTP specification. The loop filter, represented by the equations
given below, produces a VCO correction voltage <$Ev sub c>, which
controls the oscillator frequency and thus the phase <$Etheta sub o>.

The PLL behavior is completely determined by its open-loop, Laplace
transfer function <$EG(s)> in the s domain. Since both frequency and
phase corrections are required, an appropriate design consists of a
type-II PLL, which is defined by the function

<$EG(s)~=~{omega sub c sup 2} over {tau sup 2 s sup 2}~( 1 ~+~{tau s}
over omega sub z )> ,

where <$Eomega sub c> is the crossover frequency (also called loop
gain), <$Eomega sub z> is the corner frequency (required for loop
stability) and <$Etau> determines the PLL time constant and thus the
bandwidth. While this is a first-order function and some improvement in
phase noise might be gained from a higher-order function, in practice
the improvement is lost due to the effects of the clock-filter delay, as
described below.

The open-loop transfer function <$EG(s)> is constructed by breaking the
loop at point a on Figure 12 and computing the ratio of the output phase
<$Etheta sub o (s)> to the reference phase <$Etheta sub r (s)>. This
function is the product of the individual transfer functions for the
phase detector, clock filter, loop filter and VCO. The phase detector
delivers a voltage <$Ev sub d (t)~=~ theta sub r (t)>, so its transfer
function is simply <$EF sub d (s)~=~1>, expressed in V/rad. The VCO
delivers a frequency change <$EDELTA omega ~=~{roman d~theta sub o (t)}
over {roman dt}~=~alpha {v sub c (t)}>, where <$Ealpha> is the VCO gain
in rad/V-sec and <$Etheta sub o (t)~=~alpha~int v sub c (t)~dt>. Its
transfer function is the Laplace transform of the integral, <$EF sub o
(s)~=~alpha over s>, expressed in rad/V. The clock filter contributes a
stochastic delay due to the clock-filter algorithm; but, for present
purposes, this delay will be assumed a constant T, so its transfer
function is the Laplace transform of the delay, <$EF sub s (s)~=~e sup
{- Ts}>. Let <$EF(s)> be the transfer function of the loop filter, which
has yet to be determined. The open-loop transfer function <$EG(s)> is
the product of these four individual transfer functions:

<$EG(s)~=~{omega sub c sup 2} over {tau sup 2 s sup 2}~( 1 ~+~{tau s}
over omega sub z )~=~F sub d (s) F sub s (s) F(s) F sub o (s)~=~1e sup
{-Ts}~F(s)~alpha over s> .

For the moment, assume that the product <$ETs> is small, so that <$Ee
sup {-Ts}~approx ~1>. Making the following substitutions,

<$Eomega sub c sup 2~=~alpha over { K sub f}~~~~> and <$E~~~~omega sub
z~=~K sub g over {K sub f}>

and rearranging yields

<$EF(s)~=~1 over {K sub g~tau}~+~1 over {K sub f~tau sup 2 s }> ,

which corresponds to a constant term plus an integrating term scaled by
the PLL time constant <$Etau>. This form is convenient for
implementation as a sampled-data system, as described later.

With the parameter values given in Table 10, the Bode plot of the open-
loop transfer function <$EG(s)> consists of a <196>12 dB/octave line
which intersects the 0-dB baseline at <$Eomega sub c~=~2 sup -12> rad/s,
together with a +6 dB/octave line at the corner frequency <$Eomega sub
z~=~2 sup -14> rad/s. The damping factor <$Ezeta~=~omega sub c over {2
omega sub z}~=~2> suggests the PLL will be stable and have a large phase
margin together with a low overshoot. However, if the clock-filter delay
T is not small compared to the loop delay, which is approximately equal
to <$E1 over omega sub c>, the above analysis becomes unreliable and the
loop can become unstable. With the values determined as above, T is
ordinarily small enough to be neglected.

Assuming the output is taken at <$Ev sub s>, the closed-loop transfer
function <$EH(s)> is

<$EH(s)~==~{v sub s (s)} over {theta sub r (s)}~=~{F sub d (s) e sup {-
Ts}} over {1~+~G(s)}> .

If only the relative response is needed and the clock-filter delay can
be neglected, <$EH(s)> can be written

<$EH(s)~=~1 over {1~+~G(s)}~=~s sup 2 over {s sup 2~+~omega sub c sup 2
over {omega sub z~tau} s~+~omega sub c sup 2 over tau sup 2}> .

For some input function <$EI(s)> the output function <$EI(s)H(s)> can be
inverted to find the time response. Using a unit-step input <$EI(s)~=~1
over s> and the values determined as above, This yields a PLL risetime
of about 52 minutes, a maximum overshoot of about 4.8 percent in about
1.7 hours and a settling time to within one percent of the initial
offset in about 8.7 hours.
Parameter Management

A very important feature of the NTP PLL design is the ability to adapt
its behavior to match the prevailing stability of the local oscillator
and transmission conditions in the network. This is done using the
<$Ealpha> and <$Etau> parameters shown in Table 10. Mechanisms for doing
this are described in following sections.

Adjusting VCO Gain (<$Ebold alpha>)

The <$Ealpha> parameter is determined by the maximum frequency tolerance
of the local oscillator and the maximum jitter requirements of the
timekeeping system. This parameter is usually an architecture constant
and fixed during system operation. In the implementation model described
below, the reciprocal of <$Ealpha>, called the adjustment interval
<$Esigma>, determines the time between corrections of the local clock,
and thus the value of <$Ealpha>. The value of <$Esigma> can be
determined by the following procedure.

The maximum frequency tolerance for board-mounted, uncompensated quartz-
crystal oscillators is probably in the range of 10-4 (100 ppm). Many if
not most Internet timekeeping systems can tolerate jitter to at least
the order of the intrinsic local-clock resolution, called precision in
the NTP specification, which is commonly in the range from one to 20 ms.
Assuming 10-3 s peak-to-peak as the most demanding case, the interval
between clock corrections must be no more than <$Esigma~=~10 sup -3 over
{2 roman~x~10 sup -4}~=~5> sec. For the NTP reference model
<$Esigma~=~4> sec in order to allow for known features of the Unix
operating-system kernel. However, in order to support future anticipated
improvements in accuracy possible with faster workstations, it may be
useful to decrease <$Esigma> to as little as one-tenth the present
value.

Note that if <$Esigma> is changed, it is necessary to adjust the
parameters <$EK sub f> and <$EK sub g> in order to retain the same loop
bandwidth; in particular, the same <$Eomega sub c> and <$Eomega sub z>.
Since <$Ealpha> varies as the reciprocal of <$Esigma>, if <$Esigma> is
changed to something other than 22, as in Table 10, it is necessary to
divide both <$EK sub f> and <$EK sub g> by <$Esigma over 4> to obtain
the new values.

Adjusting PLL Bandwidth (<$Ebold tau>)

A key feature of the type-II PLL design is its capability to compensate
for the intrinsic frequency errors of the local oscillator. This
requires a initial period of adaptation in order to refine the frequency
estimate (see later sections of this appendix). The <$Etau> parameter
determines the PLL time constant and thus the loop bandwidth, which is
approximately equal to <$E{omega sub c} over tau>. When operated with a
relatively large bandwidth (small <$Etau>), as in the analysis above,
the PLL adapts quickly to changes in the input reference signal, but has
poor long term stability. Thus, it is possible to accumulate substantial
errors if the system is deprived of the reference signal for an extended
period. When operated with a relatively small bandwidth (large <$Etau>),
the PLL adapts slowly to changes in the input reference signal, and may
even fail to lock onto it. Assuming the frequency estimate has
stabilized, it is possible for the PLL to coast for an extended period
without external corrections and without accumulating significant error.

In order to achieve the best performance without requiring individual
tailoring of the loop bandwidth, it is necessary to compute each value
of <$Etau> based on the measured values of offset, delay and dispersion,
as produced by the NTP protocol itself. The traditional way of doing
this in precision timekeeping systems based on cesium clocks, is to
relate <$Etau> to the Allan variance, which is defined as the mean of
the first-order differences of sequential samples measured during a
specified interval <$Etau>,

<$Esigma sub y sup 2 ( tau )~=~1 over {2(N~-~1)}sum from { i = 1 } to
{N-1 }~[y(i~+~1)~-~y(i)] sup 2> ,

where y is the fractional frequency measured with respect to the local
timescale and N is the number of samples.

In the NTP local-clock model the Allan variance (called the compliance,
h in Table 11) is approximated on a continuous basis by exponentially
averaging the first-order differences of the offset samples using an
empirically determined averaging constant. Using somewhat ad-hoc mapping
functions determined from simulation and experience, the compliance is
manipulated to produce the loop time constant and update interval.

The NTP Clock Model

The PLL behavior can also be described by a set of recurrence equations,
which depend upon several variables and constants. The variables and
parameters used in these equations are shown in Tables 9, 10 and
11<$&tab11>. Note the use of powers of two, which facilitates
implementation using arithmetic shifts and avoids the requirement for a
multiply/divide capability.

A capsule overview of the design may be helpful in understanding how it
operates. The logical clock is continuously adjusted in small increments
at fixed intervals of <$Esigma>. The increments are determined while
updating the variables shown in Tables 9 and 11, which are computed from
received NTP messages as described in the NTP specification. Updates
computed from these messages occur at discrete times as each is
received. The intervals <$Emu> between updates are variable and can
range up to about 17 minutes. As part of update processing the
compliance h is computed and used to adjust the PLL time constant
<$Etau>. Finally, the update interval <$Erho> for transmitted NTP
messages is determined as a fixed multiple of <$Etau>.

Updates are numbered from zero, with those in the neighborhood of the
ith update shown in Figure 13<$&fig13>. All variables are initialized at
<$Ei~=~0> to zero, except the time constant <$Etau (0)~=~tau>, poll
interval <$Emu (0)~=~tau> (from Table 10) and compliance <$Eh (0)~=~K
sub s>. After an interval <$Emu (i)> (<$Ei~>>~0>) from the previous
update the ith update arrives at time <$Et(i)> including the time
offset <$Ev sub s (i)>. Then, after an interval <$Emu (i~+~1)> the

<$Ei+1 roman th> update arrives at time <$Et(i~+~1)> including the time
offset <$Ev sub s (i~+~1)>. When the update <$Ev sub s (i)> is received,
the frequency error <$Ef(i~+~1)> and phase error <$Eg(i~+~1)> are
computed:

<$Ef(i~+~1)~=~f(i)~+~{mu (i) v sub s (i)} over {tau (i) sup 2 }>
,<$E~~~~~g(i~+~1)~=~{v sub s (i)} over {tau (i)}> .

Note that these computations depend on the value of the time constant
<$Etau (i)> and poll interval <$Emu (i)> previously computed from the
<$Ei-1 roman th> update. Then, the time constant for the next interval
is computed from the current value of the compliance <$Eh(i)>

<$Etau (i~+~1)~=~roman max [K sub s~-~|~h(i)|,~1]> .

Next, using the new value of <$Etau>, called <$Etau prime> to avoid
confusion, the poll interval is computed

<$Erho (i~+~1)~=~K sub u~tau prime> .

Finally, the compliance <$Eh(i~+~1)> is recomputed for use in the <$Ei+1
roman th> update:
<$Eh(i~+~1)~=~h(i)~+~{K sub t~tau prime v sub s (i)~-~h(i) }over K sub
h> .

The factor <$Etau prime> in the above has the effect of adjusting the
bandwidth of the PLL as a function of compliance. When the compliance
has been low over some relatively long period, <$Etau prime> is
increased and the bandwidth is decreased. In this mode small timing
fluctuations due to jitter in the network are suppressed and the PLL
attains the most accurate frequency estimate. On the other hand, if the
compliance becomes high due to greatly increased jitter or a systematic
frequency offset, <$Etau prime> is decreased and the bandwidth is
increased. In this mode the PLL is most adaptive to transients which can
occur due to reboot of the system or a major timing error. In order to
maintain optimum stability, the poll interval <$Erho> is varied directly
with <$Etau>.

A model suitable for simulation and parameter refinement can be
constructed from the above recurrence relations. It is convenient to set
the temporary variable <$Ea~=~g(i~+~1)>. At each adjustment interval
<$Esigma> the quantity <$Ea over K sub g~+~{f(i~+~1)} over K sub f> is
added to the local-clock phase and the quantity <$Ea over K sub g> is
subtracted from a. For convenience, let n be the greatest integer in
<$E{mu (i)} over sigma>; that is, the number of adjustments that occur
in the ith interval. Thus, at the end of the ith interval just before
the <$Ei+1 roman th> update, the VCO control voltage is:

<$Ev sub c (i~+~1)~=~v sub c (i)~+~{[1~-~(1~-~1 over K sub g ) sup n
]}~{g(i~+~1)} ~+~n over {K sub f }~{ f(i~+~1)}~.>

Detailed simulation of the NTP PLL with the values specified in Tables
9, 10 and 11 and the clock filter described in the NTP specification
results in the following characteristics: For a 100-ms phase change the
loop reaches zero error in 39 minutes, overshoots 7 ms at 54 minutes and
settles to less than 1 ms in about six hours. For a 50-ppm frequency
change the loop reaches 1 ppm in about 16 hours and 0.1 ppm in about 26
hours. When the magnitude of correction exceeds a few milliseconds or a
few ppm for more than a few updates, the compliance begins to increase,
which causes the loop time constant and update interval to decrease.
When the magnitude of correction falls below about 0.1 ppm for a few
hours, the compliance begins to decrease, which causes the loop time
constant and update interval to increase. The effect is to provide a
broad capture range exceeding 4 s per day, yet the capability to resolve
oscillator skew well below 1 ms per day. These characteristics are
appropriate for typical crystal-controlled oscillators with or without
temperature compensation or oven control.

Appendix H. Analysis of Errors and Correctness Principles

Introduction

This appendix contains an analysis of errors arising in the generation
and processing of NTP timestamps and the determination of delays and
offsets. It establishes error bounds as a function of measured roundtrip
delay and dispersion to the root (primary reference source) of the
synchronization subnet. It also discusses correctness assertions about
these error bounds and the time-transfer, filtering and selection
algorithms used in NTP.

The notation <$Ew~=~[u,~v]> in the following describes the interval in
which u is the lower limit and v the upper limit, inclusive. Thus,
<$Eu~=~min (w)~<<=~v~=~max (w)>, and for scalar a,
<$Ew~+~a~=~[u~+~a,~v~+~a]>. Table 12<$&tab12> shows a summary of other
notation used in the analysis. The notation <$E<<~x~>>> designates the
(infinite) average of x, which is usually approximated by an exponential
average, while the notation <$Ex hat> designates an estimator for x. The
lower-case Greek letters <$Etheta>, <$Edelta> and <$Eepsilon> are used
to designate measurement data for the local clock to a peer clock, while
the upper-case Greek letters <$ETHETA>, <$EDELTA> and <$EEPSILON> are
used to designate measurement data for the local clock relative to the
primary reference source at the root of the synchronization subnet.
Exceptions will be noted as they arise.

Timestamp Errors

The standard second (1 s) is defined as <169>9,192,631,770 periods of
the radiation corresponding to the transition between the two hyperfine
levels of the ground state of the cesium-133 atom<170> [ALL74b], which
implies a granularity of about 1.1x10-10 s. Other intervals can be

determined as rational multiples of 1 s. While NTP time has an inherent
resolution of about 2.3x10-10 s, local clocks ordinarily have
resolutions much worse than this, so the inherent error in resolving NTP
time relative to the 1 s can be neglected.

In this analysis the local clock is represented by a counter/divider
which increments at intervals of s seconds and is driven by an
oscillator which operates at frequency <$Ef sub c~=~n over s> for some
integer n. A timestamp <$ET(t)> is determined by reading the clock at an
arbitrary time t (the argument t will be usually omitted for
conciseness). Strictly speaking, s is not known exactly, but can be
assumed bounded from above by the maximum reading error <$Erho>. The
reading error itself is represented by the random variable r bounded by
the interval <$E[-~rho ,~0]>, where <$Erho> depends on the particular
clock implementation. Since the intervals between reading the same clock
are almost always independent of and much larger than s, successive
readings can be considered independent and identically distributed. The
frequency error of the clock oscillator is represented by the random
variable f bounded by the interval <$E[-~phi ,~phi ]>, where <$Ephi>
represents the maximum frequency tolerance of the oscillator throughout
its service life. While f for a particular clock is a random variable
with respect to the population of all clocks, for any one clock it
ordinarily changes only slowly with time and can usually be assumed a
constant for that clock. Thus, an NTP timestamp can be represented by
the random variable T:

<$ET~=~t~+~r~+~f tau> ,

where t represents a clock reading, <$Etau> represents the time interval
since this reading and minor approximations inherent in the measurement
of <$Etau> are neglected.

In order to assess the nature and expected magnitude of timestamp errors
and the calculations based on them, it is useful to examine the
characteristics of the probability density functions (pdf) <$Ep sub r
(x)> and <$Ep sub f (x)> for r and f respectively. Assuming the clock
reading and counting processes are independent, the pdf for r is uniform
over the interval <$E[-~rho ,~0]>. With conventional manufacturing
processes and temperature variations the pdf for f can be approximated
by a truncated, zero-mean Gaussian distribution with standard deviation
<$Esigma>. In conventional manufacturing processes <$Esigma> is
maneuvered so that the fraction of samples rejected outside the interval
<$E[-~phi ,~phi ]> is acceptable. The pdf for the total timestamp error
<$Eepsilon (x)> is thus the sum of the r and f contributions, computed
as

<$Eepsilon (x)~ =~ int~from {- inf } to inf p sub r (t) p sub f (x~-~t)
d t> ,

which appears as a bell-shaped curve, symmetric about <$E-~rho over 2>

and bounded by the interval

<$E[ min (r)~+~min (f tau ),~max (r)~+~max (f tau )]~=~[-~rho ~-~phi tau
,~phi tau ]> .
Since f changes only slowly over time for any single clock,

<$Eepsilon~==~[ min (r)~+~f tau ,~max (r)~+~f tau ]~=~ [-~ rho ,~0]~+~f
tau> ,

where <$Eepsilon> without argument designates the interval and
<$Eepsilon (x)> designates the pdf. In the following development
subscripts will be used on various quantities to indicate to which
entity or timestamp the quantity applies. Occasionally, <$Eepsilon> will
be used to designate an absolute maximum error, rather than the
interval, but the distinction will be clear from context.

Measurement Errors

In NTP the roundtrip delay and clock offset between two peers A and B
are determined by a procedure in which timestamps are exchanged via the
network paths between them. The procedure involves the four most recent
timestamps numbered as shown in Figure 14<$&fig14>, where the <$Etheta
sub 0> represents the true clock offset of peer B relative to peer A.
The <$ET sub 1> and <$ET sub 4> timestamps are determined relative to
the A clock, while the <$ET sub 2> and <$ET sub 3> timestamps are
determined relative to the B clock. The measured roundtrip delay
<$Edelta> and clock offset <$Etheta> of B relative to A are given by

<$Edelta~=~(T sub 4~-~T sub 1 )~-~(T sub 3~-~T sub 2
)~~~~and~~~~theta~=~{(T sub 2~-~T sub 1 )~+~(T sub 3~-~T sub 4 )} over
2> .

The errors inherent in determining the timestamps T1, T2, T3 and T4 are,
respectively,

<$Eepsilon sub 1~=~[-~rho sub A ,~0]>, <$E~epsilon sub 2~=~[-~rho sub B
,~0]>, <$E~epsilon sub 3~=~[-~rho sub B ,~0]~+~f sub B (T sub 3 ~-~T sub
2 )>, <$E~epsilon sub 4~=~[-~rho sub A ,~0]~+~f sub A (T sub 4 ~-~T sub
1 )> .

For specific peers A and B, where <$Ef sub A> and <$Ef sub B> can be
considered constants, the interval containing the maximum error inherent
in determining <$Edelta> is given by

<$E[ min ( epsilon sub 4 )~-~max ( epsilon sub 1 )~-~max ( epsilon sub 3
)~+~min ( epsilon sub 2 ),~ max ( epsilon sub 4 )~-~min ( epsilon sub 1
)~-~min ( epsilon sub 3 )~+~max ( epsilon sub 2 )]>
<$E=~[-~rho sub A~-~rho sub B ,~rho sub A ~+~rho sub B ]~+~f sub A (T
sub 4~-~T sub 1 )~-~f sub B (T sub 3~-~T sub 2 )> .

In the NTP local clock model the residual frequency errors <$Ef sub A>
and <$Ef sub B> are minimized through the use of a type-II phase-lock
loop (PLL). Under most conditions these errors will be small and can be
ignored. The pdf for the remaining errors is symmetric, so that <$Edelta
hat~=~<< delta >>> is an unbiased maximum-likelihood estimator for the
true roundtrip delay, independent of the particular values of <$Erho sub
A> and <$Erho sub B>.

However, in order to reliably bound the errors under all conditions of
component variation and operational regimes, the design of the PLL and
the tolerance of its intrinsic oscillator must be controlled so that it
is not possible under any circumstances for <$Ef sub A> or <$Ef sub B>
to exceed the bounds <$E[-~phi sub A ,~phi sub A ]> or <$E[-~phi sub B
,~phi sub B ]>, respectively. Setting <$Erho~=~max ( rho sub A ,~rho sub
B )> for convenience, the absolute maximum error <$Eepsilon sub delta>
inherent in determining roundtrip delay <$Edelta> is given by

<$Eepsilon sub delta~==~rho~+~phi sub A (T sub 4~-~T sub 1 )~+~phi sub B
(T sub 3~-~T sub 2 )> ,
neglecting residuals.

As in the case for <$Edelta>, where <$Ef sub A> and <$Ef sub B> can be
considered constants, the interval containing the maximum error inherent
in determining <$Etheta> is given by

<$E{[ min ( epsilon sub 2 )~-~max ( epsilon sub 1 )~+~min ( epsilon sub
3 )~-~max ( epsilon sub 4 ),~ max ( epsilon sub 2 )~-~min ( epsilon sub
1 )~+~max ( epsilon sub 3 )~-~min ( epsilon sub 4 )]} over 2>
<$E=~[ -~rho sub B ,~rho sub A ]~+~{f sub B (T sub 3~-~T sub 2 )~-~f sub
A (T sub 4 ~-~T sub 1 )} over 2> .

Under most conditions the errors due to <$Ef sub A> and <$Ef sub B> will
be small and can be ignored. If <$Erho sub A~=~rho sub B~=~rho>; that
is, if both the A and B clocks have the same resolution, the pdf for the
remaining errors is symmetric, so that <$Etheta hat~=~<< theta >>> is an
unbiased maximum-likelihood estimator for the true clock offset <$Etheta
sub 0>, independent of the particular value of <$Erho>. If <$Erho sub
A~!=~rho sub B>, <$E<< theta >>> is not an unbiased estimator; however,
the bias error is in the order of

<$E{rho sub A~-~rho sub B } over 2> .

and can usually be neglected.

Again setting <$Erho~=~max ( rho sub A ,~rho sub B )> for convenience,
the absolute maximum error <$Eepsilon sub theta> inherent in determining
clock offset <$Etheta> is given by

<$Eepsilon sub theta~==~{rho~+~phi sub A (T sub 4~-~T sub 1 )~+~phi sub
B (T sub 3~-~T sub 2 )} over 2 > .

Network Errors

In practice, errors due to stochastic network delays usually dominate.
In general, it is not possible to characterize network delays as a
stationary random process, since network queues can grow and shrink in
chaotic fashion and arriving customer traffic is frequently bursty.
However, it is a simple exercise to calculate bounds on clock offset
errors as a function of measured delay. Let <$ET sub 2~-~T sub 1~=~a>
and <$ET sub 3~-~T sub 4~=~b>. Then,

<$Edelta~=~a~-~b~~~~ and ~~~~theta~=~{a~+~b} over 2> .

The true offset of B relative to A is called <$Etheta sub 0> in Figure
14. Let x denote the actual delay between the departure of a message
from A and its arrival at B. Therefore, <$Ex~+~theta sub 0~=~T sub 2~-~T
sub 1~==~a>. Since x must be positive in our universe, <$Ex~=~a~-~theta
sub 0~>>=~0>, which requires <$Etheta sub 0~<<=~a>. A similar argument
requires that <$Eb~<<=~theta sub 0>, so surely <$Eb~<<=~theta sub
0~<<=~a>. This inequality can also be expressed

<$Eb~=~{a~+~b} over 2~-~{a~-~b} over 2~<<=~theta sub 0~<<=~{a~+~b} over
2~+~{a~-~b} over 2~=~a> ,

which is equivalent to

<$Etheta~-~delta over 2~<<=~theta sub 0~<<=~theta~+~delta over 2> .

In the previous section bounds on delay and offset errors were
determined. Thus, the inequality can be written

<$Etheta~-~epsilon sub theta~-~{delta~+~epsilon sub delta} over
2~<<=~theta sub 0~<<=~theta~+~epsilon sub theta~+~{delta~+~ epsilon sub
delta } over 2> ,
where <$Eepsilon sub theta> is the maximum offset error and <$Eepsilon
sub delta> is the maximum delay error derived previously. The quantity

<$Eepsilon~=~epsilon sub theta~+~epsilon sub delta over 2~=~rho~+~phi
sub A (T sub 4~-~T sub 1 )~+~phi sub B (T sub 3~-~T sub 2 )> ,

called the peer dispersion, defines the maximum error in the inequality.
Thus, the correctness interval I can be defined as the interval

<$EI~=~[ theta~-~delta over 2~-~epsilon ,~theta~+~delta over 2~+~epsilon
]> ,

in which the clock offset <$EC~=~theta> is the midpoint. By
construction, the true offset <$Etheta sub 0> must lie somewhere in this
interval.

Inherited Errors

As described in the NTP specification, the NTP time server maintains the
local clock <$ETHETA>, together with the root roundtrip delay <$EDELTA>
and root dispersion <$EEPSILON> relative to the primary reference source
at the root of the synchronization subnet. The values of these variables
are either included in each update message or can be derived as
described in the NTP specification. In addition, the protocol exchange
and clock-filter algorithm provide the clock offset <$Etheta> and
roundtrip delay <$Edelta> of the local clock relative to the peer clock,
as well as various error accumulations as described below. The following
discussion establishes how errors inherent in the time-transfer process
accumulate within the subnet and contribute to the overall error budget
at each server.

An NTP measurement update includes three parts: clock offset <$Etheta>,
roundtrip delay <$Edelta> and maximum error or dispersion <$Eepsilon> of
the local clock relative to a peer clock. In case of a primary clock
update, these values are usually all zero, although <$Eepsilon> can be
tailored to reflect the specified maximum error of the primary reference
source itself. In other cases <$Etheta> and <$Edelta> are calculated
directly from the four most recent timestamps, as described in the NTP
specification. The dispersion <$Eepsilon> includes the following
contributions:

1.

Each time the local clock is read a reading error is incurred due to the
finite granularity or precision of the implementation. This is called
the measurement dispersion <$Erho>.

2.

Once an offset is determined, an error due to frequency offset or skew
accumulates with time. This is called the skew dispersion <$Ephi tau>,
where <$Ephi> represents the skew-rate constant (<$Eroman NTP.MAXSKEW
over NTP.MAXAGE> in the NTP specification) and <$Etau> is the interval
since the dispersion was last updated.

3

When a series of offsets are determined at regular intervals and
accumulated in a window of samples, as in the NTP clock-filter
algorithm, the (estimated) additional error due to offset sample
variance is called the filter dispersion <$Eepsilon sub sigma>.

4.

When a number of peers are considered for synchronization and two or
more are determined to be correctly synchronized to a primary reference

source, as in the NTP clock-selection algorithm, the (estimated)
additional error due to offset sample variance is called the selection
dispersion <$Eepsilon sub xi>.

Figure 15<$&fig15> shows how these errors accumulate in the ordinary
course of NTP processing. Received messages from a single peer are
represented by the packet variables. From the four most recent
timestamps T1, T2, T3 and T4 the clock offset and roundtrip delay sample
for the local clock relative to the peer clock are calculated directly.
Included in the message are the root roundtrip delay <$EDELTA prime> and
root dispersion <$EEPSILON prime> of the peer itself; however, before
sending, the peer adds the measurement dispersion <$Erho> and skew
dispersion <$Ephi tau>, where these quantities are determined by the
peer and <$Etau> is the interval according to the peer clock since its
clock was last updated.

The NTP clock-filter procedure saves the most recent samples <$Etheta
sub i> and <$Edelta sub i> in the clock filter as described in the NTP
specification. The quantities <$Erho> and <$Ephi> characterize the local
clock maximum reading error and frequency error, respectively. Each
sample includes the dispersion <$Eepsilon sub i~=~rho~+~phi (T sub 4~-~T
sub 1 )>, which is set upon arrival. Each time a new sample arrives all
samples in the filter are updated with the skew dispersion <$Ephi tau
sub i>, where <$Etau sub i> is the interval since the last sample
arrived, as recorded in the variable peer.update. The clock-filter
algorithm determines the selected clock offset <$Etheta> (peer.offset),
together with the associated roundtrip delay <$Edelta> (peer.delay) and
filter dispersion <$Eepsilon sub sigma>, which is added to the
associated sample dispersion <$Eepsilon sub i> to form the peer
dispersion <$Eepsilon> (peer.dispersion).

The NTP clock-selection procedure selects a single peer to become the
synchronization source as described in the NTP specification. The
operation of the algorithm determines the final clock offset <$ETHETA>
(local clock), roundtrip delay <$EDELTA> (sys.rootdelay) and dispersion
<$EEPSILON> (sys.rootdispersion) relative to the root of the
synchronization subnet, as shown in Figure 15. Note the inclusion of the
selected peer dispersion and skew accumulation since the dispersion was
last updated, as well as the select dispersion <$Eepsilon sub xi>
computed by the clock-select algorithm itself. Also, note that, in order
to preserve overall synchronization subnet stability, the final clock
offset <$ETHETA> is in fact determined from the offset of the local
clock relative to the peer clock, rather than the root of the subnet.
Finally, note that the packet variables <$EDELTA prime> and <$EEPSILON
prime> are in fact determined from the latest message received, not at
the precise time the offset selected by the clock-filter algorithm was
determined. Minor errors arising due to these simplifications will be
ignored. Thus, the total dispersion accumulation relative to the root of
the synchronization subnet is

<$EEPSILON~=~epsilon~+~phi tau~+~epsilon sub xi~+~| THETA |~+~EPSILON
prime > ,

where <$Etau> is the time since the peer variables were last updated and
<$E| THETA |> is the initial absolute error in setting the local clock.

The three values of clock offset, roundtrip delay and dispersion are all
additive; that is, if <$ETHETA sub i>, <$EDELTA sub i> and <$EEPSILON
sub i> represent the values at peer i relative to the root of the
synchronization subnet, the values

<$ETHETA sub j (t)~==~THETA sub i~+~theta sub j (t)> ,   <$EDELTA sub j
(t)~==~DELTA sub i~+~delta sub j> ,   <$EEPSILON sub j (t)~==~EPSILON
sub i~+~epsilon sub i~+~epsilon sub j (t)> ,
represent the clock offset, roundtrip delay and dispersion of peer j at
time t. The time dependence of <$Etheta sub j (t)> and <$Eepsilon sub j
(t)> represents the local-clock correction and dispersion accumulated
since the last update was received from peer i, while the term
<$Eepsilon sub i> represents the dispersion accumulated by peer i from
the time its clock was last set until the latest update was sent to peer
j. Note that, while the offset of the local clock relative to the peer
clock can be determined directly, the offset relative to the root of the
synchronization subnet is not directly determinable, except on a
probabilistic basis and within the bounds established in this and the
previous section.

The NTP synchronization subnet topology is that of a tree rooted at the
primary server(s). Thus, there is an unbroken path from every time
server to the primary reference source. Accuracy and stability are
proportional to synchronization distance <$ELAMBDA>, defined as

<$ELAMBDA~==~EPSILON~+~DELTA over 2> .

The selection algorithm favors the minimum-distance paths and thus
maximizes accuracy and stability. Since <$ETHETA sub 0>, <$EDELTA sub 0>
and <$EEPSILON sub 0> are all zero, the sum of the clock offsets,
roundtrip delays and dispersions of each server along the minimum-
distance path from the root of the synchronization subnet to a given
server i are the clock offset <$ETHETA sub i>, roundtrip delay <$EDELTA
sub i> and dispersion <$EEPSILON sub i> inherited by and characteristic
of that server.

Correctness Principles

In order to minimize the occurrence of errors due to incorrect clocks
and maximize the reliability of the service, NTP relies on multiple
peers and disjoint peer paths whenever possible. In the previous
development it was shown that, if the primary reference source at the
root of the synchronization subnet is in fact a correct clock, then the
true offset <$Etheta sub 0> relative to that clock must be contained in

the interval

<$E[ THETA~-~LAMBDA ,~THETA~+~LAMBDA ]~==~[ THETA~-~EPSILON~-~DELTA over
2 ,~THETA~+~EPSILON~+~DELTA over 2 ]> .

When a number of clocks are involved, it is not clear beforehand which
are correct and which are not; however, as cited previously, there are a
number of techniques based on clustering and filtering principles which
yield a high probability of detecting and discarding incorrect clocks.
Marzullo and Owicki [MAR85] devised an algorithm designed to find an
appropriate interval containing the correct time given the confidence
intervals of m clocks, of which no more than f are considered incorrect.
The algorithm finds the smallest single intersection containing all
points in at least <$Em~-~f> of the given confidence intervals.

Figure 16<$&fig16> illustrates the operation of this algorithm with a
scenario involving four clocks A, B, C and D, with the calculated time
(shown by the <F128>ì<F255> symbol) and confidence interval shown for
each. These intervals are computed as described in previous sections of
this appendix. For instance, any point in the A interval may possibly
represent the actual time associated with that clock. If all clocks are
correct, there must exist a nonempty intersection including all four
intervals; but, clearly this is not the case in this scenario. However,
if it is assumed that one of the clocks is incorrect (e.g., D), it might
be possible to find a nonempty intersection including all but one of the
intervals. If not, it might be possible to find a nonempty intersection
including all but two of the intervals and so on.

The algorithm proposed by DEC for use in the Digital Time Service
[DEC89] is based on these principles. For the scenario illustrated in
Figure 16, it computes the interval for <$Em~=~4> clocks, three of which
turn out to be correct and one not. The low endpoint of the intersection
is found as follows. A variable f is initialized with the number of
presumed incorrect clocks, in this case zero, and a counter i is
initialized at zero. Starting from the lowest endpoint, the algorithm
increments i at each low endpoint, decrements i at each high endpoint,
and stops when <$Ei~>>=~m~-~f>. The counter records the number of
intersections and thus the number of presumed correct clocks. In the
example the counter never reaches four, so f is increased by one and the
procedure is repeated. This time the counter reaches three and stops at
the low endpoint of the intersection marked DTS. The upper endpoint of
this intersection is found using a similar procedure.

This algorithm will always find the smallest single intersection
containing points in at least one of the original <$Em~-~f> confidence
intervals as long as the number of incorrect clocks is less than half
the total <$Ef~<<~m over 2>. However, some points in the intersection
may not be contained in all <$Em~-~f> of the original intervals;
moreover, some or all of the calculated times (such as for C in Figure
16) may lie outside the intersection. In the NTP clock-selection

procedure the above algorithm is modified so as to include at least
<$Em~-~f> of the calculated times. In the modified algorithm a counter c
is initialized at zero. When starting from either endpoint, c is
incremented at each calculated time; however, neither f nor c are reset
between finding the low and high endpoints of the intersection. If after
both endpoints have been found <$Ec~>>~f>, f is increased by one and the
entire procedure is repeated. The revised algorithm finds the smallest
intersection of <$Em~-~f> intervals containing at least <$Em~-~f>
calculated times. As shown in Figure 16, the modified algorithm produces
the intersection marked NTP and including the calculated time for C.

In the NTP clock-selection procedure the peers represented by the clocks
in the final intersection, called the survivors, are placed on a
candidate list. In the remaining steps of the procedure one or more
survivors may be discarded from the list as outlyers. Finally, the
clock-combining algorithm described in Appendix F provides a weighted
average of the remaining survivors based on synchronization distance.
The resulting estimates represent a synthetic peer with offset between
the maximum and minimum offsets of the remaining survivors. This defines
the clock offset <$ETHETA>, total roundtrip total delay <$EDELTA> and
total dispersion <$EEPSILON> which the local clock inherits. In
principle, these values could be included in the time interface provided
by the operating system to the user, so that the user could evaluate the
quality of indications directly.

Appendix I. Selected C-Language Program Listings

Following are C-language program listings of selected algorithms
described in the NTP specification. While these have been tested as part
of a software simulator using data collected in regular operation, they
do not necessarily represent a standard implementation, since many other
implementations could in principle conform to the NTP specification.

Common Definitions and Variables

The following definitions are common to all procedures and peers.

#define NMAX 40                                 /* max clocks */
#define FMAX 8                                  /* max filter size */
#define HZ 1000                                 /* clock rate */
#define MAXSTRAT 15                             /* max stratum */
#define MAXSKEW 1                               /* max skew error per
MAXAGE */
#define MAXAGE 86400                            /* max clock age */
#define MAXDISP 16                              /* max dispersion */
#define MINCLOCK 3                              /* min survivor clocks
*/
#define MAXCLOCK 10                             /* min candidate clocks
*/
#define FILTER .5                                       /* filter weight

*/
#define SELECT .75                              /* select weight */

The folowing are peer state variables (one set for each peer).

double filtp[NMAX][FMAX];                               /* offset
samples */
double fildp[NMAX][FMAX];                       /* delay samples */
double filep[NMAX][FMAX];                       /* dispersion samples */
double tp[NMAX];                                        /* offset */
double dp[NMAX];                                        /* delay */
double ep[NMAX];                                        /* dispersion */
double rp[NMAX];                                        /* last offset
*/
double utc[NMAX];                                       /* update tstamp
*/
int st[NMAX];                                           /* stratum */

The following are system state variables and constants.

double rho = 1./HZ;                                     /* max reading
error */
double phi = MAXSKEW/MAXAGE;            /* max skew rate */
double bot, top;                                        /* confidence
interval limits */
double theta;                                           /* clock offset
*/
double delta;                                           /* roundtrip
delay */
double epsil;                                           /* dispersion */
double tstamp;                                  /* current time */
int source;                                             /* clock source
*/
int n1, n2;                                             /* min/max clock
ids */

The folowing are temporary lists shared by all peers and procedures.

double list[3*NMAX];                            /* temporary list*/
int index[3*NMAX];                                      /* index list */

Clock<196>Filter Algorithm

/*
   clock filter algorithm

   n = peer id, offset = sample offset, delay = sample delay, disp =
sample dispersion;
   computes tp[n] = peer offset, dp[n] = peer delay, ep[n] = peer
dispersion

 */

void filter(int n, double offset, double delay, double disp) {

        int i, j, k, m;                                 /* int temps */
        double x;                                       /* double temps
*/

        for (i = FMAX<196>1; i >> 0; i<196> <196>) {            /*
update/shift filter */
                filtp[n][i] = filtp[n][i<196>1]; fildp[n][i] =
fildp[n][i<196>1];
                filep[n][i] = filep[n][i<196>1]+phi*(tstamp<196>utc[n]);
                }
        utc[n] = tstamp; filtp[n][0] = offset<196>tp[0]; fildp[n][0] =
delay; filep[n][0] = disp;
        m = 0;                                          /*
construct/sort temp list */
        for (i = 0; i << FMAX; i++) {
                if (filep[n][i] >>= MAXDISP) continue;
                list[m] = filep[n][i]+fildp[n][i]/2.; index[m] = i;
                for (j = 0; j << m; j++) {
                        if (list[j] >> list[m]) {
                                x = list[j]; k = index[j]; list[j] =
list[m]; index[j] = index[m];
                                list[m] = x; index[m] = k;
                                }
                        }
                m = m+1;
                }

        if (m <<= 0) ep[n] = MAXDISP;           /* compute filter
dispersion */           
        else {
                ep[n] = 0;
                for (i = FMAX<196>1; i >>= 0; i<196> <196>) {
                        if (i << m) x =
fabs(filtp[n][index[0]]<196>filtp[n][index[i]]);
                        else x = MAXDISP;
                        ep[n] = FILTER*(ep[n]+x);
                        }
                i = index[0]; ep[n] = ep[n]+filep[n][i]; tp[n] =
filtp[n][i]; dp[n] = fildp[n][i];
                }
        return;
        }

Interval Intersection Algorithm

/*

   compute interval intersection

   computes bot = lowpoint, top = highpoint (bot >> top if no
intersection)
*/

void dts() {

        int f;                                          /* intersection
ceiling */
        int end;                                        /* endpoint
counter */
        int clk;                                                /*
falseticker counter */
        int i, j, k, m, n;                              /* int temps */
        double x, y;                                    /* double temps
*/

        m = 0; i = 0;
        for (n = n1; n <<= n2; n++) {   /* construct endpoint list */
                if (ep[n] >>= MAXDISP) continue;
                m = m+1;
                list[i] = tp[n]<196>dist(n); index[i] = <196>1; /*
lowpoint */
                for (j = 0; j << i; j++) {
                        if ((list[j] >> list[i]) || ((list[j] ==
list[i]) && (index[j] >> index[i]))) {
                                x = list[j]; k = index[j]; list[j] =
list[i]; index[j] = index[i];
                                list[i] = x; index[i] = k;
                                }
                        }
                i = i+1;

                list[i] = tp[n]; index[i] = 0;          /* midpoint */
                for (j = 0; j << i; j++) {
                        if ((list[j] >> list[i]) || ((list[j] ==
list[i]) && (index[j] >> index[i]))) {
                                x = list[j]; k = index[j]; list[j] =
list[i]; index[j] = index[i];
                                list[i] = x; index[i] = k;
                                }
                        }
                i = i+1;

                list[i] = tp[n]+dist(n); index[i] = 1;  /* highpoint */
                for (j = 0; j << i; j++) {
                        if ((list[j] >> list[i]) || ((list[j] ==
list[i]) && (index[j] >> index[i]))) {
                                x = list[j]; k = index[j]; list[j] =

list[i]; index[j] = index[i];
                                list[i] = x; index[i] = k;
                                }
                        }
                i = i+1;
                }

        if (m <<= 0) return;
        for (f = 0; f << m/2; f++) {                    /* find
intersection */
                clk = 0; end = 0;                       /* lowpoint */
                for (j = 0; j << i; j++) {
                        end = end<196>index[j]; bot = list[j];
                        if (end >>= (m<196>f)) break;
                        if (index[j] == 0) clk = clk+1;
                        }
                end = 0;                                /* highpoint */
                for (j = i<196>1; j >>= 0; j<196> <196>) {
                        end = end+index[j]; top = list[j];
                        if (end >>= (m<196>f)) break;
                        if (index[j] == 0) clk = clk+1;
                        }
                if (clk <<= f) break;
                }
        return;
        }

Clock<196>Selection Algorithm

/*
   select best subset of clocks in candidate list

   bot = lowpoint, top = highpoint; constructs index = candidate index
list,
   m = number of candidates, source = clock source,
   theta = clock offset, delta = roundtrip delay, epsil = dispersion
*/

void select() {

        double xi;                                      /* max select
dispersion */
        double eps;                                     /* min peer
dispersion */
        int i, j, k, n;                                 /* int temps */
        double x, y, z;                         /* double temps */

        m = 0;
        for (n = n1; n <<= n2; n++) {   /* make/sort candidate list */
                if ((st[n] >> 0) && (st[n] << MAXSTRAT) && (tp[n] >>=

bot) && (tp[n] <<= top)) {
                        list[m] = MAXDISP*st[n]+dist(n); index[m] = n;
                        for (j = 0; j << m; j++) {
                                if (list[j] >> list[m]) {
                                        x = list[j]; k = index[j];
list[j] = list[m]; index[j] = index[m];
                                        list[m] = x; index[m] = k;
                                        }
                                }
                        m = m+1;
                        }
                }
        if (m <<= 0) {
                source = 0; return;
                }
        if (m >> MAXCLOCK) m = MAXCLOCK;

        while (1) {                                     /* cast out
falsetickers */
                xi = 0.; eps = MAXDISP;
                for (j = 0; j << m; j++) {
                        x = 0.;
                        for (k = m<196>1; k >>= 0; k<196> <196>)
                                x =
SELECT*(x+fabs(tp[index[j]]<196>tp[index[k]]));
                        if (x >> xi) {
                                xi = x; i = j;          /* max(xi) */
                                }
                        x = ep[index[j]]+phi*(tstamp<196>utc[index[j]]);
                        if (x << eps) eps = x;          /* min(eps) */
                        }
                if ((xi <<= eps) || (m <<= MINCLOCK)) break;
                if (index[i] == source) source = 0;
                for (j = i; j << m<196>1; j++) index[j] = index[j+1];
                m = m<196>1;
                }

        i = index[0];                                   /* declare
winner */
        if (source != i)
                if (source == 0) source = i;
                else if (st[i] << st[source]) source = i;
        theta = combine(); delta = dp[i]; epsil =
ep[i]+phi*(tstamp<196>utc[i])+xi;
        return;
        }

Clock<196>Combining Procedure

/*

   compute weighted ensemble average

   index = candidate index list, m = number of candidates; returns
combined clock offset
*/

double combine() {

        int i;                                          /* int temps */
        double x, y, z;                         /* double temps */
        z = 0. ; y = 0.;
        for (i = 0; i << m; i++) {                      /* compute
weighted offset */
                j = index[i]; x = dist(j)); z = z+tp[j]/x; y = y+1./x;
                }
        return z/y;                                     /* normalize */
        }

Subroutine to Compute Synchronization Distance

/*
   compute synchronization distance

   n = peer id; returns synchronization distance
 */

double dist(int n) {

        return ep[n]+phi*(tstamp<196>utc[n])+fabs(dp[n])/2.;
        }

Security considerations
see Section 3.6 and Appendix C

Author's address
David L. Mills
Electrical Engineering Department
University of Delaware
Newark, DE 19716
Phone (302) 451<196>8247
EMail mills@udel.edu