libit

00001 /*
00002    libit - Library for basic source and channel coding functions
00003    Copyright (C) 2005-2005 Vivien Chappelier, Herve Jegou
00004 
00005    This library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Library General Public
00007    License as published by the Free Software Foundation; either
00008    version 2 of the License, or (at your option) any later version.
00009 
00010    This library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Library General Public License for more details.
00014 
00015    You should have received a copy of the GNU Library General Public
00016    License along with this library; if not, write to the Free
00017    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00018 */
00019 
00020 /*------------------------------------------------------------*/
00021 /** @file test_math.c Test program for Math functions         */
00022 /*------------------------------------------------------------*/
00023 
00024 #include <stdio.h>
00025 #include <it/math.h>
00026 #include <it/types.h>
00027 #include <it/io.h>
00028 
00029 #define SIGMA 3.0
00030 #define INTEGRAL_STEP 100000
00031 
00032 /* the function returning one on one plus x^2 (integrates to arctan(x)) */
00033 it_function (one_on_one_plus_x_square)
00034 {
00035   return (1. / (1. + x * x));
00036 }
00037 
00038 /* a function returning x at the power of 4 */
00039 it_function (polynomial)
00040 {
00041   return (x * x * x * x);
00042 }
00043 
00044 /* the gaussian function with parameter sigma */
00045 it_function_args (gaussian_pdf)
00046 {
00047   double sigma;     /* standard deviation */
00048 };
00049 
00050 it_function (gaussian_pdf)
00051 {
00052   double sigma = it_this->sigma;
00053 
00054   return (1. / (sqrt (2. * M_PI) * sigma) *
00055     exp (-x * x / (2. * sigma * sigma)));
00056 }
00057 
00058 int main ()
00059 {
00060   it_function_args (gaussian_pdf) gaussian_args;
00061   it_function_args (itf_differentiate) differentiate_args;
00062   it_function_args (itf_differentiate) differentiate_args2;
00063   it_function_args (itf_diff2) diff2_args;
00064   it_function_args (itf_sum) sum_args;
00065   it_function_args (itf_mul) mul_args;
00066   it_function_args (itf_compose) compose_args;
00067   it_function_args (itf_integrate_romberg) integrate2_args;
00068   it_function_args (itf_integrate_trapezoid) integrate3_args;
00069 
00070   /* parametric function call */
00071   gaussian_args.sigma = SIGMA;
00072   it_printf ("gaussian pdf in 1 = %f\n", gaussian_pdf (1.0, &gaussian_args));
00073 
00074   /* function sum */
00075   sum_args.f = itf_identity;
00076   sum_args.f_args = NULL;
00077   sum_args.g = polynomial;
00078   sum_args.g_args = NULL;
00079   it_printf ("x + x^4 in 2 = %f (%f)\n",
00080        itf_sum (2, &sum_args), 2. + 2. * 2. * 2. * 2.);
00081 
00082   /* function product */
00083   mul_args.f = itf_identity;
00084   mul_args.f_args = NULL;
00085   mul_args.g = one_on_one_plus_x_square;
00086   mul_args.g_args = NULL;
00087   it_printf ("x / (1 + x^2) in 2 = %f (%f)\n",
00088        itf_mul (2, &mul_args), 2. / (1 + 2. * 2.));
00089 
00090   /* function composition */
00091   compose_args.f = gaussian_pdf;
00092   compose_args.f_args = &gaussian_args;
00093   compose_args.g = one_on_one_plus_x_square;
00094   compose_args.g_args = NULL;
00095   it_printf ("gaussian(1/(1 + x^2)) in 2 = %f (%f)\n",
00096        itf_compose (2, &compose_args),
00097        (1. / (sqrt (2. * M_PI) * SIGMA) *
00098         exp (-1. / (50. * SIGMA * SIGMA))));
00099 
00100   /* derivatives */
00101   it_printf ("derivative of arctan in 2 = %f (%f)\n",
00102        it_differentiate (IT_FUNCTION (atan), NULL, 2.0),
00103        1. / (1. + 2. * 2.));
00104   it_printf ("derivative of gaussian pdf in 2 = %.10f (%.10f)\n",
00105        it_differentiate (gaussian_pdf, &gaussian_args, 2.0),
00106        -1. * 2. / (SIGMA * SIGMA) * gaussian_pdf (2.0, &gaussian_args));
00107 
00108   /* composing operators */
00109   differentiate_args.function = polynomial;
00110   differentiate_args.args = NULL;
00111   differentiate_args2.function = itf_differentiate;
00112   differentiate_args2.args = &differentiate_args;
00113   it_printf ("d of x^4 in 2 = %f (%f)\n",
00114        itf_differentiate (2.0, &differentiate_args), 4. * 2. * 2. * 2.);
00115   it_printf ("d^2 of x^4 in 2 = %.03f (%.03f)\n",
00116        itf_differentiate (2, &differentiate_args2), 12. * 2. * 2.);
00117   diff2_args.function = polynomial;
00118   diff2_args.args = NULL;
00119   it_printf ("d^2 of x^4 in 2 = %.04f (%.04f)\n", itf_diff2 (2, &diff2_args), 12. * 2. * 2.); /* smaller error */
00120 
00121   /* integrals */ 
00122     it_printf ("integral (default) of 1/(1+x^2) in [0,1] = %.12f (%.12f)\n",
00123          it_integrate (one_on_one_plus_x_square, NULL, 0.0, 1.0),
00124          atan (1));
00125 
00126   integrate2_args.function = one_on_one_plus_x_square;
00127   integrate2_args.args = NULL;
00128   integrate2_args.a = 0.0;
00129   integrate2_args.N = 10;
00130   integrate2_args.p = 14;
00131   it_printf ("integral (romberg) of 1/(1+x^2) in [0,1] = %.12f (%.12f)\n",
00132        itf_integrate_romberg (1.0, &integrate2_args), atan (1));
00133 
00134   integrate3_args.function = one_on_one_plus_x_square;
00135   integrate3_args.args = NULL;
00136   integrate3_args.a = 0.0;
00137   integrate3_args.N = 100000;
00138   it_printf ("integral (trapezoid) of 1/(1+x^2) in [0,1] = %.12f (%.12f)\n",
00139        itf_integrate_trapezoid (1.0, &integrate3_args), atan (1));
00140 
00141 
00142   it_printf ("erfinv(-0.2) = %f\n", erfinv (-0.2));
00143 
00144   return (0);
00145 }
examples/test_math/test_math.c
