htswanalysis/MACS/lib/gsl/gsl-1.11/dht/test.c

   1 /* dht/test_dht.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman
  21  */
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <stdio.h>
  25 #include <math.h>
  26 #include <gsl/gsl_ieee_utils.h>
  27 #include <gsl/gsl_test.h>
  28 #include <gsl/gsl_dht.h>
  29
  30
  31 /* Test exact small transform.
  32  */
  33 int
  34 test_dht_exact(void)
  35 {
  36   int stat = 0;
  37   double f_in[3] = { 1.0, 2.0, 3.0 };
  38   double f_out[3];
  39   gsl_dht * t = gsl_dht_new(3, 1.0, 1.0);
  40   gsl_dht_apply(t, f_in, f_out);
  41
  42   /* Check values. */
  43   if(fabs( f_out[0]-( 0.375254649407520))/0.375254649407520 > 1.0e-14) stat++;
  44   if(fabs( f_out[1]-(-0.133507872695560))/0.133507872695560 > 1.0e-14) stat++;
  45   if(fabs( f_out[2]-( 0.044679925143840))/0.044679925143840 > 1.0e-14) stat++;
  46
  47
  48   /* Check inverse.
  49    * We have to adjust the normalization
  50    * so we can use the same precalculated transform.
  51    */
  52   gsl_dht_apply(t, f_out, f_in);
  53   f_in[0] *= 13.323691936314223*13.323691936314223;  /* jzero[1,4]^2 */
  54   f_in[1] *= 13.323691936314223*13.323691936314223;
  55   f_in[2] *= 13.323691936314223*13.323691936314223;
  56
  57   /* The loss of precision on the inverse
  58    * is a little surprising. However, this
  59    * thing is quite tricky since the band-limited
  60    * function represented by the samples {1,2,3}
  61    * need not be very nice. Like in any spectral
  62    * application, you really have to have some
  63    * a-priori knowledge of the underlying function.
  64    */
  65   if(fabs( f_in[0]-1.0)/1.0 > 2.0e-05) stat++;
  66   if(fabs( f_in[1]-2.0)/2.0 > 2.0e-05) stat++;
  67   if(fabs( f_in[2]-3.0)/3.0 > 2.0e-05) stat++;
  68
  69   gsl_dht_free(t);
  70
  71   return stat;
  72 }
  73
  74
  75
  76 /* Test the transform
  77  * Integrate[x J_0(a x) / (x^2 + 1), {x,0,Inf}] = K_0(a)
  78  */
  79 int
  80 test_dht_simple(void)
  81 {
  82   int stat = 0;
  83   int n;
  84   double f_in[128];
  85   double f_out[128];
  86   gsl_dht * t = gsl_dht_new(128, 0.0, 100.0);
  87
  88   for(n=0; n<128; n++) {
  89     const double x = gsl_dht_x_sample(t, n);
  90     f_in[n] = 1.0/(1.0+x*x);
  91   }
  92
  93   gsl_dht_apply(t, f_in, f_out);
  94
  95   /* This is a difficult transform to calculate this way,
  96    * since it does not satisfy the boundary condition and
  97    * it dies quite slowly. So it is not meaningful to
  98    * compare this to high accuracy. We only check
  99    * that it seems to be working.
 100    */
 101   if(fabs( f_out[0]-4.00)/4.00 > 0.02) stat++;
 102   if(fabs( f_out[5]-1.84)/1.84 > 0.02) stat++;
 103   if(fabs(f_out[10]-1.27)/1.27 > 0.02) stat++;
 104   if(fabs(f_out[35]-0.352)/0.352 > 0.02) stat++;
 105   if(fabs(f_out[100]-0.0237)/0.0237 > 0.02) stat++;
 106
 107   gsl_dht_free(t);
 108
 109   return stat;
 110 }
 111
 112
 113 /* Test the transform
 114  * Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
 115  */
 116 int
 117 test_dht_exp1(void)
 118 {
 119   int stat = 0;
 120   int n;
 121   double f_in[128];
 122   double f_out[128];
 123   gsl_dht * t = gsl_dht_new(128, 1.0, 20.0);
 124
 125   for(n=0; n<128; n++) {
 126     const double x = gsl_dht_x_sample(t, n);
 127     f_in[n] = exp(-x);
 128   }
 129
 130   gsl_dht_apply(t, f_in, f_out);
 131
 132   /* Spot check.
 133    * Note that the systematic errors in the calculation
 134    * are quite large, so it is meaningless to compare
 135    * to a high accuracy.
 136    */
 137   if(fabs( f_out[0]-0.181)/0.181 > 0.02) stat++;
 138   if(fabs( f_out[5]-0.357)/0.357 > 0.02) stat++;
 139   if(fabs(f_out[10]-0.211)/0.211 > 0.02) stat++;
 140   if(fabs(f_out[35]-0.0289)/0.0289 > 0.02) stat++;
 141   if(fabs(f_out[100]-0.00221)/0.00211 > 0.02) stat++;
 142
 143   gsl_dht_free(t);
 144
 145   return stat;
 146 }
 147
 148
 149 /* Test the transform
 150  * Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
 151  */
 152 int
 153 test_dht_poly1(void)
 154 {
 155   int stat = 0;
 156   int n;
 157   double f_in[128];
 158   double f_out[128];
 159   gsl_dht * t = gsl_dht_new(128, 1.0, 1.0);
 160
 161   for(n=0; n<128; n++) {
 162     const double x = gsl_dht_x_sample(t, n);
 163     f_in[n] = x * (1.0 - x*x);
 164   }
 165
 166   gsl_dht_apply(t, f_in, f_out);
 167
 168   /* Spot check. This function satisfies the boundary condition,
 169    * so the accuracy should be ok.
 170    */
 171   if(fabs( f_out[0]-0.057274214)/0.057274214    > 1.0e-07) stat++;
 172   if(fabs( f_out[5]-(-0.000190850))/0.000190850 > 1.0e-05) stat++;
 173   if(fabs(f_out[10]-0.000024342)/0.000024342    > 1.0e-04) stat++;
 174   if(fabs(f_out[35]-(-4.04e-07))/4.04e-07       > 1.0e-03) stat++;
 175   if(fabs(f_out[100]-1.0e-08)/1.0e-08           > 0.25)    stat++;
 176
 177   gsl_dht_free(t);
 178
 179   return stat;
 180 }
 181
 182
 183 int main()
 184 {
 185   gsl_ieee_env_setup ();
 186
 187   gsl_test( test_dht_exact(),   "Small Exact DHT");
 188   gsl_test( test_dht_simple(),  "Simple  DHT");
 189   gsl_test( test_dht_exp1(),    "Exp  J1 DHT");
 190   gsl_test( test_dht_poly1(),   "Poly J1 DHT");
 191
 192   exit (gsl_test_summary());
 193 }