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