import { BigFloatBase } from './BigFloatBase' ; import { BigFloat32 } from './BigFloat32' ; export const dekkerSplitter = ( 1 << 27 ) + 1 ; export const limbsPerDigit53 = Math . log ( 10 ) / ( 53 * Math . log ( 2 ) ) ; /** See Shewchuk page 7. */ /** Error-free addition of two floating point numbers.
  * See Shewchuk page 8. Note that output order is swapped! */ function twoSum ( a , b , sum ) { const estimate = a + b ; const b2 = estimate - a ; const a2 = estimate - b2 ; sum [ 0 ] = ( a - a2 ) + ( b - b2 ) ; sum [ 1 ] = estimate ; return ( sum ) ; } return function twoProduct ( a , b ) { tempProduct = a * b ; const a2 = a * dekkerSplitter ; const aHi = a2 - ( a2 - a ) ; const aLo = a - aHi ; const b2 = b * dekkerSplitter ; const bHi = b2 - ( b2 - b ) ; const bLo = b - bHi ; return ( aLo * bLo - ( tempProduct - aHi * bHi - aLo * bHi - aHi * bLo ) ) ; } /** Arbitrary precision floating point number. Based on a multiple-component
  * expansion format and error free transformations.
  *
  * Maximum exponent is the same as for plain JavaScript numbers,
  * least significant representable binary digit is 2^-1074. */ export class BigFloat53 implements BigFloatBase < BigFloat53 > { param value constructor ( value ? , base ? ) { if ( value ) this . setValue ( value , base ) ; } clone ( ) { return ( new BigFloat53 ( ) . setBig ( this ) ) ; } return setZero ( ) { this . len = 0 ; return ( this ) ; } setValue ( other , base ? ) { if ( typeof ( other ) == 'number' ) { return ( this . setNumber ( other ) ) ; } if ( other instanceof BigFloat53 ) { return ( this . setBig ( other ) ) ; } return ( this . setString ( other . toString ( ) , base || 10 ) ) ; } private setBig ( other ) { const len = other . len ; this . len = len ; for ( let pos = 0 ; pos < len ; ++ pos ) { this . limbList [ pos ] = other . limbList [ pos ] ; } return ( this ) ; } param value return private setNumber ( value ) { this . limbList [ 0 ] = value ; this . len = value && 1 ; return ( this ) ; } private setString ( value , base ) { temp32 [ 0 ] . setValue ( value , base ) ; this . len = temp32 [ 0 ] . getExpansion ( this . limbList ) ; this . normalize ( ) ; return ( this ) ; } param a param b return setSum ( a , b ) { this . len = 2 ; twoSum ( a , b , this . limbList ) ; return ( this ) ; } param a param b return setProduct ( a , b ) { this . len = 2 ; this . limbList [ 0 ] = twoProduct ( a , b ) ; this . limbList [ 1 ] = tempProduct ; return ( this ) ; } /** See Compress from Shewchuk page 25. */ normalize ( ) { const limbList = this . limbList ; let len = this . len ; let limb ; if ( len ) { let a = len - 1 ; let b = len - 1 ; let q = limbList [ a ] ; let err ; while ( a ) { limb = limbList [ -- a ] ; err = q ; q += limb ; err = limb - ( q - err ) ; limbList [ b ] = q ; b -= err && 1 ; q = err || q ; } limbList [ b ] = q ; while ( ++ b < len ) { limb = limbList [ b ] ; err = q ; q += limb ; err -= q - limb ; limbList [ a ] = err ; a += err && 1 ; } limbList [ a ] = q ; this . len = a + ( q && 1 ) ; } return ( this ) ; } param b param product return mulSmall ( b , product ) { const limbList = this . limbList ; const productLimbs = product . limbList ; const count = this . len ; let t1 , t2 , t3 ; let srcPos = 0 , dstPos = 0 ; /** Write output limb and move to next, unless a zero was written. */ function writeLimb ( limb ) { productLimbs [ dstPos ] = limb ; dstPos += limb && 1 ; } writeLimb ( twoProduct ( limbList [ srcPos ++ ] , b ) ) ; let q = tempProduct ; while ( srcPos < count ) { t1 = twoProduct ( limbList [ srcPos ++ ] , b ) ; t2 = q + t1 ; t3 = t2 - q ; writeLimb ( q - ( t2 - t3 ) + ( t1 - t3 ) ) ; q = tempProduct + t2 ; writeLimb ( t2 - ( q - tempProduct ) ) ; } productLimbs [ dstPos ] = q ; product . len = dstPos + ( q && 1 ) ; return ( product ) ; } param multiplier param product return private mulBig ( multiplier , product ) { const limbList = multiplier . limbList ; let pos = multiplier . len ; if ( ! pos ) return ( product . setZero ( ) ) ; -- pos ; this . mulSmall ( limbList [ pos ] , pos ? temp53 [ pos & 1 ] : product ) ; while ( pos ) { -- pos ; this . mulSmall ( limbList [ pos ] , product ) . addBig ( temp53 [ ~ pos & 1 ] , 1 , pos ? temp53 [ pos & 1 ] : product ) ; } return ( product ) ; } param multiplier param product return mul ( multiplier , product ? ) { product = product || new BigFloat53 ( ) ; if ( typeof ( multiplier ) == 'number' ) { return ( this . mulSmall ( multiplier , product ) ) ; } if ( product == this ) throw ( new Error ( 'Cannot multiply in place' ) ) ; return ( this . mulBig ( multiplier , product ) ) ; } cmp : ( other : number | BigFloat53 ) => number ; isZero ( ) { const limbList = this . limbList ; let pos = this . len ; while ( pos -- ) { if ( limbList [ pos ] ) return ( false ) ; } return ( true ) ; } getSign ( ) { let t = this . len ; return ( t && ( t = this . limbList [ t - 1 ] ) && ( t > 0 ? 1 : - 1 ) ) ; } /** Return an arbitrary number with sign matching the result of this - other. */ deltaFrom ( other ) { let t = this . len ; let sign = this . getSign ( ) ; let diff = sign ; if ( typeof ( other ) != 'number' ) { t = other . len ; diff -= t && ( t = other . limbList [ t - 1 ] ) && ( t > 0 ? 1 : - 1 ) ; if ( diff || ! sign ) return ( diff ) ; this . addBig ( other , - 1 , temp53 [ 0 ] ) ; } else { diff -= other && ( other > 0 ? 1 : - 1 ) ; if ( diff || ! sign ) return ( diff ) ; this . addSmall ( - other , temp53 [ 0 ] ) ; } t = temp53 [ 0 ] . len ; return ( t && temp53 [ 0 ] . limbList [ t - 1 ] ) ; } param b param sum return addSmall ( b , sum ) { const limbList = this . limbList ; const sumLimbs = sum . limbList ; const count = this . len ; let estimate ; let a , a2 , b2 , err ; let srcPos = 0 , dstPos = 0 ; while ( srcPos < count ) { a = limbList [ srcPos ++ ] ; estimate = a + b ; b2 = estimate - a ; a -= estimate - b2 ; err = a + ( b - b2 ) ; sumLimbs [ dstPos ] = err ; dstPos += err && 1 ; b = estimate ; } sumLimbs [ dstPos ] = b ; sum . len = dstPos + ( b && 1 ) ; return ( sum ) ; } param sign param sum return private addBig ( addend , sign , sum ) { const augendLimbs = this . limbList ; const addendLimbs = addend . limbList ; const sumLimbs = sum . limbList ; let count = this . len + addend . len ; let nextAugendPos = 0 ; let nextAddendPos = 0 ; let nextSumPos = 0 ; /** Latest limb of augend. */ let a = augendLimbs [ nextAugendPos ++ ] ; /** Latest limb of addend. */ let b = addendLimbs [ nextAddendPos ++ ] * sign ; /** Magnitude of latest augend limb. */ let a2 = a < 0 ? - a : a ; /** Magnitude of latest addend limb. */ let b2 = b < 0 ? - b : b ; let nextLimb , nextLimb2 , prevLimb ; let err ; if ( ! count ) return ( sum . setZero ( ) ) ; augendLimbs [ this . len ] = Infinity ; addendLimbs [ addend . len ] = Infinity ; /** Get next smallest limb from either augend or addend.
          * This avoids merging the two limb lists. */ function getNextLimb ( ) { let result ; if ( a2 < b2 ) { result = a ; a = augendLimbs [ nextAugendPos ++ ] ; a2 = a < 0 ? - a : a ; } else { result = b ; b = addendLimbs [ nextAddendPos ++ ] * sign ; b2 = b < 0 ? - b : b ; } return ( result ) ; } let limb = getNextLimb ( ) ; while ( -- count ) { nextLimb = getNextLimb ( ) ; prevLimb = limb ; limb += nextLimb ; nextLimb2 = limb - prevLimb ; err = ( prevLimb - ( limb - nextLimb2 ) ) + ( nextLimb - nextLimb2 ) ; sumLimbs [ nextSumPos ] = err ; nextSumPos += err && 1 ; } sumLimbs [ nextSumPos ] = limb ; sum . len = nextSumPos + ( limb && 1 ) ; return ( sum ) ; } private addSub ( addend , sign , result ? ) { result = result || new BigFloat53 ( ) ; if ( typeof ( addend ) == 'number' ) return ( this . addSmall ( sign * addend , result ) ) ; return ( this . addBig ( addend , sign , result ) ) ; } param addend param sum return add ( addend , sum ? ) { return ( this . addSub ( addend , 1 , sum ) ) ; } param subtrahend param difference return sub ( subtrahend , difference ? ) { return ( this . addSub ( subtrahend , - 1 , difference ) ) ; } /** Round towards zero, to (at least) given number of base 2^53 fractional digits. */ truncate ( fractionLimbCount ) { this . normalize ( ) ; const limbList = this . limbList ; let len = this . len ; let lo = 0 ; let hi = len ; let mid = 0 ; let limb = 0 ; while ( lo < hi ) { mid = ( lo + hi ) >> 1 ; limb = limbList [ mid ] ; if ( limb > - 1 && limb < 1 ) { lo = mid + 1 ; } else { hi = mid ; } } if ( mid && ( limb <= - 1 || limb >= 1 ) ) { limb = limbList [ -- mid ] ; } mid -= fractionLimbCount - 1 ; if ( mid > 0 ) { this . len -= mid ; len = this . len ; for ( let pos = 0 ; pos < len ; ++ pos ) { limbList [ pos ] = limbList [ pos + mid ] ; } } return ( this ) ; } round ( decimalCount ) { return ( this . truncate ( 1 + ~ ~ ( decimalCount * limbsPerDigit53 ) ) ) ; } valueOf ( ) { const limbList = this . limbList ; const len = this . len ; let result = 0 ; for ( let pos = 0 ; pos < len ; ++ pos ) { result += limbList [ pos ] ; } return ( result ) ; } return toString ( base ? ) { const limbList = this . limbList ; let pos = this . len ; temp32 [ pos & 1 ] . setZero ( ) ; while ( pos -- ) { temp32 [ ~ pos & 1 ] . add ( limbList [ pos ] , temp32 [ pos & 1 ] ) ; } return ( temp32 [ ~ pos & 1 ] . toString ( base ) ) ; } /** List of components ordered by increasing exponent. */ private limbList : number [ ] = [ ] ; private len : number ; } BigFloat53 . prototype . cmp = BigFloat53 . prototype . deltaFrom ; /** Latest approximate product from twoProduct. */ let tempProduct = 0 ; /** Temporary values for internal calculations. */ const temp32 = [ new BigFloat32 ( ) , new BigFloat32 ( ) ] ; /** Temporary values for internal calculations. */ const temp53 = [ new BigFloat53 ( ) , new BigFloat53 ( ) ] ;