version 1.9, 2015/08/06 00:43:20 |
version 1.11, 2015/08/06 23:41:52 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.8 2015/08/05 03:46:35 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.10 2015/08/06 09:12:29 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
#include <math.h> |
#include <math.h> |
|
|
|
extern int mpfr_roundmode; |
|
|
Num tobf(Num,int); |
Num tobf(Num,int); |
|
|
#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) |
#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) |
Line 14 void strtobf(char *s,BF *p) |
|
Line 16 void strtobf(char *s,BF *p) |
|
BF r; |
BF r; |
NEWBF(r); |
NEWBF(r); |
mpfr_init(r->body); |
mpfr_init(r->body); |
mpfr_set_str(r->body,s,10,MPFR_RNDN); |
mpfr_set_str(r->body,s,10,mpfr_roundmode); |
*p = r; |
*p = r; |
} |
} |
|
|
double mpfrtodbl(mpfr_t a) |
double mpfrtodbl(mpfr_t a) |
{ |
{ |
return mpfr_get_d(a,MPFR_RNDN); |
return mpfr_get_d(a,mpfr_roundmode); |
} |
} |
|
|
Num tobf(Num a,int prec) |
Num tobf(Num a,int prec) |
Line 44 Num tobf(Num a,int prec) |
|
Line 46 Num tobf(Num a,int prec) |
|
break; |
break; |
case N_R: |
case N_R: |
prec ? mpfr_init2(r,prec) : mpfr_init(r); |
prec ? mpfr_init2(r,prec) : mpfr_init(r); |
mpfr_init_set_d(r,((Real)a)->body,MPFR_RNDN); |
mpfr_init_set_d(r,((Real)a)->body,mpfr_roundmode); |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
return (Num)d; |
return (Num)d; |
break; |
break; |
Line 54 Num tobf(Num a,int prec) |
|
Line 56 Num tobf(Num a,int prec) |
|
mpz_init(z); |
mpz_init(z); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
if ( sgn < 0 ) mpz_neg(z,z); |
if ( sgn < 0 ) mpz_neg(z,z); |
mpfr_init_set_z(r,z,MPFR_RNDN); |
mpfr_init_set_z(r,z,mpfr_roundmode); |
} else { |
} else { |
mpq_init(q); |
mpq_init(q); |
mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn)); |
mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn)); |
if ( sgn < 0 ) mpq_neg(q,q); |
if ( sgn < 0 ) mpq_neg(q,q); |
mpfr_init_set_q(r,q,MPFR_RNDN); |
mpfr_init_set_q(r,q,mpfr_roundmode); |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
return (Num)d; |
return (Num)d; |
Line 91 void addbf(Num a,Num b,Num *c) |
|
Line 93 void addbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
z = ztogz((Q)b); |
z = ztogz((Q)b); |
mpfr_add_z(r,((BF)a)->body,z->body,MPFR_RNDN); |
mpfr_add_z(r,((BF)a)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)b); |
q = qtogq((Q)b); |
mpfr_add_q(r,((BF)a)->body,q->body,MPFR_RNDN); |
mpfr_add_q(r,((BF)a)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_add_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN); |
mpfr_add_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode); |
break; |
break; |
case N_B: |
case N_B: |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 115 void addbf(Num a,Num b,Num *c) |
|
Line 117 void addbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
if ( INT((Q)a) ) { |
if ( INT((Q)a) ) { |
z = ztogz((Q)a); |
z = ztogz((Q)a); |
mpfr_add_z(r,((BF)b)->body,z->body,MPFR_RNDN); |
mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)a); |
q = qtogq((Q)a); |
mpfr_add_q(r,((BF)b)->body,q->body,MPFR_RNDN); |
mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_add_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN); |
mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 151 void subbf(Num a,Num b,Num *c) |
|
Line 153 void subbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
z = ztogz((Q)b); |
z = ztogz((Q)b); |
mpfr_sub_z(r,((BF)a)->body,z->body,MPFR_RNDN); |
mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)b); |
q = qtogq((Q)b); |
mpfr_sub_q(r,((BF)a)->body,q->body,MPFR_RNDN); |
mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN); |
mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode); |
break; |
break; |
case N_B: |
case N_B: |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 172 void subbf(Num a,Num b,Num *c) |
|
Line 174 void subbf(Num a,Num b,Num *c) |
|
} else { /* NID(b)==N_B */ |
} else { /* NID(b)==N_B */ |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
mpfr_init2(s,BFPREC(b)); |
|
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
if ( INT((Q)a) ) { |
if ( INT((Q)a) ) { |
z = ztogz((Q)a); |
z = ztogz((Q)a); |
mpfr_sub_z(s,((BF)b)->body,z->body,MPFR_RNDN); |
mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)a); |
q = qtogq((Q)a); |
mpfr_sub_q(s,((BF)b)->body,q->body,MPFR_RNDN); |
mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode); |
} |
} |
mpfr_neg(r,s,MPFR_RNDN); |
mpfr_neg(r,r,mpfr_roundmode); |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 212 void mulbf(Num a,Num b,Num *c) |
|
Line 213 void mulbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
z = ztogz((Q)b); |
z = ztogz((Q)b); |
mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN); |
mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)b); |
q = qtogq((Q)b); |
mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN); |
mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN); |
mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode); |
break; |
break; |
case N_B: |
case N_B: |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 236 void mulbf(Num a,Num b,Num *c) |
|
Line 237 void mulbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
if ( INT((Q)a) ) { |
if ( INT((Q)a) ) { |
z = ztogz((Q)a); |
z = ztogz((Q)a); |
mpfr_mul_z(r,((BF)b)->body,z->body,MPFR_RNDN); |
mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)a); |
q = qtogq((Q)a); |
mpfr_mul_q(r,((BF)b)->body,q->body,MPFR_RNDN); |
mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN); |
mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 272 void divbf(Num a,Num b,Num *c) |
|
Line 273 void divbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
z = ztogz((Q)b); |
z = ztogz((Q)b); |
mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN); |
mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
q = qtogq((Q)b); |
q = qtogq((Q)b); |
mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN); |
mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_init2(r,MAX(BFPREC(a),53)); |
mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN); |
mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode); |
break; |
break; |
case N_B: |
case N_B: |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 296 void divbf(Num a,Num b,Num *c) |
|
Line 297 void divbf(Num a,Num b,Num *c) |
|
/* XXX : mpfr_z_div and mpfr_q_div are not implemented */ |
/* XXX : mpfr_z_div and mpfr_q_div are not implemented */ |
a = tobf(a,BFPREC(b)); |
a = tobf(a,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_init2(r,MAX(BFPREC(b),53)); |
mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 328 void pwrbf(Num a,Num b,Num *c) |
|
Line 329 void pwrbf(Num a,Num b,Num *c) |
|
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
z = ztogz((Q)b); |
z = ztogz((Q)b); |
mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN); |
mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode); |
} else { |
} else { |
b = tobf(b,BFPREC(a)); |
b = tobf(b,BFPREC(a)); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
} |
} |
break; |
break; |
case N_R: |
case N_R: |
Line 339 void pwrbf(Num a,Num b,Num *c) |
|
Line 340 void pwrbf(Num a,Num b,Num *c) |
|
prec = MAX(BFPREC(a),53); |
prec = MAX(BFPREC(a),53); |
mpfr_init2(r,prec); |
mpfr_init2(r,prec); |
b = tobf(b,prec); |
b = tobf(b,prec); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
case N_B: |
case N_B: |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 353 void pwrbf(Num a,Num b,Num *c) |
|
Line 354 void pwrbf(Num a,Num b,Num *c) |
|
case N_Q: |
case N_Q: |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
a = tobf(a,BFPREC(b)); |
a = tobf(a,BFPREC(b)); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
case N_R: |
case N_R: |
/* double precision = 53 */ |
/* double precision = 53 */ |
prec = MAX(BFPREC(a),53); |
prec = MAX(BFPREC(a),53); |
mpfr_init2(r,prec); |
mpfr_init2(r,prec); |
a = tobf(a,prec); |
a = tobf(a,prec); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
Line 379 void chsgnbf(Num a,Num *c) |
|
Line 380 void chsgnbf(Num a,Num *c) |
|
(*chsgnnumt[NID(a)])(a,c); |
(*chsgnnumt[NID(a)])(a,c); |
else { |
else { |
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
mpfr_neg(r,((BF)a)->body,MPFR_RNDN); |
mpfr_neg(r,((BF)a)->body,mpfr_roundmode); |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} |
} |