=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/bf.c,v retrieving revision 1.1.1.1 retrieving revision 1.7 diff -u -p -r1.1.1.1 -r1.7 --- OpenXM_contrib2/asir2000/engine/bf.c 1999/12/03 07:39:08 1.1.1.1 +++ OpenXM_contrib2/asir2000/engine/bf.c 2015/08/04 06:20:45 1.7 @@ -1,151 +1,213 @@ -/* $OpenXM: OpenXM/src/asir99/engine/bf.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */ +/* + * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.6 2009/07/15 02:19:47 noro Exp $ + */ #include "ca.h" -#if PARI #include "base.h" #include -#include "genpari.h" -extern long prec; +Num tobf(Num,int); -void ritopa(Obj,GEN *); -void patori(GEN,Obj *); +#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) -void addbf(a,b,c) -Num a,b; -Num *c; +void strtobf(char *s,BF *p) { - GEN pa,pb,z; - long ltop,lbot; + BF r; + NEWBF(r); + mpfr_init(r->body); + mpfr_set_str(r->body,s,10,MPFR_RNDN); + *p = r; +} - if ( !a ) - *c = b; - else if ( !b ) - *c = a; - else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) - (*addnumt[MIN(NID(a),NID(b))])(a,b,c); - else { - ltop = avma; ritopa((Obj)a,&pa); - ritopa((Obj)b,&pb); lbot = avma; - z = gerepile(ltop,lbot,gadd(pa,pb)); - patori(z,(Obj *)c); cgiv(z); - } +double mpfrtodbl(mpfr_t a) +{ + return mpfr_get_d(a,MPFR_RNDN); } -void subbf(a,b,c) -Num a,b; -Num *c; +Num tobf(Num a,int prec) { - GEN pa,pb,z; - long ltop,lbot; + mpfr_t r; + mpz_t z; + mpq_t q; + BF d; + N nm,dn; + int sgn; - if ( !a ) - (*chsgnnumt[NID(b)])(b,c); - else if ( !b ) - *c = a; - else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) - (*subnumt[MIN(NID(a),NID(b))])(a,b,c); - else { - ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; - z = gerepile(ltop,lbot,gsub(pa,pb)); - patori(z,(Obj *)c); cgiv(z); - } + 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,((Real)a)->body,MPFR_RNDN); + MPFRTOBF(r,d); + return (Num)d; + break; + case N_Q: + nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a); + if ( INT((Q)a) ) { + mpz_init(z); + mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); + if ( sgn < 0 ) mpz_neg(z,z); + mpfr_init_set_z(r,z,MPFR_RNDN); + } else { + mpq_init(q); + 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)); + if ( sgn < 0 ) mpq_neg(q,q); + mpfr_init_set_q(r,q,MPFR_RNDN); + } + MPFRTOBF(r,d); + return (Num)d; + break; + default: + error("tobf : invalid argument"); + break; + } + } } -void mulbf(a,b,c) -Num a,b; -Num *c; +void addbf(Num a,Num b,Num *c) { - GEN pa,pb,z; - long ltop,lbot; + mpfr_t r; + BF d; - if ( !a || !b ) - *c = 0; - else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) - (*mulnumt[MIN(NID(a),NID(b))])(a,b,c); - else { - ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; - z = gerepile(ltop,lbot,gmul(pa,pb)); - patori(z,(Obj *)c); cgiv(z); - } + if ( !a ) + *c = b; + else if ( !b ) + *c = a; + else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) + (*addnumt[MIN(NID(a),NID(b))])(a,b,c); + else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); + mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } } -void divbf(a,b,c) -Num a,b; -Num *c; +void subbf(Num a,Num b,Num *c) { - GEN pa,pb,z; - long ltop,lbot; + mpfr_t r; + BF d; - if ( !b ) - error("divbf : division by 0"); - else if ( !a ) - *c = 0; - else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) - (*divnumt[MIN(NID(a),NID(b))])(a,b,c); - else { - ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; - z = gerepile(ltop,lbot,gdiv(pa,pb)); - patori(z,(Obj *)c); cgiv(z); - } + if ( !a ) + (*chsgnnumt[NID(b)])(b,c); + else if ( !b ) + *c = a; + else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) + (*subnumt[MIN(NID(a),NID(b))])(a,b,c); + else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); + mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } } -void pwrbf(a,e,c) -Num a,e; -Num *c; +void mulbf(Num a,Num b,Num *c) { - GEN pa,pe,z; - long ltop,lbot; + mpfr_t r; + BF d; + int prec; - if ( !e ) - *c = (Num)ONE; - else if ( !a ) - *c = 0; - else { - ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma; - z = gerepile(ltop,lbot,gpui(pa,pe,prec)); - patori(z,(Obj *)c); cgiv(z); - } + if ( !a || !b ) + *c = 0; + else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) + (*mulnumt[MIN(NID(a),NID(b))])(a,b,c); + else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); + mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } } -void chsgnbf(a,c) -Num a,*c; +void divbf(Num a,Num b,Num *c) { - BF t; - GEN z; - int s; + mpfr_t r; + BF d; - if ( !a ) - *c = 0; - else if ( NID(a) <= N_A ) - (*chsgnnumt[NID(a)])(a,c); - else { - z = (GEN)((BF)a)->body; s = lg(z); NEWBF(t,s); - bcopy((char *)a,(char *)t,sizeof(struct oBF)+((s-1)*sizeof(long))); - z = (GEN)((BF)t)->body; setsigne(z,-signe(z)); - *c = (Num)t; - } + if ( !b ) + error("divbf : division by 0"); + else if ( !a ) + *c = 0; + else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) + (*divnumt[MIN(NID(a),NID(b))])(a,b,c); + else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); + mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } } -int cmpbf(a,b) -Num a,b; +void pwrbf(Num a,Num b,Num *c) { - GEN pa,pb; - int s; + mpfr_t r; + BF d; - if ( !a ) { - if ( !b || (NID(b)<=N_A) ) - return (*cmpnumt[NID(b)])(a,b); - else - return -signe(((BF)b)->body); - } else if ( !b ) { - if ( !a || (NID(a)<=N_A) ) - return (*cmpnumt[NID(a)])(a,b); - else - return signe(((BF)a)->body); - } else { - ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); - s = gcmp(pa,pb); cgiv(pb); cgiv(pa); - return s; - } + if ( !b ) + *c = (Num)ONE; + else if ( !a ) + *c = 0; + else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) + (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c); + else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); + mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } } -#endif /* PARI */ + +void chsgnbf(Num a,Num *c) +{ + mpfr_t r; + BF d; + + if ( !a ) + *c = 0; + else if ( NID(a) <= N_A ) + (*chsgnnumt[NID(a)])(a,c); + else { + mpfr_init2(r,BFPREC(a)); + mpfr_neg(r,((BF)a)->body,MPFR_RNDN); + MPFRTOBF(r,d); + *c = (Num)d; + } +} + +int cmpbf(Num a,Num b) +{ + if ( !a ) { + if ( !b || (NID(b)<=N_A) ) + return (*cmpnumt[NID(b)])(a,b); + else + return -mpfr_sgn(((BF)a)->body); + } else if ( !b ) { + if ( !a || (NID(a)<=N_A) ) + return (*cmpnumt[NID(a)])(a,b); + else + return mpfr_sgn(((BF)a)->body); + } else { + if ( NID(a) != N_B ) a = tobf(a,0); + if ( NID(b) != N_B ) b = tobf(b,0); + return mpfr_cmp(((BF)a)->body,((BF)b)->body); + } +}