1// Copyright 2018 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 5package math 6 7import ( 8 "math/bits" 9) 10 11// reduceThreshold is the maximum value where the reduction using Pi/4 12// in 3 float64 parts still gives accurate results. Above this 13// threshold Payne-Hanek range reduction must be used. 14const reduceThreshold = (1 << 52) / (4 / Pi) 15 16// trigReduce implements Payne-Hanek range reduction by Pi/4 17// for x > 0. It returns the integer part mod 8 (j) and 18// the fractional part (z) of x / (Pi/4). 19// The implementation is based on: 20// "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" 21// K. C. Ng et al, March 24, 1992 22// The simulated multi-precision calculation of x*B uses 64-bit integer arithmetic. 23func trigReduce(x float64) (j uint64, z float64) { 24 const PI4 = Pi / 4 25 if x < PI4 { 26 return 0, x 27 } 28 // Extract out the integer and exponent such that, 29 // x = ix * 2 ** exp. 30 ix := Float64bits(x) 31 exp := int(ix>>shift&mask) - bias - shift 32 ix &^= mask << shift 33 ix |= 1 << shift 34 // Use the exponent to extract the 3 appropriate uint64 digits from mPi4, 35 // B ~ (z0, z1, z2), such that the product leading digit has the exponent -61. 36 // Note, exp >= -53 since x >= PI4 and exp < 971 for maximum float64. 37 digit, bitshift := uint(exp+61)/64, uint(exp+61)%64 38 z0 := (mPi4[digit] << bitshift) | (mPi4[digit+1] >> (64 - bitshift)) 39 z1 := (mPi4[digit+1] << bitshift) | (mPi4[digit+2] >> (64 - bitshift)) 40 z2 := (mPi4[digit+2] << bitshift) | (mPi4[digit+3] >> (64 - bitshift)) 41 // Multiply mantissa by the digits and extract the upper two digits (hi, lo). 42 z2hi, _ := bits.Mul64(z2, ix) 43 z1hi, z1lo := bits.Mul64(z1, ix) 44 z0lo := z0 * ix 45 lo, c := bits.Add64(z1lo, z2hi, 0) 46 hi, _ := bits.Add64(z0lo, z1hi, c) 47 // The top 3 bits are j. 48 j = hi >> 61 49 // Extract the fraction and find its magnitude. 50 hi = hi<<3 | lo>>61 51 lz := uint(bits.LeadingZeros64(hi)) 52 e := uint64(bias - (lz + 1)) 53 // Clear implicit mantissa bit and shift into place. 54 hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1))) 55 hi >>= 64 - shift 56 // Include the exponent and convert to a float. 57 hi |= e << shift 58 z = Float64frombits(hi) 59 // Map zeros to origin. 60 if j&1 == 1 { 61 j++ 62 j &= 7 63 z-- 64 } 65 // Multiply the fractional part by pi/4. 66 return j, z * PI4 67} 68 69// mPi4 is the binary digits of 4/pi as a uint64 array, 70// that is, 4/pi = Sum mPi4[i]*2^(-64*i) 71// 19 64-bit digits and the leading one bit give 1217 bits 72// of precision to handle the largest possible float64 exponent. 73var mPi4 = [...]uint64{ 74 0x0000000000000001, 75 0x45f306dc9c882a53, 76 0xf84eafa3ea69bb81, 77 0xb6c52b3278872083, 78 0xfca2c757bd778ac3, 79 0x6e48dc74849ba5c0, 80 0x0c925dd413a32439, 81 0xfc3bd63962534e7d, 82 0xd1046bea5d768909, 83 0xd338e04d68befc82, 84 0x7323ac7306a673e9, 85 0x3908bf177bf25076, 86 0x3ff12fffbc0b301f, 87 0xde5e2316b414da3e, 88 0xda6cfd9e4f96136e, 89 0x9e8c7ecd3cbfd45a, 90 0xea4f758fd7cbe2f6, 91 0x7a0e73ef14a525d4, 92 0xd7f6bf623f1aba10, 93 0xac06608df8f6d757, 94} 95