=================================================================== RCS file: /home/cvs/OpenXM/src/ox_pari/ox_pari.c,v retrieving revision 1.6 retrieving revision 1.14 diff -u -p -r1.6 -r1.14 --- OpenXM/src/ox_pari/ox_pari.c 2015/08/18 02:24:04 1.6 +++ OpenXM/src/ox_pari/ox_pari.c 2016/08/23 03:03:26 1.14 @@ -1,44 +1,15 @@ -/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.5 2015/08/17 07:19:16 noro Exp $ */ +/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.13 2016/08/01 01:35:01 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_bf *GEN_to_cmo_bf(GEN z); -cmo_list *GEN_to_cmo_list(GEN z); -cmo_complex *GEN_to_cmo_cmo_complex(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 @@ -114,7 +85,8 @@ void pops(int n) int sm_mathcap() { - mathcap_init(OX_PARI_VERSION, ID_STRING, "ox_pari", NULL, NULL); + char *opts[] = {"no_ox_reset", NULL}; + mathcap_init2(OX_PARI_VERSION, ID_STRING, "ox_pari", NULL, NULL, opts); push((cmo*)oxf_cmo_mathcap(fd_rw)); return 0; } @@ -130,9 +102,9 @@ int sm_popCMO() return SM_popCMO; } -cmo_error2 *make_error2(int code) +cmo_error2 *make_error2(char *message) { - return (cmo_error2 *) new_cmo_int32(code); + return new_cmo_error2((cmo *)new_cmo_string(message)); } int get_i() @@ -143,7 +115,7 @@ int get_i() }else if(c->tag == CMO_ZZ) { return mpz_get_si(((cmo_zz *)c)->mpz); } - make_error2(-1); + make_error2("get_i : invalid object"); return 0; } @@ -153,353 +125,22 @@ char *get_str() if(c->tag == CMO_STRING) { return ((cmo_string *)c)->s; } - make_error2(-1); + make_error2("get_str : invalid object"); return ""; } -int cmo2int(cmo *c) +int ismatrix(GEN z) { - if(c->tag == CMO_INT32) { - return ((cmo_int32 *)c)->i; - }else if(c->tag == CMO_ZZ) { - return mpz_get_si(((cmo_zz *)c)->mpz); - } else if(c->tag == CMO_NULL){ - return 0; - } + int len,col,i; - return 0; + if ( typ(z) != t_VEC ) return 0; + if ( typ((GEN)z[1]) != t_VEC ) return 0; + len = lg(z); col = lg((GEN)z[1]); + for ( i = 2; i < len; i++ ) + if ( lg((GEN)z[i]) != col ) return 0; + return 1; } -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,4); - 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; - cell *cell; - - z = cgetg(c->length+1,17); - for ( cell = c->head->next; cell != c->head; cell = cell->next ) { - z[cell->exp] = (long)cmo_to_GEN(cell->cmo); - } - return z; -} - -GEN cmo_complex_to_GEN(cmo_complex *c) -{ - GEN z; - - z = cgetg(3,6); - 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,10); - 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: - 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_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; -} - - -GEN cmo_to_GEN(cmo *c) -{ - switch ( c->tag ) { - case CMO_ZERO: - return gen_0; - case CMO_ZZ: /* int */ - return cmo_zz_to_GEN((cmo_zz *)c); - 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) -{ - if ( gcmp0(z) ) - return new_cmo_zero(); - switch ( typ(z) ) { - case 1: /* int */ - return (cmo *)GEN_to_cmo_zz(z); - case 2: /* bigfloat */ - return (cmo *)GEN_to_cmo_bf(z); - case 6: /* complex */ - return (cmo *)GEN_to_cmo_complex(z); - case 17: case 18: /* vector */ - return (cmo *)GEN_to_cmo_list(z); - case 19: /* matrix */ - return (cmo *)GEN_to_cmo_list(shallowtrans(z)); - default: - return (cmo *)make_error2(typ(z)); - } -} - -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}, - - /* mat/mat */ - {"adj",adj,1}, - {"lll",lll,1}, - {"lllgen",lllgen,1}, - {"lllgram",lllgram,1}, - {"lllgramgen",lllgramgen,1}, - {"lllgramint",lllgramint,1}, - {"lllgramkerim",lllgramkerim,1}, - {"lllgramkerimgen",lllgramkerimgen,1}, - {"lllint",lllint,1}, - {"lllkerim",lllkerim,1}, - {"lllkerimgen",lllkerimgen,1}, - {"trans",gtrans,1}, - {"eigen",eigen,1}, - {"hermite",hnf,1}, - {"mat",gtomat,1}, - {"matrixqz2",matrixqz2,1}, - {"matrixqz3",matrixqz3,1}, - {"hess",hess,1}, - - /* mat/poly */ - {"det",det,1}, - {"det2",det2,1}, - - /* poly/poly */ - {"centerlift",centerlift,1}, - {"content",content,1}, - - /* poly/array */ - {"galois",galois,1}, - {"roots",roots,1}, - -}; - -#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 sm_executeFunction() { pari_sp av0; @@ -509,37 +150,40 @@ int sm_executeFunction() cmo *ret; GEN z,m; struct parif *parif; + unsigned long prec; + char buf[BUFSIZ]; if ( setjmp(GP_DATA->env) ) { - printf("sm_executeFunction : an error occured.\n");fflush(stdout); - push((cmo*)make_error2(0)); + sprintf(buf,"sm_executeFunction : an error occured in PARI."); + push((cmo*)make_error2(buf)); return -1; } cmo_string *func = (cmo_string *)pop(); if(func->tag != CMO_STRING) { - printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout); - push((cmo*)make_error2(0)); + sprintf(buf,"sm_executeFunction : func->tag=%d is not CMO_STRING",func->tag); + push((cmo*)make_error2(buf)); return -1; } c = (cmo_int32 *)pop(); ac = c->i; if ( ac > PARI_MAX_AC ) { - push((cmo*)make_error2(0)); + push((cmo*)make_error2("sm_executeFunction : too many arguments")); return -1; } for ( i = 0; i < ac; i++ ) { av[i] = (cmo *)pop(); - fprintf(stderr,"arg%d:",i); - print_cmo(av[i]); - fprintf(stderr,"\n"); +// fprintf(stderr,"arg%d:",i); +// print_cmo(av[i]); +// fprintf(stderr,"\n"); } if( strcmp( func->s, "exit" ) == 0 ) exit(0); parif =search_parif(func->s); if ( !parif ) { - push((cmo*)make_error2(0)); + sprintf(buf,"%s : not implemented",func->s); + push((cmo*)make_error2(buf)); return -1; } else if ( parif->type == 0 ) { /* one long int variable */ @@ -549,24 +193,45 @@ int sm_executeFunction() push(ret); return 0; } else if ( parif->type == 1 ) { - /* one variable possibly with prec */ - unsigned long prec; - + /* one number/poly/matrix argument possibly with prec */ av0 = avma; z = cmo_to_GEN(av[0]); - if ( ac == 2 ) { - prec = cmo_to_int(av[1])*3.32193/32+3; - } else - prec = precreal; - printf("input : "); - output(z); + prec = ac==2 ? cmo_to_int(av[1])*3.32193/32+3 : precreal; + if ( ismatrix(z) ) { + int i,len; + len = lg(z); + for ( i = 1; i < len; i++ ) + settyp(z[i],t_COL); + settyp(z,t_MAT); + z = shallowtrans(z); + } + printf("input : "); output(z); m = (*parif->f)(z,prec); ret = GEN_to_cmo(m); avma = av0; push(ret); return 0; + } else if ( parif->type == 2 ) { + /* one number/poly/matrix argument with flag=0 */ + av0 = avma; + z = cmo_to_GEN(av[0]); + if ( ismatrix(z) ) { + int i,len; + len = lg(z); + for ( i = 1; i < len; i++ ) + settyp(z[i],t_COL); + settyp(z,t_MAT); + z = shallowtrans(z); + } + printf("input : "); output(z); + m = (*parif->f)(z,0); + ret = GEN_to_cmo(m); + avma = av0; + push(ret); + return 0; } else { - push((cmo*)make_error2(0)); + sprintf(buf,"%s : not implemented",func->s); + push((cmo*)make_error2(buf)); return -1; } } @@ -586,6 +251,9 @@ int receive_and_execute_sm_command() break; case SM_setMathCap: pop(); + break; + case SM_shutdown: + exit(0); break; default: printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);