version 1.11, 2015/08/06 23:41:52 |
version 1.14, 2017/09/04 01:57:53 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.10 2015/08/06 09:12:29 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.13 2015/09/14 05:47:09 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
Line 32 Num tobf(Num a,int prec) |
|
Line 32 Num tobf(Num a,int prec) |
|
mpq_t q; |
mpq_t q; |
BF d; |
BF d; |
N nm,dn; |
N nm,dn; |
|
C c; |
|
Num re,im; |
int sgn; |
int sgn; |
|
|
if ( !a ) { |
if ( !a ) { |
Line 67 Num tobf(Num a,int prec) |
|
Line 69 Num tobf(Num a,int prec) |
|
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
return (Num)d; |
return (Num)d; |
break; |
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: |
default: |
error("tobf : invalid argument"); |
error("tobf : invalid argument"); |
break; |
break; |
Line 85 void addbf(Num a,Num b,Num *c) |
|
Line 92 void addbf(Num a,Num b,Num *c) |
|
*c = b; |
*c = b; |
else if ( !b ) |
else if ( !b ) |
*c = a; |
*c = a; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) ) |
(*addnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*addnumt[MAX(NID(a),NID(b))])(a,b,c); |
else if ( NID(a) == N_B ) { |
else if ( NID(a) == N_B ) { |
switch ( NID(b) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
Line 108 void addbf(Num a,Num b,Num *c) |
|
Line 115 void addbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
|
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b) == N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
Line 128 void addbf(Num a,Num b,Num *c) |
|
Line 138 void addbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
|
break; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)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) |
void subbf(Num a,Num b,Num *c) |
Line 145 void subbf(Num a,Num b,Num *c) |
|
Line 163 void subbf(Num a,Num b,Num *c) |
|
(*chsgnnumt[NID(b)])(b,c); |
(*chsgnnumt[NID(b)])(b,c); |
else if ( !b ) |
else if ( !b ) |
*c = a; |
*c = a; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) ) |
(*subnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*subnumt[MAX(NID(a),NID(b))])(a,b,c); |
else if ( NID(a) == N_B ) { |
else if ( NID(a) == N_B ) { |
switch ( NID(b) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
Line 168 void subbf(Num a,Num b,Num *c) |
|
Line 186 void subbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b)==N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
Line 189 void subbf(Num a,Num b,Num *c) |
|
Line 209 void subbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
|
|
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)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) |
void mulbf(Num a,Num b,Num *c) |
Line 205 void mulbf(Num a,Num b,Num *c) |
|
Line 233 void mulbf(Num a,Num b,Num *c) |
|
|
|
if ( !a || !b ) |
if ( !a || !b ) |
*c = 0; |
*c = 0; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) ) |
(*mulnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*mulnumt[MAX(NID(a),NID(b))])(a,b,c); |
else if ( NID(a) == N_B ) { |
else if ( NID(a) == N_B ) { |
switch ( NID(b) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
Line 228 void mulbf(Num a,Num b,Num *c) |
|
Line 256 void mulbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b) == N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
Line 248 void mulbf(Num a,Num b,Num *c) |
|
Line 278 void mulbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)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) |
void divbf(Num a,Num b,Num *c) |
Line 265 void divbf(Num a,Num b,Num *c) |
|
Line 303 void divbf(Num a,Num b,Num *c) |
|
error("divbf : division by 0"); |
error("divbf : division by 0"); |
else if ( !a ) |
else if ( !a ) |
*c = 0; |
*c = 0; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) ) |
(*divnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*divnumt[MAX(NID(a),NID(b))])(a,b,c); |
else if ( NID(a) == N_B ) { |
else if ( NID(a) == N_B ) { |
switch ( NID(b) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
Line 288 void divbf(Num a,Num b,Num *c) |
|
Line 326 void divbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b)==N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
/* XXX : mpfr_z_div and mpfr_q_div are not implemented */ |
/* XXX : mpfr_z_div and mpfr_q_div are not implemented */ |
Line 304 void divbf(Num a,Num b,Num *c) |
|
Line 344 void divbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)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) |
void pwrbf(Num a,Num b,Num *c) |
Line 321 void pwrbf(Num a,Num b,Num *c) |
|
Line 369 void pwrbf(Num a,Num b,Num *c) |
|
*c = (Num)ONE; |
*c = (Num)ONE; |
else if ( !a ) |
else if ( !a ) |
*c = 0; |
*c = 0; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) ) |
(*pwrnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*pwrnumt[MAX(NID(a),NID(b))])(a,b,c); |
else if ( NID(a) == N_B ) { |
else if ( NID(a) == N_B ) { |
switch ( NID(b) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
Line 346 void pwrbf(Num a,Num b,Num *c) |
|
Line 394 void pwrbf(Num a,Num b,Num *c) |
|
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_roundmode); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b)==N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
mpfr_init2(r,BFPREC(b)); |
mpfr_init2(r,BFPREC(b)); |
Line 363 void pwrbf(Num a,Num b,Num *c) |
|
Line 413 void pwrbf(Num a,Num b,Num *c) |
|
a = tobf(a,prec); |
a = tobf(a,prec); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)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) |
void chsgnbf(Num a,Num *c) |
Line 376 void chsgnbf(Num a,Num *c) |
|
Line 434 void chsgnbf(Num a,Num *c) |
|
|
|
if ( !a ) |
if ( !a ) |
*c = 0; |
*c = 0; |
else if ( NID(a) <= N_A ) |
else if ( NID(a) <= N_R ) |
(*chsgnnumt[NID(a)])(a,c); |
(*chsgnnumt[NID(a)])(a,c); |
else { |
else if ( NID(a) == N_B ) { |
mpfr_init2(r,BFPREC(a)); |
mpfr_init2(r,BFPREC(a)); |
mpfr_neg(r,((BF)a)->body,mpfr_roundmode); |
mpfr_neg(r,((BF)a)->body,mpfr_roundmode); |
MPFRTOBF(r,d); |
MPFRTOBF(r,d); |
*c = (Num)d; |
*c = (Num)d; |
} |
} else |
|
error("chsgnbf : invalid argument"); |
} |
} |
|
|
int cmpbf(Num a,Num b) |
int cmpbf(Num a,Num b) |
Line 393 int cmpbf(Num a,Num b) |
|
Line 452 int cmpbf(Num a,Num b) |
|
GQ q; |
GQ q; |
|
|
if ( !a ) { |
if ( !a ) { |
if ( !b || (NID(b)<=N_A) ) |
if ( !b ) return 0; |
|
else if ( NID(b)<=N_R ) |
return (*cmpnumt[NID(b)])(a,b); |
return (*cmpnumt[NID(b)])(a,b); |
|
else if ( NID(b)==N_B ) |
|
return -mpfr_sgn(((BF)b)->body); |
else |
else |
return -mpfr_sgn(((BF)a)->body); |
goto err; |
} else if ( !b ) { |
} else if ( !b ) { |
if ( !a || (NID(a)<=N_A) ) |
if ( NID(a)<=N_R ) |
return (*cmpnumt[NID(a)])(a,b); |
return (*cmpnumt[NID(a)])(a,b); |
else |
else if ( NID(a)==N_B ) |
return mpfr_sgn(((BF)a)->body); |
return mpfr_sgn(((BF)a)->body); |
} else if ( NID(a) == N_B ) { |
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) ) { |
switch ( NID(b) ) { |
case N_Q: |
case N_Q: |
if ( INT((Q)b) ) { |
if ( INT((Q)b) ) { |
Line 420 int cmpbf(Num a,Num b) |
|
Line 486 int cmpbf(Num a,Num b) |
|
case N_B: |
case N_B: |
ret = mpfr_cmp(((BF)a)->body,((BF)b)->body); |
ret = mpfr_cmp(((BF)a)->body,((BF)b)->body); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
return ret; |
return ret; |
} else { /* NID(b)==N_B */ |
} else if ( NID(b)==N_B ) { |
switch ( NID(a) ) { |
switch ( NID(a) ) { |
case N_Q: |
case N_Q: |
if ( INT((Q)a) ) { |
if ( INT((Q)a) ) { |
Line 437 int cmpbf(Num a,Num b) |
|
Line 505 int cmpbf(Num a,Num b) |
|
/* double precision = 53 */ |
/* double precision = 53 */ |
ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body); |
ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body); |
break; |
break; |
|
default: |
|
goto err; |
} |
} |
return -ret; |
return -ret; |
} |
} |
|
err: error("cmpbf : cannot compare"); |
} |
} |