154359Sroberto/* 254359Sroberto * systime -- routines to fiddle a UNIX clock. 3132451Sroberto * 4132451Sroberto * ATTENTION: Get approval from Dave Mills on all changes to this file! 5132451Sroberto * 654359Sroberto */ 782498Sroberto#include "ntp_machine.h" 882498Sroberto#include "ntp_fp.h" 982498Sroberto#include "ntp_syslog.h" 1082498Sroberto#include "ntp_unixtime.h" 1182498Sroberto#include "ntp_stdlib.h" 12182007Sroberto#include "ntp_random.h" 13182007Sroberto#include "ntpd.h" /* for sys_precision */ 1454359Sroberto 15132451Sroberto#ifdef SIM 16182007Sroberto# include "ntpsim.h" 17132451Sroberto#endif /*SIM */ 18132451Sroberto 1954359Sroberto#ifdef HAVE_SYS_PARAM_H 2054359Sroberto# include <sys/param.h> 2154359Sroberto#endif 2254359Sroberto#ifdef HAVE_UTMP_H 2354359Sroberto# include <utmp.h> 2454359Sroberto#endif /* HAVE_UTMP_H */ 2554359Sroberto#ifdef HAVE_UTMPX_H 2654359Sroberto# include <utmpx.h> 2754359Sroberto#endif /* HAVE_UTMPX_H */ 2854359Sroberto 2954359Sroberto/* 30132451Sroberto * These routines (get_systime, step_systime, adj_systime) implement an 31132451Sroberto * interface between the system independent NTP clock and the Unix 32132451Sroberto * system clock in various architectures and operating systems. 33132451Sroberto * 34132451Sroberto * Time is a precious quantity in these routines and every effort is 35132451Sroberto * made to minimize errors by always rounding toward zero and amortizing 36132451Sroberto * adjustment residues. By default the adjustment quantum is 1 us for 37132451Sroberto * the usual Unix tickadj() system call, but this can be increased if 38182007Sroberto * necessary by the tick configuration command. For instance, when the 39132451Sroberto * adjtime() quantum is a clock tick for a 100-Hz clock, the quantum 40132451Sroberto * should be 10 ms. 4154359Sroberto */ 42182007Sroberto#if defined RELIANTUNIX_CLOCK || defined SCO5_CLOCK 43182007Srobertodouble sys_tick = 10e-3; /* 10 ms tickadj() */ 44182007Sroberto#else 45182007Srobertodouble sys_tick = 1e-6; /* 1 us tickadj() */ 46182007Sroberto#endif 47132451Srobertodouble sys_residual = 0; /* adjustment residue (s) */ 4854359Sroberto 49132451Sroberto#ifndef SIM 5054359Sroberto 5154359Sroberto/* 52132451Sroberto * get_systime - return system time in NTP timestamp format. 5354359Sroberto */ 5454359Srobertovoid 5554359Srobertoget_systime( 56132451Sroberto l_fp *now /* system time */ 5754359Sroberto ) 5854359Sroberto{ 5954359Sroberto double dtemp; 6054359Sroberto 61132451Sroberto#if defined(HAVE_CLOCK_GETTIME) || defined(HAVE_GETCLOCK) 62132451Sroberto struct timespec ts; /* seconds and nanoseconds */ 63132451Sroberto 6454359Sroberto /* 65132451Sroberto * Convert Unix clock from seconds and nanoseconds to seconds. 66182007Sroberto * The bottom is only two bits down, so no need for fuzz. 67182007Sroberto * Some systems don't have that level of precision, however... 6854359Sroberto */ 6954359Sroberto# ifdef HAVE_CLOCK_GETTIME 70132451Sroberto clock_gettime(CLOCK_REALTIME, &ts); 7154359Sroberto# else 72132451Sroberto getclock(TIMEOFDAY, &ts); 7354359Sroberto# endif 7454359Sroberto now->l_i = ts.tv_sec + JAN_1970; 75132451Sroberto dtemp = ts.tv_nsec / 1e9; 76132451Sroberto 77132451Sroberto#else /* HAVE_CLOCK_GETTIME || HAVE_GETCLOCK */ 78132451Sroberto struct timeval tv; /* seconds and microseconds */ 79132451Sroberto 80132451Sroberto /* 81132451Sroberto * Convert Unix clock from seconds and microseconds to seconds. 82182007Sroberto * Add in unbiased random fuzz beneath the microsecond. 83132451Sroberto */ 84132451Sroberto GETTIMEOFDAY(&tv, NULL); 8554359Sroberto now->l_i = tv.tv_sec + JAN_1970; 86132451Sroberto dtemp = tv.tv_usec / 1e6; 8754359Sroberto 88132451Sroberto#endif /* HAVE_CLOCK_GETTIME || HAVE_GETCLOCK */ 8954359Sroberto 90132451Sroberto /* 91182007Sroberto * ntp_random() produces 31 bits (always nonnegative). 92182007Sroberto * This bit is done only after the precision has been 93182007Sroberto * determined. 94182007Sroberto */ 95182007Sroberto if (sys_precision != 0) 96182007Sroberto dtemp += (ntp_random() / FRAC - .5) / (1 << 97182007Sroberto -sys_precision); 98182007Sroberto 99182007Sroberto /* 100132451Sroberto * Renormalize to seconds past 1900 and fraction. 101132451Sroberto */ 102132451Sroberto dtemp += sys_residual; 103132451Sroberto if (dtemp >= 1) { 104132451Sroberto dtemp -= 1; 10554359Sroberto now->l_i++; 106182007Sroberto } else if (dtemp < 0) { 107132451Sroberto dtemp += 1; 108132451Sroberto now->l_i--; 109132451Sroberto } 110132451Sroberto dtemp *= FRAC; 11154359Sroberto now->l_uf = (u_int32)dtemp; 11254359Sroberto} 11354359Sroberto 11454359Sroberto 11554359Sroberto/* 116132451Sroberto * adj_systime - adjust system time by the argument. 11754359Sroberto */ 11856746Sroberto#if !defined SYS_WINNT 119132451Srobertoint /* 0 okay, 1 error */ 12054359Srobertoadj_systime( 121132451Sroberto double now /* adjustment (s) */ 12254359Sroberto ) 12354359Sroberto{ 124132451Sroberto struct timeval adjtv; /* new adjustment */ 125132451Sroberto struct timeval oadjtv; /* residual adjustment */ 126132451Sroberto double dtemp; 127132451Sroberto long ticks; 128132451Sroberto int isneg = 0; 12954359Sroberto 13054359Sroberto /* 131132451Sroberto * Most Unix adjtime() implementations adjust the system clock 132132451Sroberto * in microsecond quanta, but some adjust in 10-ms quanta. We 133132451Sroberto * carefully round the adjustment to the nearest quantum, then 134132451Sroberto * adjust in quanta and keep the residue for later. 13554359Sroberto */ 136132451Sroberto dtemp = now + sys_residual; 13754359Sroberto if (dtemp < 0) { 13854359Sroberto isneg = 1; 13954359Sroberto dtemp = -dtemp; 14054359Sroberto } 141132451Sroberto adjtv.tv_sec = (long)dtemp; 142132451Sroberto dtemp -= adjtv.tv_sec; 143132451Sroberto ticks = (long)(dtemp / sys_tick + .5); 144132451Sroberto adjtv.tv_usec = (long)(ticks * sys_tick * 1e6); 145132451Sroberto dtemp -= adjtv.tv_usec / 1e6; 146132451Sroberto sys_residual = dtemp; 14754359Sroberto 14854359Sroberto /* 149132451Sroberto * Convert to signed seconds and microseconds for the Unix 150132451Sroberto * adjtime() system call. Note we purposely lose the adjtime() 151132451Sroberto * leftover. 15254359Sroberto */ 153132451Sroberto if (isneg) { 154132451Sroberto adjtv.tv_sec = -adjtv.tv_sec; 155132451Sroberto adjtv.tv_usec = -adjtv.tv_usec; 156182007Sroberto sys_residual = -sys_residual; 15754359Sroberto } 158182007Sroberto if (adjtv.tv_sec != 0 || adjtv.tv_usec != 0) { 159182007Sroberto if (adjtime(&adjtv, &oadjtv) < 0) { 160182007Sroberto msyslog(LOG_ERR, "adj_systime: %m"); 161182007Sroberto return (0); 162182007Sroberto } 163132451Sroberto } 164132451Sroberto return (1); 16554359Sroberto} 16656746Sroberto#endif 16754359Sroberto 16854359Sroberto 16954359Sroberto/* 17054359Sroberto * step_systime - step the system clock. 17154359Sroberto */ 17254359Srobertoint 17354359Srobertostep_systime( 17454359Sroberto double now 17554359Sroberto ) 17654359Sroberto{ 17754359Sroberto struct timeval timetv, adjtv, oldtimetv; 17854359Sroberto int isneg = 0; 17954359Sroberto double dtemp; 18054359Sroberto#if defined(HAVE_CLOCK_GETTIME) || defined(HAVE_GETCLOCK) 18154359Sroberto struct timespec ts; 18254359Sroberto#endif 18354359Sroberto 18454359Sroberto dtemp = sys_residual + now; 18554359Sroberto if (dtemp < 0) { 18654359Sroberto isneg = 1; 18754359Sroberto dtemp = - dtemp; 18854359Sroberto adjtv.tv_sec = (int32)dtemp; 189132451Sroberto adjtv.tv_usec = (u_int32)((dtemp - 190132451Sroberto (double)adjtv.tv_sec) * 1e6 + .5); 19154359Sroberto } else { 19254359Sroberto adjtv.tv_sec = (int32)dtemp; 193132451Sroberto adjtv.tv_usec = (u_int32)((dtemp - 194132451Sroberto (double)adjtv.tv_sec) * 1e6 + .5); 19554359Sroberto } 19654359Sroberto#if defined(HAVE_CLOCK_GETTIME) || defined(HAVE_GETCLOCK) 197182007Sroberto# ifdef HAVE_CLOCK_GETTIME 19854359Sroberto (void) clock_gettime(CLOCK_REALTIME, &ts); 199182007Sroberto# else 20054359Sroberto (void) getclock(TIMEOFDAY, &ts); 201182007Sroberto# endif 20254359Sroberto timetv.tv_sec = ts.tv_sec; 20354359Sroberto timetv.tv_usec = ts.tv_nsec / 1000; 20454359Sroberto#else /* not HAVE_GETCLOCK */ 20554359Sroberto (void) GETTIMEOFDAY(&timetv, (struct timezone *)0); 20654359Sroberto#endif /* not HAVE_GETCLOCK */ 20754359Sroberto 20854359Sroberto oldtimetv = timetv; 20954359Sroberto 21054359Sroberto#ifdef DEBUG 21154359Sroberto if (debug) 21254359Sroberto printf("step_systime: step %.6f residual %.6f\n", now, sys_residual); 21354359Sroberto#endif 21454359Sroberto if (isneg) { 21554359Sroberto timetv.tv_sec -= adjtv.tv_sec; 21654359Sroberto timetv.tv_usec -= adjtv.tv_usec; 21754359Sroberto if (timetv.tv_usec < 0) { 21854359Sroberto timetv.tv_sec--; 21954359Sroberto timetv.tv_usec += 1000000; 22054359Sroberto } 22154359Sroberto } else { 22254359Sroberto timetv.tv_sec += adjtv.tv_sec; 22354359Sroberto timetv.tv_usec += adjtv.tv_usec; 22454359Sroberto if (timetv.tv_usec >= 1000000) { 22554359Sroberto timetv.tv_sec++; 22654359Sroberto timetv.tv_usec -= 1000000; 22754359Sroberto } 22854359Sroberto } 229132451Sroberto if (ntp_set_tod(&timetv, NULL) != 0) { 230132451Sroberto msyslog(LOG_ERR, "step-systime: %m"); 23154359Sroberto return (0); 23254359Sroberto } 23354359Sroberto sys_residual = 0; 23454359Sroberto 23554359Sroberto#ifdef NEED_HPUX_ADJTIME 23654359Sroberto /* 23754359Sroberto * CHECKME: is this correct when called by ntpdate????? 23854359Sroberto */ 23954359Sroberto _clear_adjtime(); 24054359Sroberto#endif 24154359Sroberto 24254359Sroberto /* 24354359Sroberto * FreeBSD, for example, has: 24454359Sroberto * struct utmp { 24554359Sroberto * char ut_line[UT_LINESIZE]; 24654359Sroberto * char ut_name[UT_NAMESIZE]; 24754359Sroberto * char ut_host[UT_HOSTSIZE]; 24854359Sroberto * long ut_time; 24954359Sroberto * }; 25054359Sroberto * and appends line="|", name="date", host="", time for the OLD 25154359Sroberto * and appends line="{", name="date", host="", time for the NEW 25254359Sroberto * to _PATH_WTMP . 25354359Sroberto * 25454359Sroberto * Some OSes have utmp, some have utmpx. 25554359Sroberto */ 25654359Sroberto 25754359Sroberto /* 258132451Sroberto * Write old and new time entries in utmp and wtmp if step 259132451Sroberto * adjustment is greater than one second. 26054359Sroberto * 26154359Sroberto * This might become even Uglier... 26254359Sroberto */ 26354359Sroberto if (oldtimetv.tv_sec != timetv.tv_sec) 26454359Sroberto { 26554359Sroberto#ifdef HAVE_UTMP_H 26654359Sroberto struct utmp ut; 26754359Sroberto#endif 26854359Sroberto#ifdef HAVE_UTMPX_H 26954359Sroberto struct utmpx utx; 27054359Sroberto#endif 27154359Sroberto 27254359Sroberto#ifdef HAVE_UTMP_H 27354359Sroberto memset((char *)&ut, 0, sizeof(ut)); 27454359Sroberto#endif 27554359Sroberto#ifdef HAVE_UTMPX_H 27654359Sroberto memset((char *)&utx, 0, sizeof(utx)); 27754359Sroberto#endif 27854359Sroberto 27954359Sroberto /* UTMP */ 28054359Sroberto 28154359Sroberto#ifdef UPDATE_UTMP 28254359Sroberto# ifdef HAVE_PUTUTLINE 28354359Sroberto ut.ut_type = OLD_TIME; 28454359Sroberto (void)strcpy(ut.ut_line, OTIME_MSG); 28554359Sroberto ut.ut_time = oldtimetv.tv_sec; 28654359Sroberto pututline(&ut); 28754359Sroberto setutent(); 28854359Sroberto ut.ut_type = NEW_TIME; 28954359Sroberto (void)strcpy(ut.ut_line, NTIME_MSG); 29054359Sroberto ut.ut_time = timetv.tv_sec; 29154359Sroberto pututline(&ut); 29254359Sroberto endutent(); 29354359Sroberto# else /* not HAVE_PUTUTLINE */ 29454359Sroberto# endif /* not HAVE_PUTUTLINE */ 29554359Sroberto#endif /* UPDATE_UTMP */ 29654359Sroberto 29754359Sroberto /* UTMPX */ 29854359Sroberto 29954359Sroberto#ifdef UPDATE_UTMPX 30054359Sroberto# ifdef HAVE_PUTUTXLINE 30154359Sroberto utx.ut_type = OLD_TIME; 30254359Sroberto (void)strcpy(utx.ut_line, OTIME_MSG); 30354359Sroberto utx.ut_tv = oldtimetv; 30454359Sroberto pututxline(&utx); 30554359Sroberto setutxent(); 30654359Sroberto utx.ut_type = NEW_TIME; 30754359Sroberto (void)strcpy(utx.ut_line, NTIME_MSG); 30854359Sroberto utx.ut_tv = timetv; 30954359Sroberto pututxline(&utx); 31054359Sroberto endutxent(); 31154359Sroberto# else /* not HAVE_PUTUTXLINE */ 31254359Sroberto# endif /* not HAVE_PUTUTXLINE */ 31354359Sroberto#endif /* UPDATE_UTMPX */ 31454359Sroberto 31554359Sroberto /* WTMP */ 31654359Sroberto 31754359Sroberto#ifdef UPDATE_WTMP 31854359Sroberto# ifdef HAVE_PUTUTLINE 31954359Sroberto utmpname(WTMP_FILE); 32054359Sroberto ut.ut_type = OLD_TIME; 32154359Sroberto (void)strcpy(ut.ut_line, OTIME_MSG); 32254359Sroberto ut.ut_time = oldtimetv.tv_sec; 32354359Sroberto pututline(&ut); 32454359Sroberto ut.ut_type = NEW_TIME; 32554359Sroberto (void)strcpy(ut.ut_line, NTIME_MSG); 32654359Sroberto ut.ut_time = timetv.tv_sec; 32754359Sroberto pututline(&ut); 32854359Sroberto endutent(); 32954359Sroberto# else /* not HAVE_PUTUTLINE */ 33054359Sroberto# endif /* not HAVE_PUTUTLINE */ 33154359Sroberto#endif /* UPDATE_WTMP */ 33254359Sroberto 33354359Sroberto /* WTMPX */ 33454359Sroberto 33554359Sroberto#ifdef UPDATE_WTMPX 33654359Sroberto# ifdef HAVE_PUTUTXLINE 33754359Sroberto utx.ut_type = OLD_TIME; 33854359Sroberto utx.ut_tv = oldtimetv; 33954359Sroberto (void)strcpy(utx.ut_line, OTIME_MSG); 34054359Sroberto# ifdef HAVE_UPDWTMPX 34154359Sroberto updwtmpx(WTMPX_FILE, &utx); 34254359Sroberto# else /* not HAVE_UPDWTMPX */ 34354359Sroberto# endif /* not HAVE_UPDWTMPX */ 34454359Sroberto# else /* not HAVE_PUTUTXLINE */ 34554359Sroberto# endif /* not HAVE_PUTUTXLINE */ 34654359Sroberto# ifdef HAVE_PUTUTXLINE 34754359Sroberto utx.ut_type = NEW_TIME; 34854359Sroberto utx.ut_tv = timetv; 34954359Sroberto (void)strcpy(utx.ut_line, NTIME_MSG); 35054359Sroberto# ifdef HAVE_UPDWTMPX 35154359Sroberto updwtmpx(WTMPX_FILE, &utx); 35254359Sroberto# else /* not HAVE_UPDWTMPX */ 35354359Sroberto# endif /* not HAVE_UPDWTMPX */ 35454359Sroberto# else /* not HAVE_PUTUTXLINE */ 35554359Sroberto# endif /* not HAVE_PUTUTXLINE */ 35654359Sroberto#endif /* UPDATE_WTMPX */ 35754359Sroberto 35854359Sroberto } 35954359Sroberto return (1); 36054359Sroberto} 361132451Sroberto 362132451Sroberto#else /* SIM */ 363132451Sroberto/* 364132451Sroberto * Clock routines for the simulator - Harish Nair, with help 365132451Sroberto */ 366132451Sroberto/* 367132451Sroberto * get_systime - return the system time in NTP timestamp format 368132451Sroberto */ 369132451Srobertovoid 370132451Srobertoget_systime( 371132451Sroberto l_fp *now /* current system time in l_fp */ ) 372132451Sroberto{ 373132451Sroberto /* 374132451Sroberto * To fool the code that determines the local clock precision, 375132451Sroberto * we advance the clock a minimum of 200 nanoseconds on every 376132451Sroberto * clock read. This is appropriate for a typical modern machine 377132451Sroberto * with nanosecond clocks. Note we make no attempt here to 378132451Sroberto * simulate reading error, since the error is so small. This may 379132451Sroberto * change when the need comes to implement picosecond clocks. 380132451Sroberto */ 381132451Sroberto if (ntp_node.ntp_time == ntp_node.last_time) 382132451Sroberto ntp_node.ntp_time += 200e-9; 383132451Sroberto ntp_node.last_time = ntp_node.ntp_time; 384132451Sroberto DTOLFP(ntp_node.ntp_time, now); 385132451Sroberto} 386132451Sroberto 387132451Sroberto 388132451Sroberto/* 389132451Sroberto * adj_systime - advance or retard the system clock exactly like the 390132451Sroberto * real thng. 391132451Sroberto */ 392132451Srobertoint /* always succeeds */ 393132451Srobertoadj_systime( 394132451Sroberto double now /* time adjustment (s) */ 395132451Sroberto ) 396132451Sroberto{ 397132451Sroberto struct timeval adjtv; /* new adjustment */ 398132451Sroberto double dtemp; 399132451Sroberto long ticks; 400132451Sroberto int isneg = 0; 401132451Sroberto 402132451Sroberto /* 403132451Sroberto * Most Unix adjtime() implementations adjust the system clock 404132451Sroberto * in microsecond quanta, but some adjust in 10-ms quanta. We 405132451Sroberto * carefully round the adjustment to the nearest quantum, then 406132451Sroberto * adjust in quanta and keep the residue for later. 407132451Sroberto */ 408132451Sroberto dtemp = now + sys_residual; 409132451Sroberto if (dtemp < 0) { 410132451Sroberto isneg = 1; 411132451Sroberto dtemp = -dtemp; 412132451Sroberto } 413132451Sroberto adjtv.tv_sec = (long)dtemp; 414132451Sroberto dtemp -= adjtv.tv_sec; 415132451Sroberto ticks = (long)(dtemp / sys_tick + .5); 416132451Sroberto adjtv.tv_usec = (long)(ticks * sys_tick * 1e6); 417132451Sroberto dtemp -= adjtv.tv_usec / 1e6; 418132451Sroberto sys_residual = dtemp; 419132451Sroberto 420132451Sroberto /* 421132451Sroberto * Convert to signed seconds and microseconds for the Unix 422132451Sroberto * adjtime() system call. Note we purposely lose the adjtime() 423132451Sroberto * leftover. 424132451Sroberto */ 425132451Sroberto if (isneg) { 426132451Sroberto adjtv.tv_sec = -adjtv.tv_sec; 427132451Sroberto adjtv.tv_usec = -adjtv.tv_usec; 428132451Sroberto sys_residual = -sys_residual; 429132451Sroberto } 430182007Sroberto ntp_node.adj = now; 431182007Sroberto return (1); 432132451Sroberto} 433132451Sroberto 434132451Sroberto 435132451Sroberto/* 436132451Sroberto * step_systime - step the system clock. We are religious here. 437132451Sroberto */ 438132451Srobertoint /* always succeeds */ 439132451Srobertostep_systime( 440132451Sroberto double now /* step adjustment (s) */ 441132451Sroberto ) 442132451Sroberto{ 443182007Sroberto#ifdef DEBUG 444182007Sroberto if (debug) 445182007Sroberto printf("step_systime: time %.6f adj %.6f\n", 446182007Sroberto ntp_node.ntp_time, now); 447182007Sroberto#endif 448182007Sroberto ntp_node.ntp_time += now; 449132451Sroberto return (1); 450132451Sroberto} 451132451Sroberto 452132451Sroberto/* 453132451Sroberto * node_clock - update the clocks 454132451Sroberto */ 455132451Srobertoint /* always succeeds */ 456132451Srobertonode_clock( 457132451Sroberto Node *n, /* global node pointer */ 458132451Sroberto double t /* node time */ 459132451Sroberto ) 460132451Sroberto{ 461132451Sroberto double dtemp; 462132451Sroberto 463132451Sroberto /* 464132451Sroberto * Advance client clock (ntp_time). Advance server clock 465132451Sroberto * (clk_time) adjusted for systematic and random frequency 466132451Sroberto * errors. The random error is a random walk computed as the 467132451Sroberto * integral of samples from a Gaussian distribution. 468132451Sroberto */ 469132451Sroberto dtemp = t - n->ntp_time; 470132451Sroberto n->time = t; 471132451Sroberto n->ntp_time += dtemp; 472132451Sroberto n->ferr += gauss(0, dtemp * n->fnse); 473132451Sroberto n->clk_time += dtemp * (1 + n->ferr); 474132451Sroberto 475132451Sroberto /* 476132451Sroberto * Perform the adjtime() function. If the adjustment completed 477132451Sroberto * in the previous interval, amortize the entire amount; if not, 478132451Sroberto * carry the leftover to the next interval. 479132451Sroberto */ 480132451Sroberto dtemp *= n->slew; 481132451Sroberto if (dtemp < fabs(n->adj)) { 482132451Sroberto if (n->adj < 0) { 483132451Sroberto n->adj += dtemp; 484132451Sroberto n->ntp_time -= dtemp; 485132451Sroberto } else { 486132451Sroberto n->adj -= dtemp; 487132451Sroberto n->ntp_time += dtemp; 488132451Sroberto } 489132451Sroberto } else { 490132451Sroberto n->ntp_time += n->adj; 491132451Sroberto n->adj = 0; 492132451Sroberto } 493132451Sroberto return (0); 494132451Sroberto} 495132451Sroberto 496132451Sroberto 497132451Sroberto/* 498132451Sroberto * gauss() - returns samples from a gaussion distribution 499132451Sroberto */ 500132451Srobertodouble /* Gaussian sample */ 501132451Srobertogauss( 502132451Sroberto double m, /* sample mean */ 503132451Sroberto double s /* sample standard deviation (sigma) */ 504132451Sroberto ) 505132451Sroberto{ 506132451Sroberto double q1, q2; 507132451Sroberto 508132451Sroberto /* 509132451Sroberto * Roll a sample from a Gaussian distribution with mean m and 510132451Sroberto * standard deviation s. For m = 0, s = 1, mean(y) = 0, 511132451Sroberto * std(y) = 1. 512132451Sroberto */ 513132451Sroberto if (s == 0) 514132451Sroberto return (m); 515132451Sroberto while ((q1 = drand48()) == 0); 516132451Sroberto q2 = drand48(); 517132451Sroberto return (m + s * sqrt(-2. * log(q1)) * cos(2. * PI * q2)); 518132451Sroberto} 519132451Sroberto 520132451Sroberto 521132451Sroberto/* 522132451Sroberto * poisson() - returns samples from a network delay distribution 523132451Sroberto */ 524132451Srobertodouble /* delay sample (s) */ 525132451Srobertopoisson( 526132451Sroberto double m, /* fixed propagation delay (s) */ 527132451Sroberto double s /* exponential parameter (mu) */ 528132451Sroberto ) 529132451Sroberto{ 530132451Sroberto double q1; 531132451Sroberto 532132451Sroberto /* 533132451Sroberto * Roll a sample from a composite distribution with propagation 534132451Sroberto * delay m and exponential distribution time with parameter s. 535132451Sroberto * For m = 0, s = 1, mean(y) = std(y) = 1. 536132451Sroberto */ 537132451Sroberto if (s == 0) 538132451Sroberto return (m); 539132451Sroberto while ((q1 = drand48()) == 0); 540132451Sroberto return (m - s * log(q1 * s)); 541132451Sroberto} 542132451Sroberto#endif /* SIM */ 543