Annotation of OpenXM_contrib2/asir2018/builtin/bfaux.c, Revision 1.1
1.1 ! noro 1: /* $OpenXM$ */
! 2: #include "ca.h"
! 3: #include "parse.h"
! 4:
! 5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
! 6: void Pmpfr_ai();
! 7: void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
! 8: void Pmpfr_zeta();
! 9: void Pmpfr_j0(), Pmpfr_j1();
! 10: void Pmpfr_y0(), Pmpfr_y1();
! 11: void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
! 12: void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
! 13: void Prk_ratmat();
! 14: void mp_sin(),mp_cos(),mp_tan(),mp_asin(),mp_acos(),mp_atan();
! 15: void mp_sinh(),mp_cosh(),mp_tanh(),mp_asinh(),mp_acosh(),mp_atanh();
! 16: void mp_exp(),mp_log(),mp_pow();
! 17: void mp_factorial(),mp_abs();
! 18:
! 19: struct ftab bf_tab[] = {
! 20: {"eval",Peval,-2},
! 21: {"setprec",Psetprec,-1},
! 22: {"setbprec",Psetbprec,-1},
! 23: {"setround",Psetround,-1},
! 24: {"todouble",Ptodouble,1},
! 25: {"mpfr_sin",mp_sin,-2},
! 26: {"mpfr_cos",mp_cos,-2},
! 27: {"mpfr_tan",mp_tan,-2},
! 28: {"mpfr_asin",mp_asin,-2},
! 29: {"mpfr_acos",mp_acos,-2},
! 30: {"mpfr_atan",mp_atan,-2},
! 31: {"mpfr_sinh",mp_sinh,-2},
! 32: {"mpfr_cosh",mp_cosh,-2},
! 33: {"mpfr_tanh",mp_tanh,-2},
! 34: {"mpfr_asinh",mp_asinh,-2},
! 35: {"mpfr_acosh",mp_acosh,-2},
! 36: {"mpfr_atanh",mp_atanh,-2},
! 37: {"mpfr_exp",mp_exp,-2},
! 38: {"mpfr_log",mp_log,-2},
! 39: {"mpfr_pow",mp_pow,-3},
! 40: {"mpfr_ai",Pmpfr_ai,-2},
! 41: {"mpfr_zeta",Pmpfr_zeta,-2},
! 42: {"mpfr_j0",Pmpfr_j0,-2},
! 43: {"mpfr_j1",Pmpfr_j1,-2},
! 44: {"mpfr_y0",Pmpfr_y0,-2},
! 45: {"mpfr_y1",Pmpfr_y1,-2},
! 46: {"mpfr_eint",Pmpfr_eint,-2},
! 47: {"mpfr_erf",Pmpfr_erf,-2},
! 48: {"mpfr_erfc",Pmpfr_erfc,-2},
! 49: {"mpfr_li2",Pmpfr_li2,-2},
! 50: {"mpfr_gamma",Pmpfr_gamma,-2},
! 51: {"mpfr_lngamma",Pmpfr_gamma,-2},
! 52: {"mpfr_digamma",Pmpfr_gamma,-2},
! 53: {"mpfr_floor",Pmpfr_floor,-2},
! 54: {"mpfr_ceil",Pmpfr_ceil,-2},
! 55: {"mpfr_round",Pmpfr_round,-2},
! 56: {"rk_ratmat",Prk_ratmat,7},
! 57: {0,0,0},
! 58: };
! 59:
! 60: int mpfr_roundmode = MPFR_RNDN;
! 61:
! 62: void todoublen(Num a,Num *rp)
! 63: {
! 64: double r,i;
! 65: Real real,imag;
! 66:
! 67: if ( !a ) {
! 68: *rp = 0;
! 69: return;
! 70: }
! 71: switch ( NID(a) ) {
! 72: case N_R: case N_Q: case N_B:
! 73: r = ToReal(a);
! 74: MKReal(r,real);
! 75: *rp = (Num)real;
! 76: break;
! 77: case N_C:
! 78: r = ToReal(((C)a)->r);
! 79: i = ToReal(((C)a)->i);
! 80: MKReal(r,real);
! 81: MKReal(i,imag);
! 82: reimtocplx((Num)real,(Num)imag,rp);
! 83: break;
! 84: default:
! 85: *rp = a;
! 86: break;
! 87: }
! 88: }
! 89:
! 90: void todoublep(P a,P *rp)
! 91: {
! 92: DCP dc,dcr,dcr0;
! 93:
! 94: if ( !a ) *rp = 0;
! 95: else if ( OID(a) == O_N ) todoublen((Num)a,(Num *)rp);
! 96: else {
! 97: for ( dcr0 = 0, dc = DC(a); dc; dc = NEXT(dc) ) {
! 98: NEXTDC(dcr0,dcr);
! 99: DEG(dcr) = DEG(dc);
! 100: todoublep(COEF(dc),&COEF(dcr));
! 101: }
! 102: NEXT(dcr) = 0;
! 103: MKP(VR(a),dcr0,*rp);
! 104: }
! 105: }
! 106:
! 107: void todoubler(R a,R *rp)
! 108: {
! 109: R b;
! 110:
! 111: if ( !a ) *rp = 0;
! 112: else if ( OID(a) <= O_P ) todoublep((P)a,(P *)rp);
! 113: else {
! 114: NEWR(b);
! 115: todoublep(a->nm,&b->nm);
! 116: todoublep(a->dn,&b->dn);
! 117: *rp = b;
! 118: }
! 119: }
! 120:
! 121: void todouble(Obj a,Obj *b)
! 122: {
! 123: Obj t;
! 124: LIST l;
! 125: V v;
! 126: int row,col,len;
! 127: VECT vect;
! 128: MAT mat;
! 129: int i,j;
! 130: NODE n0,n,nd;
! 131: MP m,mp,mp0;
! 132: DP d;
! 133:
! 134: if ( !a ) {
! 135: *b = 0;
! 136: return;
! 137: }
! 138: switch ( OID(a) ) {
! 139: case O_N:
! 140: todoublen((Num)a,(Num *)b);
! 141: break;
! 142: case O_P:
! 143: todoublep((P)a,(P *)b);
! 144: break;
! 145: case O_R:
! 146: todoubler((R)a,(R *)b);
! 147: break;
! 148: case O_LIST:
! 149: n0 = 0;
! 150: for ( nd = BDY((LIST)a); nd; nd = NEXT(nd) ) {
! 151: NEXTNODE(n0,n);
! 152: todouble((Obj)BDY(nd),(Obj *)&BDY(n));
! 153: }
! 154: if ( n0 )
! 155: NEXT(n) = 0;
! 156: MKLIST(l,n0);
! 157: *b = (Obj)l;
! 158: break;
! 159: case O_VECT:
! 160: len = ((VECT)a)->len;
! 161: MKVECT(vect,len);
! 162: for ( i = 0; i < len; i++ ) {
! 163: todouble((Obj)BDY((VECT)a)[i],(Obj *)&BDY(vect)[i]);
! 164: }
! 165: *b = (Obj)vect;
! 166: break;
! 167: case O_MAT:
! 168: row = ((MAT)a)->row;
! 169: col = ((MAT)a)->col;
! 170: MKMAT(mat,row,col);
! 171: for ( i = 0; i < row; i++ )
! 172: for ( j = 0; j < col; j++ ) {
! 173: todouble((Obj)BDY((MAT)a)[i][j],(Obj *)&BDY(mat)[i][j]);
! 174: }
! 175: *b = (Obj)mat;
! 176: break;
! 177: case O_DP:
! 178: mp0 = 0;
! 179: for ( m = BDY((DP)a); m; m = NEXT(m) ) {
! 180: todouble(C(m),&t);
! 181: if ( t ) {
! 182: NEXTMP(mp0,mp);
! 183: C(mp) = t;
! 184: mp->dl = m->dl;
! 185: }
! 186: }
! 187: if ( mp0 ) {
! 188: MKDP(NV((DP)a),mp0,d);
! 189: d->sugar = ((DP)a)->sugar;
! 190: *b = (Obj)d;
! 191: } else
! 192: *b = 0;
! 193:
! 194: break;
! 195: default:
! 196: error("todouble : invalid argument");
! 197: }
! 198: }
! 199:
! 200: void Ptodouble(NODE arg,Obj *rp)
! 201: {
! 202: todouble((Obj)ARG0(arg),rp);
! 203: }
! 204:
! 205: void Peval(NODE arg,Obj *rp)
! 206: {
! 207: long prec;
! 208:
! 209: asir_assert(ARG0(arg),O_R,"eval");
! 210: if ( argc(arg) == 2 ) {
! 211: prec = QTOS((Q)ARG1(arg))*3.32193;
! 212: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
! 213: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
! 214: } else
! 215: prec = 0;
! 216: evalr(CO,(Obj)ARG0(arg),prec,rp);
! 217: }
! 218:
! 219: /* set/get decimal precision */
! 220:
! 221: void Psetprec(NODE arg,Obj *rp)
! 222: {
! 223: long p;
! 224: Z q;
! 225: long prec,dprec;
! 226:
! 227: prec = mpfr_get_default_prec();
! 228: /* decimal precision */
! 229: dprec = prec*0.30103;
! 230: STOQ(dprec,q); *rp = (Obj)q;
! 231: if ( arg ) {
! 232: asir_assert(ARG0(arg),O_N,"setprec");
! 233: p = QTOS((Q)ARG0(arg))*3.32193;
! 234: if ( p > 0 )
! 235: prec = p;
! 236: }
! 237: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
! 238: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
! 239: mpfr_set_default_prec(prec);
! 240: }
! 241:
! 242: /* set/get bit precision */
! 243:
! 244: void Psetbprec(NODE arg,Obj *rp)
! 245: {
! 246: long p;
! 247: Z q;
! 248: long prec;
! 249:
! 250: prec = mpfr_get_default_prec();
! 251: STOQ(prec,q); *rp = (Obj)q;
! 252: if ( arg ) {
! 253: asir_assert(ARG0(arg),O_N,"setbprec");
! 254: p = QTOS((Q)ARG0(arg));
! 255: if ( p > 0 )
! 256: prec = p;
! 257: }
! 258: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
! 259: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
! 260: mpfr_set_default_prec(prec);
! 261: }
! 262:
! 263: void Psetround(NODE arg,Z *rp)
! 264: {
! 265: int round;
! 266:
! 267: STOQ(mpfr_roundmode,*rp);
! 268: if ( arg ) {
! 269: asir_assert(ARG0(arg),O_N,"setround");
! 270: round = QTOS((Q)ARG0(arg));
! 271: switch ( round ) {
! 272: case 0:
! 273: mpfr_roundmode = MPFR_RNDN;
! 274: break;
! 275: case 1:
! 276: mpfr_roundmode = MPFR_RNDZ;
! 277: break;
! 278: case 2:
! 279: mpfr_roundmode = MPFR_RNDU;
! 280: break;
! 281: case 3:
! 282: mpfr_roundmode = MPFR_RNDD;
! 283: break;
! 284: case 4:
! 285: mpfr_roundmode = MPFR_RNDA;
! 286: break;
! 287: case 5:
! 288: mpfr_roundmode = MPFR_RNDF;
! 289: break;
! 290: case 6:
! 291: mpfr_roundmode = MPFR_RNDNA;
! 292: break;
! 293: default:
! 294: error("setround : invalid rounding mode");
! 295: break;
! 296: }
! 297: }
! 298: }
! 299:
! 300: Num tobf(Num a,int prec);
! 301:
! 302: void mp_pi(NODE arg,BF *rp)
! 303: {
! 304: int prec;
! 305: BF r;
! 306:
! 307: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
! 308: NEWBF(r);
! 309: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
! 310: mpfr_const_pi(r->body,mpfr_roundmode);
! 311: if ( !cmpbf((Num)r,0) ) r = 0;
! 312: *rp = r;
! 313: }
! 314:
! 315: void mp_e(NODE arg,BF *rp)
! 316: {
! 317: int prec;
! 318: mpfr_t one;
! 319: BF r;
! 320:
! 321: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
! 322: NEWBF(r);
! 323: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
! 324: mpfr_init(one);
! 325: mpfr_set_ui(one,1,mpfr_roundmode);
! 326: mpfr_exp(r->body,one,mpfr_roundmode);
! 327: if ( !cmpbf((Num)r,0) ) r = 0;
! 328: *rp = r;
! 329: }
! 330:
! 331: void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
! 332: {
! 333: Num a;
! 334: int prec;
! 335: BF r,re,im;
! 336: C c;
! 337: mpc_t mpc,a1;
! 338:
! 339: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
! 340: a = tobf(ARG0(arg),prec);
! 341: if ( a && NID(a)==N_C ) {
! 342: mpc_init2(mpc,prec); mpc_init2(a1,prec);
! 343: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
! 344: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
! 345: (*mpc_f)(mpc,a1,mpfr_roundmode);
! 346: MPFRTOBF(mpc_realref(mpc),re);
! 347: MPFRTOBF(mpc_imagref(mpc),im);
! 348: if ( !cmpbf((Num)re,0) ) re = 0;
! 349: if ( !cmpbf((Num)im,0) ) im = 0;
! 350: if ( !im )
! 351: *rp = (Num)re;
! 352: else {
! 353: NEWC(c); c->r = (Num)re; c->i = (Num)im;
! 354: *rp = (Num)c;
! 355: }
! 356: } else {
! 357: NEWBF(r);
! 358: mpfr_init2(r->body,prec);
! 359: (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode);
! 360: if ( !cmpbf((Num)r,0) ) r = 0;
! 361: *rp = (Num)r;
! 362: }
! 363: }
! 364:
! 365: void mp_sin(NODE arg,Num *rp)
! 366: {
! 367: mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
! 368: }
! 369:
! 370: void mp_cos(NODE arg,Num *rp)
! 371: {
! 372: mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
! 373: }
! 374:
! 375: void mp_tan(NODE arg,Num *rp)
! 376: {
! 377: mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
! 378: }
! 379:
! 380: void mp_asin(NODE arg,Num *rp)
! 381: {
! 382: mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
! 383: }
! 384:
! 385: void mp_acos(NODE arg,Num *rp)
! 386: {
! 387: mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
! 388: }
! 389:
! 390: void mp_atan(NODE arg,Num *rp)
! 391: {
! 392: mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
! 393: }
! 394:
! 395: void mp_sinh(NODE arg,Num *rp)
! 396: {
! 397: mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
! 398: }
! 399:
! 400: void mp_cosh(NODE arg,Num *rp)
! 401: {
! 402: mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
! 403: }
! 404:
! 405: void mp_tanh(NODE arg,Num *rp)
! 406: {
! 407: mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
! 408: }
! 409:
! 410: void mp_asinh(NODE arg,Num *rp)
! 411: {
! 412: mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
! 413: }
! 414:
! 415: void mp_acosh(NODE arg,Num *rp)
! 416: {
! 417: mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
! 418: }
! 419:
! 420: void mp_atanh(NODE arg,Num *rp)
! 421: {
! 422: mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
! 423: }
! 424:
! 425: void mp_exp(NODE arg,Num *rp)
! 426: {
! 427: mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
! 428: }
! 429:
! 430: void mp_log(NODE arg,Num *rp)
! 431: {
! 432: mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
! 433: }
! 434:
! 435: void mp_abs(NODE arg,Num *rp)
! 436: {
! 437: mpfr_or_mpc(arg,mpfr_abs,mpc_abs,rp);
! 438: }
! 439:
! 440: void Pfac(NODE arg,Num *rp);
! 441:
! 442: void mp_factorial(NODE arg,Num *rp)
! 443: {
! 444: struct oNODE arg0;
! 445: Num a,a1;
! 446:
! 447: a = (Num)ARG0(arg);
! 448: if ( !a ) *rp = (Num)ONE;
! 449: else if ( INT(a) ) Pfac(arg,rp);
! 450: else {
! 451: addnum(0,a,(Num)ONE,&a1);
! 452: arg0.body = (pointer)a1;
! 453: arg0.next = arg->next;
! 454: Pmpfr_gamma(&arg0,rp);
! 455: }
! 456: }
! 457:
! 458: void mp_pow(NODE arg,Num *rp)
! 459: {
! 460: Num a,e;
! 461: int prec;
! 462: BF r,re,im;
! 463: C c;
! 464: mpc_t mpc,a1,e1;
! 465:
! 466: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
! 467: a = tobf(ARG0(arg),prec);
! 468: e = tobf(ARG1(arg),prec);
! 469: if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
! 470: mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
! 471: if ( NID(a) == N_C ) {
! 472: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
! 473: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
! 474: } else {
! 475: re = (BF)a;
! 476: mpc_set_fr(a1,re->body,mpfr_roundmode);
! 477: }
! 478: if ( NID(e) == N_C ) {
! 479: re = (BF)((C)e)->r; im = (BF)((C)e)->i;
! 480: mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
! 481: } else {
! 482: re = (BF)e;
! 483: mpc_set_fr(e1,re->body,mpfr_roundmode);
! 484: }
! 485: mpc_pow(mpc,a1,e1,mpfr_roundmode);
! 486: MPFRTOBF(mpc_realref(mpc),re);
! 487: MPFRTOBF(mpc_imagref(mpc),im);
! 488: if ( !cmpbf((Num)re,0) ) re = 0;
! 489: if ( !cmpbf((Num)im,0) ) im = 0;
! 490: if ( !im )
! 491: *rp = (Num)re;
! 492: else {
! 493: NEWC(c); c->r = (Num)re; c->i = (Num)im;
! 494: *rp = (Num)c;
! 495: }
! 496: } else {
! 497: NEWBF(r);
! 498: mpfr_init2(r->body,prec);
! 499: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
! 500: *rp = (Num)r;
! 501: }
! 502: }
! 503:
! 504: #define SETPREC \
! 505: (prec)=NEXT(arg)?QTOS((Q)ARG1(arg)):0;\
! 506: (prec)*=3.32193;\
! 507: (a)=tobf(ARG0(arg),prec);\
! 508: NEWBF(r);\
! 509: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
! 510:
! 511:
! 512: void Pmpfr_gamma(NODE arg,BF *rp)
! 513: {
! 514: Num a;
! 515: int prec;
! 516: BF r;
! 517:
! 518: SETPREC
! 519: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
! 520: *rp = r;
! 521: }
! 522:
! 523: void Pmpfr_lngamma(NODE arg,BF *rp)
! 524: {
! 525: Num a;
! 526: int prec;
! 527: BF r;
! 528:
! 529: SETPREC
! 530: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
! 531: *rp = r;
! 532: }
! 533:
! 534: void Pmpfr_digamma(NODE arg,BF *rp)
! 535: {
! 536: Num a;
! 537: int prec;
! 538: BF r;
! 539:
! 540: SETPREC
! 541: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
! 542: *rp = r;
! 543: }
! 544:
! 545: void Pmpfr_zeta(NODE arg,BF *rp)
! 546: {
! 547: Num a;
! 548: int prec;
! 549: BF r;
! 550:
! 551: SETPREC
! 552: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
! 553: *rp = r;
! 554: }
! 555:
! 556: void Pmpfr_eint(NODE arg,BF *rp)
! 557: {
! 558: Num a;
! 559: int prec;
! 560: BF r;
! 561:
! 562: SETPREC
! 563: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
! 564: *rp = r;
! 565: }
! 566:
! 567: void Pmpfr_erf(NODE arg,BF *rp)
! 568: {
! 569: Num a;
! 570: int prec;
! 571: BF r;
! 572:
! 573: SETPREC
! 574: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
! 575: *rp = r;
! 576: }
! 577:
! 578: void Pmpfr_erfc(NODE arg,BF *rp)
! 579: {
! 580: Num a;
! 581: int prec;
! 582: BF r;
! 583:
! 584: SETPREC
! 585: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
! 586: *rp = r;
! 587: }
! 588:
! 589: void Pmpfr_j0(NODE arg,BF *rp)
! 590: {
! 591: Num a;
! 592: int prec;
! 593: BF r;
! 594:
! 595: SETPREC
! 596: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
! 597: *rp = r;
! 598: }
! 599:
! 600: void Pmpfr_j1(NODE arg,BF *rp)
! 601: {
! 602: Num a;
! 603: int prec;
! 604: BF r;
! 605:
! 606: SETPREC
! 607: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
! 608: *rp = r;
! 609: }
! 610:
! 611: void Pmpfr_y0(NODE arg,BF *rp)
! 612: {
! 613: Num a;
! 614: int prec;
! 615: BF r;
! 616:
! 617: SETPREC
! 618: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
! 619: *rp = r;
! 620: }
! 621:
! 622: void Pmpfr_y1(NODE arg,BF *rp)
! 623: {
! 624: Num a;
! 625: int prec;
! 626: BF r;
! 627:
! 628: SETPREC
! 629: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
! 630: *rp = r;
! 631: }
! 632:
! 633: void Pmpfr_li2(NODE arg,BF *rp)
! 634: {
! 635: Num a;
! 636: int prec;
! 637: BF r;
! 638:
! 639: SETPREC
! 640: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
! 641: *rp = r;
! 642: }
! 643:
! 644: void Pmpfr_ai(NODE arg,BF *rp)
! 645: {
! 646: Num a;
! 647: int prec;
! 648: BF r;
! 649:
! 650: SETPREC
! 651: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
! 652: *rp = r;
! 653: }
! 654:
! 655: void Pmpfr_floor(NODE arg,Z *rp)
! 656: {
! 657: Num a;
! 658: long prec;
! 659: BF r;
! 660: mpz_t t;
! 661:
! 662: SETPREC
! 663: mpfr_floor(r->body,((BF)a)->body);
! 664: mpz_init(t);
! 665: mpfr_get_z(t,r->body,mpfr_roundmode);
! 666: MPZTOZ(t,*rp);
! 667: }
! 668:
! 669: void Pmpfr_ceil(NODE arg,Z *rp)
! 670: {
! 671: Num a;
! 672: long prec;
! 673: BF r;
! 674: mpz_t t;
! 675: Z rz;
! 676:
! 677: SETPREC
! 678: mpfr_ceil(r->body,((BF)a)->body);
! 679: mpz_init(t);
! 680: mpfr_get_z(t,r->body,mpfr_roundmode);
! 681: MPZTOZ(t,*rp);
! 682: }
! 683:
! 684: void Pmpfr_round(NODE arg,Z *rp)
! 685: {
! 686: Num a;
! 687: long prec;
! 688: BF r;
! 689: mpz_t t;
! 690:
! 691: SETPREC
! 692: mpfr_round(r->body,((BF)a)->body);
! 693: mpz_init(t);
! 694: mpfr_get_z(t,r->body,mpfr_roundmode);
! 695: MPZTOZ(t,*rp);
! 696: }
! 697:
! 698: double **almat_double(int n)
! 699: {
! 700: int i;
! 701: double **a;
! 702:
! 703: a = (double **)MALLOC(n*sizeof(double *));
! 704: for ( i = 0; i < n; i++ )
! 705: a[i] = (double *)MALLOC(n*sizeof(double));
! 706: return a;
! 707: }
! 708:
! 709: /*
! 710: * k <- (A(xi)-(sbeta-mn2/xi))f
! 711: * A(t) = (num[0]+num[1]t+...+num[d-1]*t^(d-1))/den(t)
! 712: */
! 713:
! 714: struct jv {
! 715: int j;
! 716: double v;
! 717: };
! 718:
! 719: struct smat {
! 720: int *rlen;
! 721: struct jv **row;
! 722: };
! 723:
! 724: void eval_pfaffian2(double *k,int n,int d,struct smat *num,P den,double xi,double *f)
! 725: {
! 726: struct smat ma;
! 727: struct jv *maj;
! 728: int i,j,l,s;
! 729: double t,dn;
! 730: P r;
! 731: Real u;
! 732:
! 733: memset(k,0,n*sizeof(double));
! 734: for ( i = d-1; i >= 0; i-- ) {
! 735: ma = num[i];
! 736: for ( j = 0; j < n; j++ ) {
! 737: maj = ma.row[j];
! 738: l = ma.rlen[j];
! 739: for ( t = 0, s = 0; s < l; s++, maj++ ) t += maj->v*f[maj->j];
! 740: k[j] = k[j]*xi+t;
! 741: }
! 742: }
! 743: MKReal(xi,u);
! 744: substp(CO,den,den->v,(P)u,&r); dn = ToReal(r);
! 745: for ( j = 0; j < n; j++ )
! 746: k[j] /= dn;
! 747: }
! 748:
! 749: void Prk_ratmat(NODE arg,LIST *rp)
! 750: {
! 751: VECT mat;
! 752: P den;
! 753: int ord;
! 754: double sbeta,x0,x1,xi,h,mn2,hd;
! 755: double a2,a3,a4,a5,a6;
! 756: double b21,b31,b32,b41,b42,b43,b51,b52,b53,b54,b61,b62,b63,b64,b65;
! 757: double c1,c2,c3,c4,c5,c6,c7;
! 758: VECT fv;
! 759: int step,j,i,k,n,d,len,s;
! 760: struct smat *num;
! 761: Obj **b;
! 762: MAT mati;
! 763: double *f,*w,*k1,*k2,*k3,*k4,*k5,*k6;
! 764: NODE nd,nd1;
! 765: Real x,t;
! 766: LIST l;
! 767:
! 768: ord = QTOS((Q)ARG0(arg));
! 769: mat = (VECT)ARG1(arg); den = (P)ARG2(arg);
! 770: x0 = ToReal((Num)ARG3(arg)); x1 = ToReal((Num)ARG4(arg));
! 771: step = QTOS((Q)ARG5(arg)); fv = (VECT)ARG6(arg);
! 772: h = (x1-x0)/step;
! 773:
! 774: n = fv->len;
! 775: d = mat->len;
! 776: num = (struct smat *)MALLOC(d*sizeof(struct smat));
! 777: for ( i = 0; i < d; i++ ) {
! 778: num[i].rlen = (int *)MALLOC_ATOMIC(n*sizeof(int));
! 779: num[i].row = (struct jv **)MALLOC(n*sizeof(struct jv *));
! 780: mati = (MAT)mat->body[i];
! 781: b = (Obj **)mati->body;
! 782: for ( j = 0; j < n; j++ ) {
! 783: for ( len = k = 0; k < n; k++ )
! 784: if ( b[j][k] ) len++;
! 785: num[i].rlen[j] = len;
! 786: if ( !len )
! 787: num[i].row[j] = 0;
! 788: else {
! 789: num[i].row[j] = (struct jv *)MALLOC_ATOMIC((len)*sizeof(struct jv));
! 790: for ( s = k = 0; k < n; k++ )
! 791: if ( b[j][k] ) {
! 792: num[i].row[j][s].j = k;
! 793: num[i].row[j][s].v = ToReal((Num)b[j][k]);
! 794: s++;
! 795: }
! 796: }
! 797: }
! 798: }
! 799: f = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 800: for ( j = 0; j < n; j++ )
! 801: f[j] = ToReal((Num)fv->body[j]);
! 802: w = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 803: k1 = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 804: k2 = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 805: k3 = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 806: k4 = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 807: k5 = (double *)MALLOC_ATOMIC(n*sizeof(double));
! 808: k6 = (double *)MALLOC(n*sizeof(double));
! 809: nd = 0;
! 810: switch ( ord ) {
! 811: case 4:
! 812: a2 = 1/2.0*h; b21 = 1/2.0*h;
! 813: a3 = 1/2.0*h; b31 = 0.0; b32 = 1/2.0*h;
! 814: a4 = 1.0*h; b41 = 0.0; b42 = 0.0; b43 = 1.0*h;
! 815: c1 = 1/6.0*h; c2 = 1/3.0*h; c3 = 1/3.0*h; c4 = 1/6.0*h;
! 816: for ( i = 0; i < step; i++ ) {
! 817: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
! 818: xi = x0+i*h;
! 819: eval_pfaffian2(k1,n,d,num,den,xi,f);
! 820: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
! 821: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
! 822: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
! 823: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
! 824: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
! 825: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
! 826: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j];
! 827: memcpy(f,w,n*sizeof(double));
! 828: MKReal(f[0],t);
! 829: MKReal(xi+h,x);
! 830: nd1 = mknode(2,x,t);
! 831: MKLIST(l,nd1);
! 832: MKNODE(nd1,l,nd);
! 833: nd = nd1;
! 834: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
! 835: }
! 836: MKLIST(*rp,nd);
! 837: break;
! 838: case 5:
! 839: default:
! 840: a2 = 1/4.0*h; b21 = 1/4.0*h;
! 841: a3 = 1/4.0*h; b31 = 1/8.0*h; b32 = 1/8.0*h;
! 842: a4 = 1/2.0*h; b41 = 0.0; b42 = 0.0; b43 = 1/2.0*h;
! 843: a5 = 3/4.0*h; b51 = 3/16.0*h;b52 = -3/8.0*h; b53 = 3/8.0*h; b54 = 9/16.0*h;
! 844: a6 = 1.0*h; b61 = -3/7.0*h;b62 = 8/7.0*h; b63 = 6/7.0*h; b64 = -12/7.0*h; b65 = 8/7.0*h;
! 845: c1 = 7/90.0*h; c2 = 0.0; c3 = 16/45.0*h; c4 = 2/15.0*h; c5 = 16/45.0*h; c6 = 7/90.0*h;
! 846: for ( i = 0; i < step; i++ ) {
! 847: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
! 848: xi = x0+i*h;
! 849: eval_pfaffian2(k1,n,d,num,den,xi,f);
! 850: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
! 851: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
! 852: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
! 853: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
! 854: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
! 855: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
! 856: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b51*k1[j]+b52*k2[j]+b53*k3[j]+b54*k4[j];
! 857: eval_pfaffian2(k5,n,d,num,den,xi+a5,w);
! 858: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b61*k1[j]+b62*k2[j]+b63*k3[j]+b64*k4[j]+b65*k5[j];
! 859: eval_pfaffian2(k6,n,d,num,den,xi+a6,w);
! 860: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j]+c5*k5[j]+c6*k6[j];
! 861: memcpy(f,w,n*sizeof(double));
! 862: MKReal(f[0],t);
! 863: MKReal(xi+h,x);
! 864: nd1 = mknode(2,x,t);
! 865: MKLIST(l,nd1);
! 866: MKNODE(nd1,l,nd);
! 867: nd = nd1;
! 868: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
! 869: }
! 870: MKLIST(*rp,nd);
! 871: break;
! 872: }
! 873: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bfaux.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin