2008-03-28 14:56:47 -06:00
|
|
|
// Copyright 2009 The Go Authors. All rights reserved.
|
|
|
|
// Use of this source code is governed by a BSD-style
|
|
|
|
// license that can be found in the LICENSE file.
|
|
|
|
|
2008-06-27 18:06:23 -06:00
|
|
|
package math
|
2008-03-28 14:56:47 -06:00
|
|
|
|
|
|
|
/*
|
2010-01-11 13:38:31 -07:00
|
|
|
Floating-point arctangent.
|
|
|
|
*/
|
|
|
|
|
2012-06-24 17:39:07 -06:00
|
|
|
// The original C code, the long comment, and the constants below were
|
|
|
|
// from http://netlib.sandia.gov/cephes/cmath/atan.c, available from
|
|
|
|
// http://www.netlib.org/cephes/cmath.tgz.
|
|
|
|
// The go code is a version of the original C.
|
|
|
|
//
|
|
|
|
// atan.c
|
|
|
|
// Inverse circular tangent (arctangent)
|
|
|
|
//
|
|
|
|
// SYNOPSIS:
|
|
|
|
// double x, y, atan();
|
|
|
|
// y = atan( x );
|
|
|
|
//
|
|
|
|
// DESCRIPTION:
|
|
|
|
// Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
|
|
|
|
//
|
|
|
|
// Range reduction is from three intervals into the interval from zero to 0.66.
|
|
|
|
// The approximant uses a rational function of degree 4/5 of the form
|
|
|
|
// x + x**3 P(x)/Q(x).
|
|
|
|
//
|
|
|
|
// ACCURACY:
|
|
|
|
// Relative error:
|
|
|
|
// arithmetic domain # trials peak rms
|
|
|
|
// DEC -10, 10 50000 2.4e-17 8.3e-18
|
|
|
|
// IEEE -10, 10 10^6 1.8e-16 5.0e-17
|
|
|
|
//
|
|
|
|
// Cephes Math Library Release 2.8: June, 2000
|
|
|
|
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
|
|
|
|
//
|
|
|
|
// The readme file at http://netlib.sandia.gov/cephes/ says:
|
|
|
|
// Some software in this archive may be from the book _Methods and
|
|
|
|
// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
|
|
|
|
// International, 1989) or from the Cephes Mathematical Library, a
|
|
|
|
// commercial product. In either event, it is copyrighted by the author.
|
|
|
|
// What you see here may be used freely but it comes with no support or
|
|
|
|
// guarantee.
|
|
|
|
//
|
|
|
|
// The two known misprints in the book are repaired here in the
|
|
|
|
// source listings for the gamma function and the incomplete beta
|
|
|
|
// integral.
|
|
|
|
//
|
|
|
|
// Stephen L. Moshier
|
|
|
|
// moshier@na-net.ornl.gov
|
|
|
|
|
|
|
|
// xatan evaluates a series valid in the range [0, 0.66].
|
|
|
|
func xatan(x float64) float64 {
|
2009-10-06 20:40:35 -06:00
|
|
|
const (
|
2012-06-24 17:39:07 -06:00
|
|
|
P0 = -8.750608600031904122785e-01
|
|
|
|
P1 = -1.615753718733365076637e+01
|
|
|
|
P2 = -7.500855792314704667340e+01
|
|
|
|
P3 = -1.228866684490136173410e+02
|
|
|
|
P4 = -6.485021904942025371773e+01
|
|
|
|
Q0 = +2.485846490142306297962e+01
|
|
|
|
Q1 = +1.650270098316988542046e+02
|
|
|
|
Q2 = +4.328810604912902668951e+02
|
|
|
|
Q3 = +4.853903996359136964868e+02
|
|
|
|
Q4 = +1.945506571482613964425e+02
|
2009-01-15 20:11:32 -07:00
|
|
|
)
|
2012-06-24 17:39:07 -06:00
|
|
|
z := x * x
|
|
|
|
z = z * ((((P0*z+P1)*z+P2)*z+P3)*z + P4) / (((((z+Q0)*z+Q1)*z+Q2)*z+Q3)*z + Q4)
|
|
|
|
z = x*z + x
|
|
|
|
return z
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
|
|
|
|
2010-01-11 13:38:31 -07:00
|
|
|
// satan reduces its argument (known to be positive)
|
2012-06-24 17:39:07 -06:00
|
|
|
// to the range [0, 0.66] and calls xatan.
|
|
|
|
func satan(x float64) float64 {
|
|
|
|
const (
|
|
|
|
Morebits = 6.123233995736765886130e-17 // pi/2 = PIO2 + Morebits
|
|
|
|
Tan3pio8 = 2.41421356237309504880 // tan(3*pi/8)
|
|
|
|
)
|
|
|
|
if x <= 0.66 {
|
|
|
|
return xatan(x)
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
2012-06-24 17:39:07 -06:00
|
|
|
if x > Tan3pio8 {
|
|
|
|
return Pi/2 - xatan(1/x) + Morebits
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
2012-06-24 17:39:07 -06:00
|
|
|
return Pi/4 + xatan((x-1)/(x+1)) + 0.5*Morebits
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
|
|
|
|
2010-01-11 13:38:31 -07:00
|
|
|
// Atan returns the arctangent of x.
|
2010-04-26 23:44:39 -06:00
|
|
|
//
|
|
|
|
// Special cases are:
|
2012-06-24 17:39:07 -06:00
|
|
|
// Atan(±0) = ±0
|
|
|
|
// Atan(±Inf) = ±Pi/2
|
math: regularize build
This will be nicer to the automatic tools.
It requires a few more assembly stubs
but fewer Go files.
There are a few instances where it looks like
there are new blobs of code, but they are just
being copied out of deleted files.
There is no new code here.
Suppose you have a portable implementation for Sin
and a 386-specific assembly one. The old way to
do this was to write three files
sin_decl.go
func Sin(x float64) float64 // declaration only
sin_386.s
assembly implementation
sin_port.go
func Sin(x float64) float64 { ... } // pure-Go impl
and then link in either sin_decl.go+sin_386.s or
just sin_port.go. The Makefile actually did the magic
of linking in only the _port.go files for those without
assembly and only the _decl.go files for those with
assembly, or at least some of that magic.
The biggest problem with this, beyond being hard
to explain to the build system, is that once you do
explain it to the build system, godoc knows which
of sin_port.go or sin_decl.go are involved on a given
architecture, and it (correctly) ignores the other.
That means you have to put identical doc comments
in both files.
The new approach, which is more like what we did
in the later packages math/big and sync/atomic,
is to have
sin.go
func Sin(x float64) float64 // decl only
func sin(x float64) float64 {...} // pure-Go impl
sin_386.s
// assembly for Sin (ignores sin)
sin_amd64.s
// assembly for Sin: jmp sin
sin_arm.s
// assembly for Sin: jmp sin
Once we abandon Makefiles we can put all the assembly
stubs in one source file, so the number of files will
actually go down.
Chris asked whether the branches cost anything.
Given that they are branching to pure-Go implementations
that are not typically known for their speed, the single
direct branch is not going to be noticeable. That is,
it's on the slow path.
An alternative would have been to preserve the old
"only write assembly files when there's an implementation"
and still have just one copy of the declaration of Sin
(and thus one doc comment) by doing:
sin.go
func Sin(x float64) float64 { return sin(x) }
sin_decl.go
func sin(x float64) float64 // declaration only
sin_386.s
// assembly for sin
sin_port.go
func sin(x float64) float64 { portable code }
In this version everyone would link in sin.go and
then either sin_decl.go+sin_386.s or sin_port.go.
This has an extra function call on all paths, including
the "fast path" to get to assembly, and it triples the
number of Go files involved compared to what I did
in this CL. On the other hand you don't have to
write assembly stubs. After starting down this path
I decided that the assembly stubs were the easier
approach.
As for generating the assembly stubs on the fly, much
of the goal here is to eliminate magic from the build
process, so that zero-configuration tools like goinstall
or the new go tool can handle this package.
R=golang-dev, r, cw, iant, r
CC=golang-dev
https://golang.org/cl/5488057
2011-12-13 13:20:12 -07:00
|
|
|
func Atan(x float64) float64
|
|
|
|
|
|
|
|
func atan(x float64) float64 {
|
2010-04-26 23:44:39 -06:00
|
|
|
if x == 0 {
|
|
|
|
return x
|
|
|
|
}
|
2009-03-05 14:31:01 -07:00
|
|
|
if x > 0 {
|
2009-11-09 13:07:39 -07:00
|
|
|
return satan(x)
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
2009-12-15 16:35:38 -07:00
|
|
|
return -satan(-x)
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|