octave_packages/gsl-1.0.8/test_hyperg.m

   1 ## Copyright (C) 2008   Raymond E. Rogers
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ##  any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, see .
  15
  16 #
  17 #
  18 # Testing initialization/header for special functions
  19 # Octave format
  20 # It is intended to read test_hyperg.c from the GNU GSL hyperg test file
  21 # and use the test lines.  If you don't want to duplicate thier
  22 # errors, comment out the bad tests with //
  23 # Compare test_hyperg.c vs. test_hyperg_corr.c
  24 #
  25 # R. Rogers v1, testing only hyperg_1F1 and hyperg_U in double
  26 #       KummerM and KummerU
  27 #       Some array verification is done; since I messed it up once.
  28
  29
  30 global GSL_DBL_EPSILON=        2.2204460492503131e-16
  31 global GSL_SQRT_DBL_EPSILON=   1.4901161193847656e-08
  32 global TEST_TOL0= (2.0*GSL_DBL_EPSILON)
  33 global TEST_TOL1=  (16.0*GSL_DBL_EPSILON)
  34 global TEST_TOL2=  (256.0*GSL_DBL_EPSILON)
  35 global TEST_TOL3=  (2048.0*GSL_DBL_EPSILON)
  36 global TEST_TOL4=  (16384.0*GSL_DBL_EPSILON)
  37 global TEST_TOL5=  (131072.0*GSL_DBL_EPSILON)
  38 global TEST_TOL6=  (1048576.0*GSL_DBL_EPSILON)
  39 global TEST_SQRT_TOL0= (2.0*GSL_SQRT_DBL_EPSILON)
  40 global TEST_SNGL=  (1.0e-06)
  41
  42 global TEST_SF_INCONS=  1
  43 global TEST_SF_ERRNEG=  2
  44 global TEST_SF_TOLBAD=  4
  45 global TEST_SF_RETBAD=  8
  46 global TEST_SF_ERRBAD=  16
  47 global TEST_SF_ERRBIG=  32
  48 global TEST_SF_EXPBAD=  64
  49 global TEST_FACTOR=     100
  50
  51
  52 function res=READ_TEST_SF_hyperg_Kummer(input_name)
  53 global GSL_DBL_EPSILON
  54 global GSL_SQRT_DBL_EPSILON
  55 global TEST_TOL0
  56 global TEST_TOL1
  57 global TEST_TOL2
  58 global TEST_TOL3
  59 global TEST_TOL4
  60 global TEST_TOL5
  61 global TEST_TOL6
  62 global TEST_SQRT_TOL0
  63 global TEST_SNGL
  64
  65 global TEST_SF_INCONS
  66 global TEST_SF_ERRNEG
  67 global TEST_SF_TOLBAD
  68 global TEST_SF_RETBAD
  69 global TEST_SF_ERRBAD
  70 global TEST_SF_ERRBIG
  71 global TEST_SF_EXPBAD
  72
  73 global TEST_FACTOR
  74         [source_id,source_msg]=fopen(input_name,"r","native")
  75
  76
  77                 while (! feof(source_id))
  78                         do
  79                                 input_line=fgetl(source_id);
  80                         until( index(input_line,"//") == 0);
  81
  82                         str_p=index(input_line,"hyperg_1F1_e");
  83                         if (str_p != 0)
  84                                 # Take it apart
  85                                 string_split=char(strsplit(input_line,","));
  86                                 arg1=str2double(substr(string_split(3,:),3));
  87                                 arg2=str2double(string_split(4,:));
  88                                 arg3=str2double(string_split(5,:));
  89                                 val=str2double(string_split(7,:));
  90                                 tol=eval(string_split(8,:));
  91                                 [ret,err]=hyperg_1F1(arg1,arg2,arg3);
  92                                 # This is to prevent extanious stops on some errors
  93                                 if ret==NaN
  94                                         ret=Inf
  95                                 endif
  96                                 if (abs((ret-val)/val)<tol*TEST_FACTOR)
  97                                         printf("\n Passed KummerM: %e %e %e ",arg1,arg2,arg3)
  98                                 else
  99                                         printf("\n Failed KummerM: %s\n value=%e, computed=%e, tol=%e, returned error=%e ",input_line,val,ret,tol,err)
 100                                         printf("\n error %e", abs((ret-val)/val))
 101                                 endif
 102                         endif
 103                         str_p=index(input_line,"hyperg_U_e");
 104                         if (str_p != 0)
 105                                 # Take it apart
 106                                 string_split=char(strsplit(input_line,","));
 107                                 arg1=str2double(substr(string_split(3,:),3));
 108                                 arg2=str2double(string_split(4,:));
 109                                 arg3=str2double(string_split(5,:));
 110                                 val=str2double(string_split(7,:));
 111                                 tol=eval(string_split(8,:));
 112                                 [ret,err]=hyperg_U(arg1,arg2,arg3);
 113                                 # This is to prevent extanious stops on some errors
 114                                 if ret==NaN
 115                                         ret=Inf
 116                                 endif
 117                                 if (abs((ret-val)/val)<tol*TEST_FACTOR)
 118                                         printf("\n Passed KummerU: %e %e %e ",arg1,arg2,arg3)
 119                                 else
 120                                         printf("\n Failed KummerU: %s\n value=%e, computed=%e, tol=%e, returned error=%e ",input_line,val,ret,tol,err)
 121                                         printf("\n error %e", abs((ret-val)/val))
 122                                 endif
 123                         endif
 124                 endwhile
 125         fclose(source_id);
 126         #lets do some array tests
 127         disp("\n array tests")
 128         mat1=[1,2]
 129         disp("\n hyperg_1F1(1,1,mat1)")
 130         hyperg_1F1(1,1,mat1)
 131         disp("\n hyperg_1F1(1,mat1,1) ")
 132         hyperg_1F1(1,mat1,1)
 133         disp("\n hyperg_1F1(mat1,1,1)")
 134         hyperg_1F1(mat1,1,1)
 135         disp("\n hyperg_1F1(mat1,mat1,mat1) ")
 136         hyperg_1F1(mat1,mat1,mat1)
 137         disp("\n hyperg_U(1,1,mat1)")
 138         hyperg_U(1,1,mat1)
 139         disp("\n hyperg_U(1,mat1,1)")
 140         hyperg_U(1,mat1,1)
 141         disp("\n hyperg_U(mat1,1,1)")
 142         hyperg_U(mat1,1,1)
 143         disp("\n hyperg_U(mat1,mat1,mat1)")
 144         hyperg_U(mat1,mat1,mat1)
 145
 146         res="";
 147 end
 148 #file_name=input("file name: ","s")
 149 file_name="test_hyperg_corr.c"
 150 READ_TEST_SF_hyperg_Kummer(file_name);
 151
 152