=================================================================== RCS file: /home/cvs/OpenXM/src/ox_pari/ox_pari.c,v retrieving revision 1.7 retrieving revision 1.12 diff -u -p -r1.7 -r1.12 --- OpenXM/src/ox_pari/ox_pari.c 2015/08/18 05:04:35 1.7 +++ OpenXM/src/ox_pari/ox_pari.c 2015/08/21 01:20:30 1.12 @@ -1,45 +1,15 @@ -/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.6 2015/08/18 02:24:04 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); -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 @@ -131,9 +101,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() @@ -144,7 +114,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; } @@ -154,379 +124,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; - int i; - cell *cell; - - z = cgetg(c->length+1,17); - 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,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: - 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; -} - - -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_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 4: /* rational number */ - return (cmo *)GEN_to_cmo_qq(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)); - } -} -/* 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 sm_executeFunction() { pari_sp av0; @@ -537,38 +150,39 @@ int sm_executeFunction() 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 */ @@ -577,19 +191,13 @@ 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]); prec = ac==2 ? cmo_to_int(av[1])*3.32193/32+3 : precreal; - if ( parif->type == 2 ) { - /* matrix argument */ + if ( ismatrix(z) ) { int i,len; - - if ( typ(z) != t_VEC ) { - push((cmo*)make_error2(0)); - return -1; - } len = lg(z); for ( i = 1; i < len; i++ ) settyp(z[i],t_COL); @@ -603,7 +211,8 @@ int sm_executeFunction() push(ret); return 0; } else { - push((cmo*)make_error2(0)); + sprintf(buf,"%s : not implemented",func->s); + push((cmo*)make_error2(buf)); return -1; } } @@ -623,6 +232,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);