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