# include # include # include # include # include "toms886.h" /******************************************************************************/ void cheb ( int deg, double pt, double tcheb[] ) /******************************************************************************/ /* Purpose: CHEB computes normalized Chebyshev polynomials. Discussion: This subroutine computes the array TCHEB of normalized Chebyshev polynomials from degree 0 to DEG: T_0(x)=1, T_j(x) = sqrt(2) * cos ( j * acos(x) ) at the point x = PT. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 February 2014 Author: Original FORTRAN77 version by Marco Caliari, Stefano De Marchi, Marco Vianello. This C version by John Burkardt. Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Algorithm 886: Padua2D: Lagrange Interpolation at Padua Points on Bivariate Domains, ACM Transactions on Mathematical Software, Volume 35, Number 3, October 2008, Article 21, 11 pages. Parameters: Input, int DEG, the degree. 0 <= DEG. Input, double PT, the evaluation point. Output, double TCHEB[DEG+1], the value of the normalized Chebyshev polynomials of degrees 0 through DEG at the point PT. */ { int j; const double sqrt2 = 1.4142135623730951; if ( deg < 0 ) { return; } tcheb[0] = 1.0; if ( deg < 1 ) { return; } tcheb[1] = sqrt2 * pt; if ( deg < 2 ) { return; } tcheb[2] = 2.0 * pt * tcheb[1] - sqrt2 * tcheb[0]; /* Chebyshev recurrence. */ for ( j = 3; j <= deg; j++ ) { tcheb[j] = 2.0 * pt * tcheb[j-1] - tcheb[j-2]; } return; } /******************************************************************************/ void dgemm ( char transa, char transb, int m, int n, int k, double alpha, double a[], int lda, double b[], int ldb, double beta, double c[], int ldc ) /******************************************************************************/ /* Purpose: DGEMM computes C = alpha * A * B and related operations. Discussion: DGEMM performs one of the matrix-matrix operations C := alpha * op ( A ) * op ( B ) + beta * C, where op ( X ) is one of op ( X ) = X or op ( X ) = X', ALPHA and BETA are scalars, and A, B and C are matrices, with op ( A ) an M by K matrix, op ( B ) a K by N matrix and C an N by N matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 10 February 2014 Author: Original FORTRAN77 version by Jack Dongarra. C version by John Burkardt. Parameters: Input, char TRANSA, specifies the form of op( A ) to be used in the matrix multiplication as follows: 'N' or 'n', op ( A ) = A. 'T' or 't', op ( A ) = A'. 'C' or 'c', op ( A ) = A'. Input, char TRANSB, specifies the form of op ( B ) to be used in the matrix multiplication as follows: 'N' or 'n', op ( B ) = B. 'T' or 't', op ( B ) = B'. 'C' or 'c', op ( B ) = B'. Input, int M, the number of rows of the matrix op ( A ) and of the matrix C. 0 <= M. Input, int N, the number of columns of the matrix op ( B ) and the number of columns of the matrix C. 0 <= N. Input, int K, the number of columns of the matrix op ( A ) and the number of rows of the matrix op ( B ). 0 <= K. Input, double ALPHA, the scalar multiplier for op ( A ) * op ( B ). Input, double A(LDA,KA), where: if TRANSA is 'N' or 'n', KA is equal to K, and the leading M by K part of the array contains A; if TRANSA is not 'N' or 'n', then KA is equal to M, and the leading K by M part of the array must contain the matrix A. Input, int LDA, the first dimension of A as declared in the calling routine. When TRANSA = 'N' or 'n' then LDA must be at least max ( 1, M ), otherwise LDA must be at least max ( 1, K ). Input, double B(LDB,KB), where: if TRANSB is 'N' or 'n', kB is N, and the leading K by N part of the array contains B; if TRANSB is not 'N' or 'n', then KB is equal to K, and the leading N by K part of the array must contain the matrix B. Input, int LDB, the first dimension of B as declared in the calling routine. When TRANSB = 'N' or 'n' then LDB must be at least max ( 1, K ), otherwise LDB must be at least max ( 1, N ). Input, double BETA, the scalar multiplier for C. Input, double C(LDC,N). Before entry, the leading M by N part of this array must contain the matrix C, except when BETA is 0.0, in which case C need not be set on entry. On exit, the array C is overwritten by the M by N matrix alpha * op ( A ) * op ( B ) + beta * C. Input, int LDC, the first dimension of C as declared in the calling routine. max ( 1, M ) <= LDC. */ { int i; int j; int l; int nrowa; int nrowb; int nota; int notb; double temp; /* Set NOTA and NOTB as true if A and B respectively are not transposed and set NROWA, NCOLA and NROWB as the number of rows and columns of A and the number of rows of B respectively. */ nota = ( ( transa == 'N' ) || ( transa == 'n' ) ); if ( nota ) { nrowa = m; } else { nrowa = k; } notb = ( ( transb == 'N' ) || ( transb == 'n' ) ); if ( notb ) { nrowb = k; } else { nrowb = n; } /* Test the input parameters. */ if ( ! ( transa == 'N' || transa == 'n' || transa == 'C' || transa == 'c' || transa == 'T' || transa == 't' ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input TRANSA had illegal value.\n" ); exit ( 1 ); } if ( ! ( transb == 'N' || transb == 'n' || transb == 'C' || transb == 'c' || transb == 'T' || transb == 't' ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input TRANSB had illegal value.\n" ); exit ( 1 ); } if ( m < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input M had illegal value.\n" ); exit ( 1 ); } if ( n < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input N had illegal value.\n" ); exit ( 1 ); } if ( k < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input K had illegal value.\n" ); exit ( 1 ); } if ( lda < i4_max ( 1, nrowa ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input LDA had illegal value.\n" ); exit ( 1 ); } if ( ldb < i4_max ( 1, nrowb ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input LDB had illegal value.\n" ); exit ( 1 ); } if ( ldc < i4_max ( 1, m ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DGEMM - Fatal error!\n" ); fprintf ( stderr, " Input LDC had illegal value.\n" ); exit ( 1 ); } /* Quick return if possible. */ if ( m == 0 ) { return; } if ( n == 0 ) { return; } if ( ( alpha == 0.0 || k == 0 ) && ( beta == 1.0 ) ) { return; } /* And if alpha is 0.0. */ if ( alpha == 0.0 ) { if ( beta == 0.0 ) { for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = 0.0; } } } else { for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = beta * c[i+j*ldc]; } } } return; } /* Start the operations. */ if ( notb ) { /* Form C := alpha*A*B + beta*C. */ if ( nota ) { for ( j = 0; j < n; j++ ) { if ( beta == 0.0 ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = 0.0; } } else if ( beta != 1.0 ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = beta * c[i+j*ldc]; } } for ( l = 0; l < k; l++ ) { if ( b[l+j*ldb] != 0.0 ) { temp = alpha * b[l+j*ldb]; for ( i = 0; i < m; i++ ) { c[i+j*ldc] = c[i+j*ldc] + temp * a[i+l*lda]; } } } } } /* Form C := alpha*A'*B + beta*C */ else { for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { temp = 0.0; for ( l = 0; l < k; l++ ) { temp = temp + a[l+i*lda] * b[l+j*ldb]; } if ( beta == 0.0 ) { c[i+j*ldc] = alpha * temp; } else { c[i+j*ldc] = alpha * temp + beta * c[i+j*ldc]; } } } } } /* Form C := alpha*A*B' + beta*C */ else { if ( nota ) { for ( j = 0; j < n; j++ ) { if ( beta == 0.0 ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = 0.0; } } else if ( beta != 1.0 ) { for ( i = 0; i < m; i++ ) { c[i+j*ldc] = beta * c[i+j*ldc]; } } for ( l = 0; l < k; l++ ) { if ( b[j+l*ldb] != 0.0 ) { temp = alpha * b[j+l*ldb]; for ( i = 0; i < m; i++ ) { c[i+j*ldc] = c[i+j*ldc] + temp * a[i+l*lda]; } } } } } /* Form C := alpha*A'*B' + beta*C */ else { for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { temp = 0.0; for ( l = 0; l < k; l++ ) { temp = temp + a[l+i*lda] * b[j+l*ldb]; } if ( beta == 0.0 ) { c[i+j*ldc] = alpha * temp; } else { c[i+j*ldc] = alpha * temp + beta * c[i+j*ldc]; } } } } } return; } /******************************************************************************/ double franke ( double x, double y ) /******************************************************************************/ /* Purpose: FRANKE returns the value of the Franke function #1. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 February 2014 Author: John Burkardt Reference: Richard Franke, Scattered Data Interpolation: Tests of Some Methods, Mathematics of Computation, Volume 38, Number 157, January 1982, pages 181-200. Parameters: Input, double X, Y, the evalution points. Output, double FRANKE, the function values. */ { double f; f = 0.75 * exp ( - ( pow ( 9.0 * x - 2.0, 2 ) + pow ( 9.0 * y - 2.0, 2 ) ) / 4.0 ) + 0.75 * exp ( - ( pow ( 9.0 * x + 1.0, 2 ) ) / 49.0 - ( 9.0 * y + 1.0 ) / 10.0 ) + 0.5 * exp ( - ( pow ( 9.0 * x - 7.0, 2 ) + pow ( 9.0 * y - 3.0, 2 ) ) / 4.0 ) - 0.2 * exp ( - pow ( 9.0 * x - 4.0, 2 ) - pow ( 9.0 * y - 7.0, 2 ) ); return f; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ void padua2 ( int deg, int degmax, int npd, double wpd[], double fpd[], double raux1[], double raux2[], double c0[], double *esterr ) /******************************************************************************/ /* Purpose: PADUA2 computes the Padua interpolation coefficient matrix. Discussion: This function computes the coefficient matrix C0, in the orthonormal Chebyshev basis T_j(x)T_{k-j}(y), 0 <= j <= k <= DEG, T_0(x)=1, T_j(x) = sqrt(2) * cos(j * acos(x)), of the interpolation polynomial of degree DEG of the function values FPD at the set of NPD Padua points (PD1,PD2) in the square [-1,1]^2. The interpolant may be evaluated at an arbitrary point by the function PD2VAL. PD1, PD2 and WPD are the Padua points and weights computed by PDPTS. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 February 2014 Author: Original FORTRAN77 version by Marco Caliari, Stefano De Marchi, Marco Vianello. This C version by John Burkardt. Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Algorithm 886: Padua2D: Lagrange Interpolation at Padua Points on Bivariate Domains, ACM Transactions on Mathematical Software, Volume 35, Number 3, October 2008, Article 21, 11 pages. Parameters: Input, int DEG, the degree of approximation. Input, int DEGMAX, the maximum degree allowed. Input, int NPD, the number of Padua points. Input, double WPD[NPD], the weights. Input, double FPD[NPD], the value at the Padua points of the function to be interpolated. Workspace, double RAUX1[(DEGMAX+1)*(DEG+2)]. Workspace, double RAUX2[(DEGMAX+1)*(DEG+2)]. Output, double C0[(DEGMAX+2)*(DEG+1)], the coefficient matrix. Output, double *ESTERR, the estimated error. */ { double angle; int i; int j; int k; const double pi = 3.1415926535897931; double pt; /* Build the matrix P_2 and store it in RAUX2. */ for ( i = 0; i <= deg + 1; i++ ) { angle = ( double ) ( i ) * pi / ( double ) ( deg + 1 ); pt = - cos ( angle ); cheb ( deg, pt, raux2 + i * ( degmax + 1 ) ); } /* Build the matrix G(f) and store it in C0. */ for ( j = 0; j <= deg + 1; j++ ) { for ( i = 0; i <= degmax + 1; i++ ) { c0[i+j*(degmax+2)] = 0.0; } } k = 0; for ( j = 0; j <= deg + 1; j++ ) { for ( i = 0; i <= deg; i++ ) { if ( ( i + j ) % 2 == 0 ) { c0[i+j*(degmax+2)] = fpd[k] * wpd[k]; k = k + 1; } else { c0[i+j*(degmax+2)] = 0.0; } } } /* Compute the matrix-matrix product G(f)*P_2' and store it in RAUX1. */ dgemm ( 'n', 't', deg + 1, deg + 1, deg + 2, 1.0, c0, degmax + 2, raux2, degmax + 1, 0.0, raux1, degmax + 1 ); /* Build the matrix P_1 and store it in RAUX2. */ for ( i = 0; i <= deg; i++ ) { angle = ( double ) ( i ) * pi / ( double ) ( deg ); pt = - cos ( angle ); cheb ( deg, pt, raux2 + i * ( degmax + 1 ) ); } /* Compute the matrix-matrix product C(f) = P_1 * ( G(f) * P_2' ) and store it in C0. */ dgemm ( 'n', 'n', deg + 1, deg + 1, deg + 1, 1.0, raux2, degmax + 1, raux1, degmax + 1, 0.0, c0, degmax + 2 ); c0[deg+0*(degmax+2)] = c0[deg+0*(degmax+2)] / 2.0; /* Estimate the error. */ *esterr = 0.0; for ( j = 0; j <= 2; j++ ) { for ( i = 0; i <= deg - j; i++ ) { *esterr = *esterr + fabs ( c0[i+(deg-i-j)*(degmax+2)] ); } } *esterr = 2.0 * *esterr; return; } /******************************************************************************/ double pd2val ( int deg, int degmax, double c0[], double tg1, double tg2 ) /******************************************************************************/ /* Purpose: PD2VAL evaluates the Padua2 interpolant. Discussion: This function returns the value of the interpolant at (TG1,TG2). C0 is the matrix of the coefficients computed by PADUA2. Licensing: This code is distributed under the GNU LGPL license. Modified: 16 February 2014 Author: Original FORTRAN77 version by Marco Caliari, Stefano De Marchi, Marco Vianello. This C version by John Burkardt. Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Algorithm 886: Padua2D: Lagrange Interpolation at Padua Points on Bivariate Domains, ACM Transactions on Mathematical Software, Volume 35, Number 3, October 2008, Article 21, 11 pages. Parameters: Input, int DEG, the degree of approximation. Input, int DEGMAX, the maximum degree allowed. Input, double C0[(0:DEGMAX+1)*(0:DEG)], the coefficient matrix. Input, double TG1, TG2, the first and second coordinates of the target point. Output, double PD2VAL, the value of the interpolant at the target point. */ { int i; int j; double t; double *ttg1; double *ttg2; double value; /* Compute the normalized Chebyshev polynomials at the target point. */ ttg1 = ( double * ) malloc ( ( deg + 1 ) * sizeof ( double ) ); cheb ( deg, tg1, ttg1 ); ttg2 = ( double * ) malloc ( ( deg + 1 ) * sizeof ( double ) ); cheb ( deg, tg2, ttg2 ); /* Evaluate the interpolant */ value = 0.0; for ( i = deg; 0 <= i; i-- ) { t = 0.0; for ( j = 0; j <= i; j++ ) { t = t + ttg1[j] * c0[j+(deg-i)*(degmax+2)]; } value = value + ttg2[deg-i] * t; } free ( ttg1 ); free ( ttg2 ); return value; } /******************************************************************************/ void pdpts ( int deg, double pd1[], double pd2[], double wpd[], int *npd ) /******************************************************************************/ /* Purpose: PDPTS returns the points and weights for Padua interpolation. Discussion: This subroutine computes the first family of Padua points and weights corresponding to degree DEG. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 February 2014 Author: Original FORTRAN77 version by Marco Caliari, Stefano De Marchi, Marco Vianello. This C version by John Burkardt. Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Algorithm 886: Padua2D: Lagrange Interpolation at Padua Points on Bivariate Domains, ACM Transactions on Mathematical Software, Volume 35, Number 3, October 2008, Article 21, 11 pages. Parameters: Input, int DEG, the degree of approximation. Output, double PD1[NPD], PD2[NPD], the first and second coordinates of the Padua points Output, double WPD[NPD], the weights. Output, int *NPD, the number of Padua points. NPD = ( DEG + 1 ) * ( DEG + 2 ) / 2. */ { int itemp0; int j; int k; const double pi = 3.1415926535897931; double rtemp0; /* Compute the Padua points of the first family at degree DEG. */ if ( deg == 0 ) { pd1[0] = -1.0; pd2[0] = -1.0; wpd[0] = 2.0; *npd = 1; return; } *npd = 0; itemp0 = deg * ( deg + 1 ); rtemp0 = pi / ( double ) ( itemp0 ); for ( j = 0; j <= deg + 1; j++ ) { for ( k = ( j % 2 ); k <= deg; k = k + 2 ) { pd1[*npd] = - cos ( ( double ) ( ( deg + 1 ) * k ) * rtemp0 ); pd2[*npd] = - cos ( ( double ) ( deg * j ) * rtemp0 ); wpd[*npd] = 2.0 / ( double ) ( itemp0 ); if ( k == 0 || k == deg ) { wpd[*npd] = wpd[*npd] / 2.0; } if ( j == 0 || j == deg + 1 ) { wpd[*npd] = wpd[*npd] / 2.0; } *npd = *npd + 1; } } return; } /******************************************************************************/ double r8_huge ( ) /******************************************************************************/ /* Purpose: R8_HUGE returns a "huge" R8. Discussion: The value returned by this function is NOT required to be the maximum representable R8. This value varies from machine to machine, from compiler to compiler, and may cause problems when being printed. We simply want a "very large" but non-infinite number. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 October 2007 Author: John Burkardt Parameters: Output, double R8_HUGE, a "huge" R8 value. */ { double value; value = 1.0E+30; return value; } /******************************************************************************/ double r8_max ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MAX returns the maximum of two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, Y, the quantities to compare. Output, double R8_MAX, the maximum of X and Y. */ { double value; if ( y < x ) { value = x; } else { value = y; } return value; } /******************************************************************************/ double r8_min ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MIN returns the minimum of two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, Y, the quantities to compare. Output, double R8_MIN, the minimum of X and Y. */ { double value; if ( y < x ) { value = y; } else { value = x; } return value; } /******************************************************************************/ double r8_sign ( double x ) /******************************************************************************/ /* Purpose: R8_SIGN returns the sign of an R8. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 May 2006 Author: John Burkardt Parameters: Input, double X, the number whose sign is desired. Output, double R8_SIGN, the sign of X. */ { double value; if ( x < 0.0 ) { value = - 1.0; } else { value = + 1.0; } return value; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 June 2013 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }