tools/import-R

   1 #!/usr/local/bin/perl -w
   2 # -----------------------------------------------------------------------------
   3
   4 use strict;
   5
   6 my @files =
   7     ("dpq.h",
   8      "ftrunc.c",
   9      "dnorm.c",
  10      "pnorm.c",
  11      "qnorm.c",
  12      "plnorm.c",
  13      "qlnorm.c",
  14      "ppois.c",
  15      # "qpois.c",
  16      "stirlerr.c",
  17      "bd0.c",
  18      "dpois.c",
  19      "dgamma.c",
  20      "pgamma.c",
  21      # "qgamma.c", # Not good enough.
  22      "chebyshev.c",
  23      "lgammacor.c",
  24      "lbeta.c",
  25      "dt.c",
  26      "pt.c",
  27      "qt.c",
  28      # "pf.c",
  29      # "qf.c",
  30      "pchisq.c",
  31      "qchisq.c",
  32      "dweibull.c",
  33      "pweibull.c",
  34      "pbinom.c",
  35      "dbinom.c",
  36      "qbinom.c",
  37      "dnbinom.c",
  38      "pnbinom.c",
  39      "qnbinom.c",
  40      "dbeta.c",
  41      # "pbeta.c",
  42      # "qbeta.c",
  43      "dhyper.c",
  44      "phyper.c",
  45      "dexp.c",
  46      "pexp.c",
  47      "dgeom.c",
  48      "pgeom.c",
  49      "dcauchy.c",
  50      "pcauchy.c",
  51      "bessel.h",
  52      "bessel_i.c",
  53      "bessel_j.c",
  54      "bessel_k.c",
  55      "bessel_y.c",
  56      );
  57
  58 # These are for plugin fn-R.  They are so small it makes no sense to place
  59 # them elsewhere.
  60 push @files,
  61     (
  62      "dlnorm.c",
  63      "df.c",
  64      "dchisq.c",
  65      "qweibull.c",
  66      "qexp.c",
  67      "qgeom.c",
  68      );
  69
  70
  71 my $mathfunc = $ARGV[0];
  72 my $dir = $ARGV[1];
  73 my $subdir = "src/nmath";
  74
  75 unless (defined ($mathfunc) && -r $mathfunc &&
  76         defined ($dir) && -d "$dir/$subdir") {
  77     print STDERR "Usage: $0 mathfunc.c R-directory\n";
  78     exit 1;
  79 }
  80
  81
  82 my ($prefix,$postfix) = &read_mathfunc ($mathfunc);
  83
  84 print $prefix;
  85 print "\n";
  86 print "/* The following source code was imported from the R project.  */\n";
  87 print "/* It was automatically transformed by $0.  */\n";
  88 print "\n";
  89
  90 foreach my $file (@files) {
  91     my $cleandefs = ($file =~ /\.c$/);
  92     print "/* Imported $subdir/$file from R.  */\n";
  93     &import_file ("$dir/$subdir/$file", $cleandefs);
  94 }
  95 print $postfix;
  96
  97 # -----------------------------------------------------------------------------
  98
  99 sub import_file {
 100     my ($filename,$cleandefs) = @_;
 101
 102     my %defines = ();
 103     my $incomment = 0; # Stupid.
 104
 105     local (*FIL);
 106     open (*FIL, "<$filename") or
 107         die "$0: Cannot read $filename: $!\n";
 108   LINE:
 109     while (<FIL>) {
 110         if (/^\s*\#\s*ifndef\s+GNUMERIC_VERSION\b/ ... /^\s*\#\s*endif\b.*\bGNUMERIC_VERSION\b/) {
 111             next;
 112         }
 113
 114         if (/^\s*\#\s*define\s+([a-zA-Z0-9_]+)/) {
 115             $defines{$1} = 1;
 116         } elsif (/^\s*\#\s*undef\s+([a-zA-Z0-9_]+)/) {
 117             delete $defines{$1};
 118         } elsif (/^\s*\#\s*include\s+/) {
 119             # Ignore for now.
 120             next LINE;
 121         }
 122
 123         $_ .= <FIL> if /^static\s+double\s*$/;
 124
 125         s/\s+$//;
 126         if ($incomment) {
 127             $incomment = 0 if m|\*/|;
 128         } else {
 129             s/\bdouble\b/gnm_float/g;
 130             s/\bRboolean\b/gboolean/g;
 131
 132             s/^(\s*)(const\s+gnm_float\s+([a-zA-Z_][a-zA-Z0-9_]*)\s*\[\s*\d+\s*\]\s*=)/$1static $2/;
 133
 134             # Improve precision of "log(1-x)".
 135             s/\blog\s*\(\s*1\s*-\s*([a-zA-Z0-9_]+)\s*\)/gnm_log1p (-$1)/g;
 136
 137             # Improve precision of "log(1+x)".
 138             s/\blog\s*\(\s*1\s*\+\s*/gnm_log1p \(/g;
 139
 140             s/\bISNAN\b/gnm_isnan/g;
 141             s/\bR_(finite|FINITE)\b/gnm_finite/g;
 142             s/\b(sqrt|exp|log|pow|log1p|expm1|floor|ceil|sin|cos|sinh|tan)(\s|$|\()/gnm_$1$2/g;
 143             s/\bfabs\b/gnm_abs/g;
 144
 145             s/\bgnm_floor\s*\(\s*([a-z]+)\s*\+\s*1e-7\s*\)/gnm_fake_floor($1)/;
 146
 147             # Constants.
 148             s/\b(M_LN2|M_PI|M_PI_2|M_SQRT2|M_2PI)\b/$1gnum/g;
 149             s/\bDBL_(EPSILON|MIN|MAX)/GNM_$1/g;
 150             s/([-+]?\d*\.(\d{5,})([eE][-+]?\d+)?)/GNM_const\($1\)/g;
 151             s@(\d)\s*/\s*(\d+\.\d*)@$1 / GNM_const\($2\)@g;
 152
 153             # Fix constant quotients.
 154             s~([-+]?\d+\.\d*)(\s*/\s*)([-+e0-9.]+)~GNM_const\($1\)$2$3~;
 155
 156             # These are made static.
 157             s/^gnm_float\s+(pbeta_both|bd0|dpois_raw|chebyshev_eval|lgammacor|lbeta|pbeta_raw|dbinom_raw|bessel_[ijky]|bessel_[jy]_ex)\b/static gnm_float $1/;
 158             s/^int\s+(chebyshev_init)/static int $1/;
 159
 160             # Fix a bug...
 161             s/\(\(-37\.5193 < x\) \|\| \(x < 8\.2924\)\)/((-37.5193 < x) && (x < 8.2924))/;
 162
 163             # Optimization given our stupid gammafn.
 164             s|> 10|< 10| if /p and q are small: p <= q > 10/;
 165             s|gnm_log\(gammafn\(p\) \* \(gammafn\(q\) / gammafn\(p \+ q\)\)\)|gnm_lgamma (p) + gnm_lgamma (q) - gnm_lgamma (p + q)|;
 166
 167             s/dnorm4/dnorm/g;
 168             s/pnorm5/pnorm/g;
 169             s/qnorm5/qnorm/g;
 170
 171             s/\b(trunc|ftrunc)\b/gnm_trunc/g;
 172             s/\b(lgammafn|lgamma)\b/gnm_lgamma/g;
 173             s/\bML_NAN\b/gnm_nan/g;
 174             s/\bML_VALID\b/\!gnm_isnan/g;
 175             s/\bML_NEGINF\b/gnm_ninf/g;
 176             s/\bML_POSINF\b/gnm_pinf/g;
 177
 178             if ($filename !~ /\bpgamma\.c$/) {
 179                 # Improve precision of "lgammagnum(x+1)".
 180                 s/\bgnm_lgamma\s*\(([^()]+)\s*\+\s*1(\.0*)?\s*\)/lgamma1p ($1)/;
 181                 s/\bgamma_cody\s*\(1\.\s*\+\s*([^()]+)\s*\)/gnm_exp(lgamma1p ($1))/;
 182             }
 183             s/\bR_Log1_Exp\b/swap_log_tail/g;
 184
 185             if ($filename =~ m|/bessel_i\.c$|) {
 186                 s/\bgamma_cody\(empal\)/gnm_exp(lgamma1p(nu))/;
 187             }
 188
 189             s/\bgamma_cody\b/gnm_gamma/g;
 190
 191             if (/^(static\s+)?(gnm_float|int)\s+(chebyshev_init)\s*\(/ .. /^\}/) {
 192                 next unless s/^(static\s+)?(gnm_float|int)\s+([a-zA-Z0-9_]+)\s*\(.*//;
 193                 $_ = "/* Definition of function $3 removed.  */";
 194             }
 195
 196             if ($filename =~ m|/bessel_j\.c$|) {
 197                 if (/^\s*const static gnm_float fact\[/) {
 198                     while (!/;$/) {
 199                         $_ .= <FIL>;
 200                     }
 201                     $_ = '/* removed array fact */';
 202                 } else {
 203                     s/\bfact\s*\[([^][]*)\]/fact($1)/g;
 204                 }
 205             }
 206
 207             # Somewhat risky.
 208             s/\%([-0-9.]*)([efgEFG])/\%$1\" GNM_FORMAT_$2 \"/g;
 209
 210             s/int give_log/gboolean give_log/g;
 211             s/int log_p/gboolean log_p/g;
 212             s/int lower_tail/gboolean lower_tail/g;
 213
 214             # Fix pbinom
 215             s/\bpbeta\s*\(1\s*-\s*([^,]*),\s*([^,]*),\s*([^,]*),\s*([^,]*),\s*([^,]*)\)/pbeta ($1, $3, $2, !$4, $5)/;
 216
 217             # Fix pt.
 218             if ($filename =~ m|/pt\.c$|) {
 219                 s/(n > 4e5)/0 && $1/;
 220                 if (/(^.*\s*=\s*)pbeta\s*(\(.*\);)/) {
 221                     $_ = "$1 (n > x * x)\n" .
 222                         "\t? pbeta (x * x / (n + x * x), 0.5, n / 2, /*lower_tail*/0, log_p)\n" .
 223                         "\t: pbeta $2";
 224                 }
 225             }
 226
 227
 228             if ($filename =~ m|/qbeta\.c$| && /xinbta = 0\.5;/) {
 229                 s/0\.5/(xinbta < lower) ? gnm_sqrt (lower) : 1 - gnm_sqrt (lower)/;
 230             }
 231
 232             $incomment = 1 if m|/\*$|;
 233         }
 234
 235         print "$_\n";
 236     }
 237     close (*FIL);
 238
 239     if ($cleandefs && keys %defines) {
 240         print "/* Cleaning up done by $0:  */\n";
 241         foreach my $def (sort keys %defines) {
 242             print "#undef $def\n";
 243         }
 244     }
 245
 246     print "\n";
 247     print "/* ", ('-' x 72), " */\n";
 248 }
 249
 250 # -----------------------------------------------------------------------------
 251
 252 sub read_mathfunc {
 253     my ($filename) = @_;
 254
 255     my $prefix = '';
 256     my $postfix = '';
 257
 258     local (*FIL);
 259     open (*FIL, "<$filename") or
 260         die "$0: Cannot read $filename: $!\n";
 261     my $state = 'pre';
 262     while (<FIL>) {
 263         if ($state eq 'pre') {
 264             $prefix .= $_;
 265             $state = 'mid' if m"--- BEGIN MAGIC R SOURCE MARKER ---";
 266         }
 267         if ($state eq 'mid') {
 268             $state = 'post' if m"--- END MAGIC R SOURCE MARKER ---";
 269         }
 270         if ($state eq 'post') {
 271             $postfix .= $_;
 272         }
 273     }
 274     close (*FIL);
 275
 276     die "$0: wrong set of magic markers in $filename.\n" if $state ne 'post';
 277
 278     return ($prefix,$postfix);
 279 }
 280
 281 # -----------------------------------------------------------------------------