Annotation of OpenXM/src/ox_pari/ox_pari.c, Revision 1.6
1.6 ! noro 1: /* $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.5 2015/08/17 07:19:16 noro Exp $ */
1.1 noro 2:
3: #include <stdio.h>
4: #include <stdlib.h>
5: #include <string.h>
1.2 noro 6: #include "pari/pari.h"
1.3 noro 7: #include "pari/paripriv.h"
1.1 noro 8: #include "gmp.h"
1.2 noro 9: #include "gmp-impl.h"
1.3 noro 10: #include "mpfr.h"
1.1 noro 11: #include "ox_toolkit.h"
12: OXFILE *fd_rw;
13:
1.3 noro 14: #define MPFR_PREC(x) ((x)->_mpfr_prec)
15: #define MPFR_EXP(x) ((x)->_mpfr_exp)
16: #define MPFR_MANT(x) ((x)->_mpfr_d)
17: #define MPFR_LAST_LIMB(x) ((MPFR_PREC (x) - 1) / GMP_NUMB_BITS)
18: #define MPFR_LIMB_SIZE(x) (MPFR_LAST_LIMB (x) + 1)
19:
1.1 noro 20: static int stack_size = 0;
21: static int stack_pointer = 0;
22: static cmo **stack = NULL;
23: extern int debug_print;
1.2 noro 24: long paristack=10000000;
1.1 noro 25:
26: void init_pari(void);
1.2 noro 27: cmo *GEN_to_cmo(GEN z);
28: cmo_zz *GEN_to_cmo_zz(GEN z);
1.3 noro 29: cmo_bf *GEN_to_cmo_bf(GEN z);
1.2 noro 30: cmo_list *GEN_to_cmo_list(GEN z);
1.6 ! noro 31: cmo_complex *GEN_to_cmo_cmo_complex(GEN z);
1.2 noro 32: GEN cmo_to_GEN(cmo *c);
1.6 ! noro 33: GEN cmo_int32_to_GEN(cmo_int32 *c);
1.2 noro 34: GEN cmo_zz_to_GEN(cmo_zz *c);
1.6 ! noro 35: GEN cmo_qq_to_GEN(cmo_qq *c);
1.3 noro 36: GEN cmo_bf_to_GEN(cmo_bf *c);
1.6 ! noro 37: GEN cmo_list_to_GEN(cmo_list *c);
! 38: GEN cmo_rp_to_GEN(cmo_recursive_polynomial *c);
! 39: GEN cmo_up_to_GEN(cmo_polynomial_in_one_variable *c);
! 40: GEN cmo_complex_to_GEN(cmo_complex *c);
1.1 noro 41:
42: #define INIT_S_SIZE 2048
43: #define EXT_S_SIZE 2048
44:
1.3 noro 45: void *gc_realloc(void *p,size_t osize,size_t nsize)
46: {
47: return (void *)GC_realloc(p,nsize);
48: }
49:
50: void gc_free(void *p,size_t size)
51: {
52: GC_free(p);
53: }
54:
55: void init_gc()
56: {
1.4 noro 57: GC_INIT();
1.3 noro 58: mp_set_memory_functions(GC_malloc,gc_realloc,gc_free);
59: }
60:
1.1 noro 61: void init_pari()
62: {
1.2 noro 63: pari_init(paristack,2);
1.1 noro 64: }
65:
66: int initialize_stack()
67: {
1.4 noro 68: stack_pointer = 0;
69: stack_size = INIT_S_SIZE;
70: stack = MALLOC(stack_size*sizeof(cmo*));
71: return 0;
1.1 noro 72: }
73:
74: static int extend_stack()
75: {
1.4 noro 76: int size2 = stack_size + EXT_S_SIZE;
77: cmo **stack2 = MALLOC(size2*sizeof(cmo*));
78: memcpy(stack2, stack, stack_size*sizeof(cmo *));
79: free(stack);
80: stack = stack2;
81: stack_size = size2;
82: return 0;
1.1 noro 83: }
84:
85: int push(cmo* m)
86: {
1.4 noro 87: stack[stack_pointer] = m;
88: stack_pointer++;
89: if(stack_pointer >= stack_size) {
90: extend_stack();
91: }
92: return 0;
1.1 noro 93: }
94:
95: cmo* pop()
96: {
1.4 noro 97: if(stack_pointer > 0) {
98: stack_pointer--;
99: return stack[stack_pointer];
100: }
101: return new_cmo_null();
1.1 noro 102: }
103:
104: void pops(int n)
105: {
1.4 noro 106: stack_pointer -= n;
107: if(stack_pointer < 0) {
108: stack_pointer = 0;
109: }
1.1 noro 110: }
111:
112: #define OX_PARI_VERSION 20150731
113: #define ID_STRING "2015/07/31 15:00:00"
114:
115: int sm_mathcap()
116: {
1.4 noro 117: mathcap_init(OX_PARI_VERSION, ID_STRING, "ox_pari", NULL, NULL);
118: push((cmo*)oxf_cmo_mathcap(fd_rw));
119: return 0;
1.1 noro 120: }
121:
122: int sm_popCMO()
123: {
1.4 noro 124: cmo* m = pop();
1.1 noro 125:
1.4 noro 126: if(m != NULL) {
127: send_ox_cmo(fd_rw, m);
128: return 0;
129: }
130: return SM_popCMO;
1.1 noro 131: }
132:
133: cmo_error2 *make_error2(int code)
134: {
1.4 noro 135: return (cmo_error2 *) new_cmo_int32(code);
1.1 noro 136: }
137:
138: int get_i()
139: {
1.4 noro 140: cmo *c = pop();
141: if(c->tag == CMO_INT32) {
142: return ((cmo_int32 *)c)->i;
143: }else if(c->tag == CMO_ZZ) {
144: return mpz_get_si(((cmo_zz *)c)->mpz);
145: }
146: make_error2(-1);
147: return 0;
1.1 noro 148: }
149:
150: char *get_str()
151: {
1.4 noro 152: cmo *c = pop();
153: if(c->tag == CMO_STRING) {
154: return ((cmo_string *)c)->s;
155: }
156: make_error2(-1);
157: return "";
1.1 noro 158: }
159:
160: int cmo2int(cmo *c)
161: {
1.4 noro 162: if(c->tag == CMO_INT32) {
163: return ((cmo_int32 *)c)->i;
164: }else if(c->tag == CMO_ZZ) {
165: return mpz_get_si(((cmo_zz *)c)->mpz);
166: } else if(c->tag == CMO_NULL){
167: return 0;
168: }
1.1 noro 169:
1.4 noro 170: return 0;
1.1 noro 171: }
172:
1.6 ! noro 173: GEN cmo_int32_to_GEN(cmo_int32 *c)
! 174: {
! 175: GEN z;
! 176: int i,sgn;
! 177:
! 178: i = c->i;
! 179: if ( !i ) return gen_0;
! 180: z = cgeti(3);
! 181: sgn = 1;
! 182: if ( i < 0 ) {
! 183: i = -i;
! 184: sgn = -1;
! 185: }
! 186: z[2] = i;
! 187: setsigne(z,sgn);
! 188: setlgefint(z,lg(z));
! 189: return z;
! 190: }
! 191:
1.2 noro 192: GEN cmo_zz_to_GEN(cmo_zz *c)
193: {
194: mpz_ptr mpz;
195: GEN z;
196: long *ptr;
197: int j,sgn,len;
198:
199: mpz = c->mpz;
200: sgn = mpz_sgn(mpz);
201: len = ABSIZ(mpz);
202: ptr = (long *)PTR(mpz);
203: z = cgeti(len+2);
204: for ( j = 0; j < len; j++ )
205: z[len-j+1] = ptr[j];
206: setsigne(z,sgn);
207: setlgefint(z,lg(z));
208: return z;
209: }
210:
1.6 ! noro 211: GEN cmo_qq_to_GEN(cmo_qq *c)
! 212: {
! 213: GEN z,nm,den;
! 214:
! 215: z = cgetg(3,4);
! 216: nm = cmo_zz_to_GEN(new_cmo_zz_set_mpz(mpq_numref(c->mpq)));
! 217: den = cmo_zz_to_GEN(new_cmo_zz_set_mpz(mpq_denref(c->mpq)));
! 218: z[1] = (long)nm;
! 219: z[2] = (long)den;
! 220: return z;
! 221: }
! 222:
1.3 noro 223: GEN cmo_bf_to_GEN(cmo_bf *c)
224: {
225: mpfr_ptr mpfr;
226: GEN z;
227: int sgn,len,j;
228: long exp;
229: long *ptr;
230:
231: mpfr = c->mpfr;
232: sgn = MPFR_SIGN(mpfr);
233: exp = MPFR_EXP(mpfr)-1;
234: len = MPFR_LIMB_SIZE(mpfr);
235: ptr = (long *)MPFR_MANT(mpfr);
236: z = cgetr(len+2);
237: for ( j = 0; j < len; j++ )
238: z[len-j+1] = ptr[j];
239: z[1] = evalsigne(sgn)|evalexpo(exp);
240: setsigne(z,sgn);
241: return z;
242: }
243:
1.6 ! noro 244: /* list->vector */
! 245:
! 246: GEN cmo_list_to_GEN(cmo_list *c)
! 247: {
! 248: GEN z;
! 249: cell *cell;
! 250:
! 251: z = cgetg(c->length+1,17);
! 252: for ( cell = c->head->next; cell != c->head; cell = cell->next ) {
! 253: z[cell->exp] = (long)cmo_to_GEN(cell->cmo);
! 254: }
! 255: return z;
! 256: }
! 257:
! 258: GEN cmo_complex_to_GEN(cmo_complex *c)
! 259: {
! 260: GEN z;
! 261:
! 262: z = cgetg(3,6);
! 263: z[1] = (long)cmo_to_GEN(c->re);
! 264: z[2] = (long)cmo_to_GEN(c->im);
! 265: return z;
! 266: }
! 267:
! 268: GEN cmo_up_to_GEN(cmo_polynomial_in_one_variable *c)
! 269: {
! 270: GEN z;
! 271: int d,i;
! 272: cell *cell;
! 273:
! 274: d = c->head->next->exp;
! 275: z = cgetg(d+3,10);
! 276: setsigne(z,1);
! 277: setvarn(z,c->var);
! 278: setlgef(z,d+3);
! 279: for ( i = 2; i <= d+2; i++ )
! 280: z[i] = (long)gen_0;
! 281: for ( cell = c->head->next; cell != c->head; cell = cell->next ) {
! 282: z[2+cell->exp] = (long)cmo_to_GEN(cell->cmo);
! 283: }
! 284: return z;
! 285: }
! 286:
! 287: cmo_list *current_ringdef;
! 288:
! 289: void register_variables(cmo_list *ringdef)
! 290: {
! 291: current_ringdef = ringdef;
! 292: }
! 293:
! 294: GEN cmo_rp_to_GEN(cmo_recursive_polynomial *c)
! 295: {
! 296: register_variables(c->ringdef);
! 297: switch ( c->coef->tag ) {
! 298: case CMO_ZERO:
! 299: return gen_0;
! 300: case CMO_INT32:
! 301: return cmo_int32_to_GEN((cmo_int32 *)c->coef);
! 302: case CMO_ZZ:
! 303: return cmo_zz_to_GEN((cmo_zz *)c->coef);
! 304: case CMO_QQ:
! 305: return cmo_qq_to_GEN((cmo_qq *)c->coef);
! 306: case CMO_POLYNOMIAL_IN_ONE_VARIABLE:
! 307: return cmo_up_to_GEN((cmo_polynomial_in_one_variable *)c->coef);
! 308: default:
! 309: return 0;
! 310: }
! 311: }
! 312:
1.2 noro 313: cmo_zz *GEN_to_cmo_zz(GEN z)
314: {
315: cmo_zz *c;
316:
317: c = new_cmo_zz();
318: mpz_import(c->mpz,lgef(z)-2,1,sizeof(long),0,0,&z[2]);
319: if ( signe(z) < 0 )
320: mpz_neg(c->mpz,c->mpz);
321: return c;
322: }
323:
1.3 noro 324: cmo_bf *GEN_to_cmo_bf(GEN z)
325: {
326: cmo_bf *c;
327: int len,prec,j;
328: long *ptr;
329:
330: c = new_cmo_bf();
331: len = lg(z)-2;
332: prec = len*sizeof(long)*8;
333: mpfr_init2(c->mpfr,prec);
334: ptr = (long *)MPFR_MANT(c->mpfr);
335: for ( j = 0; j < len; j++ )
336: ptr[j] = z[len-j+1];
337: MPFR_EXP(c->mpfr) = (long long)(expo(z)+1);
338: MPFR_SIGN(c->mpfr) = gsigne(z);
339: return c;
340: }
341:
342:
1.2 noro 343: cmo_list *GEN_to_cmo_list(GEN z)
344: {
345: cmo_list *c;
346: cmo *ob;
347: int i,len;
348:
349: c = new_cmo_list();
350: len = lg(z)-1;
351: for ( i = 1; i <= len; i++ ) {
352: ob = GEN_to_cmo((GEN)z[i]);
353: c = list_append(c,ob);
354: }
355: return c;
356: }
357:
1.6 ! noro 358: cmo_complex *GEN_to_cmo_complex(GEN z)
! 359: {
! 360: cmo_complex *c;
! 361:
! 362: c = new_cmo_complex();
! 363: c->re = GEN_to_cmo((GEN)z[1]);
! 364: c->im = GEN_to_cmo((GEN)z[2]);
! 365: return c;
! 366: }
! 367:
1.2 noro 368:
369: GEN cmo_to_GEN(cmo *c)
370: {
371: switch ( c->tag ) {
372: case CMO_ZERO:
1.3 noro 373: return gen_0;
1.2 noro 374: case CMO_ZZ: /* int */
375: return cmo_zz_to_GEN((cmo_zz *)c);
1.3 noro 376: case CMO_BIGFLOAT: /* bigfloat */
377: return cmo_bf_to_GEN((cmo_bf *)c);
1.6 ! noro 378: case CMO_LIST:
! 379: return cmo_list_to_GEN((cmo_list *)c);
! 380: case CMO_RECURSIVE_POLYNOMIAL:
! 381: return cmo_rp_to_GEN((cmo_recursive_polynomial *)c);
! 382: case CMO_POLYNOMIAL_IN_ONE_VARIABLE:
! 383: return cmo_up_to_GEN((cmo_polynomial_in_one_variable *)c);
1.2 noro 384: default:
385: return 0;
386: }
387: }
388:
389: cmo *GEN_to_cmo(GEN z)
390: {
391: if ( gcmp0(z) )
392: return new_cmo_zero();
393: switch ( typ(z) ) {
394: case 1: /* int */
395: return (cmo *)GEN_to_cmo_zz(z);
1.3 noro 396: case 2: /* bigfloat */
397: return (cmo *)GEN_to_cmo_bf(z);
1.6 ! noro 398: case 6: /* complex */
! 399: return (cmo *)GEN_to_cmo_complex(z);
1.2 noro 400: case 17: case 18: /* vector */
401: return (cmo *)GEN_to_cmo_list(z);
402: case 19: /* matrix */
403: return (cmo *)GEN_to_cmo_list(shallowtrans(z));
404: default:
405: return (cmo *)make_error2(typ(z));
406: }
407: }
408:
1.3 noro 409: struct parif {
410: char *name;
1.4 noro 411: GEN (*f)();
1.3 noro 412: int type;
413: } parif_tab[] = {
1.4 noro 414: /* (ulong)allocatemoremem(ulong) */
415: {"allocatemem",(GEN (*)())allocatemoremem,0},
416: /* num/num */
417: {"abs",gabs,1},
418: {"erfc",gerfc,1},
419: {"arg",garg,1},
420: {"isqrt",racine,1},
421: {"gamma",ggamma,1},
422: {"zeta",gzeta,1},
423: {"floor",gfloor,1},
424: {"frac",gfrac,1},
425: {"imag",gimag,1},
426: {"conj",gconj,1},
427: {"ceil",gceil,1},
428: {"isprime",gisprime,2},
429: {"bigomega",gbigomega,1},
430: {"denom",denom,1},
431: {"numer",numer,1},
432: {"lngamma",glngamma,1},
433: {"logagm",glogagm,1},
434: {"classno",classno,1},
435: {"classno2",classno2,1},
436: {"dilog",dilog,1},
437: {"disc",discsr,1},
438: {"discf",discf,1},
439: {"nextprime",nextprime,1},
440: {"eintg1",eint1,1},
441: {"eta",eta,1},
442: {"issqfree",gissquarefree,1},
443: {"issquare",gcarreparfait,1},
444: {"gamh",ggamd,1},
445: {"hclassno",classno3,1},
446:
447: /* num/array */
448: {"binary",binaire,1},
449: {"factorint",factorint,2},
450: {"factor",Z_factor,1},
451: {"cf",gcf,1},
452: {"divisors",divisors,1},
453: {"smallfact",smallfact,1},
454:
455: /* mat/mat */
456: {"adj",adj,1},
457: {"lll",lll,1},
458: {"lllgen",lllgen,1},
459: {"lllgram",lllgram,1},
460: {"lllgramgen",lllgramgen,1},
461: {"lllgramint",lllgramint,1},
462: {"lllgramkerim",lllgramkerim,1},
463: {"lllgramkerimgen",lllgramkerimgen,1},
464: {"lllint",lllint,1},
465: {"lllkerim",lllkerim,1},
466: {"lllkerimgen",lllkerimgen,1},
467: {"trans",gtrans,1},
468: {"eigen",eigen,1},
469: {"hermite",hnf,1},
470: {"mat",gtomat,1},
471: {"matrixqz2",matrixqz2,1},
472: {"matrixqz3",matrixqz3,1},
473: {"hess",hess,1},
474:
475: /* mat/poly */
476: {"det",det,1},
477: {"det2",det2,1},
478:
479: /* poly/poly */
480: {"centerlift",centerlift,1},
481: {"content",content,1},
482:
483: /* poly/array */
484: {"galois",galois,1},
485: {"roots",roots,1},
486:
1.3 noro 487: };
1.2 noro 488:
489: #define PARI_MAX_AC 64
490:
1.3 noro 491: struct parif *search_parif(char *name)
492: {
493: int tablen,i;
494:
495: tablen = sizeof(parif_tab)/sizeof(struct parif);
496: for ( i = 0; i < tablen; i++ ) {
497: if ( !strcmp(parif_tab[i].name,name) )
498: return &parif_tab[i];
499: }
500: return 0;
501: }
502:
1.1 noro 503: int sm_executeFunction()
504: {
1.5 noro 505: pari_sp av0;
1.2 noro 506: int ac,i;
507: cmo_int32 *c;
508: cmo *av[PARI_MAX_AC];
509: cmo *ret;
510: GEN z,m;
1.3 noro 511: struct parif *parif;
1.2 noro 512:
1.3 noro 513: if ( setjmp(GP_DATA->env) ) {
1.4 noro 514: printf("sm_executeFunction : an error occured.\n");fflush(stdout);
515: push((cmo*)make_error2(0));
516: return -1;
517: }
518: cmo_string *func = (cmo_string *)pop();
519: if(func->tag != CMO_STRING) {
520: printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout);
521: push((cmo*)make_error2(0));
522: return -1;
523: }
1.1 noro 524:
1.4 noro 525: c = (cmo_int32 *)pop();
1.2 noro 526: ac = c->i;
527: if ( ac > PARI_MAX_AC ) {
1.4 noro 528: push((cmo*)make_error2(0));
529: return -1;
1.2 noro 530: }
531: for ( i = 0; i < ac; i++ ) {
532: av[i] = (cmo *)pop();
533: fprintf(stderr,"arg%d:",i);
534: print_cmo(av[i]);
535: fprintf(stderr,"\n");
536: }
1.4 noro 537: if( strcmp( func->s, "exit" ) == 0 )
538: exit(0);
1.3 noro 539:
540: parif =search_parif(func->s);
541: if ( !parif ) {
1.4 noro 542: push((cmo*)make_error2(0));
543: return -1;
1.3 noro 544: } else if ( parif->type == 0 ) {
545: /* one long int variable */
546: int a = cmo_to_int(av[0]);
1.4 noro 547: a = (int)(parif->f)(a);
1.3 noro 548: ret = (cmo *)new_cmo_int32(a);
1.2 noro 549: push(ret);
1.4 noro 550: return 0;
1.3 noro 551: } else if ( parif->type == 1 ) {
552: /* one variable possibly with prec */
553: unsigned long prec;
554:
1.5 noro 555: av0 = avma;
1.2 noro 556: z = cmo_to_GEN(av[0]);
1.3 noro 557: if ( ac == 2 ) {
558: prec = cmo_to_int(av[1])*3.32193/32+3;
559: } else
560: prec = precreal;
1.6 ! noro 561: printf("input : ");
! 562: output(z);
1.3 noro 563: m = (*parif->f)(z,prec);
1.2 noro 564: ret = GEN_to_cmo(m);
1.5 noro 565: avma = av0;
1.2 noro 566: push(ret);
1.4 noro 567: return 0;
1.3 noro 568: } else {
1.4 noro 569: push((cmo*)make_error2(0));
570: return -1;
1.3 noro 571: }
1.1 noro 572: }
573:
574: int receive_and_execute_sm_command()
575: {
1.4 noro 576: int code = receive_int32(fd_rw);
577: switch(code) {
578: case SM_popCMO:
579: sm_popCMO();
580: break;
581: case SM_executeFunction:
582: sm_executeFunction();
583: break;
584: case SM_mathcap:
585: sm_mathcap();
586: break;
587: case SM_setMathCap:
588: pop();
589: break;
590: default:
591: printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);
592: break;
593: }
594: return 0;
1.1 noro 595: }
596:
597: int receive()
598: {
1.4 noro 599: int tag;
1.1 noro 600:
1.4 noro 601: tag = receive_ox_tag(fd_rw);
602: switch(tag) {
603: case OX_DATA:
604: printf("receive : ox_data %d\n",tag);fflush(stdout);
605: push(receive_cmo(fd_rw));
606: break;
607: case OX_COMMAND:
608: printf("receive : ox_command %d\n",tag);fflush(stdout);
609: receive_and_execute_sm_command();
610: break;
611: default:
612: printf("receive : tag=%d\n",tag);fflush(stdout);
613: }
614: return 0;
1.1 noro 615: }
616:
617: int main()
618: {
1.3 noro 619: init_gc();
1.4 noro 620: ox_stderr_init(stderr);
621: initialize_stack();
622: init_pari();
623:
624: fprintf(stderr,"ox_pari\n");
625:
626: fd_rw = oxf_open(3);
627: oxf_determine_byteorder_server(fd_rw);
628:
629: while(1){
630: receive();
631: }
1.1 noro 632: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ox_pari.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / ox_pari