=================================================================== RCS file: /home/cvs/OpenXM/src/ox_pari/ox_pari.c,v retrieving revision 1.4 retrieving revision 1.7 diff -u -p -r1.4 -r1.7 --- OpenXM/src/ox_pari/ox_pari.c 2015/08/17 06:14:37 1.4 +++ OpenXM/src/ox_pari/ox_pari.c 2015/08/18 05:04:35 1.7 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.3 2015/08/17 05:18:35 noro Exp $ */ +/* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.6 2015/08/18 02:24:04 noro Exp $ */ #include #include @@ -26,11 +26,19 @@ 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 @@ -163,6 +171,25 @@ int cmo2int(cmo *c) return 0; } +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; @@ -182,6 +209,18 @@ GEN cmo_zz_to_GEN(cmo_zz *c) 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; @@ -203,6 +242,77 @@ GEN cmo_bf_to_GEN(cmo_bf *c) 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; @@ -214,6 +324,22 @@ cmo_zz *GEN_to_cmo_zz(GEN z) 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; @@ -248,16 +374,33 @@ cmo_list *GEN_to_cmo_list(GEN z) 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; } @@ -272,6 +415,10 @@ cmo *GEN_to_cmo(GEN z) 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 */ @@ -280,6 +427,7 @@ cmo *GEN_to_cmo(GEN z) return (cmo *)make_error2(typ(z)); } } +/* type=1 : num/poly arg, type=2 : matrix arg */ struct parif { char *name; @@ -327,30 +475,6 @@ struct parif { {"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}, @@ -359,6 +483,34 @@ struct parif { {"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 @@ -377,14 +529,16 @@ struct parif *search_parif(char *name) int sm_executeFunction() { - long ltop,lbot; + pari_sp av0; int ac,i; cmo_int32 *c; cmo *av[PARI_MAX_AC]; cmo *ret; GEN z,m; struct parif *parif; + unsigned long prec; + if ( setjmp(GP_DATA->env) ) { printf("sm_executeFunction : an error occured.\n");fflush(stdout); push((cmo*)make_error2(0)); @@ -423,20 +577,29 @@ int sm_executeFunction() ret = (cmo *)new_cmo_int32(a); push(ret); return 0; - } else if ( parif->type == 1 ) { - /* one variable possibly with prec */ - unsigned long prec; - - ltop = avma; + } else if ( parif->type == 1 || parif->type == 2 ) { + /* 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; + prec = ac==2 ? cmo_to_int(av[1])*3.32193/32+3 : precreal; + if ( parif->type == 2 ) { + /* matrix argument */ + 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); + settyp(z,t_MAT); + z = shallowtrans(z); + } + printf("input : "); output(z); m = (*parif->f)(z,prec); - lbot = avma; ret = GEN_to_cmo(m); -// gerepile(ltop,lbot,0); + avma = av0; push(ret); return 0; } else {