This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
update Math-Complex to CPAN version 1.58
[perl5.git] / cpan / Math-Complex / lib / Math / Trig.pm
1 #
2 # Trigonometric functions, mostly inherited from Math::Complex.
3 # -- Jarkko Hietaniemi, since April 1997
4 # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
5 #
6
7 package Math::Trig;
8
9 { use 5.006; }
10 use strict;
11
12 use Math::Complex 1.57;
13 use Math::Complex qw(:trig :pi);
14 require Exporter;
15
16 our @ISA = qw(Exporter);
17
18 our $VERSION = 1.22;
19
20 my @angcnv = qw(rad2deg rad2grad
21                 deg2rad deg2grad
22                 grad2rad grad2deg);
23
24 my @areal = qw(asin_real acos_real);
25
26 our @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
27            @angcnv, @areal);
28
29 my @rdlcnv = qw(cartesian_to_cylindrical
30                 cartesian_to_spherical
31                 cylindrical_to_cartesian
32                 cylindrical_to_spherical
33                 spherical_to_cartesian
34                 spherical_to_cylindrical);
35
36 my @greatcircle = qw(
37                      great_circle_distance
38                      great_circle_direction
39                      great_circle_bearing
40                      great_circle_waypoint
41                      great_circle_midpoint
42                      great_circle_destination
43                     );
44
45 my @pi = qw(pi pi2 pi4 pip2 pip4);
46
47 our @EXPORT_OK = (@rdlcnv, @greatcircle, @pi, 'Inf');
48
49 # See e.g. the following pages:
50 # http://www.movable-type.co.uk/scripts/LatLong.html
51 # http://williams.best.vwh.net/avform.htm
52
53 our %EXPORT_TAGS = ('radial' => [ @rdlcnv ],
54                 'great_circle' => [ @greatcircle ],
55                 'pi'     => [ @pi ]);
56
57 sub _DR  () { pi2/360 }
58 sub _RD  () { 360/pi2 }
59 sub _DG  () { 400/360 }
60 sub _GD  () { 360/400 }
61 sub _RG  () { 400/pi2 }
62 sub _GR  () { pi2/400 }
63
64 #
65 # Truncating remainder.
66 #
67
68 sub _remt ($$) {
69     # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
70     $_[0] - $_[1] * int($_[0] / $_[1]);
71 }
72
73 #
74 # Angle conversions.
75 #
76
77 sub rad2rad($)     { _remt($_[0], pi2) }
78
79 sub deg2deg($)     { _remt($_[0], 360) }
80
81 sub grad2grad($)   { _remt($_[0], 400) }
82
83 sub rad2deg ($;$)  { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) }
84
85 sub deg2rad ($;$)  { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) }
86
87 sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) }
88
89 sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) }
90
91 sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) }
92
93 sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) }
94
95 #
96 # acos and asin functions which always return a real number
97 #
98
99 sub acos_real {
100     return 0  if $_[0] >=  1;
101     return pi if $_[0] <= -1;
102     return acos($_[0]);
103 }
104
105 sub asin_real {
106     return  &pip2 if $_[0] >=  1;
107     return -&pip2 if $_[0] <= -1;
108     return asin($_[0]);
109 }
110
111 sub cartesian_to_spherical {
112     my ( $x, $y, $z ) = @_;
113
114     my $rho = sqrt( $x * $x + $y * $y + $z * $z );
115
116     return ( $rho,
117              atan2( $y, $x ),
118              $rho ? acos_real( $z / $rho ) : 0 );
119 }
120
121 sub spherical_to_cartesian {
122     my ( $rho, $theta, $phi ) = @_;
123
124     return ( $rho * cos( $theta ) * sin( $phi ),
125              $rho * sin( $theta ) * sin( $phi ),
126              $rho * cos( $phi   ) );
127 }
128
129 sub spherical_to_cylindrical {
130     my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
131
132     return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
133 }
134
135 sub cartesian_to_cylindrical {
136     my ( $x, $y, $z ) = @_;
137
138     return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
139 }
140
141 sub cylindrical_to_cartesian {
142     my ( $rho, $theta, $z ) = @_;
143
144     return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
145 }
146
147 sub cylindrical_to_spherical {
148     return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
149 }
150
151 sub great_circle_distance {
152     my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
153
154     $rho = 1 unless defined $rho; # Default to the unit sphere.
155
156     my $lat0 = pip2 - $phi0;
157     my $lat1 = pip2 - $phi1;
158
159     return $rho *
160         acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
161                    sin( $lat0 ) * sin( $lat1 ) );
162 }
163
164 sub great_circle_direction {
165     my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
166
167     my $lat0 = pip2 - $phi0;
168     my $lat1 = pip2 - $phi1;
169
170     return rad2rad(pi2 -
171         atan2(sin($theta0-$theta1) * cos($lat1),
172                 cos($lat0) * sin($lat1) -
173                     sin($lat0) * cos($lat1) * cos($theta0-$theta1)));
174 }
175
176 *great_circle_bearing         = \&great_circle_direction;
177
178 sub great_circle_waypoint {
179     my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
180
181     $point = 0.5 unless defined $point;
182
183     my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
184
185     return undef if $d == pi;
186
187     my $sd = sin($d);
188
189     return ($theta0, $phi0) if $sd == 0;
190
191     my $A = sin((1 - $point) * $d) / $sd;
192     my $B = sin(     $point  * $d) / $sd;
193
194     my $lat0 = pip2 - $phi0;
195     my $lat1 = pip2 - $phi1;
196
197     my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
198     my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
199     my $z = $A * sin($lat0)                + $B * sin($lat1);
200
201     my $theta = atan2($y, $x);
202     my $phi   = acos_real($z);
203
204     return ($theta, $phi);
205 }
206
207 sub great_circle_midpoint {
208     great_circle_waypoint(@_[0..3], 0.5);
209 }
210
211 sub great_circle_destination {
212     my ( $theta0, $phi0, $dir0, $dst ) = @_;
213
214     my $lat0 = pip2 - $phi0;
215
216     my $phi1   = asin_real(sin($lat0)*cos($dst) +
217                            cos($lat0)*sin($dst)*cos($dir0));
218
219     my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
220                                  cos($dst)-sin($lat0)*sin($phi1));
221
222     my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
223
224     $dir1 -= pi2 if $dir1 > pi2;
225
226     return ($theta1, $phi1, $dir1);
227 }
228
229 1;
230
231 __END__
232 =pod
233
234 =head1 NAME
235
236 Math::Trig - trigonometric functions
237
238 =head1 SYNOPSIS
239
240     use Math::Trig;
241
242     $x = tan(0.9);
243     $y = acos(3.7);
244     $z = asin(2.4);
245
246     $halfpi = pi/2;
247
248     $rad = deg2rad(120);
249
250     # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
251     use Math::Trig ':pi';
252
253     # Import the conversions between cartesian/spherical/cylindrical.
254     use Math::Trig ':radial';
255
256         # Import the great circle formulas.
257     use Math::Trig ':great_circle';
258
259 =head1 DESCRIPTION
260
261 C<Math::Trig> defines many trigonometric functions not defined by the
262 core Perl which defines only the C<sin()> and C<cos()>.  The constant
263 B<pi> is also defined as are a few convenience functions for angle
264 conversions, and I<great circle formulas> for spherical movement.
265
266 =head1 TRIGONOMETRIC FUNCTIONS
267
268 The tangent
269
270 =over 4
271
272 =item B<tan>
273
274 =back
275
276 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
277 are aliases)
278
279 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
280
281 The arcus (also known as the inverse) functions of the sine, cosine,
282 and tangent
283
284 B<asin>, B<acos>, B<atan>
285
286 The principal value of the arc tangent of y/x
287
288 B<atan2>(y, x)
289
290 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
291 and acotan/acot are aliases).  Note that atan2(0, 0) is not well-defined.
292
293 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
294
295 The hyperbolic sine, cosine, and tangent
296
297 B<sinh>, B<cosh>, B<tanh>
298
299 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
300 and cotanh/coth are aliases)
301
302 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
303
304 The area (also known as the inverse) functions of the hyperbolic
305 sine, cosine, and tangent
306
307 B<asinh>, B<acosh>, B<atanh>
308
309 The area cofunctions of the hyperbolic sine, cosine, and tangent
310 (acsch/acosech and acoth/acotanh are aliases)
311
312 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
313
314 The trigonometric constant B<pi> and some of handy multiples
315 of it are also defined.
316
317 B<pi, pi2, pi4, pip2, pip4>
318
319 =head2 ERRORS DUE TO DIVISION BY ZERO
320
321 The following functions
322
323     acoth
324     acsc
325     acsch
326     asec
327     asech
328     atanh
329     cot
330     coth
331     csc
332     csch
333     sec
334     sech
335     tan
336     tanh
337
338 cannot be computed for all arguments because that would mean dividing
339 by zero or taking logarithm of zero. These situations cause fatal
340 runtime errors looking like this
341
342     cot(0): Division by zero.
343     (Because in the definition of cot(0), the divisor sin(0) is 0)
344     Died at ...
345
346 or
347
348     atanh(-1): Logarithm of zero.
349     Died at...
350
351 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
352 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
353 C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
354 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
355 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
356 pi>, where I<k> is any integer.
357
358 Note that atan2(0, 0) is not well-defined.
359
360 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
361
362 Please note that some of the trigonometric functions can break out
363 from the B<real axis> into the B<complex plane>. For example
364 C<asin(2)> has no definition for plain real numbers but it has
365 definition for complex numbers.
366
367 In Perl terms this means that supplying the usual Perl numbers (also
368 known as scalars, please see L<perldata>) as input for the
369 trigonometric functions might produce as output results that no more
370 are simple real numbers: instead they are complex numbers.
371
372 The C<Math::Trig> handles this by using the C<Math::Complex> package
373 which knows how to handle complex numbers, please see L<Math::Complex>
374 for more information. In practice you need not to worry about getting
375 complex numbers as results because the C<Math::Complex> takes care of
376 details like for example how to display complex numbers. For example:
377
378     print asin(2), "\n";
379
380 should produce something like this (take or leave few last decimals):
381
382     1.5707963267949-1.31695789692482i
383
384 That is, a complex number with the real part of approximately C<1.571>
385 and the imaginary part of approximately C<-1.317>.
386
387 =head1 PLANE ANGLE CONVERSIONS
388
389 (Plane, 2-dimensional) angles may be converted with the following functions.
390
391 =over
392
393 =item deg2rad
394
395     $radians  = deg2rad($degrees);
396
397 =item grad2rad
398
399     $radians  = grad2rad($gradians);
400
401 =item rad2deg
402
403     $degrees  = rad2deg($radians);
404
405 =item grad2deg
406
407     $degrees  = grad2deg($gradians);
408
409 =item deg2grad
410
411     $gradians = deg2grad($degrees);
412
413 =item rad2grad
414
415     $gradians = rad2grad($radians);
416
417 =back
418
419 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
420 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
421 If you don't want this, supply a true second argument:
422
423     $zillions_of_radians  = deg2rad($zillions_of_degrees, 1);
424     $negative_degrees     = rad2deg($negative_radians, 1);
425
426 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
427 grad2grad().
428
429 =over 4
430
431 =item rad2rad
432
433     $radians_wrapped_by_2pi = rad2rad($radians);
434
435 =item deg2deg
436
437     $degrees_wrapped_by_360 = deg2deg($degrees);
438
439 =item grad2grad
440
441     $gradians_wrapped_by_400 = grad2grad($gradians);
442
443 =back
444
445 =head1 RADIAL COORDINATE CONVERSIONS
446
447 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
448 systems, explained shortly in more detail.
449
450 You can import radial coordinate conversion functions by using the
451 C<:radial> tag:
452
453     use Math::Trig ':radial';
454
455     ($rho, $theta, $z)     = cartesian_to_cylindrical($x, $y, $z);
456     ($rho, $theta, $phi)   = cartesian_to_spherical($x, $y, $z);
457     ($x, $y, $z)           = cylindrical_to_cartesian($rho, $theta, $z);
458     ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
459     ($x, $y, $z)           = spherical_to_cartesian($rho, $theta, $phi);
460     ($rho_c, $theta, $z)   = spherical_to_cylindrical($rho_s, $theta, $phi);
461
462 B<All angles are in radians>.
463
464 =head2 COORDINATE SYSTEMS
465
466 B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
467
468 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
469 coordinates which define a point in three-dimensional space.  They are
470 based on a sphere surface.  The radius of the sphere is B<rho>, also
471 known as the I<radial> coordinate.  The angle in the I<xy>-plane
472 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
473 coordinate.  The angle from the I<z>-axis is B<phi>, also known as the
474 I<polar> coordinate.  The North Pole is therefore I<0, 0, rho>, and
475 the Gulf of Guinea (think of the missing big chunk of Africa) I<0,
476 pi/2, rho>.  In geographical terms I<phi> is latitude (northward
477 positive, southward negative) and I<theta> is longitude (eastward
478 positive, westward negative).
479
480 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
481 some texts define the I<phi> to start from the horizontal plane, some
482 texts use I<r> in place of I<rho>.
483
484 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
485 coordinates which define a point in three-dimensional space.  They are
486 based on a cylinder surface.  The radius of the cylinder is B<rho>,
487 also known as the I<radial> coordinate.  The angle in the I<xy>-plane
488 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
489 coordinate.  The third coordinate is the I<z>, pointing up from the
490 B<theta>-plane.
491
492 =head2 3-D ANGLE CONVERSIONS
493
494 Conversions to and from spherical and cylindrical coordinates are
495 available.  Please notice that the conversions are not necessarily
496 reversible because of the equalities like I<pi> angles being equal to
497 I<-pi> angles.
498
499 =over 4
500
501 =item cartesian_to_cylindrical
502
503     ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
504
505 =item cartesian_to_spherical
506
507     ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
508
509 =item cylindrical_to_cartesian
510
511     ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
512
513 =item cylindrical_to_spherical
514
515     ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
516
517 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
518
519 =item spherical_to_cartesian
520
521     ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
522
523 =item spherical_to_cylindrical
524
525     ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
526
527 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
528
529 =back
530
531 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
532
533 A great circle is section of a circle that contains the circle
534 diameter: the shortest distance between two (non-antipodal) points on
535 the spherical surface goes along the great circle connecting those two
536 points.
537
538 =head2 great_circle_distance
539
540 You can compute spherical distances, called B<great circle distances>,
541 by importing the great_circle_distance() function:
542
543   use Math::Trig 'great_circle_distance';
544
545   $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
546
547 The I<great circle distance> is the shortest distance between two
548 points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
549 optional, it defaults to 1 (the unit sphere), therefore the distance
550 defaults to radians.
551
552 If you think geographically the I<theta> are longitudes: zero at the
553 Greenwhich meridian, eastward positive, westward negative -- and the
554 I<phi> are latitudes: zero at the North Pole, northward positive,
555 southward negative.  B<NOTE>: this formula thinks in mathematics, not
556 geographically: the I<phi> zero is at the North Pole, not at the
557 Equator on the west coast of Africa (Bay of Guinea).  You need to
558 subtract your geographical coordinates from I<pi/2> (also known as 90
559 degrees).
560
561   $distance = great_circle_distance($lon0, pi/2 - $lat0,
562                                     $lon1, pi/2 - $lat1, $rho);
563
564 =head2 great_circle_direction
565
566 The direction you must follow the great circle (also known as I<bearing>)
567 can be computed by the great_circle_direction() function:
568
569   use Math::Trig 'great_circle_direction';
570
571   $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
572
573 =head2 great_circle_bearing
574
575 Alias 'great_circle_bearing' for 'great_circle_direction' is also available.
576
577   use Math::Trig 'great_circle_bearing';
578
579   $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1);
580
581 The result of great_circle_direction is in radians, zero indicating
582 straight north, pi or -pi straight south, pi/2 straight west, and
583 -pi/2 straight east.
584
585 =head2 great_circle_destination
586
587 You can inversely compute the destination if you know the
588 starting point, direction, and distance:
589
590   use Math::Trig 'great_circle_destination';
591
592   # $diro is the original direction,
593   # for example from great_circle_bearing().
594   # $distance is the angular distance in radians,
595   # for example from great_circle_distance().
596   # $thetad and $phid are the destination coordinates,
597   # $dird is the final direction at the destination.
598
599   ($thetad, $phid, $dird) =
600     great_circle_destination($theta, $phi, $diro, $distance);
601
602 or the midpoint if you know the end points:
603
604 =head2 great_circle_midpoint
605
606   use Math::Trig 'great_circle_midpoint';
607
608   ($thetam, $phim) =
609     great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
610
611 The great_circle_midpoint() is just a special case of
612
613 =head2 great_circle_waypoint
614
615   use Math::Trig 'great_circle_waypoint';
616
617   ($thetai, $phii) =
618     great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
619
620 Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
621 $phi1).  Note that antipodal points (where their distance is I<pi>
622 radians) do not have waypoints between them (they would have an an
623 "equator" between them), and therefore C<undef> is returned for
624 antipodal points.  If the points are the same and the distance
625 therefore zero and all waypoints therefore identical, the first point
626 (either point) is returned.
627
628 The thetas, phis, direction, and distance in the above are all in radians.
629
630 You can import all the great circle formulas by
631
632   use Math::Trig ':great_circle';
633
634 Notice that the resulting directions might be somewhat surprising if
635 you are looking at a flat worldmap: in such map projections the great
636 circles quite often do not look like the shortest routes --  but for
637 example the shortest possible routes from Europe or North America to
638 Asia do often cross the polar regions.  (The common Mercator projection
639 does B<not> show great circles as straight lines: straight lines in the
640 Mercator projection are lines of constant bearing.)
641
642 =head1 EXAMPLES
643
644 To calculate the distance between London (51.3N 0.5W) and Tokyo
645 (35.7N 139.8E) in kilometers:
646
647     use Math::Trig qw(great_circle_distance deg2rad);
648
649     # Notice the 90 - latitude: phi zero is at the North Pole.
650     sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
651     my @L = NESW( -0.5, 51.3);
652     my @T = NESW(139.8, 35.7);
653     my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
654
655 The direction you would have to go from London to Tokyo (in radians,
656 straight north being zero, straight east being pi/2).
657
658     use Math::Trig qw(great_circle_direction);
659
660     my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
661
662 The midpoint between London and Tokyo being
663
664     use Math::Trig qw(great_circle_midpoint);
665
666     my @M = great_circle_midpoint(@L, @T);
667
668 or about 69 N 89 E, in the frozen wastes of Siberia.
669
670 B<NOTE>: you B<cannot> get from A to B like this:
671
672    Dist = great_circle_distance(A, B)
673    Dir  = great_circle_direction(A, B)
674    C    = great_circle_destination(A, Dist, Dir)
675
676 and expect C to be B, because the bearing constantly changes when
677 going from A to B (except in some special case like the meridians or
678 the circles of latitudes) and in great_circle_destination() one gives
679 a B<constant> bearing to follow.
680
681 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
682
683 The answers may be off by few percentages because of the irregular
684 (slightly aspherical) form of the Earth.  The errors are at worst
685 about 0.55%, but generally below 0.3%.
686
687 =head2 Real-valued asin and acos
688
689 For small inputs asin() and acos() may return complex numbers even
690 when real numbers would be enough and correct, this happens because of
691 floating-point inaccuracies.  You can see these inaccuracies for
692 example by trying theses:
693
694   print cos(1e-6)**2+sin(1e-6)**2 - 1,"\n";
695   printf "%.20f", cos(1e-6)**2+sin(1e-6)**2,"\n";
696
697 which will print something like this
698
699   -1.11022302462516e-16
700   0.99999999999999988898
701
702 even though the expected results are of course exactly zero and one.
703 The formulas used to compute asin() and acos() are quite sensitive to
704 this, and therefore they might accidentally slip into the complex
705 plane even when they should not.  To counter this there are two
706 interfaces that are guaranteed to return a real-valued output.
707
708 =over 4
709
710 =item asin_real
711
712     use Math::Trig qw(asin_real);
713
714     $real_angle = asin_real($input_sin);
715
716 Return a real-valued arcus sine if the input is between [-1, 1],
717 B<inclusive> the endpoints.  For inputs greater than one, pi/2
718 is returned.  For inputs less than minus one, -pi/2 is returned.
719
720 =item acos_real
721
722     use Math::Trig qw(acos_real);
723
724     $real_angle = acos_real($input_cos);
725
726 Return a real-valued arcus cosine if the input is between [-1, 1],
727 B<inclusive> the endpoints.  For inputs greater than one, zero
728 is returned.  For inputs less than minus one, pi is returned.
729
730 =back
731
732 =head1 BUGS
733
734 Saying C<use Math::Trig;> exports many mathematical routines in the
735 caller environment and even overrides some (C<sin>, C<cos>).  This is
736 construed as a feature by the Authors, actually... ;-)
737
738 The code is not optimized for speed, especially because we use
739 C<Math::Complex> and thus go quite near complex numbers while doing
740 the computations even when the arguments are not. This, however,
741 cannot be completely avoided if we want things like C<asin(2)> to give
742 an answer instead of giving a fatal runtime error.
743
744 Do not attempt navigation using these formulas.
745
746 L<Math::Complex>
747
748 =head1 AUTHORS
749
750 Jarkko Hietaniemi <F<jhi!at!iki.fi>>,
751 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>,
752 Zefram <zefram@fysh.org>
753
754 =head1 LICENSE
755
756 This library is free software; you can redistribute it and/or modify
757 it under the same terms as Perl itself. 
758
759 =cut
760
761 # eof