Annotation of OpenXM_contrib/gnuplot/standard.c, Revision 1.1.1.3
1.1 maekawa 1: #ifndef lint
1.1.1.3 ! ohara 2: static char *RCSid = "$Id: standard.c,v 1.5.2.1 2000/05/07 16:44:57 lhecking Exp $";
1.1 maekawa 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 {
1.1.1.3 ! ohara 575: double alpha = sqrt((x + 1) * (x + 1) + y * y) / 2
! 576: + sqrt((x - 1) * (x - 1) + y * y) / 2;
! 577: double beta = sqrt((x + 1) * (x + 1) + y * y) / 2
! 578: - sqrt((x - 1) * (x - 1) + y * y) / 2;
1.1 maekawa 579: if (beta > 1)
580: beta = 1; /* Avoid rounding error problems */
581: else if (beta < -1)
582: beta = -1;
1.1.1.3 ! ohara 583: push( Gcomplex(&a, (y > 0? -1: 1)*acos(beta) / ang2rad,
! 584: log(alpha + sqrt(alpha*alpha-1)) / ang2rad));
1.1 maekawa 585: }
586: }
587:
588: void f_atan()
589: {
590: struct value a;
591: register double x, y, u, v, w, z;
592: (void) pop(&a);
593: x = real(&a); y = imag(&a);
594: if (y == 0.0)
595: push( Gcomplex(&a,atan(x)/ang2rad, 0.0) );
596: else if (x == 0.0 && fabs(y) >= 1.0) {
597: undefined = TRUE;
598: push(Gcomplex(&a,0.0, 0.0));
599: } else {
600: if (x >= 0) {
601: u = x;
602: v = y;
603: } else {
604: u = -x;
605: v = -y;
606: }
607:
608: z = atan(2*u/(1-u*u-v*v));
609: w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
610: if (z < 0)
611: z = z + 2*PI_ON_TWO;
612: if (x < 0) {
613: z = -z;
614: w = -w;
615: }
616: push( Gcomplex(&a,0.5*z/ang2rad, w) );
617: }
618: }
619:
620: /* real parts only */
621: void f_atan2()
622: {
623: struct value a;
624: double x;
625: double y;
626:
627: x = real(pop(&a));
628: y = real(pop(&a));
629:
630: if (x == 0.0 && y == 0.0) {
631: undefined = TRUE;
632: push(Ginteger(&a,0));
633: }
634: push(Gcomplex(&a,atan2(y,x)/ang2rad,0.0));
635: }
636:
637:
638: void f_sinh()
639: {
640: struct value a;
641: (void) pop(&a);
642: push( Gcomplex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
643: }
644:
645: void f_cosh()
646: {
647: struct value a;
648: (void) pop(&a);
649: push( Gcomplex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
650: }
651:
652: void f_tanh()
653: {
654: struct value a;
655: register double den;
656: (void) pop(&a);
657: den = cosh(2*real(&a)) + cos(2*imag(&a));
658: push( Gcomplex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
659: }
660:
661: void f_asinh()
662: {
663: struct value a; /* asinh(z) = -I*asin(I*z) */
664: register double alpha, beta, x, y;
665: (void) pop(&a);
666: x = -imag(&a); y = real(&a);
667: if (y == 0.0 && fabs(x) <= 1.0) {
668: push( Gcomplex(&a, 0.0, -asin(x)/ang2rad) );
669: } else if (y == 0.0) {
670: push( Gcomplex(&a, 0.0, 0.0) );
671: undefined = TRUE;
672: } else if (x == 0.0) {
673: push( Gcomplex(&a, log(y+sqrt(y*y+1))/ang2rad, 0.0) );
674: } else {
675: beta = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
676: alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
677: push( Gcomplex(&a, log(alpha + sqrt(alpha*alpha-1))/ang2rad, -asin(beta)/ang2rad) );
678: }
679: }
680:
681: void f_acosh()
682: {
683: struct value a;
684: register double alpha, beta, x, y; /* acosh(z) = I*acos(z) */
685: (void) pop(&a);
686: x = real(&a); y = imag(&a);
687: if (y == 0.0 && fabs(x) <= 1.0) {
688: push( Gcomplex(&a, 0.0, acos(x)/ang2rad) );
689: } else if (y == 0) {
690: push( Gcomplex(&a, log(x+sqrt(x*x-1))/ang2rad, 0.0) );
691: } else {
1.1.1.3 ! ohara 692: alpha = sqrt((x + 1) * (x + 1) + y * y) / 2
! 693: + sqrt((x - 1) * (x - 1) + y * y) / 2;
! 694: beta = sqrt((x + 1) * (x + 1) + y * y) / 2
! 695: - sqrt((x - 1) * (x - 1) + y * y) / 2;
! 696: push( Gcomplex(&a, log(alpha + sqrt(alpha * alpha - 1)) / ang2rad, (y<0 ? -1 : 1) * acos(beta) / ang2rad));
1.1 maekawa 697: }
698: }
699:
700: void f_atanh()
701: {
702: struct value a;
703: register double x, y, u, v, w, z; /* atanh(z) = -I*atan(I*z) */
704: (void) pop(&a);
705: x = -imag(&a); y = real(&a);
706: if (y == 0.0)
707: push( Gcomplex(&a, 0.0, -atan(x)/ang2rad) );
708: else if (x == 0.0 && fabs(y) >= 1.0) {
709: undefined = TRUE;
710: push(Gcomplex(&a,0.0, 0.0));
711: } else {
712: if (x >= 0) {
713: u = x;
714: v = y;
715: } else {
716: u = -x;
717: v = -y;
718: }
719:
720: z = atan(2*u/(1-u*u-v*v));
721: w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
722: if (z < 0)
723: z = z + 2*PI_ON_TWO;
724: if (x < 0) {
725: z = -z;
726: w = -w;
727: }
728: push( Gcomplex(&a, w, -0.5*z/ang2rad) );
729: }
730: }
731:
732: void f_int()
733: {
734: struct value a;
735: push( Ginteger(&a,(int)real(pop(&a))) );
736: }
737:
738:
739: void f_abs()
740: {
741: struct value a;
742: (void) pop(&a);
743: switch (a.type) {
744: case INTGR:
745: push( Ginteger(&a,abs(a.v.int_val)) );
746: break;
747: case CMPLX:
748: push( Gcomplex(&a,magnitude(&a), 0.0) );
749: }
750: }
751:
752: void f_sgn()
753: {
754: struct value a;
755: (void) pop(&a);
756: switch(a.type) {
757: case INTGR:
758: push( Ginteger(&a,(a.v.int_val > 0) ? 1 :
759: (a.v.int_val < 0) ? -1 : 0) );
760: break;
761: case CMPLX:
762: push( Ginteger(&a,(a.v.cmplx_val.real > 0.0) ? 1 :
763: (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
764: break;
765: }
766: }
767:
768:
769: void f_sqrt()
770: {
771: struct value a;
772: register double mag;
773: (void) pop(&a);
774: mag = sqrt(magnitude(&a));
775: if (imag(&a) == 0.0) {
776: if (real(&a) < 0.0)
777: push( Gcomplex(&a,0.0,mag) );
778: else
779: push( Gcomplex(&a,mag, 0.0) );
780: } else {
781: /* -pi < ang < pi, so real(sqrt(z)) >= 0 */
782: double ang = angle(&a) / 2.0;
783: push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
784: }
785: }
786:
787:
788: void f_exp()
789: {
790: struct value a;
791: register double mag, ang;
792: (void) pop(&a);
793: mag = gp_exp(real(&a));
794: ang = imag(&a);
795: push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
796: }
797:
798:
799: void f_log10()
800: {
801: struct value a;
802: (void) pop(&a);
803: push( Gcomplex(&a,log(magnitude(&a))/M_LN10, angle(&a)/M_LN10) );
804: }
805:
806:
807: void f_log()
808: {
809: struct value a;
810: (void) pop(&a);
811: push( Gcomplex(&a,log(magnitude(&a)), angle(&a)) );
812: }
813:
814:
815: void f_floor()
816: {
817: struct value a;
818:
819: (void) pop(&a);
820: switch (a.type) {
821: case INTGR:
822: push( Ginteger(&a,(int)floor((double)a.v.int_val)));
823: break;
824: case CMPLX:
825: push( Ginteger(&a,(int)floor(a.v.cmplx_val.real)));
826: }
827: }
828:
829:
830: void f_ceil()
831: {
832: struct value a;
833:
834: (void) pop(&a);
835: switch (a.type) {
836: case INTGR:
837: push( Ginteger(&a,(int)ceil((double)a.v.int_val)));
838: break;
839: case CMPLX:
840: push( Ginteger(&a,(int)ceil(a.v.cmplx_val.real)));
841: }
842: }
843:
844: /* bessel function approximations */
845: static double jzero(x)
846: double x;
847: {
848: double p, q, x2;
849: int n;
850:
851: x2 = x * x;
852: p = pjzero[8];
853: q = qjzero[8];
854: for (n=7; n>=0; n--) {
855: p = p*x2 + pjzero[n];
856: q = q*x2 + qjzero[n];
857: }
858: return(p/q);
859: }
860:
861: static double pzero(x)
862: double x;
863: {
864: double p, q, z, z2;
865: int n;
866:
867: z = 8.0 / x;
868: z2 = z * z;
869: p = ppzero[5];
870: q = qpzero[5];
871: for (n=4; n>=0; n--) {
872: p = p*z2 + ppzero[n];
873: q = q*z2 + qpzero[n];
874: }
875: return(p/q);
876: }
877:
878: static double qzero(x)
879: double x;
880: {
881: double p, q, z, z2;
882: int n;
883:
884: z = 8.0 / x;
885: z2 = z * z;
886: p = pqzero[5];
887: q = qqzero[5];
888: for (n=4; n>=0; n--) {
889: p = p*z2 + pqzero[n];
890: q = q*z2 + qqzero[n];
891: }
892: return(p/q);
893: }
894:
895: static double yzero(x)
896: double x;
897: {
898: double p, q, x2;
899: int n;
900:
901: x2 = x * x;
902: p = pyzero[8];
903: q = qyzero[8];
904: for (n=7; n>=0; n--) {
905: p = p*x2 + pyzero[n];
906: q = q*x2 + qyzero[n];
907: }
908: return(p/q);
909: }
910:
911: static double rj0(x)
912: double x;
913: {
914: if ( x <= 0.0 )
915: x = -x;
916: if ( x < 8.0 )
917: return(jzero(x));
918: else
919: return( sqrt(TWO_ON_PI/x) *
920: (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
921:
922: }
923:
924: static double ry0(x)
925: double x;
926: {
927: if ( x < 0.0 )
928: return(dzero/dzero); /* error */
929: if ( x < 8.0 )
930: return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
931: else
932: return( sqrt(TWO_ON_PI/x) *
933: (pzero(x)*sin(x-PI_ON_FOUR) +
934: (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
935:
936: }
937:
938:
939: static double jone(x)
940: double x;
941: {
942: double p, q, x2;
943: int n;
944:
945: x2 = x * x;
946: p = pjone[8];
947: q = qjone[8];
948: for (n=7; n>=0; n--) {
949: p = p*x2 + pjone[n];
950: q = q*x2 + qjone[n];
951: }
952: return(p/q);
953: }
954:
955: static double pone(x)
956: double x;
957: {
958: double p, q, z, z2;
959: int n;
960:
961: z = 8.0 / x;
962: z2 = z * z;
963: p = ppone[5];
964: q = qpone[5];
965: for (n=4; n>=0; n--) {
966: p = p*z2 + ppone[n];
967: q = q*z2 + qpone[n];
968: }
969: return(p/q);
970: }
971:
972: static double qone(x)
973: double x;
974: {
975: double p, q, z, z2;
976: int n;
977:
978: z = 8.0 / x;
979: z2 = z * z;
980: p = pqone[5];
981: q = qqone[5];
982: for (n=4; n>=0; n--) {
983: p = p*z2 + pqone[n];
984: q = q*z2 + qqone[n];
985: }
986: return(p/q);
987: }
988:
989: static double yone(x)
990: double x;
991: {
992: double p, q, x2;
993: int n;
994:
995: x2 = x * x;
996: p = 0.0;
997: q = qyone[8];
998: for (n=7; n>=0; n--) {
999: p = p*x2 + pyone[n];
1000: q = q*x2 + qyone[n];
1001: }
1002: return(p/q);
1003: }
1004:
1005: static double rj1(x)
1006: double x;
1007: {
1008: double v,w;
1009: v = x;
1010: if ( x < 0.0 )
1011: x = -x;
1012: if ( x < 8.0 )
1013: return(v*jone(x));
1014: else {
1015: w = sqrt(TWO_ON_PI/x) *
1016: (pone(x)*cos(x-THREE_PI_ON_FOUR) -
1017: 8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
1018: if (v < 0.0)
1019: w = -w;
1020: return( w );
1021: }
1022: }
1023:
1024: static double ry1(x)
1025: double x;
1026: {
1027: if ( x <= 0.0 )
1028: return(dzero/dzero); /* error */
1029: if ( x < 8.0 )
1030: return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
1031: else
1032: return( sqrt(TWO_ON_PI/x) *
1033: (pone(x)*sin(x-THREE_PI_ON_FOUR) +
1034: (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
1035: }
1036:
1037:
1038: void f_besj0()
1039: {
1040: struct value a;
1041: (void) pop(&a);
1042: if (fabs(imag(&a)) > zero)
1043: int_error("can only do bessel functions of reals",NO_CARET);
1044: push( Gcomplex(&a,rj0(real(&a)),0.0) );
1045: }
1046:
1047:
1048: void f_besj1()
1049: {
1050: struct value a;
1051: (void) pop(&a);
1052: if (fabs(imag(&a)) > zero)
1053: int_error("can only do bessel functions of reals",NO_CARET);
1054: push( Gcomplex(&a,rj1(real(&a)),0.0) );
1055: }
1056:
1057:
1058: void f_besy0()
1059: {
1060: struct value a;
1061: (void) pop(&a);
1062: if (fabs(imag(&a)) > zero)
1063: int_error("can only do bessel functions of reals",NO_CARET);
1064: if (real(&a) > 0.0)
1065: push( Gcomplex(&a,ry0(real(&a)),0.0) );
1066: else {
1067: push( Gcomplex(&a,0.0,0.0) );
1068: undefined = TRUE ;
1069: }
1070: }
1071:
1072:
1073: void f_besy1()
1074: {
1075: struct value a;
1076: (void) pop(&a);
1077: if (fabs(imag(&a)) > zero)
1078: int_error("can only do bessel functions of reals",NO_CARET);
1079: if (real(&a) > 0.0)
1080: push( Gcomplex(&a,ry1(real(&a)),0.0) );
1081: else {
1082: push( Gcomplex(&a,0.0,0.0) );
1083: undefined = TRUE ;
1084: }
1085: }
1086:
1087:
1088: /* functions for accessing fields from tm structure, for time series
1089: * they are all the same, so define a macro
1090: */
1091:
1092: #define TIMEFUNC(name, field) \
1093: void name() { \
1094: struct value a; struct tm tm; \
1095: (void) pop(&a); ggmtime(&tm, real(&a)); \
1096: push(Gcomplex(&a, (double)tm.field, 0.0)); \
1097: }
1098:
1099: TIMEFUNC(f_tmsec, tm_sec)
1100: TIMEFUNC(f_tmmin, tm_min)
1101: TIMEFUNC(f_tmhour, tm_hour)
1102: TIMEFUNC(f_tmmday, tm_mday)
1103: TIMEFUNC(f_tmmon, tm_mon)
1104: TIMEFUNC(f_tmyear, tm_year)
1105: TIMEFUNC(f_tmwday, tm_wday)
1106: TIMEFUNC(f_tmyday, tm_yday)
1107:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>