# rounding mode by wrapping (in the rndna XSub) 
# the mpfr_round_nearest_away() macro which is
# provided by the current devel (svn)
# version of the mpfr library. Avoid functions
# that take other than mpfr_t arguments (apart
# from the rounding mode argument).
# But allow mpfr_strtofr.

use strict;
use warnings;
use Math::MPFR qw(:mpfr);

use Inline C => Config =>
  LIBS => '-lmpfr -lgmp',
  BUILD_NOISY => 1;

use Inline C => <<'EOC';

#include <mpfr.h>


/*
 Create the type quad_func_t for pointer
 to a function that takes 3 mpfr_t arguments
 and 1 mp_rnd_t argument, and returns int
*/

typedef int
(*quad_func_t)(mpfr_t, mpfr_t, mpfr_t, mp_rnd_t);

/*
 Create the type trip_func_t for pointer
 to a function that takes 2 mpfr_t arguments
 and 1 mp_rnd_t argument, and returns int
*/

typedef int
(*trip_func_t)(mpfr_t, mpfr_t, mp_rnd_t);

/*
 Create the type str_func_t for pointer
 to a function that takes 1 mpfr_t argument,
 1 char * argument, 1 char** argument, 1 int
 argument and 1 mp_rnd_t argument, and
 returns int
*/

typedef int
(*str_func_t)(mpfr_t, char*, char**, int, mp_rnd_t);

SV* get_mpfr_mul_func_ptr() {
    return newSViv(PTR2IV(mpfr_mul));
}

SV* get_mpfr_add_func_ptr() {
    return newSViv(PTR2IV(mpfr_add));
}

SV* get_mpfr_div_func_ptr() {
    return newSViv(PTR2IV(mpfr_div));
}

SV* get_mpfr_sub_func_ptr() {
    return newSViv(PTR2IV(mpfr_sub));
}

SV* get_mpfr_pow_func_ptr() {
    return newSViv(PTR2IV(mpfr_pow));
}

SV* get_mpfr_sin_func_ptr() {
    return newSViv(PTR2IV(mpfr_sin));
}

SV* get_mpfr_cos_func_ptr() {
    return newSViv(PTR2IV(mpfr_cos));
}

SV* get_mpfr_log_func_ptr() {
    return newSViv(PTR2IV(mpfr_log));
}

SV* get_mpfr_set_func_ptr() {
    return newSViv(PTR2IV(mpfr_set));
}

SV* get_mpfr_strtofr_func_ptr() {
    return newSViv(PTR2IV(mpfr_strtofr));
}

void rndna(SV * func_ptr, ...) {
  dXSARGS;
  int ret;
  IV temp = SvIV(ST(0));

  if(items != 3 && items != 4)
    croak("items == %d, not 3 or 4", items);

  /* minimum exponent needs to be set to at
   least 1 more than minimum allowable */
  if(mpfr_get_emin_min() >= mpfr_get_emin())
    mpfr_set_emin(mpfr_get_emin_min() + 1);

  if(items == 3) {
    trip_func_t foo = INT2PTR(trip_func_t, temp);
    ret = mpfr_round_nearest_away(
      foo,
      *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(1))))),
      *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(2)))))
                                  );
  }

  if(items == 4) {
    if(SvPOK(ST(2))) { /* mpfr_strtofr */
      str_func_t foo = INT2PTR(str_func_t, temp);
      ret = mpfr_round_nearest_away(
        foo,
        *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(1))))),
        SvPV_nolen(ST(2)),
        NULL,
        SvIV(ST(3))
                                    );
    }
    else {
      quad_func_t foo = INT2PTR(quad_func_t, temp);
      ret = mpfr_round_nearest_away(
        foo,
        *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(1))))),
        *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(2))))),
        *(INT2PTR(mpfr_t *, SvIVX(SvRV(ST(3)))))
                                    );
    }
  }

  sp = mark;
  XPUSHs(sv_2mortal(newSViv(ret)));
  XSRETURN(1);
}
                  
EOC

my $rop = Math::MPFR->new('1.3');
my $op = Math::MPFR->new('1.2');
my $ret;

my $func_ptr = get_mpfr_mul_func_ptr();

# $rop = $rop * $op;
$ret = rndna($func_ptr, $rop, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_add_func_ptr();

# $rop = $rop + $op;
$ret = rndna($func_ptr, $rop, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_div_func_ptr();

# $rop = $rop / $op;
$ret = rndna($func_ptr, $rop, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_sub_func_ptr();

# $rop = $rop - $op;
$ret = rndna($func_ptr, $rop, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_pow_func_ptr();

# $rop = $rop ** $op;
$ret = rndna($func_ptr, $rop, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_sin_func_ptr();

# $rop = sin($op);
$ret = rndna($func_ptr, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_cos_func_ptr();

# $rop = cos($op);
$ret = rndna($func_ptr, $rop, $op);
print "$ret $rop $op\n";

$func_ptr = get_mpfr_log_func_ptr();

# $rop = log($op);
$ret = rndna($func_ptr, $rop, $op);
print "$ret $rop $op\n";

eval {$ret = rndna($rop, $op);};

if($@ =~ / not 3 or 4/) {
  print "Trapped error is as expected\n";
}
else {
  print "Unexpected error: $@";
}

my $check = Rmpfr_init2(109);
my $p_109 = Rmpfr_init2(109);
my $p_110 = Rmpfr_init2(110);

#This string can be expressed exactly in 110 bits, with last bit = 1
my $str = "0.1387778780781445675529539585113525390625e31";

Rmpfr_set_str($check, $str, 10, MPFR_RNDN);
$ret = Rmpfr_strtofr($p_110, $str, 10, MPFR_RNDN);

my $broken;
if($ret) {
 $broken = 1; # mpfr_strtofr bug, since fixed.
              # Result is exact, should have returned 0
 print "WARNING: mpfr_strtofr() returned incorrect ternary value\n";
}

$func_ptr = get_mpfr_set_func_ptr();

$ret = rndna($func_ptr, $p_109, $p_110);
print $ret, "\n";

#$p_109 should be rounded away from zero, but
# $check should not. Therefore:

if($p_109 > $check) {print "ok 1\n"}
else {print "not ok 1\n"}

$check *= -1;
$p_110 *= -1;

$ret = rndna($func_ptr, $p_109, $p_110);
print $ret, "\n";

#$p_109 should be rounded away from zero, but
# $check should not. Therefore:

if($p_109 < $check) {print "ok 2\n"}
else {print "not ok 2\n"}

$func_ptr = get_mpfr_strtofr_func_ptr();

# clear $p_109:
Rmpfr_set_ui($p_109, 42, MPFR_RNDN);

$ret = rndna($func_ptr, $p_109, $str, 10);
print $ret, "\n";

$check *= -1; # return to its original +ve value.

if($broken) {
  if($p_109 == $check) {print "ok 3 - given presence of mpfr_strtofr b
+ug\n"}
  else {print "not ok 3 - given presence of mpfr_strtofr bug\n"}
}
else {
  if($p_109 > $check) {print "ok 3 - given absence of mpfr_strtofr bug
+\n"}
  else {print "not ok 3 - given absence of mpfr_strtofr bug\n"}
}
[download]

1 1.5600000000000001 1.2
1 2.7600000000000002 1.2
-1 2.3000000000000003 1.2
0 1.1000000000000003 1.2
-1 1.1211693641406026 1.2
-1 9.3203908596722629e-1 1.2
1 3.6235775447667362e-1 1.2
1 1.8232155679395459e-1 1.2
Trapped error is as expected
WARNING: mpfr_strtofr() returned incorrect ternary value
1
ok 1
-1
ok 2
-1
ok 3 - given presence of mpfr_strtofr bug
[download]