go-hep.org/x/hep@v0.38.1/fastjet/internal/predicates/float64Pred.go (about)

     1  // Copyright ©2017 The go-hep 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 predicates
     6  
     7  import (
     8  	"math"
     9  )
    10  
    11  const (
    12  	// macheps is the machine epsilon aka unit roundoff
    13  	// The machine epsilon is an upper bound on the absolute relative true error in
    14  	// representing a number.
    15  	// If y is the machine representation of x then |(x-y)/x| <= macheps
    16  	// https://en.wikipedia.org/wiki/Machine_epsilon
    17  	// Go's float64 type has a 52-bit fractional mantissa,
    18  	// therefore the value 2^-52
    19  	macheps = 1.0 / (1 << 52)
    20  )
    21  
    22  // float64Pred dynamically updates the potential error.
    23  //
    24  // If y is the machine representation of x then |(x-y)/x| <= macheps and |x-y| = e.
    25  // Since we want the max possible error we assume |(x-y)/x| = macheps
    26  // macheps*|x| = |x - y|
    27  // if (x-y)>=0 then macheps*|x| = x - y  ->  y = x- macheps*|x|
    28  // else macheps*|x| = -(x - y)  ->  macheps*|x|+x = y
    29  // Each one of these has two cases again. x can be positive or negative, resulting in the same two possible equations:
    30  // y/(1-macheps)=x or y/(1+macheps)=x.
    31  // Since x is unknown, we will use the larger number as factor, to avoid having an error greater than
    32  // the maxError value we have.
    33  // |y-x| = macheps*|x| -> e = macheps*|x| -> e = macheps* |y|/(1-macheps)
    34  // A Special case is when y = 0. Then we use the smallest nonzero float, because that is the max
    35  // possible error in this case.
    36  type float64Pred struct {
    37  	// n is the number
    38  	n float64
    39  	// e is the max rounding error possible
    40  	e float64
    41  }
    42  
    43  // newFloat64Pred returns a new float64Pred e set to 0.
    44  func newFloat64Pred(n float64) float64Pred {
    45  	return float64Pred{
    46  		n: n,
    47  		e: 0,
    48  	}
    49  }
    50  
    51  // // addFloat64 adds f to p and updates the potential error
    52  // func (p float64Pred) addFloat64(f float64) float64Pred {
    53  // 	p.n += f
    54  // 	if p.n == 0 {
    55  // 		p.e += math.SmallestNonzeroFloat64
    56  // 	} else {
    57  // 		p.e += macheps * math.Abs(p.n) / (1 - macheps)
    58  // 	}
    59  // 	return p
    60  // }
    61  
    62  // addFloat64Pred adds b to a and updates the potential error
    63  func (a float64Pred) addFloat64Pred(b float64Pred) float64Pred {
    64  	a.n += b.n
    65  	if a.n == 0 {
    66  		a.e += math.SmallestNonzeroFloat64 + b.e
    67  	} else {
    68  		a.e += macheps*math.Abs(a.n)/(1-macheps) + b.e
    69  	}
    70  	return a
    71  }
    72  
    73  // // subFloat64 subtracts f from p and updates the potential error
    74  // func (p float64Pred) subFloat64(f float64) float64Pred {
    75  // 	p.n -= f
    76  // 	if p.n == 0 {
    77  // 		p.e += math.SmallestNonzeroFloat64
    78  // 	} else {
    79  // 		p.e += macheps * math.Abs(p.n) / (1 - macheps)
    80  // 	}
    81  // 	return p
    82  // }
    83  
    84  // subFloat64Pred subtracts f from p and updates the potential error
    85  func (a float64Pred) subFloat64Pred(b float64Pred) float64Pred {
    86  	a.n -= b.n
    87  	if a.n == 0 {
    88  		a.e += math.SmallestNonzeroFloat64 + b.e
    89  	} else {
    90  		a.e += macheps*math.Abs(a.n)/(1-macheps) + b.e
    91  	}
    92  	return a
    93  }
    94  
    95  // mulFloat64 multiplies p with f and updates the potential error.
    96  //
    97  // mul(mul(a,b),c) = mul(a*b+error,c) = a*b*c + error*c + newError
    98  // sum(mul(a,b),c) = sum(a*b+error,c) = a*b+c + error + newError
    99  // Conclusively, when multiplications are chained, the error also depends on the value
   100  // of the number, but this does not apply to sums or subtractions.
   101  //
   102  // If this is not a chained multiplication p.e will be 0, making that part irrelevant.
   103  func (p float64Pred) mulFloat64(f float64) float64Pred {
   104  	p.n *= f
   105  	if p.n == 0 {
   106  		p.e += math.SmallestNonzeroFloat64 + p.e*math.Abs(f)
   107  	} else {
   108  		p.e += macheps*math.Abs(p.n)/(1-macheps) + p.e*math.Abs(f)
   109  	}
   110  	return p
   111  }