/*
* $OpenXM: OpenXM_contrib2/asir2018/engine/bf.c,v 1.1 2018/09/19 05:45:07 noro Exp $
*/
#include "ca.h"
#include "base.h"
#include <math.h>
extern int mpfr_roundmode;
Num tobf(Num,int);
#define BFPREC(a) (BDY((BF)(a))->_mpfr_prec)
void strtobf(char *s,BF *p)
{
BF r;
NEWBF(r);
mpfr_init(BDY(r));
mpfr_set_str(BDY(r),s,10,mpfr_roundmode);
*p = r;
}
double mpfrtodbl(mpfr_t a)
{
return mpfr_get_d(a,mpfr_roundmode);
}
Num tobf(Num a,int prec)
{
mpfr_t r;
mpz_t z;
mpq_t q;
BF d;
C c;
Num re,im;
int sgn;
if ( !a ) {
prec ? mpfr_init2(r,prec) : mpfr_init(r);
mpfr_set_zero(r,1);
MPFRTOBF(r,d);
return (Num)d;
} else {
switch ( NID(a) ) {
case N_B:
return a;
break;
case N_R:
prec ? mpfr_init2(r,prec) : mpfr_init(r);
mpfr_init_set_d(r,BDY((Real)a),mpfr_roundmode);
MPFRTOBF(r,d);
return (Num)d;
break;
case N_Q:
if ( INT(a) )
mpfr_init_set_z(r,BDY((Z)a),mpfr_roundmode);
else
mpfr_init_set_q(r,BDY((Q)a),mpfr_roundmode);
MPFRTOBF(r,d);
return (Num)d;
break;
case N_C:
re = tobf(((C)a)->r,prec); im = tobf(((C)a)->i,prec);
NEWC(c); c->r = re; c->i = im;
return (Num)c;
break;
default:
error("tobf : invalid argument");
return 0;
break;
}
}
}
void addbf(Num a,Num b,Num *c)
{
mpfr_t r;
BF d;
if ( !a )
*c = b;
else if ( !b )
*c = a;
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
(*addnumt[MAX(NID(a),NID(b))])(a,b,c);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
mpfr_init2(r,BFPREC(a));
if ( INT(b) )
mpfr_add_z(r,BDY((BF)a),BDY((Z)b),mpfr_roundmode);
else
mpfr_add_q(r,BDY((BF)a),BDY((Q)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(a),53));
mpfr_add_d(r,BDY((BF)a),BDY((Real)b),mpfr_roundmode);
break;
case N_B:
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
mpfr_add(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
break;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else if ( NID(b) == N_B ) {
switch ( NID(a) ) {
case N_Q:
mpfr_init2(r,BFPREC(b));
if ( INT(a) )
mpfr_add_z(r,BDY((BF)b),BDY((Z)a),mpfr_roundmode);
else
mpfr_add_q(r,BDY((BF)b),BDY((Q)a),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(b),53));
mpfr_add_d(r,BDY((BF)b),BDY((Real)a),mpfr_roundmode);
break;
default:
goto err;
break;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else
goto err;
if ( !cmpbf(*c,0) ) *c = 0;
return;
err: error("addbf : invalid argument");
}
void subbf(Num a,Num b,Num *c)
{
mpfr_t r,s;
BF d;
if ( !a )
(*chsgnnumt[NID(b)])(b,c);
else if ( !b )
*c = a;
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
(*subnumt[MAX(NID(a),NID(b))])(a,b,c);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
mpfr_init2(r,BFPREC(a));
if ( INT(b) )
mpfr_sub_z(r,BDY((BF)a),BDY((Z)b),mpfr_roundmode);
else
mpfr_sub_q(r,BDY((BF)a),BDY((Q)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(a),53));
mpfr_sub_d(r,BDY((BF)a),BDY((Real)b),mpfr_roundmode);
break;
case N_B:
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
mpfr_sub(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else if ( NID(b)==N_B ) {
switch ( NID(a) ) {
case N_Q:
mpfr_init2(r,BFPREC(b));
if ( INT(a) )
mpfr_sub_z(r,BDY((BF)b),BDY((Z)a),mpfr_roundmode);
else
mpfr_sub_q(r,BDY((BF)b),BDY((Q)a),mpfr_roundmode);
mpfr_neg(r,r,mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(b),53));
mpfr_d_sub(r,BDY((Real)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else
goto err;
if ( !cmpbf(*c,0) ) *c = 0;
return;
err: error("subbf : invalid argument");
}
void mulbf(Num a,Num b,Num *c)
{
mpfr_t r;
BF d;
int prec;
if ( !a || !b )
*c = 0;
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
(*mulnumt[MAX(NID(a),NID(b))])(a,b,c);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
mpfr_init2(r,BFPREC(a));
if ( INT(b) )
mpfr_mul_z(r,BDY((BF)a),BDY((Z)b),mpfr_roundmode);
else
mpfr_mul_q(r,BDY((BF)a),BDY((Q)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(a),53));
mpfr_mul_d(r,BDY((BF)a),BDY((Real)b),mpfr_roundmode);
break;
case N_B:
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
mpfr_mul(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else if ( NID(b) == N_B ) {
switch ( NID(a) ) {
case N_Q:
mpfr_init2(r,BFPREC(b));
if ( INT(a) )
mpfr_mul_z(r,BDY((BF)b),BDY((Z)a),mpfr_roundmode);
else
mpfr_mul_q(r,BDY((BF)b),BDY((Q)a),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(b),53));
mpfr_mul_d(r,BDY((BF)b),BDY((Real)a),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else
goto err;
if ( !cmpbf(*c,0) ) *c = 0;
return;
err: error("mulbf : invalid argument");
}
void divbf(Num a,Num b,Num *c)
{
mpfr_t s,r;
BF d;
if ( !b )
error("divbf : division by 0");
else if ( !a )
*c = 0;
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
(*divnumt[MAX(NID(a),NID(b))])(a,b,c);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
mpfr_init2(r,BFPREC(a));
if ( INT(b) )
mpfr_div_z(r,BDY((BF)a),BDY((Z)b),mpfr_roundmode);
else
mpfr_div_q(r,BDY((BF)a),BDY((Q)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(a),53));
mpfr_div_d(r,BDY((BF)a),BDY((Real)b),mpfr_roundmode);
break;
case N_B:
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
mpfr_div(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else if ( NID(b)==N_B ) {
switch ( NID(a) ) {
case N_Q:
/* XXX : mpfr_z_div and mpfr_q_div are not implemented */
a = tobf(a,BFPREC(b));
mpfr_init2(r,BFPREC(b));
mpfr_div(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
mpfr_init2(r,MAX(BFPREC(b),53));
mpfr_d_div(r,BDY((Real)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else
goto err;
if ( !cmpbf(*c,0) ) *c = 0;
return;
err: error("mulbf : invalid argument");
}
void pwrbf(Num a,Num b,Num *c)
{
int prec;
mpfr_t r;
BF d;
if ( !b )
*c = (Num)ONE;
else if ( !a )
*c = 0;
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
(*pwrnumt[MAX(NID(a),NID(b))])(a,b,c);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
mpfr_init2(r,BFPREC(a));
if ( INT(b) ) {
mpfr_pow_z(r,BDY((BF)a),BDY((Z)b),mpfr_roundmode);
} else {
b = tobf(b,BFPREC(a));
mpfr_pow(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
}
break;
case N_R:
/* double precision = 53 */
prec = MAX(BFPREC(a),53);
mpfr_init2(r,prec);
b = tobf(b,prec);
mpfr_pow(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
case N_B:
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
mpfr_pow(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else if ( NID(b)==N_B ) {
switch ( NID(a) ) {
case N_Q:
mpfr_init2(r,BFPREC(b));
a = tobf(a,BFPREC(b));
mpfr_pow(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
case N_R:
/* double precision = 53 */
prec = MAX(BFPREC(a),53);
mpfr_init2(r,prec);
a = tobf(a,prec);
mpfr_pow(r,BDY((BF)a),BDY((BF)b),mpfr_roundmode);
break;
default:
goto err;
}
MPFRTOBF(r,d);
*c = (Num)d;
} else
goto err;
if ( !cmpbf(*c,0) ) *c = 0;
return;
err: error("pwrbf : invalid argument");
}
void chsgnbf(Num a,Num *c)
{
mpfr_t r;
BF d;
if ( !a )
*c = 0;
else if ( NID(a) <= N_R )
(*chsgnnumt[NID(a)])(a,c);
else if ( NID(a) == N_B ) {
mpfr_init2(r,BFPREC(a));
mpfr_neg(r,BDY((BF)a),mpfr_roundmode);
MPFRTOBF(r,d);
*c = (Num)d;
} else
error("chsgnbf : invalid argument");
}
int cmpbf(Num a,Num b)
{
int ret;
if ( !a ) {
if ( !b ) return 0;
else if ( NID(b)<=N_R )
return (*cmpnumt[NID(b)])(a,b);
else if ( NID(b)==N_B )
return -mpfr_sgn(BDY((BF)b));
else
goto err;
} else if ( !b ) {
if ( NID(a)<=N_R )
return (*cmpnumt[NID(a)])(a,b);
else if ( NID(a)==N_B )
return mpfr_sgn(BDY((BF)a));
else
goto err;
} else if ( NID(a) <= N_R && NID(b) <= N_R )
return (*cmpnumt[MAX(NID(a),NID(b))])(a,b);
else if ( NID(a) == N_B ) {
switch ( NID(b) ) {
case N_Q:
if ( INT(b) )
ret = mpfr_cmp_z(BDY((BF)a),BDY((Z)b));
else
ret = mpfr_cmp_q(BDY((BF)a),BDY((Q)b));
break;
case N_R:
/* double precision = 53 */
ret = mpfr_cmp_d(BDY((BF)a),BDY((Real)b));
break;
case N_B:
ret = mpfr_cmp(BDY((BF)a),BDY((BF)b));
break;
default:
goto err;
}
return ret;
} else if ( NID(b)==N_B ) {
switch ( NID(a) ) {
case N_Q:
if ( INT(a) )
ret = mpfr_cmp_z(BDY((BF)b),BDY((Z)a));
else
ret = mpfr_cmp_q(BDY((BF)b),BDY((Q)a));
break;
case N_R:
/* double precision = 53 */
ret = mpfr_cmp_d(BDY((BF)b),BDY((Real)a));
break;
default:
goto err;
}
return -ret;
}
err: error("cmpbf : cannot compare");
return 0;
}