Annotation of OpenXM_contrib/gnuplot/standard.c, Revision 1.1
1.1 ! maekawa 1: #ifndef lint
! 2: static char *RCSid = "$Id: standard.c,v 1.24 1998/04/14 00:16:21 drd Exp $";
! 3: #endif
! 4:
! 5: /* GNUPLOT - standard.c */
! 6:
! 7: /*[
! 8: * Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
! 9: *
! 10: * Permission to use, copy, and distribute this software and its
! 11: * documentation for any purpose with or without fee is hereby granted,
! 12: * provided that the above copyright notice appear in all copies and
! 13: * that both that copyright notice and this permission notice appear
! 14: * in supporting documentation.
! 15: *
! 16: * Permission to modify the software is granted, but not the right to
! 17: * distribute the complete modified source code. Modifications are to
! 18: * be distributed as patches to the released version. Permission to
! 19: * distribute binaries produced by compiling modified sources is granted,
! 20: * provided you
! 21: * 1. distribute the corresponding source modifications from the
! 22: * released version in the form of a patch file along with the binaries,
! 23: * 2. add special version identification to distinguish your version
! 24: * in addition to the base release version number,
! 25: * 3. provide your name and address as the primary contact for the
! 26: * support of your modified version, and
! 27: * 4. retain our contact information in regard to use of the base
! 28: * software.
! 29: * Permission to distribute the released version of the source code along
! 30: * with corresponding source modifications in the form of a patch file is
! 31: * granted with same provisions 2 through 4 for binary distributions.
! 32: *
! 33: * This software is provided "as is" without express or implied warranty
! 34: * to the extent permitted by applicable law.
! 35: ]*/
! 36:
! 37: #include "plot.h"
! 38: #include "setshow.h" /* for ang2rad */
! 39: #include "fnproto.h"
! 40:
! 41: extern struct value stack[STACK_DEPTH];
! 42: extern int s_p;
! 43: extern double zero;
! 44:
! 45: static double jzero __PROTO((double x));
! 46: static double pzero __PROTO((double x));
! 47: static double qzero __PROTO((double x));
! 48: static double yzero __PROTO((double x));
! 49: static double rj0 __PROTO((double x));
! 50: static double ry0 __PROTO((double x));
! 51: static double jone __PROTO((double x));
! 52: static double pone __PROTO((double x));
! 53: static double qone __PROTO((double x));
! 54: static double yone __PROTO((double x));
! 55: static double rj1 __PROTO((double x));
! 56: static double ry1 __PROTO((double x));
! 57:
! 58: /* The bessel function approximations here are from
! 59: * "Computer Approximations"
! 60: * by Hart, Cheney et al.
! 61: * John Wiley & Sons, 1968
! 62: */
! 63:
! 64: /* There appears to be a mistake in Hart, Cheney et al. on page 149.
! 65: * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
! 66: * Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
! 67: * In the functions below, Qn(x) is implementated using the later
! 68: * equation.
! 69: * These bessel functions are accurate to about 1e-13
! 70: */
! 71:
! 72: #if (defined (ATARI) || defined (MTOS)) && defined(__PUREC__)
! 73: /* Sorry. But PUREC bugs here.
! 74: * These bessel functions are NOT accurate to about 1e-13
! 75: */
! 76:
! 77: #define PI_ON_FOUR 0.785398163397448309615661
! 78: #define PI_ON_TWO 1.570796326794896619231313
! 79: #define THREE_PI_ON_FOUR 2.356194490192344928846982
! 80: #define TWO_ON_PI 0.636619772367581343075535
! 81:
! 82: static double dzero = 0.0;
! 83:
! 84: /* jzero for x in [0,8]
! 85: * Index 5849, 19.22 digits precision
! 86: */
! 87: static double pjzero[] = {
! 88: 0.493378725179413356181681e+21,
! 89: -0.117915762910761053603844e+21,
! 90: 0.638205934107235656228943e+19,
! 91: -0.136762035308817138686542e+18,
! 92: 0.143435493914034611166432e+16,
! 93: -0.808522203485379387119947e+13,
! 94: 0.250715828553688194555516e+11,
! 95: -0.405041237183313270636066e+8,
! 96: 0.268578685698001498141585e+5
! 97: };
! 98:
! 99: static double qjzero[] = {
! 100: 0.493378725179413356211328e+21,
! 101: 0.542891838409228516020019e+19,
! 102: 0.302463561670946269862733e+17,
! 103: 0.112775673967979850705603e+15,
! 104: 0.312304311494121317257247e+12,
! 105: 0.669998767298223967181403e+9,
! 106: 0.111463609846298537818240e+7,
! 107: 0.136306365232897060444281e+4,
! 108: 0.1e+1
! 109: };
! 110:
! 111: /* pzero for x in [8,inf]
! 112: * Index 6548, 18.16 digits precision
! 113: */
! 114: static double ppzero[] = {
! 115: 0.227790901973046843022700e+5,
! 116: 0.413453866395807657967802e+5,
! 117: 0.211705233808649443219340e+5,
! 118: 0.348064864432492703474453e+4,
! 119: 0.153762019090083542957717e+3,
! 120: 0.889615484242104552360748e+0
! 121: };
! 122:
! 123: static double qpzero[] = {
! 124: 0.227790901973046843176842e+5,
! 125: 0.413704124955104166398920e+5,
! 126: 0.212153505618801157304226e+5,
! 127: 0.350287351382356082073561e+4,
! 128: 0.157111598580808936490685e+3,
! 129: 0.1e+1
! 130: };
! 131:
! 132: /* qzero for x in [8,inf]
! 133: * Index 6948, 18.33 digits precision
! 134: */
! 135: static double pqzero[] = {
! 136: -0.892266002008000940984692e+2,
! 137: -0.185919536443429938002522e+3,
! 138: -0.111834299204827376112621e+3,
! 139: -0.223002616662141984716992e+2,
! 140: -0.124410267458356384591379e+1,
! 141: -0.8803330304868075181663e-2,
! 142: };
! 143:
! 144: static double qqzero[] = {
! 145: 0.571050241285120619052476e+4,
! 146: 0.119511315434346136469526e+5,
! 147: 0.726427801692110188369134e+4,
! 148: 0.148872312322837565816135e+4,
! 149: 0.905937695949931258588188e+2,
! 150: 0.1e+1
! 151: };
! 152:
! 153:
! 154: /* yzero for x in [0,8]
! 155: * Index 6245, 18.78 digits precision
! 156: */
! 157: static double pyzero[] = {
! 158: -0.275028667862910958370193e+20,
! 159: 0.658747327571955492599940e+20,
! 160: -0.524706558111276494129735e+19,
! 161: 0.137562431639934407857134e+18,
! 162: -0.164860581718572947312208e+16,
! 163: 0.102552085968639428450917e+14,
! 164: -0.343637122297904037817103e+11,
! 165: 0.591521346568688965427383e+8,
! 166: -0.413703549793314855412524e+5
! 167: };
! 168:
! 169: static double qyzero[] = {
! 170: 0.372645883898616588198998e+21,
! 171: 0.419241704341083997390477e+19,
! 172: 0.239288304349978185743936e+17,
! 173: 0.916203803407518526248915e+14,
! 174: 0.261306575504108124956848e+12,
! 175: 0.579512264070072953738009e+9,
! 176: 0.100170264128890626566665e+7,
! 177: 0.128245277247899380417633e+4,
! 178: 0.1e+1
! 179: };
! 180:
! 181:
! 182: /* jone for x in [0,8]
! 183: * Index 6050, 20.98 digits precision
! 184: */
! 185: static double pjone[] = {
! 186: 0.581199354001606143928051e+21,
! 187: -0.667210656892491629802094e+20,
! 188: 0.231643358063400229793182e+19,
! 189: -0.358881756991010605074364e+17,
! 190: 0.290879526383477540973760e+15,
! 191: -0.132298348033212645312547e+13,
! 192: 0.341323418230170053909129e+10,
! 193: -0.469575353064299585976716e+7,
! 194: 0.270112271089232341485679e+4
! 195: };
! 196:
! 197: static double qjone[] = {
! 198: 0.116239870800321228785853e+22,
! 199: 0.118577071219032099983711e+20,
! 200: 0.609206139891752174610520e+17,
! 201: 0.208166122130760735124018e+15,
! 202: 0.524371026216764971540673e+12,
! 203: 0.101386351435867398996705e+10,
! 204: 0.150179359499858550592110e+7,
! 205: 0.160693157348148780197092e+4,
! 206: 0.1e+1
! 207: };
! 208:
! 209:
! 210: /* pone for x in [8,inf]
! 211: * Index 6749, 18.11 digits precision
! 212: */
! 213: static double ppone[] = {
! 214: 0.352246649133679798341724e+5,
! 215: 0.627588452471612812690057e+5,
! 216: 0.313539631109159574238670e+5,
! 217: 0.498548320605943384345005e+4,
! 218: 0.211152918285396238210572e+3,
! 219: 0.12571716929145341558495e+1
! 220: };
! 221:
! 222: static double qpone[] = {
! 223: 0.352246649133679798068390e+5,
! 224: 0.626943469593560511888834e+5,
! 225: 0.312404063819041039923016e+5,
! 226: 0.493039649018108897938610e+4,
! 227: 0.203077518913475932229357e+3,
! 228: 0.1e+1
! 229: };
! 230:
! 231: /* qone for x in [8,inf]
! 232: * Index 7149, 18.28 digits precision
! 233: */
! 234: static double pqone[] = {
! 235: 0.351175191430355282253332e+3,
! 236: 0.721039180490447503928086e+3,
! 237: 0.425987301165444238988699e+3,
! 238: 0.831898957673850827325226e+2,
! 239: 0.45681716295512267064405e+1,
! 240: 0.3532840052740123642735e-1
! 241: };
! 242:
! 243: static double qqone[] = {
! 244: 0.749173741718091277145195e+4,
! 245: 0.154141773392650970499848e+5,
! 246: 0.915223170151699227059047e+4,
! 247: 0.181118670055235135067242e+4,
! 248: 0.103818758546213372877664e+3,
! 249: 0.1e+1
! 250: };
! 251:
! 252:
! 253: /* yone for x in [0,8]
! 254: * Index 6444, 18.24 digits precision
! 255: */
! 256: static double pyone[] = {
! 257: -0.292382196153296254310105e+20,
! 258: 0.774852068218683964508809e+19,
! 259: -0.344104806308411444618546e+18,
! 260: 0.591516076049007061849632e+16,
! 261: -0.486331694256717507482813e+14,
! 262: 0.204969667374566218261980e+12,
! 263: -0.428947196885524880182182e+9,
! 264: 0.355692400983052605669132e+6
! 265: };
! 266:
! 267: static double qyone[] = {
! 268: 0.149131151130292035017408e+21,
! 269: 0.181866284170613498688507e+19,
! 270: 0.113163938269888452690508e+17,
! 271: 0.475517358888813771309277e+14,
! 272: 0.150022169915670898716637e+12,
! 273: 0.371666079862193028559693e+9,
! 274: 0.726914730719888456980191e+6,
! 275: 0.107269614377892552332213e+4,
! 276: 0.1e+1
! 277: };
! 278:
! 279: #else
! 280:
! 281: #define PI_ON_FOUR 0.78539816339744830961566084581987572
! 282: #define PI_ON_TWO 1.57079632679489661923131269163975144
! 283: #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
! 284: #define TWO_ON_PI 0.63661977236758134307553505349005744
! 285:
! 286: static double dzero = 0.0;
! 287:
! 288: /* jzero for x in [0,8]
! 289: * Index 5849, 19.22 digits precision
! 290: */
! 291: static double GPFAR pjzero[] = {
! 292: 0.4933787251794133561816813446e+21,
! 293: -0.11791576291076105360384408e+21,
! 294: 0.6382059341072356562289432465e+19,
! 295: -0.1367620353088171386865416609e+18,
! 296: 0.1434354939140346111664316553e+16,
! 297: -0.8085222034853793871199468171e+13,
! 298: 0.2507158285536881945555156435e+11,
! 299: -0.4050412371833132706360663322e+8,
! 300: 0.2685786856980014981415848441e+5
! 301: };
! 302:
! 303: static double GPFAR qjzero[] = {
! 304: 0.4933787251794133562113278438e+21,
! 305: 0.5428918384092285160200195092e+19,
! 306: 0.3024635616709462698627330784e+17,
! 307: 0.1127756739679798507056031594e+15,
! 308: 0.3123043114941213172572469442e+12,
! 309: 0.669998767298223967181402866e+9,
! 310: 0.1114636098462985378182402543e+7,
! 311: 0.1363063652328970604442810507e+4,
! 312: 0.1e+1
! 313: };
! 314:
! 315: /* pzero for x in [8,inf]
! 316: * Index 6548, 18.16 digits precision
! 317: */
! 318: static double GPFAR ppzero[] = {
! 319: 0.2277909019730468430227002627e+5,
! 320: 0.4134538663958076579678016384e+5,
! 321: 0.2117052338086494432193395727e+5,
! 322: 0.348064864432492703474453111e+4,
! 323: 0.15376201909008354295771715e+3,
! 324: 0.889615484242104552360748e+0
! 325: };
! 326:
! 327: static double GPFAR qpzero[] = {
! 328: 0.2277909019730468431768423768e+5,
! 329: 0.4137041249551041663989198384e+5,
! 330: 0.2121535056188011573042256764e+5,
! 331: 0.350287351382356082073561423e+4,
! 332: 0.15711159858080893649068482e+3,
! 333: 0.1e+1
! 334: };
! 335:
! 336: /* qzero for x in [8,inf]
! 337: * Index 6948, 18.33 digits precision
! 338: */
! 339: static double GPFAR pqzero[] = {
! 340: -0.8922660020080009409846916e+2,
! 341: -0.18591953644342993800252169e+3,
! 342: -0.11183429920482737611262123e+3,
! 343: -0.2230026166621419847169915e+2,
! 344: -0.124410267458356384591379e+1,
! 345: -0.8803330304868075181663e-2,
! 346: };
! 347:
! 348: static double GPFAR qqzero[] = {
! 349: 0.571050241285120619052476459e+4,
! 350: 0.1195113154343461364695265329e+5,
! 351: 0.726427801692110188369134506e+4,
! 352: 0.148872312322837565816134698e+4,
! 353: 0.9059376959499312585881878e+2,
! 354: 0.1e+1
! 355: };
! 356:
! 357:
! 358: /* yzero for x in [0,8]
! 359: * Index 6245, 18.78 digits precision
! 360: */
! 361: static double GPFAR pyzero[] = {
! 362: -0.2750286678629109583701933175e+20,
! 363: 0.6587473275719554925999402049e+20,
! 364: -0.5247065581112764941297350814e+19,
! 365: 0.1375624316399344078571335453e+18,
! 366: -0.1648605817185729473122082537e+16,
! 367: 0.1025520859686394284509167421e+14,
! 368: -0.3436371222979040378171030138e+11,
! 369: 0.5915213465686889654273830069e+8,
! 370: -0.4137035497933148554125235152e+5
! 371: };
! 372:
! 373: static double GPFAR qyzero[] = {
! 374: 0.3726458838986165881989980739e+21,
! 375: 0.4192417043410839973904769661e+19,
! 376: 0.2392883043499781857439356652e+17,
! 377: 0.9162038034075185262489147968e+14,
! 378: 0.2613065755041081249568482092e+12,
! 379: 0.5795122640700729537380087915e+9,
! 380: 0.1001702641288906265666651753e+7,
! 381: 0.1282452772478993804176329391e+4,
! 382: 0.1e+1
! 383: };
! 384:
! 385:
! 386: /* jone for x in [0,8]
! 387: * Index 6050, 20.98 digits precision
! 388: */
! 389: static double GPFAR pjone[] = {
! 390: 0.581199354001606143928050809e+21,
! 391: -0.6672106568924916298020941484e+20,
! 392: 0.2316433580634002297931815435e+19,
! 393: -0.3588817569910106050743641413e+17,
! 394: 0.2908795263834775409737601689e+15,
! 395: -0.1322983480332126453125473247e+13,
! 396: 0.3413234182301700539091292655e+10,
! 397: -0.4695753530642995859767162166e+7,
! 398: 0.270112271089232341485679099e+4
! 399: };
! 400:
! 401: static double GPFAR qjone[] = {
! 402: 0.11623987080032122878585294e+22,
! 403: 0.1185770712190320999837113348e+20,
! 404: 0.6092061398917521746105196863e+17,
! 405: 0.2081661221307607351240184229e+15,
! 406: 0.5243710262167649715406728642e+12,
! 407: 0.1013863514358673989967045588e+10,
! 408: 0.1501793594998585505921097578e+7,
! 409: 0.1606931573481487801970916749e+4,
! 410: 0.1e+1
! 411: };
! 412:
! 413:
! 414: /* pone for x in [8,inf]
! 415: * Index 6749, 18.11 digits precision
! 416: */
! 417: static double GPFAR ppone[] = {
! 418: 0.352246649133679798341724373e+5,
! 419: 0.62758845247161281269005675e+5,
! 420: 0.313539631109159574238669888e+5,
! 421: 0.49854832060594338434500455e+4,
! 422: 0.2111529182853962382105718e+3,
! 423: 0.12571716929145341558495e+1
! 424: };
! 425:
! 426: static double GPFAR qpone[] = {
! 427: 0.352246649133679798068390431e+5,
! 428: 0.626943469593560511888833731e+5,
! 429: 0.312404063819041039923015703e+5,
! 430: 0.4930396490181088979386097e+4,
! 431: 0.2030775189134759322293574e+3,
! 432: 0.1e+1
! 433: };
! 434:
! 435: /* qone for x in [8,inf]
! 436: * Index 7149, 18.28 digits precision
! 437: */
! 438: static double GPFAR pqone[] = {
! 439: 0.3511751914303552822533318e+3,
! 440: 0.7210391804904475039280863e+3,
! 441: 0.4259873011654442389886993e+3,
! 442: 0.831898957673850827325226e+2,
! 443: 0.45681716295512267064405e+1,
! 444: 0.3532840052740123642735e-1
! 445: };
! 446:
! 447: static double GPFAR qqone[] = {
! 448: 0.74917374171809127714519505e+4,
! 449: 0.154141773392650970499848051e+5,
! 450: 0.91522317015169922705904727e+4,
! 451: 0.18111867005523513506724158e+4,
! 452: 0.1038187585462133728776636e+3,
! 453: 0.1e+1
! 454: };
! 455:
! 456:
! 457: /* yone for x in [0,8]
! 458: * Index 6444, 18.24 digits precision
! 459: */
! 460: static double GPFAR pyone[] = {
! 461: -0.2923821961532962543101048748e+20,
! 462: 0.7748520682186839645088094202e+19,
! 463: -0.3441048063084114446185461344e+18,
! 464: 0.5915160760490070618496315281e+16,
! 465: -0.4863316942567175074828129117e+14,
! 466: 0.2049696673745662182619800495e+12,
! 467: -0.4289471968855248801821819588e+9,
! 468: 0.3556924009830526056691325215e+6
! 469: };
! 470:
! 471: static double GPFAR qyone[] = {
! 472: 0.1491311511302920350174081355e+21,
! 473: 0.1818662841706134986885065935e+19,
! 474: 0.113163938269888452690508283e+17,
! 475: 0.4755173588888137713092774006e+14,
! 476: 0.1500221699156708987166369115e+12,
! 477: 0.3716660798621930285596927703e+9,
! 478: 0.726914730719888456980191315e+6,
! 479: 0.10726961437789255233221267e+4,
! 480: 0.1e+1
! 481: };
! 482:
! 483: #endif /* (ATARI || MTOS) && __PUREC__ */
! 484:
! 485: void f_real()
! 486: {
! 487: struct value a;
! 488: push( Gcomplex(&a,real(pop(&a)), 0.0) );
! 489: }
! 490:
! 491: void f_imag()
! 492: {
! 493: struct value a;
! 494: push( Gcomplex(&a,imag(pop(&a)), 0.0) );
! 495: }
! 496:
! 497:
! 498: /* ang2rad is 1 if we are in radians, or pi/180 if we are in degrees */
! 499:
! 500: void f_arg()
! 501: {
! 502: struct value a;
! 503: push( Gcomplex(&a,angle(pop(&a))/ang2rad, 0.0) );
! 504: }
! 505:
! 506: void f_conjg()
! 507: {
! 508: struct value a;
! 509: (void) pop(&a);
! 510: push( Gcomplex(&a,real(&a),-imag(&a) ));
! 511: }
! 512:
! 513: /* ang2rad is 1 if we are in radians, or pi/180 if we are in degrees */
! 514:
! 515: void f_sin()
! 516: {
! 517: struct value a;
! 518: (void) pop(&a);
! 519: push( Gcomplex(&a,sin(ang2rad*real(&a))*cosh(ang2rad*imag(&a)), cos(ang2rad*real(&a))*sinh(ang2rad*imag(&a))) );
! 520: }
! 521:
! 522: void f_cos()
! 523: {
! 524: struct value a;
! 525: (void) pop(&a);
! 526: push( Gcomplex(&a,cos(ang2rad*real(&a))*cosh(ang2rad*imag(&a)), -sin(ang2rad*real(&a))*sinh(ang2rad*imag(&a))));
! 527: }
! 528:
! 529: void f_tan()
! 530: {
! 531: struct value a;
! 532: register double den;
! 533: (void) pop(&a);
! 534: if (imag(&a) == 0.0)
! 535: push( Gcomplex(&a,tan(ang2rad*real(&a)),0.0) );
! 536: else {
! 537: den = cos(2*ang2rad*real(&a))+cosh(2*ang2rad*imag(&a));
! 538: if (den == 0.0) {
! 539: undefined = TRUE;
! 540: push( &a );
! 541: }
! 542: else
! 543: push( Gcomplex(&a,sin(2*ang2rad*real(&a))/den, sinh(2*ang2rad*imag(&a))/den) );
! 544: }
! 545: }
! 546:
! 547: void f_asin()
! 548: {
! 549: struct value a;
! 550: register double alpha, beta, x, y;
! 551: (void) pop(&a);
! 552: x = real(&a); y = imag(&a);
! 553: if (y == 0.0 && fabs(x) <= 1.0) {
! 554: push( Gcomplex(&a,asin(x)/ang2rad,0.0) );
! 555: } else if (x == 0.0) {
! 556: push( Gcomplex(&a, 0.0, -log(-y+sqrt(y*y+1))/ang2rad) );
! 557: } else {
! 558: beta = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
! 559: if (beta > 1) beta = 1; /* Avoid rounding error problems */
! 560: alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
! 561: push( Gcomplex(&a,asin(beta)/ang2rad, -log(alpha + sqrt(alpha*alpha-1))/ang2rad) );
! 562: }
! 563: }
! 564:
! 565: void f_acos()
! 566: {
! 567: struct value a;
! 568: register double x, y;
! 569: (void) pop(&a);
! 570: x = real(&a); y = imag(&a);
! 571: if (y == 0.0 && fabs(x) <= 1.0) {
! 572: /* real result */
! 573: push( Gcomplex(&a,acos(x)/ang2rad,0.0) );
! 574: } else {
! 575: double alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
! 576: double beta = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
! 577: if (beta > 1)
! 578: beta = 1; /* Avoid rounding error problems */
! 579: else if (beta < -1)
! 580: beta = -1;
! 581: push( Gcomplex(&a,acos(beta)/ang2rad, log(alpha + sqrt(alpha*alpha-1))/ang2rad));
! 582: }
! 583: }
! 584:
! 585: void f_atan()
! 586: {
! 587: struct value a;
! 588: register double x, y, u, v, w, z;
! 589: (void) pop(&a);
! 590: x = real(&a); y = imag(&a);
! 591: if (y == 0.0)
! 592: push( Gcomplex(&a,atan(x)/ang2rad, 0.0) );
! 593: else if (x == 0.0 && fabs(y) >= 1.0) {
! 594: undefined = TRUE;
! 595: push(Gcomplex(&a,0.0, 0.0));
! 596: } else {
! 597: if (x >= 0) {
! 598: u = x;
! 599: v = y;
! 600: } else {
! 601: u = -x;
! 602: v = -y;
! 603: }
! 604:
! 605: z = atan(2*u/(1-u*u-v*v));
! 606: w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
! 607: if (z < 0)
! 608: z = z + 2*PI_ON_TWO;
! 609: if (x < 0) {
! 610: z = -z;
! 611: w = -w;
! 612: }
! 613: push( Gcomplex(&a,0.5*z/ang2rad, w) );
! 614: }
! 615: }
! 616:
! 617: /* real parts only */
! 618: void f_atan2()
! 619: {
! 620: struct value a;
! 621: double x;
! 622: double y;
! 623:
! 624: x = real(pop(&a));
! 625: y = real(pop(&a));
! 626:
! 627: if (x == 0.0 && y == 0.0) {
! 628: undefined = TRUE;
! 629: push(Ginteger(&a,0));
! 630: }
! 631: push(Gcomplex(&a,atan2(y,x)/ang2rad,0.0));
! 632: }
! 633:
! 634:
! 635: void f_sinh()
! 636: {
! 637: struct value a;
! 638: (void) pop(&a);
! 639: push( Gcomplex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
! 640: }
! 641:
! 642: void f_cosh()
! 643: {
! 644: struct value a;
! 645: (void) pop(&a);
! 646: push( Gcomplex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
! 647: }
! 648:
! 649: void f_tanh()
! 650: {
! 651: struct value a;
! 652: register double den;
! 653: (void) pop(&a);
! 654: den = cosh(2*real(&a)) + cos(2*imag(&a));
! 655: push( Gcomplex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
! 656: }
! 657:
! 658: void f_asinh()
! 659: {
! 660: struct value a; /* asinh(z) = -I*asin(I*z) */
! 661: register double alpha, beta, x, y;
! 662: (void) pop(&a);
! 663: x = -imag(&a); y = real(&a);
! 664: if (y == 0.0 && fabs(x) <= 1.0) {
! 665: push( Gcomplex(&a, 0.0, -asin(x)/ang2rad) );
! 666: } else if (y == 0.0) {
! 667: push( Gcomplex(&a, 0.0, 0.0) );
! 668: undefined = TRUE;
! 669: } else if (x == 0.0) {
! 670: push( Gcomplex(&a, log(y+sqrt(y*y+1))/ang2rad, 0.0) );
! 671: } else {
! 672: beta = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
! 673: alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
! 674: push( Gcomplex(&a, log(alpha + sqrt(alpha*alpha-1))/ang2rad, -asin(beta)/ang2rad) );
! 675: }
! 676: }
! 677:
! 678: void f_acosh()
! 679: {
! 680: struct value a;
! 681: register double alpha, beta, x, y; /* acosh(z) = I*acos(z) */
! 682: (void) pop(&a);
! 683: x = real(&a); y = imag(&a);
! 684: if (y == 0.0 && fabs(x) <= 1.0) {
! 685: push( Gcomplex(&a, 0.0, acos(x)/ang2rad) );
! 686: } else if (y == 0) {
! 687: push( Gcomplex(&a, log(x+sqrt(x*x-1))/ang2rad, 0.0) );
! 688: } else {
! 689: alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
! 690: beta = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
! 691: push( Gcomplex(&a, log(alpha + sqrt(alpha*alpha-1))/ang2rad, acos(beta)/ang2rad));
! 692: }
! 693: }
! 694:
! 695: void f_atanh()
! 696: {
! 697: struct value a;
! 698: register double x, y, u, v, w, z; /* atanh(z) = -I*atan(I*z) */
! 699: (void) pop(&a);
! 700: x = -imag(&a); y = real(&a);
! 701: if (y == 0.0)
! 702: push( Gcomplex(&a, 0.0, -atan(x)/ang2rad) );
! 703: else if (x == 0.0 && fabs(y) >= 1.0) {
! 704: undefined = TRUE;
! 705: push(Gcomplex(&a,0.0, 0.0));
! 706: } else {
! 707: if (x >= 0) {
! 708: u = x;
! 709: v = y;
! 710: } else {
! 711: u = -x;
! 712: v = -y;
! 713: }
! 714:
! 715: z = atan(2*u/(1-u*u-v*v));
! 716: w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
! 717: if (z < 0)
! 718: z = z + 2*PI_ON_TWO;
! 719: if (x < 0) {
! 720: z = -z;
! 721: w = -w;
! 722: }
! 723: push( Gcomplex(&a, w, -0.5*z/ang2rad) );
! 724: }
! 725: }
! 726:
! 727: void f_int()
! 728: {
! 729: struct value a;
! 730: push( Ginteger(&a,(int)real(pop(&a))) );
! 731: }
! 732:
! 733:
! 734: void f_abs()
! 735: {
! 736: struct value a;
! 737: (void) pop(&a);
! 738: switch (a.type) {
! 739: case INTGR:
! 740: push( Ginteger(&a,abs(a.v.int_val)) );
! 741: break;
! 742: case CMPLX:
! 743: push( Gcomplex(&a,magnitude(&a), 0.0) );
! 744: }
! 745: }
! 746:
! 747: void f_sgn()
! 748: {
! 749: struct value a;
! 750: (void) pop(&a);
! 751: switch(a.type) {
! 752: case INTGR:
! 753: push( Ginteger(&a,(a.v.int_val > 0) ? 1 :
! 754: (a.v.int_val < 0) ? -1 : 0) );
! 755: break;
! 756: case CMPLX:
! 757: push( Ginteger(&a,(a.v.cmplx_val.real > 0.0) ? 1 :
! 758: (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
! 759: break;
! 760: }
! 761: }
! 762:
! 763:
! 764: void f_sqrt()
! 765: {
! 766: struct value a;
! 767: register double mag;
! 768: (void) pop(&a);
! 769: mag = sqrt(magnitude(&a));
! 770: if (imag(&a) == 0.0) {
! 771: if (real(&a) < 0.0)
! 772: push( Gcomplex(&a,0.0,mag) );
! 773: else
! 774: push( Gcomplex(&a,mag, 0.0) );
! 775: } else {
! 776: /* -pi < ang < pi, so real(sqrt(z)) >= 0 */
! 777: double ang = angle(&a) / 2.0;
! 778: push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
! 779: }
! 780: }
! 781:
! 782:
! 783: void f_exp()
! 784: {
! 785: struct value a;
! 786: register double mag, ang;
! 787: (void) pop(&a);
! 788: mag = gp_exp(real(&a));
! 789: ang = imag(&a);
! 790: push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
! 791: }
! 792:
! 793:
! 794: void f_log10()
! 795: {
! 796: struct value a;
! 797: (void) pop(&a);
! 798: push( Gcomplex(&a,log(magnitude(&a))/M_LN10, angle(&a)/M_LN10) );
! 799: }
! 800:
! 801:
! 802: void f_log()
! 803: {
! 804: struct value a;
! 805: (void) pop(&a);
! 806: push( Gcomplex(&a,log(magnitude(&a)), angle(&a)) );
! 807: }
! 808:
! 809:
! 810: void f_floor()
! 811: {
! 812: struct value a;
! 813:
! 814: (void) pop(&a);
! 815: switch (a.type) {
! 816: case INTGR:
! 817: push( Ginteger(&a,(int)floor((double)a.v.int_val)));
! 818: break;
! 819: case CMPLX:
! 820: push( Ginteger(&a,(int)floor(a.v.cmplx_val.real)));
! 821: }
! 822: }
! 823:
! 824:
! 825: void f_ceil()
! 826: {
! 827: struct value a;
! 828:
! 829: (void) pop(&a);
! 830: switch (a.type) {
! 831: case INTGR:
! 832: push( Ginteger(&a,(int)ceil((double)a.v.int_val)));
! 833: break;
! 834: case CMPLX:
! 835: push( Ginteger(&a,(int)ceil(a.v.cmplx_val.real)));
! 836: }
! 837: }
! 838:
! 839: /* bessel function approximations */
! 840: static double jzero(x)
! 841: double x;
! 842: {
! 843: double p, q, x2;
! 844: int n;
! 845:
! 846: x2 = x * x;
! 847: p = pjzero[8];
! 848: q = qjzero[8];
! 849: for (n=7; n>=0; n--) {
! 850: p = p*x2 + pjzero[n];
! 851: q = q*x2 + qjzero[n];
! 852: }
! 853: return(p/q);
! 854: }
! 855:
! 856: static double pzero(x)
! 857: double x;
! 858: {
! 859: double p, q, z, z2;
! 860: int n;
! 861:
! 862: z = 8.0 / x;
! 863: z2 = z * z;
! 864: p = ppzero[5];
! 865: q = qpzero[5];
! 866: for (n=4; n>=0; n--) {
! 867: p = p*z2 + ppzero[n];
! 868: q = q*z2 + qpzero[n];
! 869: }
! 870: return(p/q);
! 871: }
! 872:
! 873: static double qzero(x)
! 874: double x;
! 875: {
! 876: double p, q, z, z2;
! 877: int n;
! 878:
! 879: z = 8.0 / x;
! 880: z2 = z * z;
! 881: p = pqzero[5];
! 882: q = qqzero[5];
! 883: for (n=4; n>=0; n--) {
! 884: p = p*z2 + pqzero[n];
! 885: q = q*z2 + qqzero[n];
! 886: }
! 887: return(p/q);
! 888: }
! 889:
! 890: static double yzero(x)
! 891: double x;
! 892: {
! 893: double p, q, x2;
! 894: int n;
! 895:
! 896: x2 = x * x;
! 897: p = pyzero[8];
! 898: q = qyzero[8];
! 899: for (n=7; n>=0; n--) {
! 900: p = p*x2 + pyzero[n];
! 901: q = q*x2 + qyzero[n];
! 902: }
! 903: return(p/q);
! 904: }
! 905:
! 906: static double rj0(x)
! 907: double x;
! 908: {
! 909: if ( x <= 0.0 )
! 910: x = -x;
! 911: if ( x < 8.0 )
! 912: return(jzero(x));
! 913: else
! 914: return( sqrt(TWO_ON_PI/x) *
! 915: (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
! 916:
! 917: }
! 918:
! 919: static double ry0(x)
! 920: double x;
! 921: {
! 922: if ( x < 0.0 )
! 923: return(dzero/dzero); /* error */
! 924: if ( x < 8.0 )
! 925: return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
! 926: else
! 927: return( sqrt(TWO_ON_PI/x) *
! 928: (pzero(x)*sin(x-PI_ON_FOUR) +
! 929: (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
! 930:
! 931: }
! 932:
! 933:
! 934: static double jone(x)
! 935: double x;
! 936: {
! 937: double p, q, x2;
! 938: int n;
! 939:
! 940: x2 = x * x;
! 941: p = pjone[8];
! 942: q = qjone[8];
! 943: for (n=7; n>=0; n--) {
! 944: p = p*x2 + pjone[n];
! 945: q = q*x2 + qjone[n];
! 946: }
! 947: return(p/q);
! 948: }
! 949:
! 950: static double pone(x)
! 951: double x;
! 952: {
! 953: double p, q, z, z2;
! 954: int n;
! 955:
! 956: z = 8.0 / x;
! 957: z2 = z * z;
! 958: p = ppone[5];
! 959: q = qpone[5];
! 960: for (n=4; n>=0; n--) {
! 961: p = p*z2 + ppone[n];
! 962: q = q*z2 + qpone[n];
! 963: }
! 964: return(p/q);
! 965: }
! 966:
! 967: static double qone(x)
! 968: double x;
! 969: {
! 970: double p, q, z, z2;
! 971: int n;
! 972:
! 973: z = 8.0 / x;
! 974: z2 = z * z;
! 975: p = pqone[5];
! 976: q = qqone[5];
! 977: for (n=4; n>=0; n--) {
! 978: p = p*z2 + pqone[n];
! 979: q = q*z2 + qqone[n];
! 980: }
! 981: return(p/q);
! 982: }
! 983:
! 984: static double yone(x)
! 985: double x;
! 986: {
! 987: double p, q, x2;
! 988: int n;
! 989:
! 990: x2 = x * x;
! 991: p = 0.0;
! 992: q = qyone[8];
! 993: for (n=7; n>=0; n--) {
! 994: p = p*x2 + pyone[n];
! 995: q = q*x2 + qyone[n];
! 996: }
! 997: return(p/q);
! 998: }
! 999:
! 1000: static double rj1(x)
! 1001: double x;
! 1002: {
! 1003: double v,w;
! 1004: v = x;
! 1005: if ( x < 0.0 )
! 1006: x = -x;
! 1007: if ( x < 8.0 )
! 1008: return(v*jone(x));
! 1009: else {
! 1010: w = sqrt(TWO_ON_PI/x) *
! 1011: (pone(x)*cos(x-THREE_PI_ON_FOUR) -
! 1012: 8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
! 1013: if (v < 0.0)
! 1014: w = -w;
! 1015: return( w );
! 1016: }
! 1017: }
! 1018:
! 1019: static double ry1(x)
! 1020: double x;
! 1021: {
! 1022: if ( x <= 0.0 )
! 1023: return(dzero/dzero); /* error */
! 1024: if ( x < 8.0 )
! 1025: return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
! 1026: else
! 1027: return( sqrt(TWO_ON_PI/x) *
! 1028: (pone(x)*sin(x-THREE_PI_ON_FOUR) +
! 1029: (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
! 1030: }
! 1031:
! 1032:
! 1033: void f_besj0()
! 1034: {
! 1035: struct value a;
! 1036: (void) pop(&a);
! 1037: if (fabs(imag(&a)) > zero)
! 1038: int_error("can only do bessel functions of reals",NO_CARET);
! 1039: push( Gcomplex(&a,rj0(real(&a)),0.0) );
! 1040: }
! 1041:
! 1042:
! 1043: void f_besj1()
! 1044: {
! 1045: struct value a;
! 1046: (void) pop(&a);
! 1047: if (fabs(imag(&a)) > zero)
! 1048: int_error("can only do bessel functions of reals",NO_CARET);
! 1049: push( Gcomplex(&a,rj1(real(&a)),0.0) );
! 1050: }
! 1051:
! 1052:
! 1053: void f_besy0()
! 1054: {
! 1055: struct value a;
! 1056: (void) pop(&a);
! 1057: if (fabs(imag(&a)) > zero)
! 1058: int_error("can only do bessel functions of reals",NO_CARET);
! 1059: if (real(&a) > 0.0)
! 1060: push( Gcomplex(&a,ry0(real(&a)),0.0) );
! 1061: else {
! 1062: push( Gcomplex(&a,0.0,0.0) );
! 1063: undefined = TRUE ;
! 1064: }
! 1065: }
! 1066:
! 1067:
! 1068: void f_besy1()
! 1069: {
! 1070: struct value a;
! 1071: (void) pop(&a);
! 1072: if (fabs(imag(&a)) > zero)
! 1073: int_error("can only do bessel functions of reals",NO_CARET);
! 1074: if (real(&a) > 0.0)
! 1075: push( Gcomplex(&a,ry1(real(&a)),0.0) );
! 1076: else {
! 1077: push( Gcomplex(&a,0.0,0.0) );
! 1078: undefined = TRUE ;
! 1079: }
! 1080: }
! 1081:
! 1082:
! 1083: /* functions for accessing fields from tm structure, for time series
! 1084: * they are all the same, so define a macro
! 1085: */
! 1086:
! 1087: #define TIMEFUNC(name, field) \
! 1088: void name() { \
! 1089: struct value a; struct tm tm; \
! 1090: (void) pop(&a); ggmtime(&tm, real(&a)); \
! 1091: push(Gcomplex(&a, (double)tm.field, 0.0)); \
! 1092: }
! 1093:
! 1094: TIMEFUNC(f_tmsec, tm_sec)
! 1095: TIMEFUNC(f_tmmin, tm_min)
! 1096: TIMEFUNC(f_tmhour, tm_hour)
! 1097: TIMEFUNC(f_tmmday, tm_mday)
! 1098: TIMEFUNC(f_tmmon, tm_mon)
! 1099: TIMEFUNC(f_tmyear, tm_year)
! 1100: TIMEFUNC(f_tmwday, tm_wday)
! 1101: TIMEFUNC(f_tmyday, tm_yday)
! 1102:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot