git.maemo.org Git - gnuplot/blob - src/standard.c

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