version 1.1.1.1, 1999/12/03 07:39:08 |
version 1.7, 2015/08/04 06:20:45 |
|
|
/* $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" |
#include "ca.h" |
#if PARI |
|
#include "base.h" |
#include "base.h" |
#include <math.h> |
#include <math.h> |
#include "genpari.h" |
|
|
|
extern long prec; |
Num tobf(Num,int); |
|
|
void ritopa(Obj,GEN *); |
#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) |
void patori(GEN,Obj *); |
|
|
|
void addbf(a,b,c) |
void strtobf(char *s,BF *p) |
Num a,b; |
|
Num *c; |
|
{ |
{ |
GEN pa,pb,z; |
BF r; |
long ltop,lbot; |
NEWBF(r); |
|
mpfr_init(r->body); |
|
mpfr_set_str(r->body,s,10,MPFR_RNDN); |
|
*p = r; |
|
} |
|
|
if ( !a ) |
double mpfrtodbl(mpfr_t a) |
*c = b; |
{ |
else if ( !b ) |
return mpfr_get_d(a,MPFR_RNDN); |
*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); |
|
} |
|
} |
} |
|
|
void subbf(a,b,c) |
Num tobf(Num a,int prec) |
Num a,b; |
|
Num *c; |
|
{ |
{ |
GEN pa,pb,z; |
mpfr_t r; |
long ltop,lbot; |
mpz_t z; |
|
mpq_t q; |
|
BF d; |
|
N nm,dn; |
|
int sgn; |
|
|
if ( !a ) |
if ( !a ) { |
(*chsgnnumt[NID(b)])(b,c); |
prec ? mpfr_init2(r,prec) : mpfr_init(r); |
else if ( !b ) |
mpfr_set_zero(r,1); |
*c = a; |
MPFRTOBF(r,d); |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
return (Num)d; |
(*subnumt[MIN(NID(a),NID(b))])(a,b,c); |
} else { |
else { |
switch ( NID(a) ) { |
ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; |
case N_B: |
z = gerepile(ltop,lbot,gsub(pa,pb)); |
return a; |
patori(z,(Obj *)c); cgiv(z); |
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) |
void addbf(Num a,Num b,Num *c) |
Num a,b; |
|
Num *c; |
|
{ |
{ |
GEN pa,pb,z; |
mpfr_t r; |
long ltop,lbot; |
BF d; |
|
|
if ( !a || !b ) |
if ( !a ) |
*c = 0; |
*c = b; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( !b ) |
(*mulnumt[MIN(NID(a),NID(b))])(a,b,c); |
*c = a; |
else { |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; |
(*addnumt[MIN(NID(a),NID(b))])(a,b,c); |
z = gerepile(ltop,lbot,gmul(pa,pb)); |
else { |
patori(z,(Obj *)c); cgiv(z); |
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) |
void subbf(Num a,Num b,Num *c) |
Num a,b; |
|
Num *c; |
|
{ |
{ |
GEN pa,pb,z; |
mpfr_t r; |
long ltop,lbot; |
BF d; |
|
|
if ( !b ) |
if ( !a ) |
error("divbf : division by 0"); |
(*chsgnnumt[NID(b)])(b,c); |
else if ( !a ) |
else if ( !b ) |
*c = 0; |
*c = a; |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
(*divnumt[MIN(NID(a),NID(b))])(a,b,c); |
(*subnumt[MIN(NID(a),NID(b))])(a,b,c); |
else { |
else { |
ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma; |
if ( NID(a) != N_B ) a = tobf(a,0); |
z = gerepile(ltop,lbot,gdiv(pa,pb)); |
if ( NID(b) != N_B ) b = tobf(b,0); |
patori(z,(Obj *)c); cgiv(z); |
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) |
void mulbf(Num a,Num b,Num *c) |
Num a,e; |
|
Num *c; |
|
{ |
{ |
GEN pa,pe,z; |
mpfr_t r; |
long ltop,lbot; |
BF d; |
|
int prec; |
|
|
if ( !e ) |
if ( !a || !b ) |
*c = (Num)ONE; |
*c = 0; |
else if ( !a ) |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
*c = 0; |
(*mulnumt[MIN(NID(a),NID(b))])(a,b,c); |
else { |
else { |
ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma; |
if ( NID(a) != N_B ) a = tobf(a,0); |
z = gerepile(ltop,lbot,gpui(pa,pe,prec)); |
if ( NID(b) != N_B ) b = tobf(b,0); |
patori(z,(Obj *)c); cgiv(z); |
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) |
void divbf(Num a,Num b,Num *c) |
Num a,*c; |
|
{ |
{ |
BF t; |
mpfr_t r; |
GEN z; |
BF d; |
int s; |
|
|
|
if ( !a ) |
if ( !b ) |
*c = 0; |
error("divbf : division by 0"); |
else if ( NID(a) <= N_A ) |
else if ( !a ) |
(*chsgnnumt[NID(a)])(a,c); |
*c = 0; |
else { |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
z = (GEN)((BF)a)->body; s = lg(z); NEWBF(t,s); |
(*divnumt[MIN(NID(a),NID(b))])(a,b,c); |
bcopy((char *)a,(char *)t,sizeof(struct oBF)+((s-1)*sizeof(long))); |
else { |
z = (GEN)((BF)t)->body; setsigne(z,-signe(z)); |
if ( NID(a) != N_B ) a = tobf(a,0); |
*c = (Num)t; |
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) |
void pwrbf(Num a,Num b,Num *c) |
Num a,b; |
|
{ |
{ |
GEN pa,pb; |
mpfr_t r; |
int s; |
BF d; |
|
|
if ( !a ) { |
if ( !b ) |
if ( !b || (NID(b)<=N_A) ) |
*c = (Num)ONE; |
return (*cmpnumt[NID(b)])(a,b); |
else if ( !a ) |
else |
*c = 0; |
return -signe(((BF)b)->body); |
else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) ) |
} else if ( !b ) { |
(*pwrnumt[MIN(NID(a),NID(b))])(a,b,c); |
if ( !a || (NID(a)<=N_A) ) |
else { |
return (*cmpnumt[NID(a)])(a,b); |
if ( NID(a) != N_B ) a = tobf(a,0); |
else |
if ( NID(b) != N_B ) b = tobf(b,0); |
return signe(((BF)a)->body); |
mpfr_init2(r,MAX(BFPREC(a),BFPREC(b))); |
} else { |
mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN); |
ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); |
MPFRTOBF(r,d); |
s = gcmp(pa,pb); cgiv(pb); cgiv(pa); |
*c = (Num)d; |
return s; |
} |
} |
|
} |
} |
#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); |
|
} |
|
} |