← Back to branch summary

~john-koepi/ubuntu/trusty/golang/default

« back to all changes in this revision

Viewing changes to src/pkg/math/hypot_port.go

  • Committer: Bazaar Package Importer
  • Author(s): Ondřej Surý
  • Date: 2011-04-20 17:36:48 UTC
  • Revision ID: james.westby@ubuntu.com-20110420173648-ifergoxyrm832trd
Tags: upstream-2011.03.07.1
Import upstream version 2011.03.07.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pkg/math/hypot_port.go
 
1
// Copyright 2009-2010 The Go Authors. All rights reserved.
 
2
// Use of this source code is governed by a BSD-style
 
3
// license that can be found in the LICENSE file.
 
4
 
 
5
package math
 
6
 
 
7
/*
 
8
        Hypot -- sqrt(p*p + q*q), but overflows only if the result does.
 
9
        See:
 
10
                Cleve Moler and Donald Morrison,
 
11
                Replacing Square Roots by Pythagorean Sums
 
12
                IBM Journal of Research and Development,
 
13
                Vol. 27, Number 6, pp. 577-581, Nov. 1983
 
14
*/
 
15
 
 
16
// Hypot computes Sqrt(p*p + q*q), taking care to avoid
 
17
// unnecessary overflow and underflow.
 
18
//
 
19
// Special cases are:
 
20
//      Hypot(p, q) = +Inf if p or q is infinite
 
21
//      Hypot(p, q) = NaN if p or q is NaN
 
22
func hypotGo(p, q float64) float64 {
 
23
        // TODO(rsc): Remove manual inlining of IsNaN, IsInf
 
24
        // when compiler does it for us
 
25
        // special cases
 
26
        switch {
 
27
        case p < -MaxFloat64 || p > MaxFloat64 || q < -MaxFloat64 || q > MaxFloat64: // IsInf(p, 0) || IsInf(q, 0):
 
28
                return Inf(1)
 
29
        case p != p || q != q: // IsNaN(p) || IsNaN(q):
 
30
                return NaN()
 
31
        }
 
32
        if p < 0 {
 
33
                p = -p
 
34
        }
 
35
        if q < 0 {
 
36
                q = -q
 
37
        }
 
38
 
 
39
        if p < q {
 
40
                p, q = q, p
 
41
        }
 
42
 
 
43
        if p == 0 {
 
44
                return 0
 
45
        }
 
46
 
 
47
        pfac := p
 
48
        q = q / p
 
49
        r := q
 
50
        p = 1
 
51
        for {
 
52
                r = r * r
 
53
                s := r + 4
 
54
                if s == 4 {
 
55
                        return p * pfac
 
56
                }
 
57
                r = r / s
 
58
                p = p + 2*r*p
 
59
                q = q * r
 
60
                r = q / p
 
61
        }
 
62
        panic("unreachable")
 
63
}