=================================================================== RCS file: /home/cvs/OpenXM/src/ox_pari/ox_pari.c,v retrieving revision 1.11 retrieving revision 1.12 diff -u -p -r1.11 -r1.12 --- OpenXM/src/ox_pari/ox_pari.c 2015/08/21 00:53:53 1.11 +++ OpenXM/src/ox_pari/ox_pari.c 2015/08/21 01:20:30 1.12 @@ -1,49 +1,15 @@ -/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.10 2015/08/20 08:56:30 noro Exp $ */ +/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.11 2015/08/21 00:53:53 noro Exp $ */ -#include -#include -#include -#include "pari/pari.h" -#include "pari/paripriv.h" -#include "gmp.h" -#include "gmp-impl.h" -#include "mpfr.h" -#include "ox_toolkit.h" +#include "ox_pari.h" + OXFILE *fd_rw; -#define MPFR_PREC(x) ((x)->_mpfr_prec) -#define MPFR_EXP(x) ((x)->_mpfr_exp) -#define MPFR_MANT(x) ((x)->_mpfr_d) -#define MPFR_LAST_LIMB(x) ((MPFR_PREC (x) - 1) / GMP_NUMB_BITS) -#define MPFR_LIMB_SIZE(x) (MPFR_LAST_LIMB (x) + 1) - static int stack_size = 0; static int stack_pointer = 0; static cmo **stack = NULL; extern int debug_print; long paristack=10000000; -void init_pari(void); -cmo *GEN_to_cmo(GEN z); -cmo_zz *GEN_to_cmo_zz(GEN z); -cmo_qq *GEN_to_cmo_qq(GEN z); -cmo_bf *GEN_to_cmo_bf(GEN z); -cmo_list *GEN_to_cmo_list(GEN z); -cmo_complex *GEN_to_cmo_cmo_complex(GEN z); -cmo_polynomial_in_one_variable *GEN_to_cmo_up(GEN z); -cmo_recursive_polynomial *GEN_to_cmo_rp(GEN z); - -GEN cmo_to_GEN(cmo *c); -GEN cmo_int32_to_GEN(cmo_int32 *c); -GEN cmo_zz_to_GEN(cmo_zz *c); -GEN cmo_qq_to_GEN(cmo_qq *c); -GEN cmo_bf_to_GEN(cmo_bf *c); -GEN cmo_list_to_GEN(cmo_list *c); -GEN cmo_rp_to_GEN(cmo_recursive_polynomial *c); -GEN cmo_up_to_GEN(cmo_polynomial_in_one_variable *c); -GEN cmo_complex_to_GEN(cmo_complex *c); - - #define INIT_S_SIZE 2048 #define EXT_S_SIZE 2048 @@ -162,394 +128,6 @@ char *get_str() return ""; } -GEN cmo_int32_to_GEN(cmo_int32 *c) -{ - GEN z; - int i,sgn; - - i = c->i; - if ( !i ) return gen_0; - z = cgeti(3); - sgn = 1; - if ( i < 0 ) { - i = -i; - sgn = -1; - } - z[2] = i; - setsigne(z,sgn); - setlgefint(z,lg(z)); - return z; -} - -GEN cmo_zz_to_GEN(cmo_zz *c) -{ - mpz_ptr mpz; - GEN z; - long *ptr; - int j,sgn,len; - - mpz = c->mpz; - sgn = mpz_sgn(mpz); - len = ABSIZ(mpz); - ptr = (long *)PTR(mpz); - z = cgeti(len+2); - for ( j = 0; j < len; j++ ) - z[len-j+1] = ptr[j]; - setsigne(z,sgn); - setlgefint(z,lg(z)); - return z; -} - -GEN cmo_qq_to_GEN(cmo_qq *c) -{ - GEN z,nm,den; - - z = cgetg(3,t_FRAC); - nm = cmo_zz_to_GEN(new_cmo_zz_set_mpz(mpq_numref(c->mpq))); - den = cmo_zz_to_GEN(new_cmo_zz_set_mpz(mpq_denref(c->mpq))); - z[1] = (long)nm; - z[2] = (long)den; - return z; -} - -GEN cmo_bf_to_GEN(cmo_bf *c) -{ - mpfr_ptr mpfr; - GEN z; - int sgn,len,j; - long exp; - long *ptr; - - mpfr = c->mpfr; - sgn = MPFR_SIGN(mpfr); - exp = MPFR_EXP(mpfr)-1; - len = MPFR_LIMB_SIZE(mpfr); - ptr = (long *)MPFR_MANT(mpfr); - z = cgetr(len+2); - for ( j = 0; j < len; j++ ) - z[len-j+1] = ptr[j]; - z[1] = evalsigne(sgn)|evalexpo(exp); - setsigne(z,sgn); - return z; -} - -/* list->vector */ - -GEN cmo_list_to_GEN(cmo_list *c) -{ - GEN z; - int i; - cell *cell; - - z = cgetg(c->length+1,t_VEC); - for ( i = 0, cell = c->head->next; cell != c->head; cell = cell->next, i++ ) { - z[i+1] = (long)cmo_to_GEN(cell->cmo); - } - return z; -} - -GEN cmo_complex_to_GEN(cmo_complex *c) -{ - GEN z; - - z = cgetg(3,t_COMPLEX); - z[1] = (long)cmo_to_GEN(c->re); - z[2] = (long)cmo_to_GEN(c->im); - return z; -} - -GEN cmo_up_to_GEN(cmo_polynomial_in_one_variable *c) -{ - GEN z; - int d,i; - cell *cell; - - d = c->head->next->exp; - z = cgetg(d+3,t_POL); - setsigne(z,1); - setvarn(z,c->var); - setlgef(z,d+3); - for ( i = 2; i <= d+2; i++ ) - z[i] = (long)gen_0; - for ( cell = c->head->next; cell != c->head; cell = cell->next ) { - z[2+cell->exp] = (long)cmo_to_GEN(cell->cmo); - } - return z; -} - -cmo_list *current_ringdef; - -void register_variables(cmo_list *ringdef) -{ - current_ringdef = ringdef; -} - -GEN cmo_rp_to_GEN(cmo_recursive_polynomial *c) -{ - register_variables(c->ringdef); - switch ( c->coef->tag ) { - case CMO_ZERO: - case CMO_NULL: - return gen_0; - case CMO_INT32: - return cmo_int32_to_GEN((cmo_int32 *)c->coef); - case CMO_ZZ: - return cmo_zz_to_GEN((cmo_zz *)c->coef); - case CMO_QQ: - return cmo_qq_to_GEN((cmo_qq *)c->coef); - case CMO_POLYNOMIAL_IN_ONE_VARIABLE: - return cmo_up_to_GEN((cmo_polynomial_in_one_variable *)c->coef); - default: - return 0; - } -} - -cmo_zz *GEN_to_cmo_zz(GEN z) -{ - cmo_zz *c; - - c = new_cmo_zz(); - mpz_import(c->mpz,lgef(z)-2,1,sizeof(long),0,0,&z[2]); - if ( signe(z) < 0 ) - mpz_neg(c->mpz,c->mpz); - return c; -} - -cmo_qq *GEN_to_cmo_qq(GEN z) -{ - cmo_qq *c; - GEN num,den; - - num = (GEN)z[1]; - den = (GEN)z[2]; - c = new_cmo_qq(); - mpz_import(mpq_numref(c->mpq),lgef(num)-2,1,sizeof(long),0,0,&num[2]); - mpz_import(mpq_denref(c->mpq),lgef(num)-2,1,sizeof(long),0,0,&den[2]); - if ( signe(num)*signe(den) < 0 ) - mpz_neg(mpq_numref(c->mpq),mpq_numref(c->mpq)); - return c; -} - - -cmo_bf *GEN_to_cmo_bf(GEN z) -{ - cmo_bf *c; - int len,prec,j; - long *ptr; - - c = new_cmo_bf(); - len = lg(z)-2; - prec = len*sizeof(long)*8; - mpfr_init2(c->mpfr,prec); - ptr = (long *)MPFR_MANT(c->mpfr); - for ( j = 0; j < len; j++ ) - ptr[j] = z[len-j+1]; - MPFR_EXP(c->mpfr) = (long long)(expo(z)+1); - MPFR_SIGN(c->mpfr) = gsigne(z); - return c; -} - - -cmo_list *GEN_to_cmo_list(GEN z) -{ - cmo_list *c; - cmo *ob; - int i,len; - - c = new_cmo_list(); - len = lg(z)-1; - for ( i = 1; i <= len; i++ ) { - ob = GEN_to_cmo((GEN)z[i]); - c = list_append(c,ob); - } - return c; -} - -cmo_complex *GEN_to_cmo_complex(GEN z) -{ - cmo_complex *c; - - c = new_cmo_complex(); - c->re = GEN_to_cmo((GEN)z[1]); - c->im = GEN_to_cmo((GEN)z[2]); - return c; -} - -cmo_polynomial_in_one_variable *GEN_to_cmo_up(GEN z) -{ - cmo_polynomial_in_one_variable *c; - int i; - cmo *coef; - - c = new_cmo_polynomial_in_one_variable(varn(z)); - for ( i = lg(z)-1; i >= 2; i-- ) - if ( (GEN)z[i] != gen_0 ) { - coef = GEN_to_cmo((GEN)z[i]); - list_append_monomial((cmo_list *)c, coef, i-2); - } - return c; -} - -cmo_recursive_polynomial *GEN_to_cmo_rp(GEN z) -{ - cmo_recursive_polynomial *c; - - if ( !signe(z) ) return (cmo_recursive_polynomial *)new_cmo_zero(); - c = new_cmo_recursive_polynomial(current_ringdef,(cmo *)GEN_to_cmo_up(z)); - return c; -} - -GEN cmo_to_GEN(cmo *c) -{ - switch ( c->tag ) { - case CMO_ZERO: - case CMO_NULL: - return gen_0; - case CMO_ZZ: /* int */ - return cmo_zz_to_GEN((cmo_zz *)c); - case CMO_COMPLEX: /* complex */ - return cmo_complex_to_GEN((cmo_complex *)c); - case CMO_IEEE_DOUBLE_FLOAT: - return dbltor(((cmo_double *)c)->d); - case CMO_BIGFLOAT: /* bigfloat */ - return cmo_bf_to_GEN((cmo_bf *)c); - case CMO_LIST: - return cmo_list_to_GEN((cmo_list *)c); - case CMO_RECURSIVE_POLYNOMIAL: - return cmo_rp_to_GEN((cmo_recursive_polynomial *)c); - case CMO_POLYNOMIAL_IN_ONE_VARIABLE: - return cmo_up_to_GEN((cmo_polynomial_in_one_variable *)c); - default: - return 0; - } -} - -cmo *GEN_to_cmo(GEN z) -{ - char buf[BUFSIZ]; - - if ( gcmp0(z) ) - return new_cmo_zero(); - switch ( typ(z) ) { - case t_INT: /* int */ - return (cmo *)GEN_to_cmo_zz(z); - case t_REAL: /* bigfloat */ - return (cmo *)GEN_to_cmo_bf(z); - case t_FRAC: /* rational number */ - return (cmo *)GEN_to_cmo_qq(z); - case t_COMPLEX: /* complex */ - return (cmo *)GEN_to_cmo_complex(z); - case t_POL: - return (cmo *)GEN_to_cmo_rp(z); - case t_VEC: case t_COL: /* vector */ - return (cmo *)GEN_to_cmo_list(z); - case t_MAT: /* matrix */ - return (cmo *)GEN_to_cmo_list(shallowtrans(z)); - default: - sprintf(buf,"GEN_to_cmo : unsupported type=%d",(int)typ(z)); - return (cmo *)make_error2(buf); - } -} -/* type=1 : num/poly arg, type=2 : matrix arg */ - -struct parif { - char *name; - GEN (*f)(); - int type; -} parif_tab[] = { -/* (ulong)allocatemoremem(ulong) */ - {"allocatemem",(GEN (*)())allocatemoremem,0}, -/* num/num */ - {"abs",gabs,1}, - {"erfc",gerfc,1}, - {"arg",garg,1}, - {"isqrt",racine,1}, - {"gamma",ggamma,1}, - {"zeta",gzeta,1}, - {"floor",gfloor,1}, - {"frac",gfrac,1}, - {"imag",gimag,1}, - {"conj",gconj,1}, - {"ceil",gceil,1}, - {"isprime",gisprime,2}, - {"bigomega",gbigomega,1}, - {"denom",denom,1}, - {"numer",numer,1}, - {"lngamma",glngamma,1}, - {"logagm",glogagm,1}, - {"classno",classno,1}, - {"classno2",classno2,1}, - {"dilog",dilog,1}, - {"disc",discsr,1}, - {"discf",discf,1}, - {"nextprime",nextprime,1}, - {"eintg1",eint1,1}, - {"eta",eta,1}, - {"issqfree",gissquarefree,1}, - {"issquare",gcarreparfait,1}, - {"gamh",ggamd,1}, - {"hclassno",classno3,1}, - - /* num/array */ - {"binary",binaire,1}, - {"factorint",factorint,2}, - {"factor",Z_factor,1}, - {"cf",gcf,1}, - {"divisors",divisors,1}, - {"smallfact",smallfact,1}, - - /* poly/poly */ - {"centerlift",centerlift,1}, - {"content",content,1}, - - /* poly/array */ - {"galois",galois,1}, - {"roots",roots,1}, - - /* mat/mat */ - {"adj",adj,2}, - {"lll",lll,2}, - {"lllgen",lllgen,2}, - {"lllgram",lllgram,2}, - {"lllgramgen",lllgramgen,2}, - {"lllgramint",lllgramint,2}, - {"lllgramkerim",lllgramkerim,2}, - {"lllgramkerimgen",lllgramkerimgen,2}, - {"lllint",lllint,2}, - {"lllkerim",lllkerim,2}, - {"lllkerimgen",lllkerimgen,2}, - {"trans",gtrans,2}, - {"eigen",eigen,2}, - {"hermite",hnf,2}, - {"mat",gtomat,2}, - {"matrixqz2",matrixqz2,2}, - {"matrixqz3",matrixqz3,2}, - {"hess",hess,2}, - {"ker",ker,2}, - {"keri",keri,2}, - {"kerint",kerint,2}, - {"kerintg1",kerint1,2}, - - /* mat/poly */ - {"det",det,2}, - {"det2",det2,2}, - -}; - -#define PARI_MAX_AC 64 - -struct parif *search_parif(char *name) -{ - int tablen,i; - - tablen = sizeof(parif_tab)/sizeof(struct parif); - for ( i = 0; i < tablen; i++ ) { - if ( !strcmp(parif_tab[i].name,name) ) - return &parif_tab[i]; - } - return 0; -} - int ismatrix(GEN z) { int len,col,i; @@ -613,7 +191,7 @@ int sm_executeFunction() ret = (cmo *)new_cmo_int32(a); push(ret); return 0; - } else if ( parif->type == 1 || parif->type == 2 ) { + } else if ( parif->type == 1 ) { /* one number/poly/matrix argument possibly with prec */ av0 = avma; z = cmo_to_GEN(av[0]);