/*---------------------------------------------------------------------------
  Copyright 2017-2022, Microsoft Research, Daan Leijen.

  This is free software; you can redistribute it and/or modify it under the
  terms of the Apache License, Version 2.0. A copy of the License can be
  found in the LICENSE file at the root of this distribution.
---------------------------------------------------------------------------*/

/* 128-bit double-double floating point numbers.

The `:ddouble` type implements [double double][ddwiki] 128-bit floating point numbers
as a pair of IEEE `:float64` values. This extends the precision to 31 decimal digits
(versus 15 for `:float64`), but keeps the same range as
a `:float64` with a maximum  value of about 1.8·10^308^ (`dd-max`). Because
`:float64`s usually have hardware support, a `:ddouble` is usually much faster
than arbitrary precision floating point numbers.

## Precision {-}

If you add two regular 64-bit `:float64` values you can quickly notice the imprecision
due to decimal fractions that cannot be represented precisely. For example:
```
> 0.1 + 0.2
0.30000000000000004
```

This happens because the decimal `0.1` and `0.2` cannot be represented exactly
by a `:float64` which is encoded in base 2. For example, if we show `0.1` to 20 digits, we get:
```
> 0.1.show(20)
"0.10000000000000000555"
```

When we convert the constant `0.1` to a `:ddouble` we can see with even more precision how it is
approximated:

```
> 0.1.ddouble
0.1000000000000000055511151231258
```

However, if we convert the number `0.1` directly to a `:ddouble` (instead of going through
a `:float64` constant), we can represent `0.1` more precisely:

```
> ddouble("0.1")   // for convenience; using `ddouble-exp(1,-1)` is more efficient
0.1
```

This is possible because a `:ddouble` uses two `:float64`s to represent numbers; in this
case the first (_hi_) float64 is (the inexact) `0.1` while the second one (_lo_) is the
correction to that number. The `show-sum`  (or `show-hex`) function shows this internal pair of numbers:

```
> ddouble("0.1").show-sum(20)
"0.10000000000000000555 + -5.55111512312578301027e-18"

> ddouble("0.1") + ddouble("0.2")
0.3
```

Generally, a `:ddouble` _d_ is represented by a pair of `:float64`s, _hi_ and _lo_,
such that _d_ equals _hi_+_lo_, where \|_lo_\| ≤ 0.5·ulp(_hi_).
Note that despite the extra precision, underflow/overflow and machine precision are
not as well-defined as with a regular `:float64` [@lin;@shewchuk;@Dekker].
Nevertheless, using a `:ddouble` often prevents instability in practice over
`:float64` computations due to  rounding errors or combining very large and
small numbers.

Take the "thin triangle" problem for example [@Goldberg:float;@Kahan:triangle].
The challenge is to compute
the area of a very thin triangle with sides $a$, $b$, and $c$. The short sides
 $b$ and $c$ are just 3 units in the last place (ulp) shorter than
the longest side $a$. Using the value $s = (a + b + c)/2$, the area of
a triangle is $A = \sqrt{s(s-a)(s-b)(s-c)}$. This is troublesome if $s$ is
close to $a$ since $s-a$ magnifies any rounding error. Let's take $a=7$,
and $b = c = 7/2 + 3\cdot2^{-111}$; according to [@gustafson:posit],
“If the units are in light-years, then the short sides are only longer than
half the long side by $1/200$th the diameter of a proton. Yet that pops the triangle
up to about the width of a doorway at the thickest point”. Using various
128-bit numbers we get:
````
exact          : 3.14784204874900425235885265494550774498...e-16
128-bit ddouble: 3.147842048749004252358852654945\([47]{color:#F88}\)e-16
128-bit posit  : 3.147842048749004252358852654945507744\([39]{color:#F88}\)e-16
128-bit ieee   : 3.\([63481490842332134725920516158057682878]{color:#F88}\)e-16
````
For this kind of example, a `:ddouble` has better precision than a
regular 128-bit IEEE float since it can combine very large and
small values. Note that Kahan [@Kahan:triangle] shows how to rewrite the area equation
to avoid magnifying rounding errors -- in that case the result for IEEE 128-bit floats becomes:
````
128-bit ieee x : 3.147842048749004252358852654945507\([92210]{color:#F88}\)e-16
````

The implementation is based closely on the excellent [QD] C++ library [@Hida:qd;@Hida:qdlib],
and assumes proper 64-bit IEEE `:float64`s with correct rounding.
Integers can be represented precisely up to 30 decimal digits (and a bit more...
up to 2^106^ - 2).

## References {-}

~ Bibliography { caption:"0" }

~~ BibItem { #dekker; bibitem-label:"[1]"; searchterm:"Dekker+A+Floating-Point+Technique+for+Extending+the+Available+Precision" }
T.\ Dekker.
_A Floating-Point Technique for Extending the Available Precision_.
Numerische\ Mathematik, vol. 18 (3), June 1971, 224--242.
[pdf](https://gdz.sub.uni-goettingen.de/dms/resolveppn/?PPN=GDZPPN001170007).
~~

~~ BibItem { #Goldberg:float; bibitem-label:"[2]"}
David Goldberg.
_What Every Computer Scientist Should Know About Floating-point Arithmetic_.
ACM Computing Survey. vol. 23 (1), March 1991. doi: [10.1145/103162.103163](https://doi.org/10.1145/103162.103163).
[pdf](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.22.6768&rep=rep1&type=pdf).
~~

~~ BibItem { #Gustafson:posit; bibitem-label:"[3]"; searchterm:"Gustafson+John+Beating+floating+point+at+its+own+game+posit+arithmetic"; }
John\ L. Gustafson and Isaac\ Yonemoto.
_Beating Floating Point at its Own Game: Posit Arithmetic_. 2017.
[pdf](https://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf).
~~

~~ BibItem { #Hida:qd; bibitem-label:"[4]"; searchterm:"Hida+Quad-Double+Arithmetic:+Algorithms,+Implementation,+and+Application" }
Yozo Hida,\ Xiaoye S.\ Li, and David\ H.\ Bailey.
_Quad-Double Arithmetic: Algorithms, Implementation, and Application_.
Lawrence Berkeley National Laboratory Technical Report LBNL-46996. 2000. [pdf](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.4.5769).
~~

~~ BibItem { #Hida:qdlib; bibitem-label:"[5]"; searchterm:"Hida+Library+for+double-double+and+quad-double+arithmetic" }
Yozo Hida,\ Xiaoye\ S.\ Li, and David\ H.\ Bailey.
_Library for double-double and quad-double arithmetic._
(2007). [pdf](https://www.jaist.ac.jp/~s1410018/papers/qd.pdf).
~~

~~ Bibitem { #lin; bibitem-label:"[6]"; searchterm:'Linnainmaa+"Software+for+Doubled-Precision+Floating-Point+Computations"' }
Seppo\ Linnainmaa.
_Software for Doubled-Precision Floating-Point Computations_.
ACM Transactions on Mathematical Software (TOMS), vol. 7 (3), Sept. 1981, 272--283.
~~

~~ Bibitem { #Kahan:triangle; bibitem-label:"[7]"}
William\ Kahan.
_Miscalculating Area and Angles of a Needle-like Triangle_.
Lecture notes, 2004, [pdf](https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf).
~~

~~ BibItem { #shewchuk; bibitem-label:"[8]"; }
Jonathan\ Richard\ Shewchuk.
_Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates_.
Discrete & Computational Geometry, vol. 18, 305--363, 1997. [pdf](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
~~

~

[ddwiki]: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic
[qd]:     https://crd-legacy.lbl.gov/~dhbailey/mpdist

\/
*/
module std/num/ddoublestd/num/ddouble

import std/core/undivstd/core/undiv
import std/num/float64std/num/float64
import std/num/decimalstd/num/decimal
import std/text/parsestd/text/parse

/* The `:ddouble` type implements [double double][ddwiki] 128-bit floating point numbers
as a pair of IEEE `:float64` values. This extends the precision to 31 decimal digits
(versus 15 for `:float64`), but keeps the same range as
a `:float64` with a maximum  value of about 1.8·10^308^. Because
`:float64`s usually have hardware support, a `:ddouble` is usually much faster
than arbitrary precision floating point numbers.

Internally a `:ddouble` _d_ is represented as a pair of `:float64`s, _hi_ and _lo_,
such that the number represented by _d_ is _hi_+_lo_, where \|_lo_\| ≤ 0.5·ulp(_hi_).
*/
abstract value struct ddoublestd/num/ddouble/ddouble: V
  histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 : float64std/core/types/float64: V
  lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 : float64std/core/types/float64: V

// Create a `:ddouble` from a `:float64`.
pub fun float64/ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble( dd: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(dd: float64,0.0literal: float64
hex64= 0x0p+0)

val maxprecisestd/num/ddouble/maxprecise: int : intstd/core/types/int: V = 9007199254740991literal: int
dec  = 9007199254740991
hex64= 0x001FFFFFFFFFFFFF
bit64= 0b0000000000011111111111111111111111111111111111111111111111111111
val minprecisestd/num/ddouble/minprecise: int : intstd/core/types/int: V = ~std/core/int/(~): (i : int) -> intmaxprecisestd/num/ddouble/maxprecise: int

fun is-precisestd/num/ddouble/is-precise: (i : int) -> bool(ii: int : intstd/core/types/int: V )result: -> total bool : boolstd/core/types/bool: V
  (ii: int >=std/core/int/(>=): (x : int, y : int) -> bool minprecisestd/num/ddouble/minprecise: int &&std/core/types/(&&): (x : bool, y : bool) -> bool ii: int <=std/core/int/(<=): (x : int, y : int) -> bool maxprecisestd/num/ddouble/maxprecise: int)


// Create a `:ddouble` from an `:int`.
// A `:ddouble` can represent integers precisely up to 30 digits.
// If an integer is passed that is out of range an infinity is returned.
pub fun int/ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble( ii: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> ddouble(ii: int,0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000)

fun ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> ddouble( ii: int : intstd/core/types/int: V, ee: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if ii: int.is-precisestd/num/ddouble/is-precise: (i : int) -> bool then small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( ii: int, ee: int ) else
    // ensure we convert at least 31 (and at most 42) digits precisely
    val pp: int    = ii: int.count-digitsstd/core/int/count-digits: (i : int) -> int
    val xx: int    = ii: int
    val pxpx: int   = pp: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec  = 14
hex8 = 0x0E
bit8 = 0b00001110
    val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)hihi: int,yy: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = xx: int.cdivmod-exp10std/core/int/cdivmod-exp10: (i : int, n : int) -> (int, int)(pxpx: int)
    val pypy: int   = pxpx: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec  = 14
hex8 = 0x0E
bit8 = 0b00001110
    if pypy: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then
      // if 28 or less digits we can convert with two doubles
      // trace("ddouble-i: " + i.show + ", hi:" + hi.show + ",  lo: " + y.show)
      small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble(hihi: int, pxpx: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( yy: int, ee: int)

    else
      // otherwise use 2 initial 14 digits and a trailing value
      val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)midmid: int,zz: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = yy: int.cdivmod-exp10std/core/int/cdivmod-exp10: (i : int, n : int) -> (int, int)(pypy: int)
      val pzpz: int = pypy: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec  = 14
hex8 = 0x0E
bit8 = 0b00001110
      val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)lolo: int,ploplo: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = if pzpz: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)zz: int,0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) else (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)zz: int.cdiv-exp10std/core/int/cdiv-exp10: (i : int, n : int) -> int(pzpz: int), pzpz: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)
      //trace("ddouble-i: " + i.show + ", hi:" + hi.show + ", mid: " + mid.show + ", lo: " + lo.show)
      //trace("  px: " + px.show + ", py : " + py.show + ", plo: " + plo.show + ", e: " + e.show )
      small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( hihi: int, pxpx: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( midmid: int, pypy: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( lolo: int, ploplo: int +std/core/int/(+): (x : int, y : int) -> int ee: int))


fun small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( ii: int : intstd/core/types/int: V, ee: int : intstd/core/types/int: V )result: -> total ddouble
  val dddd: ddouble = ii: int.float64std/num/float64/float64: (i : int) -> float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble
  if ee: int.is-zerostd/core/int/is-zero: (x : int) -> bool then dddd: ddouble else dddd: ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(ee: int)

// &pi;
pub val dd-pistd/num/ddouble/dd-pi: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+1literal: float64
hex64= 0x1.921fb54442d18p1, 0x1.1A62633145C07p-53literal: float64
hex64= 0x1.1a62633145c07p-53) // 3.1415926535897932384626433832795028841971693993751
// pub val dd-pi-qd      = Ddouble(3.1415926535897931160,   1.2246467991473529607e-16)  // QD's pi constant is one bit off...

// 2&pi;
pub val dd-twopistd/num/ddouble/dd-twopi: ddouble   = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+2literal: float64
hex64= 0x1.921fb54442d18p2, 0x1.1A62633145C07p-52literal: float64
hex64= 0x1.1a62633145c07p-52)

// &pi;/2
pub val dd-pi2std/num/ddouble/dd-pi2: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+0literal: float64
hex64= 0x1.921fb54442d18p0, 0x1.1A62633145C07p-54literal: float64
hex64= 0x1.1a62633145c07p-54)

// &pi;/4
pub val dd-pi4std/num/ddouble/dd-pi4: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p-1literal: float64
hex64= 0x1.921fb54442d18p-1, 0x1.1A62633145C07p-55literal: float64
hex64= 0x1.1a62633145c07p-55)
// QD:                         7.853981633974482790e-01,3.061616997868383018e-17)

// &pi;/16
val dd-pi16std/num/ddouble/dd-pi16: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0x1.921FB54442D18p+1literal: float64
hex64= 0x1.921fb54442d18p1, 0x1.1A62633145C07p-53literal: float64
hex64= 0x1.1a62633145c07p-53)

// 3&pi;/4
pub val dd-pi34std/num/ddouble/dd-pi34: ddouble    = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.2D97C7F3321D2p+1literal: float64
hex64= 0x1.2d97c7f3321d2p1, 0x1.A79394C9E8A0Bp-54literal: float64
hex64= 0x1.a79394c9e8a0bp-54)
// QD:                          2.356194490192344837e+00,9.1848509936051484375e-17)  // QD is again one bit off

// The _e_ constant.
pub val dd-estd/num/ddouble/dd-e: ddouble       = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.5BF0A8B145769p+1literal: float64
hex64= 0x1.5bf0a8b145769p1, 0x1.4D57EE2B1013Ap-53literal: float64
hex64= 0x1.4d57ee2b1013ap-53) // 2.7182818284590452353602874713526624977572470937

// The natural logarithm of 2
pub val dd-ln2std/num/ddouble/dd-ln2: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.62E42FEFA39EFp-1literal: float64
hex64= 0x1.62e42fefa39efp-1, 0x1.ABC9E3B39803Fp-56literal: float64
hex64= 0x1.abc9e3b39803fp-56) // 0.69314718055994530941723212145817656807550013436026

// The natural logarithm of 10
pub val dd-ln10std/num/ddouble/dd-ln10: ddouble    = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.26BB1BBB55516p+1literal: float64
hex64= 0x1.26bb1bbb55516p1, -0x1.F48AD494EA3E9p-53literal: float64
hex64= -0x1.f48ad494ea3e9p-53) // 2.3025850929940456840179914546843642076011014886288

// The base-2 logarithm of _e_.
pub val dd-log2estd/num/ddouble/dd-log2e: ddouble   = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.71547652B82FEp+0literal: float64
hex64= 0x1.71547652b82fep0, 0x1.777D0FFDA0D23p-56literal: float64
hex64= 0x1.777d0ffda0d23p-56) // 1.442695040888963407359924681001892137426645954153

// The base-10 logarithm of _e_.
pub val dd-log10estd/num/ddouble/dd-log10e: ddouble  = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.BCB7B1526E50Ep-2literal: float64
hex64= 0x1.bcb7b1526e50ep-2, 0x1.95355BAAAFAD3p-57literal: float64
hex64= 0x1.95355baaafad3p-57) // 0.43429448190325182765112891891660508229439700580367

// The square-root of 2
pub val dd-sqrt2std/num/ddouble/dd-sqrt2: ddouble   = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.6A09E667F3BCDp+0literal: float64
hex64= 0x1.6a09e667f3bcdp0, -0x1.BDD3413B26456p-54literal: float64
hex64= -0x1.bdd3413b26456p-54) // 1.4142135623730950488016887242096980785696718753769

// `1.0 / sqrt(2.0)`
pub val dd-sqrt12std/num/ddouble/dd-sqrt12: ddouble  = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.6A09E667F3BCDp-1literal: float64
hex64= 0x1.6a09e667f3bcdp-1, -0x1.BDD3413B26456p-55literal: float64
hex64= -0x1.bdd3413b26456p-55) // 0.70710678118654752440084436210484903928483593768847

// [Euler's constant](https://en.wikipedia.org/wiki/Euler%E2%80%93Mascheroni_constant)
pub val dd-eulerstd/num/ddouble/dd-euler: ddouble   = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.2788CFC6FB619p-1literal: float64
hex64= 0x1.2788cfc6fb619p-1, -0x1.6CB90701FBFABp-58literal: float64
hex64= -0x1.6cb90701fbfabp-58) //0.57721566490153286060651209008240243104215933593992

// The 'machine epsilon': this is not well-defined for a `:ddouble` in general since
// the difference between 1.0 and the next representable `:ddouble` value is `dd-true-min`.
// Instead, we take the square of `flt-epsilon`, i.e. 2^-104^.
pub val dd-epsilonstd/num/ddouble/dd-epsilon: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1p-104literal: float64
hex64= 0x1p-104,0.0literal: float64
hex64= 0x0p+0)

// 8*dd-epsilon
val dd-epsilon8std/num/ddouble/dd-epsilon8: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.0p-101literal: float64
hex64= 0x1p-101, 0x0.0p+0literal: float64
hex64= 0x0p+0)

// The maximum representable `:ddouble`
pub val dd-maxstd/num/ddouble/dd-max: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0x1.FFFFFFFFFFFFFp+1023literal: float64
hex64= 0x1.fffffffffffffp1023, 0x1.FFFFFFFFFFFFFp+969literal: float64
hex64= 0x1.fffffffffffffp969 )

// The smallest positive `:ddouble` that is still normalized
pub val dd-minstd/num/ddouble/dd-min: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( flt-minstd/num/float64/flt-min: float64, 0.0literal: float64
hex64= 0x0p+0 )   // normalized = 2^-1022

// The smallest positive `:ddouble`  (which is subnormal).
pub val dd-true-minstd/num/ddouble/dd-true-min: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( flt-true-minstd/num/float64/flt-true-min: float64, 0.0literal: float64
hex64= 0x0p+0 )

// Not-A-Number
pub val dd-nanstd/num/ddouble/dd-nan: ddouble     = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(nanstd/num/float64/nan: float64,0.0literal: float64
hex64= 0x0p+0)

// Positive infinity
pub val dd-posinfstd/num/ddouble/dd-posinf: ddouble  = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(posinfstd/num/float64/posinf: float64,0.0literal: float64
hex64= 0x0p+0)

// Negative infinity
pub val dd-neginfstd/num/ddouble/dd-neginf: ddouble  = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(neginfstd/num/float64/neginf: float64,0.0literal: float64
hex64= 0x0p+0)

// maximal precision in decimal digits of a `:ddouble`.
pub val dd-max-precstd/num/ddouble/dd-max-prec: int = 31literal: int
dec  = 31
hex8 = 0x1F
bit8 = 0b00011111

val dd-default-precstd/num/ddouble/dd-default-prec: int = -31literal: int
dec  = -31
hex8 = 0xE1
bit8 = 0b11100001

/*------------------------------------------------------
   Compare etc.
------------------------------------------------------*/

// Zero constant
pub val zerostd/num/ddouble/zero: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(0.0literal: float64
hex64= 0x0p+0)

// One
pub val onestd/num/ddouble/one: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(1.0literal: float64
hex64= 0x1p0)

// Ten (`10.ddouble`)
pub val tenstd/num/ddouble/ten: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(10.0literal: float64
hex64= 0x1.4p3)

val twostd/num/ddouble/two: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(2.0literal: float64
hex64= 0x1p1)

// Is this `:ddouble` equal to is-zero
pub fun is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-zerostd/num/float64/is-zero: (d : float64) -> bool

// Is this `:ddouble` negative?
pub fun is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool

// Is this `:ddouble` positive?
pub fun is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-posstd/num/float64/is-pos: (d : float64) -> bool

// Return the sign of the `:ddouble`.
pub fun is-signstd/num/ddouble/is-sign: (x : ddouble) -> order( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total order : orderstd/core/types/order: V
  if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64==std/num/float64/(==): (x : float64, y : float64) -> bool0.0literal: float64
hex64= 0x0p+0 then Eqstd/core/types/Eq: order
  elif xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 <std/num/float64/(<): (x : float64, y : float64) -> bool 0.0literal: float64
hex64= 0x0p+0 then Ltstd/core/types/Lt: order
  else Gtstd/core/types/Gt: order

// Is this `:ddouble` not-a-number?
pub fun is-nanstd/num/ddouble/is-nan: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-nanstd/num/float64/is-nan: (d : float64) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-nanstd/num/float64/is-nan: (d : float64) -> bool

// Is this a finite `:ddouble`? (i.e. not `is-nan` or `is-inf`)
pub fun is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool

// Is this an infinite value.
pub fun is-infstd/num/ddouble/is-inf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-infstd/num/float64/is-inf: (d : float64) -> bool

// Does `x` equal positive infinity?
pub fun is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-posinfstd/num/float64/is-posinf: (d : float64) -> bool

// Does `x` equal negative infinity?
pub fun is-neginfstd/num/ddouble/is-neginf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-neginfstd/num/float64/is-neginf: (d : float64) -> bool

// Compare two `:ddouble` values.
pub fun cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total order : orderstd/core/types/order: V
  match cmpstd/num/float64/cmp: (x : float64, y : float64) -> order(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
    Eqstd/core/types/Eq: order  -> cmpstd/num/float64/cmp: (x : float64, y : float64) -> order(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)
    ordord: order -> ordord: order


pub fun (>)std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool (xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Gtstd/core/types/Gt: order }
pub fun (>=)std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Ltstd/core/types/Lt: order }
pub fun (==)std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Eqstd/core/types/Eq: order }
pub fun (!=)std/num/ddouble/(!=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Eqstd/core/types/Eq: order }
pub fun (<)std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool (xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Ltstd/core/types/Lt: order }
pub fun (<=)std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Gtstd/core/types/Gt: order }

// The minimum of `x` and `y`.
pub fun minstd/num/ddouble/min: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble <=std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool yy: ddouble then xx: ddouble else yy: ddouble

// The maximum of `x` and `y`
pub fun maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble >=std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> bool yy: ddouble then xx: ddouble else yy: ddouble

/*------------------------------------------------------
   Addition
------------------------------------------------------*/

value struct edoublestd/num/ddouble/edouble: V
  numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64 : float64std/core/types/float64: V
  errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 : float64std/core/types/float64: V

// often called `twosum` in literature (see [@shewchuk])
fun sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V
  val zz: float64    = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64
  val diffdiff: float64 = zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64
  val errerr: float64  = (xx: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)
  Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64)

fun dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val zz: float64    = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64
  val diffdiff: float64 = zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64
  val errerr: float64  = (xx: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64)

// As `sum` but with `x.abs >= y.abs`
fun quicksumstd/num/ddouble/quicksum: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V
  val zz: float64   = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64
  val errerr: float64 = yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64)
  Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64)

fun dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool returnreturn: ddouble ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: float64)std/core/types/Unit: ()
  val zz: float64   = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64
  val errerr: float64 = yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64)
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64)

// Add two `:ddouble`s
pub fun (+)std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val z1z1: edouble = sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  val lolo: edouble = sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)
  val e1e1: float64 = z1z1: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 lolo: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64
  val z2z2: edouble = quicksumstd/num/ddouble/quicksum: (x : float64, y : float64) -> edouble(z1z1: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,e1e1: float64)
  val e2e2: float64 = z2z2: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 lolo: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64
  dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(z2z2: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,e2e2: float64)


// Create a `:ddouble` as the sum of two `:float64`'s.
pub fun ddoublestd/num/ddouble/ddouble: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if yy: float64.is-zerostd/num/float64/is-zero: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: float64) else dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble(xx: float64,yy: float64)

// Negate a `:ddouble`.
pub fun (~)std/num/ddouble/(~): (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(~std/num/float64/(~): (f : float64) -> float64xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,~std/num/float64/(~): (f : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)

// Subtract two values.
pub fun (-)std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (~std/num/ddouble/(~): (x : ddouble) -> ddoubleyy: ddouble)

// Return the absolute value.
pub fun absstd/num/ddouble/abs: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoublexx: ddouble else xx: ddouble

// Increment by one.
pub fun incstd/num/ddouble/inc: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble :ddoublestd/num/ddouble/ddouble: V
  xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble

// Decrement by one.
pub fun decstd/num/ddouble/dec: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble :ddoublestd/num/ddouble/ddouble: V
  xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble


/*------------------------------------------------------
   Multiply and divide
------------------------------------------------------*/

/*
val splitter = 0x1.0000002p+27 // 134217729.0 = 2^27 + 1
val splitbound  =  0x1.0p996 // 6.696928794914171e+299  = 2^996   // Note, QD seems one bit off
val two28 = 0x1.0p28 // 268435456.0 = 2^28

fun split( x : float64 ) : (float64,float64)
  if x > splitbound || x < ~splitbound
    val y = x * 0x1.0p-28  // 3.7252902984619140625e-09 = 2^-28
    val t = y * splitter
    val hi = t - (t - y)
    val lo = y - hi
    (hi*two28, lo*two28)

  else
    val t = x * splitter
    val hi = t - (t - x)
    val lo = x - hi
    (hi,lo)


// often called `twoproduct` in literature (see [@shewchuk])
fun prod( x : float64, y : float64 ) : edouble
  val z = x*y
  val (xhi,xlo) = split(x)
  val (yhi,ylo) = split(y)
  val err = ((xhi*yhi - z) + (xhi*ylo + xlo*yhi)) + (xlo*ylo)
  Edouble(z,if z.is-finite then err else z)

fun prodsqr( x : float64 ) : edouble
  val z = x*x
  val (h,l) = split(x)
  val err = ((h*h - z) + (2.0*h*l)) + (l*l)
  Edouble(z,if z.is-finite then err else z)

*/

// often called `twoproduct` in literature (see [@shewchuk])
fun prodstd/num/ddouble/prod: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V
  val zz: float64 = xx: float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: float64
  val errerr: float64 = fmaddstd/num/float64/fmadd: (x : float64, y : float64, z : float64) -> float64(xx: float64,yy: float64,~std/num/float64/(~): (f : float64) -> float64zz: float64)
  Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64,errerr: float64)

fun prodsqrstd/num/ddouble/prodsqr: (x : float64) -> edouble( xx: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V
  val zz: float64 = xx: float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: float64
  val errerr: float64 = fmaddstd/num/float64/fmadd: (x : float64, y : float64, z : float64) -> float64(xx: float64,xx: float64,~std/num/float64/(~): (f : float64) -> float64zz: float64)
  Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64,errerr: float64)

// Multiply two `:ddouble`s
pub fun (*)std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val zz: edouble  = prodstd/num/ddouble/prod: (x : float64, y : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64, yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  val ee: float64  = zz: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 (xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(zz: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64, ee: float64)

// Multiply `x` with itself.
pub fun sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val zz: edouble = prodsqrstd/num/ddouble/prodsqr: (x : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  val ee: float64 = (zz: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 (2.0literal: float64
hex64= 0x1p1*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)
  dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(zz: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,ee: float64)

// Divide two `:ddouble`s
pub fun (/)std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val q1q1: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  if !std/core/types/bool/(!): (b : bool) -> boolq1q1: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool !std/core/types/bool/(!): (b : bool) -> boolyy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool returnreturn: ddouble q1q1: ddouble
  val r1r1: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble q1q1: ddouble)
  val q2q2: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(r1r1: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  val r2r2: ddouble = r1r1: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble q2q2: ddouble)
  val q3q3: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(r2r2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  val qq: ddouble  = dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(q1q1: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,q2q2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
  qq: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble q3q3: ddouble

// Remainder of two `:ddouble`s
pub fun (%)std/num/ddouble/(%): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val nn: ddouble = (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble
  xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (nn: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleyy: ddouble)

// Division and remainder of two `:ddouble`s
pub fun divremstd/num/ddouble/divrem: (x : ddouble, y : ddouble) -> (ddouble, ddouble)( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (ddouble, ddouble) : (std/core/types/tuple2: (V, V) -> Vddoublestd/num/ddouble/ddouble: V,ddoublestd/num/ddouble/ddouble: V)
  val nn: ddouble = (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble
  (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)nn: ddouble, xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (nn: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleyy: ddouble))std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)

/*------------------------------------------------------
   Rounding
------------------------------------------------------*/

// Convert a `:ddouble` to a `:float64` (losing precision)
pub fun float64std/num/ddouble/float64: (x : ddouble) -> float64( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total float64 : float64std/core/types/float64: V
  xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64

// Convert a `:ddouble` to the nearest integer (rounding to the nearest even number in case of a tie)
pub fun intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nonfinnonfin: ? int : intstd/core/types/int: V = 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 )result: -> total int : intstd/core/types/int: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then nonfinnonfin: int else xx: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.intstd/num/decimal/int: (x : decimal, rnd : ? round) -> int

// Round a `:ddouble` to the nearest integer, rounding to the nearest even number in case of a tie.
pub fun roundstd/num/ddouble/round: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.roundstd/num/float64/round: (d : float64) -> float64
  val diffdiff: float64 = rr: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64
  if diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool 0.0literal: float64
hex64= 0x0p+0 then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.roundstd/num/float64/round: (d : float64) -> float64)
  elif diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool 0.5literal: float64
hex64= 0x1p-1 &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 1.0literal: float64
hex64= 0x1p0)
  elif diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool -0.5literal: float64
hex64= -0x1p-1 &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-posstd/num/float64/is-pos: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 1.0literal: float64
hex64= 0x1p0)
  else ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64)

// Round to negative infinity.
pub fun floorstd/num/ddouble/floor: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.floorstd/num/float64/floor: (d : float64) -> float64
  if rr: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64
   then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.floorstd/num/float64/floor: (d : float64) -> float64)
   else Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(rr: float64,0.0literal: float64
hex64= 0x0p+0)

// Round to positive infinity.
pub fun ceilingstd/num/ddouble/ceiling: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.ceilingstd/num/float64/ceiling: (d : float64) -> float64
  if rr: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64
   then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.ceilingstd/num/float64/ceiling: (d : float64) -> float64)
   else Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(rr: float64,0.0literal: float64
hex64= 0x0p+0)

// Round towards zero.
pub fun truncatestd/num/ddouble/truncate: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then xx: ddouble.ceilingstd/num/ddouble/ceiling: (x : ddouble) -> ddouble else xx: ddouble.floorstd/num/ddouble/floor: (x : ddouble) -> ddouble

// The fraction of `x` such that `x.truncate + x.fraction == x`.
pub fun fractionstd/num/ddouble/fraction: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.truncatestd/num/ddouble/truncate: (x : ddouble) -> ddouble

// The _floored_ fraction of `x`. This is always positive, such that `x.floor + x.ffraction == x`.
pub fun ffractionstd/num/ddouble/ffraction: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.floorstd/num/ddouble/floor: (x : ddouble) -> ddouble

// Round a `:ddouble` to a specified precision.
// Uses `round` if the precision is smaller or equal to zero.
pub fun round-to-precstd/num/ddouble/round-to-prec: (x : ddouble, prec : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: int : intstd/core/types/int: V  )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if precprec: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then xx: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble
  elif precprec: int >std/core/int/(>): (x : int, y : int) -> bool dd-max-precstd/num/ddouble/dd-max-prec: int then xx: ddouble
  else
    val pp: ddouble : ddoublestd/num/ddouble/ddouble: V = powi10std/num/ddouble/powi10: (exp : int) -> ddouble(precprec: int)
    (xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble pp: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble pp: ddouble


/*------------------------------------------------------
   Exponents (needed for showing ddouble's)
------------------------------------------------------*/

// 'Load exponent': returns `x`&middot;2^`exp`^.
pub fun ldexpstd/num/ddouble/ldexp: (x : ddouble, exp : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(ldexpstd/num/float64/ldexp: (x : float64, e : int) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,expexp: int),ldexpstd/num/float64/ldexp: (x : float64, e : int) -> float64(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,expexp: int))

fun npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, accacc: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if nn: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then accacc: ddouble
  elif nn: int.is-evenstd/core/int/is-even: (i : int) -> bool then npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xx: ddouble), accacc: ddouble, pretend-decreasingstd/core/undiv/pretend-decreasing: (x : int) -> int(nn: int/std/core/int/(/): (x : int, y : int) -> int2literal: int
dec  = 2
hex8 = 0x02
bit8 = 0b00000010) )
  else npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( xx: ddouble, xx: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleaccacc: ddouble, pretend-decreasingstd/core/undiv/pretend-decreasing: (x : int) -> int(nn: int.decstd/core/int/dec: (i : int) -> int))

fun npwrstd/num/ddouble/npwr: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if nn: int==std/core/int/(==): (x : int, y : int) -> bool0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then (if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then dd-nanstd/num/ddouble/dd-nan: ddouble else onestd/num/ddouble/one: ddouble)
  elif nn: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec  = 1
hex8 = 0x01
bit8 = 0b00000001 then xx: ddouble
  else xx: ddouble.npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble(onestd/num/ddouble/one: ddouble,nn: int)

// Return `x` to the power of `n`.
fun powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val pp: ddouble = npwrstd/num/ddouble/npwr: (x : ddouble, n : int) -> ddouble(xx: ddouble,nn: int.absstd/core/int/abs: (i : int) -> int)
  if nn: int.is-negstd/core/int/is-neg: (i : int) -> bool then (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble pp: ddouble) else pp: ddouble

// Return 10 to the power of `exp`.
fun powi10std/num/ddouble/powi10: (exp : int) -> ddouble( expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble(tenstd/num/ddouble/ten: ddouble,expexp: int)

// Return 10 to the power of `exp`.
pub fun exp10std/num/ddouble/exp10: (exp : ddouble) -> ddouble( expexp: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble(tenstd/num/ddouble/ten: ddouble,expexp: ddouble)

// Return 2 to the power of `exp`.
pub fun exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble( expexp: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble(twostd/num/ddouble/two: ddouble,expexp: ddouble)

// Create a `:ddouble` `x` such that `x` equals `d`&middot;10^`e`^.
pub fun float64/ddouble-expstd/num/ddouble/float64/ddouble-exp: (d : float64, e : int) -> ddouble( dd: float64 : float64std/core/types/float64: V, ee: int : intstd/core/types/int: V )result: -> total ddouble
  dd: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(ee: int)

// Create a `:ddouble` `x` such that `x` equals `i`&middot;10^`e`^.
pub fun int/ddouble-expstd/num/ddouble/int/ddouble-exp: (i : int, exp : int) -> ddouble( ii: int : intstd/core/types/int: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  ii: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(expexp: int)

fun mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if expexp: int.is-zerostd/core/int/is-zero: (x : int) -> bool then xx: ddouble else xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble powi10std/num/ddouble/powi10: (exp : int) -> ddouble(expexp: int)

/*------------------------------------------------------
   Show
------------------------------------------------------*/

// Decode a `:ddouble` `d` into two doubles `(hi,lo)` such that `d` equals  `hi`+`lo`,
// where `lo` &le; 0.5&middot;ulp(`hi`).
pub fun decodestd/num/ddouble/decode: (d : ddouble) -> (float64, float64)( dd: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (float64, float64) : (std/core/types/tuple2: (V, V) -> Vfloat64std/core/types/float64: V,float64std/core/types/float64: V)
  (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)dd: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,dd: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)

// Encode a `:ddouble` `d` from two doubles `(hi,lo)` such that `d` equals  `hi`+`lo`.
pub fun encodestd/num/ddouble/encode: (hi : float64, lo : float64) -> ddouble( hihi: float64 : float64std/core/types/float64: V, lolo: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble lolo: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble

// Convert a `:ddouble` to a `:decimal` up to a given precision `prec` (= `-1`).
// A negative precision converts precisely. Returns 0 for non-finite `:ddouble`'s.
pub fun decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = -1literal: int
dec  = -1
hex8 = 0xFF
bit8 = 0b11111111 )result: -> total decimal : decimalstd/num/decimal/decimal: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000.decimalstd/num/decimal/int/decimal: (i : int, exp : ? int) -> decimal
    else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.decimalstd/num/decimal/float64/decimal: (d : float64, max-prec : ? int) -> decimal(precprec: int) +std/num/decimal/(+): (x : decimal, y : decimal) -> decimal xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.decimalstd/num/decimal/float64/decimal: (d : float64, max-prec : ? int) -> decimal(precprec: int)

/* Show a `:ddouble` `x` precisely as the sum of two `:float64`s
in [hexadecimal notation](https://books.google.com/books?id=FgMsCwAAQBAJ&pg=PA41).
Use this if you need to guarantee that you can parse back `:ddouble`s exactly,
i.e. `x == x.show-hex.ddouble`.
```
> 0.1.ddouble.show-hex
"0x1.999999999999Ap-4 + 0x0.0p+0"
> "0.1".ddouble.show-hex
"0x1.999999999999Ap-4 + -0x1.999999999999Ap-58"
> dd-pi.show-hex
"0x1.921FB54442D18p+1 + 0x1.1A62633145C07p-53"
> dd-max.show-hex
"0x1.FFFFFFFFFFFFFp+1023 + 0x1.FFFFFFFFFFFFFp+969"
```
.
*/
pub fun show-hexstd/num/ddouble/show-hex: (x : ddouble, width : ? int, use-capitals : ? bool, pre : ? string) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, widthwidth: ? int : intstd/core/types/int: V = 1literal: int
dec  = 1
hex8 = 0x01
bit8 = 0b00000001, use-capitalsuse-capitals: ? bool : boolstd/core/types/bool: V = Truestd/core/types/True: bool, prepre: ? string : stringstd/core/types/string: V = "0x"literal: string
count= 2 )result: -> total string : stringstd/core/types/string: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string
    else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.show-hexstd/num/float64/show-hex: (d : float64, width : ? int, use-capitals : ? bool, pre : ? string) -> string(widthwidth: int,use-capitalsuse-capitals: bool,prepre: string) ++std/core/types/(++): (x : string, y : string) -> string " + "literal: string
count= 3 ++std/core/types/(++): (x : string, y : string) -> string xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.show-hexstd/num/float64/show-hex: (d : float64, width : ? int, use-capitals : ? bool, pre : ? string) -> string(widthwidth: int,use-capitalsuse-capitals: bool,prepre: string)

// Show a `:ddouble` `x` with a given precision `prec` (=`-31`).
// The precision specifies the  number of digits after the dot (in either scientific of fixed-point notation).
// If the precision is negative, _at most_ `prec` digits are displayed, while for a positive
// precision, exactly `prec` digits behind the dot are displayed.
// This uses `show-fixed` when the exponent of `x` in scientific notation is larger than -5
// and smaller than the precision (or 15 in case of a negative precision), otherwise it uses `show-exp`.
pub fun showstd/num/ddouble/show: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string
    else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.showstd/num/decimal/show: (d : decimal, prec : ? int) -> string(precprec: int)

/* Show a ddouble `x` with a given precision `prec` (=`-31`) in scientific notation.
The precision specifies the  number of digits after the dot, i.e.
the number of significant digits is `prec+1`.
If the precision is negative, _at most_ `prec` digits are displayed, and if
it is positive exactly `prec` digits are used.
```
> 1.1.ddouble.show-exp
"1.1000000000000000888178419700125"
> 1.1.ddouble.show-exp(-100)
"1.100000000000000088817841970012523233890533447265625"
> 1.1.ddouble.show-exp(5)
"1.10000"
> 1.1.ddouble.show-exp(-5)
"1.1"
```
.
*/
pub fun show-expstd/num/ddouble/show-exp: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string
    else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.show-expstd/num/decimal/show-exp: (d : decimal, prec : ? int) -> string(precprec: int)

/* Show a ddouble `x` with a given precision `prec` (=`-31`) in fixed-point notation.
The precision specifies the  number of digits after the dot.
If the precision is negative, _at most_  `prec` digits after the dot are displayed,
while for a positive precision, exactly `prec` digits are used.
```
> 0.1.ddouble.show-fixed
"0.1000000000000000055511151231258"
> 0.1.ddouble.show-fixed(-100)
"0.1000000000000000055511151231257827021181583404541015625"
> 0.1.ddouble.show-fixed(5)
"0.10000"
> 0.1.ddouble.show-fixed(-5)
"0.1"
```
.
*/
pub fun show-fixedstd/num/ddouble/show-fixed: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string
    else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.show-fixedstd/num/decimal/show-fixed: (d : decimal, prec : ? int) -> string(precprec: int)

// Show a `:ddouble` as the sum of  `:float64`s with an optional precision.
// Note: use `show-hex` for reliable round-trip parsing.
pub fun show-sumstd/num/ddouble/show-sum: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = -17literal: int
dec  = -17
hex8 = 0xEF
bit8 = 0b11101111 )result: -> total string : stringstd/core/types/string: V
  if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool
    then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string
    else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string(precprec: int) ++std/core/types/(++): (x : string, y : string) -> string " + "literal: string
count= 3 ++std/core/types/(++): (x : string, y : string) -> string xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string(precprec: int)


/*------------------------------------------------------
  Parsing
------------------------------------------------------*/

// Parse a floating point number with up to 31 digits precision.
// Return `dd-nan` if the string is an invalid number.
pub fun string/ddoublestd/num/ddouble/string/ddouble: (s : string) -> ddouble( ss: string : stringstd/core/types/string: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  ss: string.parse-ddoublestd/num/ddouble/parse-ddouble: (s : string) -> maybe<ddouble>.defaultstd/core/maybe/default: (m : maybe<ddouble>, nothing : ddouble) -> ddouble(dd-nanstd/num/ddouble/dd-nan: ddouble)

/*
// Parse a floating point number with up to 31 digits precision.
// Returns `Nothing` on invalid input.
pub fun parse-ddouble( s  : string ) : maybe<ddouble>
  val t = s.trim
  match(t.find(rx-ddouble))
    Just(cap) ->
      val whole	= cap.groups[1]
      val frac	= cap.groups[2].trim-right("0")
      val exp	= cap.groups[3].parse-int.default(0)
      val w	= (whole + frac).parse-int.default(0)
      val e	= exp - frac.count
      val x	= ddouble-int-exp(w,e)
      // trace("parse: s: " + s + "\n  w: " + w.show + ", e: " + e.show + "\n  wx: " + w.ddouble.show-sum + ", x: " + x.show + "\n  " + x.show-sum)
      Just(x)

    Nothing -> match(t.find(rx-special))
      Just(cap) ->
        val special =	if !cap.groups[1].is-empty then dd-nan
                       elif cap.groups[2] == "-" then dd-neginf
                       else dd-posinf
        Just(special)

      Nothing -> match(t.find(rx-float64-sum))
        Just(cap) ->
          match cap.groups[1].parse-float64
            Just(hi) -> match cap.groups[2].parse-float64
              Just(lo) -> Just( hi.ddouble + lo.ddouble )
              Nothing -> Nothing

            Nothing -> Nothing

        Nothing -> Nothing



val rx-ddouble = regex(r"^([\-\+]?\d+)(?:\.(\d*))?(?:[eE]([\-\+]?\d+))?$")
val rx-float64-sum = regex(r"^(\S*\d) *\+ *(\S+)$")
val rx-special = regex(r"^(nan)|([\+\-])?inf(?:inity)?$",ignorecase=True)
*/

pub fun parse-ddoublestd/num/ddouble/parse-ddouble: (s : string) -> maybe<ddouble>( ss: string  : stringstd/core/types/string: V )result: -> total maybe<ddouble> : maybestd/core/types/maybe: V -> V<ddoublestd/num/ddouble/ddouble: V>
  ss: string.trimstd/core/string/trim: (s : string) -> string.to-lowerstd/core/string/to-lower: (s : string) -> string.slicestd/core/sslice/slice: (s : string) -> sslice.parse-eofstd/text/parse/parse-eof: (input : sslice, p : () -> parse ddouble) -> parse-error<ddouble>(pddoublestd/num/ddouble/pddouble: () -> parse ddouble).maybestd/text/parse/maybe: (perr : parse-error<ddouble>) -> maybe<ddouble>

pub fun pddoublestd/num/ddouble/pddouble: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V
  (pddouble-sumstd/num/ddouble/pddouble-sum: () -> parse ddouble ||std/text/parse/(||): (p1 : parser<total,ddouble>, p2 : parser<total,ddouble>) -> parse ddouble pddouble-normalstd/num/ddouble/pddouble-normal: () -> parse ddouble)

fun pddouble-sumstd/num/ddouble/pddouble-sum: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V
  val hihi: float64 = pdoublestd/num/float64/pdouble: () -> parse float64()
  if !std/core/types/bool/(!): (b : bool) -> parse boolhihi: float64.is-finitestd/num/float64/is-finite: (d : float64) -> parse bool returnreturn: ddouble hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble   // NaN or +-Infinity
  pstringstd/text/parse/pstring: (s : string) -> parse string(" + "literal: string
count= 3)
  val lolo: float64 = pdoublestd/num/float64/pdouble: () -> parse float64()
  returnreturn: ddouble (hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> parse ddouble lolo: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble)

fun pddouble-normalstd/num/ddouble/pddouble-normal: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V
  val negneg: bool   = signstd/text/parse/sign: () -> parse bool()
  val wholewhole: string = digitsstd/text/parse/digits: () -> parse string()
  val fracfrac: string  = optionalstd/text/parse/optional: (default : string, p : parser<total,string>) -> parse string(""literal: string
count= 0, { charstd/text/parse/char: (c : char) -> parse char('.'literal: char
unicode= u002E); digitsstd/text/parse/digits: () -> parse string() } ).trim-rightstd/core/sslice/trim-right: (s : string, sub : string) -> parse string("0"literal: string
count= 1)
  val expexp: int   = optionalstd/text/parse/optional: (default : int, p : parser<total,int>) -> parse int(0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000, { charstd/text/parse/char: (c : char) -> parse char('e'literal: char
unicode= u0065); pintstd/text/parse/pint: () -> parse int() } )
  val ww: int = ((if negneg: bool then "-"literal: string
count= 1 else ""literal: string
count= 0) ++std/core/types/(++): (x : string, y : string) -> parse string wholewhole: string ++std/core/types/(++): (x : string, y : string) -> parse string fracfrac: string).parse-intstd/core/int/parse-int: (s : string, hex : ? bool) -> parse maybe<int>.defaultstd/core/maybe/default: (m : maybe<int>, nothing : int) -> parse int(0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000)
  val ee: int = expexp: int -std/core/int/(-): (x : int, y : int) -> parse int fracfrac: string.countstd/core/string/count: (s : string) -> parse int
  val xx: ddouble = ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> parse ddouble(ww: int,ee: int)
  // trace("parse-normal: w: " + w.show + ", e: " ++ e.show ++ "\n  wx: " ++ w.ddouble.show-sum ++ ", x: " ++ x.show ++ "\n  " ++ x.show-sum)
  xx: ddouble

/*------------------------------------------------------
  Advanced operations
------------------------------------------------------*/

/*------------------------------------------------------
  Roots
------------------------------------------------------*/

// The square root of a non-negative `:ddouble` `x`.
// For negative `x`, `dd-nan` is returned.
pub fun sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  // From the QD library:
  // Use Karp's strategy: if a is an approximation to sqrt(x), then
  // sqrt(x) = x*a + (x - (x*a)^2) * a/2
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble
  elif xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then dd-nanstd/num/ddouble/dd-nan: ddouble
  else
    val aa: float64  = 1.0literal: float64
hex64= 0x1p0 /std/num/float64/(/): (x : float64, y : float64) -> float64 sqrtstd/num/float64/sqrt: (d : float64) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64)
    val t1t1: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 aa: float64
    val t2t2: float64 = (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(t1t1: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble)).histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 aa: float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 0.5literal: float64
hex64= 0x1p-1
    dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble(t1t1: float64,t2t2: float64)


// The `n`-th root of a `:ddouble` number `x`.
// `n` must be positive, and if `n` is even, then
// `x` must not be negative.
pub fun nrootstd/num/ddouble/nroot: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if nn: int==std/core/int/(==): (x : int, y : int) -> bool2literal: int
dec  = 2
hex8 = 0x02
bit8 = 0b00000010 then xx: ddouble.sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble
  elif nn: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec  = 1
hex8 = 0x01
bit8 = 0b00000001 then xx: ddouble
  elif (nn: int<=std/core/int/(<=): (x : int, y : int) -> bool0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 ||std/core/types/(||): (x : bool, y : bool) -> bool (nn: int.is-evenstd/core/int/is-even: (i : int) -> bool &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool)) then dd-nanstd/num/ddouble/dd-nan: ddouble
  elif xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble
  else
    val rr: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble
    val a0a0: ddouble = expstd/num/float64/exp: (p : float64) -> float64( ~std/num/float64/(~): (f : float64) -> float64lnstd/num/float64/ln: (f : float64) -> float64(rr: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) /std/num/float64/(/): (x : float64, y : float64) -> float64 nn: int.float64std/num/float64/float64: (i : int) -> float64 ).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble
    val a1a1: ddouble = a0a0: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ((a0a0: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble rr: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble(a0a0: ddouble,nn: int))) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble nn: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)
    (if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble else onestd/num/ddouble/one: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble a1a1: ddouble


// Multiply `x` by a `:float64` `p` where `p` must be a power of 2.
fun mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, pp: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 pp: float64, xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 pp: float64)

val one-halfstd/num/ddouble/one-half: ddouble : ddoublestd/num/ddouble/ddouble: V = 0.5literal: float64
hex64= 0x1p-1.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble

fun halfstd/num/ddouble/half: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble, 0.5literal: float64
hex64= 0x1p-1 )

fun twicestd/num/ddouble/twice: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble, 2.0literal: float64
hex64= 0x1p1 )


/*------------------------------------------------------
  Exponents and logarithms
------------------------------------------------------*/

// Return _e_ (`dd-e`) to the power of `x`.
pub fun expstd/num/ddouble/exp: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  /* From the QD library:
     Strategy:  We first reduce the size of x by noting that

          exp(kr + m * ln(2)) = 2^m * exp(r)^k

     where m and k are integers.  By choosing m appropriately
     we can make |kr| <= ln(2) / 2 = 0.347.  Then exp(r) is
     evaluated using the familiar Taylor series.  Reducing the
     argument substantially speeds up the convergence.
  */
  val kk: float64 : float64std/core/types/float64: V = 512.0literal: float64
hex64= 0x1p9
  val inv-kinv-k: float64 = 1.0literal: float64
hex64= 0x1p0 /std/num/float64/(/): (x : float64, y : float64) -> float64 kk: float64

  if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 <=std/num/float64/(<=): (x : float64, y : float64) -> bool -709.0literal: float64
hex64= -0x1.628p9 then zerostd/num/ddouble/zero: ddouble
  elif xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 >=std/num/float64/(>=): (x : float64, y : float64) -> bool 709.0literal: float64
hex64= 0x1.628p9 then dd-posinfstd/num/ddouble/dd-posinf: ddouble
  elif xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble
  elif xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolonestd/num/ddouble/one: ddouble then dd-estd/num/ddouble/dd-e: ddouble
  else
    val mm: int : intstd/core/types/int: V = floorstd/num/float64/floor: (d : float64) -> (local<$3613>) float64( (xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$3613>) float64 /std/num/float64/(/): (x : float64, y : float64) -> (local<$3613>) float64 dd-ln2std/num/ddouble/dd-ln2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$3613>) float64) +std/num/float64/(+): (x : float64, y : float64) -> (local<$3613>) float64 0.5literal: float64
hex64= 0x1p-1 ).intstd/num/float64/int: (f : float64) -> (local<$3613>) int
    val rr: ddouble = mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> (local<$3613>) ddouble( xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$3613>) ddouble (mm: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> (local<$3613>) ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> (local<$3613>) ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble), inv-kinv-k: float64 )
    val pp: ddouble = rr: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> (local<$3613>) ddouble
    val tt: ddouble = rr: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$3613>) ddouble pp: ddouble.halfstd/num/ddouble/half: (x : ddouble) -> (local<$3613>) ddouble
    var accacc: local-var<$3613,ddouble> := exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> (local<$3613>) ddouble( pp: ddouble, tt: ddouble, rr: ddouble, inv-kinv-k: float64 *std/num/float64/(*): (x : float64, y : float64) -> (local<$3613>) float64 dd-epsilonstd/num/ddouble/dd-epsilon: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$3613>) float64, exp-factorsstd/num/ddouble/exp-factors: list<ddouble> )
    repeatstd/core/repeat: (n : int, action : () -> (local<$3613>) ()) -> (local<$3613>) ()(9literal: int
dec  = 9
hex8 = 0x09
bit8 = 0b00001001)
      accacc: local-var<$3613,ddouble> :=std/core/types/local-set: (v : local-var<$3613,ddouble>, assigned : ddouble) -> (local<$3613>) () accacc: ddouble.twicestd/num/ddouble/twice: (x : ddouble) -> (local<$3613>) ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$3613>) ddouble accacc: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> (local<$3613>) ddouble

    accacc: local-var<$3613,ddouble> :=std/core/types/local-set: (v : local-var<$3613,ddouble>, assigned : ddouble) -> (local<$3613>) () accacc: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$3613>) ddouble onestd/num/ddouble/one: ddouble
    ldexpstd/num/ddouble/ldexp: (x : ddouble, exp : int) -> (local<$3613>) ddouble( accacc: ddouble, mm: int )


val exp-factorsstd/num/ddouble/exp-factors: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a>
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.66666666666666657e-01literal: float64
hex64= 0x1.5555555555555p-3,  9.25185853854297066e-18literal: float64
hex64= 0x1.5555555555555p-57),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(4.16666666666666644e-02literal: float64
hex64= 0x1.5555555555555p-5,  2.31296463463574266e-18literal: float64
hex64= 0x1.5555555555555p-59),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(8.33333333333333322e-03literal: float64
hex64= 0x1.1111111111111p-7,  1.15648231731787138e-19literal: float64
hex64= 0x1.1111111111111p-63),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.38888888888888894e-03literal: float64
hex64= 0x1.6c16c16c16c17p-10, -5.30054395437357706e-20literal: float64
hex64= -0x1.f49f49f49f49fp-65),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.98412698412698413e-04literal: float64
hex64= 0x1.a01a01a01a01ap-13,  1.72095582934207053e-22literal: float64
hex64= 0x1.a01a01a01a01ap-73),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(2.48015873015873016e-05literal: float64
hex64= 0x1.a01a01a01a01ap-16,  2.15119478667758816e-23literal: float64
hex64= 0x1.a01a01a01a01ap-76)
]std/core/types/Nil: forall<a> list<a>

fun exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> ddouble( pp: ddouble : ddoublestd/num/ddouble/ddouble: V, tt: ddouble : ddoublestd/num/ddouble/ddouble: V,
                rr: ddouble : ddoublestd/num/ddouble/ddouble: V,
                epseps: float64 : float64std/core/types/float64: V,
                fsfs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V>,
                ss: ? ddouble : ddoublestd/num/ddouble/ddouble: V = zerostd/num/ddouble/zero: ddouble )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  match fsfs: list<ddouble>
    Nilstd/core/types/Nil: forall<a> list<a> -> ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble
    Consstd/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a>(ff: ddouble,fs1fs1: list<ddouble>) ->
      val s1s1: ddouble = ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble
      val p1p1: ddouble = pp: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble rr: ddouble
      val t1t1: ddouble = p1p1: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ff: ddouble
      if t1t1: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 <=std/num/float64/(<=): (x : float64, y : float64) -> bool epseps: float64
       then ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble
       else exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> ddouble( p1p1: ddouble, t1t1: ddouble, rr: ddouble, epseps: float64, fs1fs1: list<ddouble>, s1s1: ddouble )



// The natural logarithm (in base _e_) of `x`.
pub fun lnstd/num/ddouble/ln: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  /* From QD Library:
     Strategy.  The Taylor series for log converges much more
     slowly than that of exp, due to the lack of the factorial
     term in the denominator.  Hence this routine instead tries
     to determine the root of the function
         f(x) = exp(x) - a
     using Newton iteration.  The iteration is given by
         x' = x - f(x)/f'(x)
            = x - (1 - a * exp(-x))
            = x + a * exp(-x) - 1.
     Only one iteration is needed, since Newton's iteration
     approximately doubles the number of digits per iteration.
  */
  if xx: ddouble <=std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool zerostd/num/ddouble/zero: ddouble then (if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolzerostd/num/ddouble/zero: ddouble then dd-neginfstd/num/ddouble/dd-neginf: ddouble else dd-nanstd/num/ddouble/dd-nan: ddouble)
  elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then zerostd/num/ddouble/zero: ddouble
  elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool dd-estd/num/ddouble/dd-e: ddouble then onestd/num/ddouble/one: ddouble
  elif xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble
  else
    val a0a0: ddouble = lnstd/num/float64/ln: (f : float64) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble // approximate
    val a1a1: ddouble = a0a0: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble expstd/num/ddouble/exp: (x : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublea0a0: ddouble) -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble)
    a1a1: ddouble


// Return the logarithm in some base `b` of a `:ddouble` `x`
pub fun logstd/num/ddouble/log: (x : ddouble, base : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, basebase: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(basebase: ddouble)

// The logarithm in base 10 of `x`.
pub fun log10std/num/ddouble/log10: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-ln10std/num/ddouble/dd-ln10: ddouble

// The logarithm in base 2 of `x`.
pub fun log2std/num/ddouble/log2: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble

// `x` to the power of `y` both as `:ddouble`
pub fun powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  expstd/num/ddouble/exp: (x : ddouble) -> ddouble(yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble))


// Return `ln(1.0 + x)`.
// Avoids potential imprecision for small `x` where adding `1.0` explicitly
// may lead to rounding errors.
pub fun ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble else
    val yy: ddouble = onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble xx: ddouble
    val zz: ddouble = yy: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble
    if zz: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then xx: ddouble else
      lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(yy: ddouble) *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble zz: ddouble)


// Return `exp(x - 1.0)`.
// Avoids rounding errors for values of `x` very close to `1.0`.
pub fun expm1std/num/ddouble/expm1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble else
    val yy: ddouble = expstd/num/ddouble/exp: (x : ddouble) -> ddouble(xx: ddouble)
    if yy: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolonestd/num/ddouble/one: ddouble then xx: ddouble else
      val ymym: ddouble = yy: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble
      if ymym: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble else ymym: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(yy: ddouble))


fun log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  dd-log2estd/num/ddouble/dd-log2e: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(xx: ddouble)

fun exp2m1std/num/ddouble/exp2m1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  expm1std/num/ddouble/expm1: (x : ddouble) -> ddouble(dd-ln2std/num/ddouble/dd-ln2: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble xx: ddouble)

// Returns `ln(exp(x) + exp(y))`.
// Avoids overlow/underflow errors.
pub fun lnaddexpstd/num/ddouble/lnaddexp: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolyy: ddouble then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble else
    val zz: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble
    if zz: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(expstd/num/ddouble/exp: (x : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublezz: ddouble))
                else yy: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(expstd/num/ddouble/exp: (x : ddouble) -> ddouble(zz: ddouble))


// Returns `log2( exp2(x) + exp2(y) )`.
// Avoids overlow/underflow errors.
pub fun logaddexp2std/num/ddouble/logaddexp2: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolyy: ddouble then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble else
    val zz: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble
    if zz: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble(exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublezz: ddouble))
                else yy: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble(exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble(zz: ddouble))


// Return if two `:ddouble`s are nearly equal with respect to some `epsilon` (=`8*dd-epsilon`).
// The epsilon is the nearest difference for numbers around 1.0. The routine automatically
// scales the epsilon for larger and smaller numbers, and for numbers close to zero.
pub fun nearly-eqstd/num/ddouble/nearly-eq: (x : ddouble, y : ddouble, epsilon : ? ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V, epsilonepsilon: ? ddouble : ddoublestd/num/ddouble/ddouble: V = dd-epsilon8std/num/ddouble/dd-epsilon8: ddouble )result: -> total bool : boolstd/core/types/bool: V
  if xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool yy: ddouble returnreturn: bool Truestd/core/types/True: bool
  val diffdiff: ddouble = (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).absstd/num/ddouble/abs: (x : ddouble) -> ddouble
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool diffdiff: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool dd-minstd/num/ddouble/dd-min: ddouble then
    // very close to zero, scale the epsilon for subnormal numbers
    (twostd/num/ddouble/two: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublediffdiff: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool (epsilonepsilon: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble dd-minstd/num/ddouble/dd-min: ddouble))
  else
    val sumsum: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble yy: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble
    ((twostd/num/ddouble/two: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublediffdiff: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (if sumsum: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool dd-maxstd/num/ddouble/dd-max: ddouble then dd-maxstd/num/ddouble/dd-max: ddouble else sumsum: ddouble)) <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool epsilonepsilon: ddouble)

// Return if two `:ddouble`s are nearly equal with respect to an `epsilon` of `10*dd-epsilon`.
// See also `nearly-eq` which takes an explicit `epsilon`.
pub fun (~=)std/num/ddouble/(~=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V
  nearly-eqstd/num/ddouble/nearly-eq: (x : ddouble, y : ddouble, epsilon : ? ddouble) -> bool(xx: ddouble,yy: ddouble)

/*------------------------------------------------------
  Summing
------------------------------------------------------*/

// Return the sum of a list of doubles.
// Uses [Kahan-Babu&scaron;kan-Neumaier summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements)
// to minimize rounding errors. This
// is more precise as Kahan summation and about as fast.\
// `[1.0e3,1.0e97,1.0e3,-1.0e97].sum == 2000.0`\
// A. Neumaier, _Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen_.
// Math. Mechanik, 54:39--51, 1974.
pub fun list/sumstd/num/ddouble/list/sum: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  var totaltotal: local-var<$1668,ddouble> : ddoublestd/num/ddouble/ddouble: V := zerostd/num/ddouble/zero: ddouble
  var compcomp: local-var<$1668,ddouble>  : ddoublestd/num/ddouble/ddouble: V := zerostd/num/ddouble/zero: ddouble
  xsxs: list<ddouble>.foreachstd/core/list/foreach: (xs : list<ddouble>, action : (ddouble) -> (local<$1668>) ()) -> (local<$1668>) () fnfn: (x : ddouble) -> (local<$1668>) ()(xx: ddouble)
    val tt: ddouble = totaltotal: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble xx: ddouble
    val cc: ddouble = if totaltotal: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> (local<$1668>) ddouble >=std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> (local<$1668>) bool xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> (local<$1668>) ddouble then (totaltotal: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble tt: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble xx: ddouble else (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble tt: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble totaltotal: ddouble
    compcomp: local-var<$1668,ddouble>  :=std/core/types/local-set: (v : local-var<$1668,ddouble>, assigned : ddouble) -> (local<$1668>) () compcomp: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble cc: ddouble
    totaltotal: local-var<$1668,ddouble> :=std/core/types/local-set: (v : local-var<$1668,ddouble>, assigned : ddouble) -> (local<$1668>) () tt: ddouble
  totaltotal: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1668>) ddouble compcomp: ddouble;

// The hypotenuse of `x` and `y`: `sqrt(x*x + y*y)`.
// Prevents overflow for large numbers.
pub fun hypotstd/num/ddouble/hypot: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val xxxx: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble)
  val yyyy: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble)
  val lolo: ddouble = minstd/num/ddouble/min: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble)
  val hihi: ddouble = maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble)
  if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else
    val zz: ddouble  = lolo: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble
    hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble zz: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublezz: ddouble )


// The square root of the sum of the squares of three doubles.
// Prevents overflow for large numbers.
pub fun xyz/hypotstd/num/ddouble/xyz/hypot: (x : ddouble, y : ddouble, z : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V, zz: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val xxxx: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble)
  val yyyy: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble)
  val zzzz: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(zz: ddouble)
  val hihi: ddouble = maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble),zzzz: ddouble)
  if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else
    hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xxxx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(yyyy: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(zzzz: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) )


// The square root of the sum of squares of a list of doubles.
// Prevents overflow for large numbers and uses Kahan-Babu&scaron;kan-Neumaier summation
// for precision.
pub fun list/hypotstd/num/ddouble/list/hypot: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val hihi: ddouble = xsxs: list<ddouble>.abs-maxstd/num/ddouble/list/abs-max: (xs : list<ddouble>) -> ddouble
  if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else
    hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble xsxs: list<ddouble>.mapstd/core/list/map: (xs : list<ddouble>, f : (ddouble) -> ddouble) -> list<ddouble>(fnfn: (x : ddouble) -> ddouble(xx: ddouble){ sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) }).sumstd/num/ddouble/list/sum: (xs : list<ddouble>) -> ddouble.sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble


// The maximum of the absolute values.
fun abs-maxstd/num/ddouble/abs-max: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble),absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble))

// The maximum of a list of absolute values.
fun list/abs-maxstd/num/ddouble/list/abs-max: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  xsxs: list<ddouble>.foldrstd/core/list/foldr: (xs : list<ddouble>, z : ddouble, f : (ddouble, ddouble) -> ddouble) -> ddouble(ddouble/zerostd/num/ddouble/zero: ddouble,fnfn: (x : ddouble, m : ddouble) -> ddouble(xx: ddouble,mm: ddouble){ maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble),mm: ddouble) })

/*------------------------------------------------------
  Trigonometry
------------------------------------------------------*/

// The sine function of a given angle in radians.
pub fun sinstd/num/ddouble/sin: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble).fststd/core/types/tuple2/fst: (tuple2 : (ddouble, ddouble)) -> ddouble

// The cosine function of a given angle in radians.
pub fun cosstd/num/ddouble/cos: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble).sndstd/core/types/tuple2/snd: (tuple2 : (ddouble, ddouble)) -> ddouble

// The tangent of a given angle in radians.
pub fun tanstd/num/ddouble/tan: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)ss: ddouble,cc: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble)
  ss: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble cc: ddouble

// Calculate sine and cosine on an angle in radians.
pub fun sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble)( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (ddouble, ddouble) : (std/core/types/tuple2: (V, V) -> Vddoublestd/num/ddouble/ddouble: V,ddoublestd/num/ddouble/ddouble: V)
  // quick approximation for small values
  if radrad: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 <std/num/float64/(<): (x : float64, y : float64) -> bool 1.0e-11literal: float64
hex64= 0x1.5fd7fe1796495p-37 then
    // return (rad, one - 0.5.ddouble*sqr(rad))
    val ss: ddouble = radrad: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (radrad: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble 3literal: int
dec  = 3
hex8 = 0x03
bit8 = 0b00000011.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble))
    val cc: ddouble = (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble).sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble
    returnreturn: (ddouble, ddouble) (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)ss: ddouble,cc: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)

  /*
  To compute sin(x), we choose integers a, b so that
       rad = s + a * (pi/2) + b * (pi/16)
  and |s| <= pi/32.  A precomputed sin(pi/16) and
  Chebyshev series, we can compute very precisely.
  */
  // find a, b, and s.
  /*
  val z  = (rad / dd-twopi).round
  val r  = rad - (dd-twopi * z)

  val qa : float64 = floor( (r.hi / dd-pi2.hi) + 0.5)
  val a  = qa.truncate.int
  val t  = r - (dd-pi2 * qa.ddouble)
  val qb : float64 = floor((t.hi / dd-pi16.hi) + 0.5)
  val b  = qb.truncate.int

  val s  = t - (dd-pi16 * qb.ddouble)
  */

  val x1x1: ddouble = radrad: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-twopistd/num/ddouble/dd-twopi: ddouble
  val x3x3: ddouble = x1x1: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble x1x1: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble
  // s = x - a*pi/2 - b*pi/16
  val x32x32: ddouble = x3x3: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble x3x3: ddouble
  val x34x34: ddouble = x32x32: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble x32x32: ddouble // 4.ddouble * x3
  val aa: int   = x34x34: ddouble.intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int
  val bb: int   = (8literal: int
dec  = 8
hex8 = 0x08
bit8 = 0b00001000.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (x34x34: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble aa: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)).intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int
  val ss: ddouble   = dd-pistd/num/ddouble/dd-pi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (x32x32: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ((8literal: int
dec  = 8
hex8 = 0x08
bit8 = 0b00001000*std/core/int/(*): (int, int) -> intaa: int +std/core/int/(+): (x : int, y : int) -> int bb: int).ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble 16literal: int
dec  = 16
hex8 = 0x10
bit8 = 0b00010000.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble))
  val s2s2: ddouble = ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble

  // Use the Chebyshev series for best precision.
  val sinssins: ddouble = ss: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ch-factorsstd/num/ddouble/ch-factors: list<ddouble>.foldlstd/core/list/foldl: (xs : list<ddouble>, z : ddouble, f : (ddouble, ddouble) -> ddouble) -> ddouble(ddouble/zerostd/num/ddouble/zero: ddouble) fnfn: (acc : ddouble, f : ddouble) -> ddouble(accacc: ddouble,ff: ddouble) { ff: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble accacc: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubles2s2: ddouble }
  val cosscoss: ddouble = (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble).sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble

  // sinb = sin(b*pi/16)
  val sinbsinb: ddouble = if bb: int >=std/core/int/(>=): (x : int, y : int) -> bool 0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then sin16std/num/ddouble/sin16: (i : int) -> ddouble(bb: int) else ~std/num/ddouble/(~): (x : ddouble) -> ddoublesin16std/num/ddouble/sin16: (i : int) -> ddouble(~std/core/int/(~): (i : int) -> intbb: int)
  val cosbcosb: ddouble = sin16std/num/ddouble/sin16: (i : int) -> ddouble(8literal: int
dec  = 8
hex8 = 0x08
bit8 = 0b00001000 -std/core/int/(-): (x : int, y : int) -> int bb: int.absstd/core/int/abs: (i : int) -> int)

  if aa: int==std/core/int/(==): (x : int, y : int) -> bool0literal: int
dec  = 0
hex8 = 0x00
bit8 = 0b00000000 then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)
  elif aa: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec  = 1
hex8 = 0x01
bit8 = 0b00000001 then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, ~std/num/ddouble/(~): (x : ddouble) -> ddoublecosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)
  elif aa: int==std/core/int/(==): (x : int, y : int) -> bool -1literal: int
dec  = -1
hex8 = 0xFF
bit8 = 0b11111111 then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble  cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble, cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)
  else (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)~std/num/ddouble/(~): (x : ddouble) -> ddoublesinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) // |a| == 2

val ch-factorsstd/num/ddouble/ch-factors: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a>
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 1.6056491947130061858e-10literal: float64
hex64= 0x1.6116039167de1p-33, 6.1925234565562595936e-27literal: float64
hex64= 0x1.ea9f4c1702653p-88 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -2.5052108052208301379e-8literal: float64
hex64= -0x1.ae64561f4804fp-26, -3.6598195022865791876e-25literal: float64
hex64= -0x1.c5104f49cf934p-82 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.0000027557319223964441402literal: float64
hex64= 0x1.71de3a556b36bp-19, -2.0315661398415506513e-22literal: float64
hex64= -0x1.eb33da7509008p-73 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -0.00019841269841269841253literal: float64
hex64= -0x1.a01a01a01a01ap-13, 6.8577289081075077176e-21literal: float64
hex64= 0x1.0313e2634850bp-67 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.0083333333333333332177literal: float64
hex64= 0x1.1111111111111p-7, 1.1563735775184918100e-19literal: float64
hex64= 0x1.110a7e6657814p-63 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -0.16666666666666665741literal: float64
hex64= -0x1.5555555555555p-3, -9.2518585321663028924e-18literal: float64
hex64= -0x1.5555555162e8fp-57 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 1.0000000000000000000literal: float64
hex64= 0x1p0, -6.0239567712403467977e-31literal: float64
hex64= -0x1.86fa265ab5ed7p-101 )
]std/core/types/Nil: forall<a> list<a>

// Return sin(i*pi/16) for 0 <= i <= 8
fun sin16std/num/ddouble/sin16: (i : int) -> ddouble( ii: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  sin16-tablestd/num/ddouble/sin16-table: vector<ddouble>.atstd/core/vector/at: (v : vector<ddouble>, index : int) -> maybe<ddouble>(ii: int).defaultstd/core/maybe/default: (m : maybe<ddouble>, nothing : ddouble) -> ddouble(dd-nanstd/num/ddouble/dd-nan: ddouble)

val sin16-tablestd/num/ddouble/sin16-table: vector<ddouble> : vectorstd/core/types/vector: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a>
  ddouble/zerostd/num/ddouble/zero: ddouble,
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.19509032201612827584literal: float64
hex64= 0x1.8f8b83c69a60bp-3, -7.9910790684617312634e-18literal: float64
hex64= -0x1.26d19b9ff8d82p-57 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.38268343236508978178literal: float64
hex64= 0x1.87de2a6aea963p-2, -1.0050772696461587612e-17literal: float64
hex64= -0x1.72cedd3d5a61p-57 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.55557023301960217765literal: float64
hex64= 0x1.1c73b39ae68c8p-1, 4.7094109405616768214e-17literal: float64
hex64= 0x1.b25dd267f66p-55 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.70710678118654757274literal: float64
hex64= 0x1.6a09e667f3bcdp-1, -4.8336466567264561092e-17literal: float64
hex64= -0x1.bdd3413b26455p-55 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.83146961230254523567literal: float64
hex64= 0x1.a9b66290ea1a3p-1, 1.4073856984728100930e-18literal: float64
hex64= 0x1.9f630e8b6dafp-60 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.92387953251128673848literal: float64
hex64= 0x1.d906bcf328d46p-1, 1.7645047084336683223e-17literal: float64
hex64= 0x1.457e610231ac4p-56 ),
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.98078528040323043058literal: float64
hex64= 0x1.f6297cff75cbp-1, 1.8546939997825014970e-17literal: float64
hex64= 0x1.562172a361fd6p-56 ),
  onestd/num/ddouble/one: ddouble
]std/core/types/Nil: forall<a> list<a>.vectorstd/core/vector/list/vector: (xs : list<ddouble>) -> vector<ddouble>

// Return `x` with the sign of `y`.
pub fun with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if yy: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddouble(xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble) else xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble

// The arc-sine of `x`. Returns the angle in radians.
pub fun asinstd/num/ddouble/asin: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val aa: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble
  if aa: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble
  elif aa: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-pi2std/num/ddouble/dd-pi2: ddouble.with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble(xx: ddouble)
  else atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( xx: ddouble, sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble) )

// The arc-cosine of `x`. Returns the angle in radians.
pub fun acosstd/num/ddouble/acos: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  val aa: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble
  if aa: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble
  elif aa: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then (if xx: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else dd-pistd/num/ddouble/dd-pi: ddouble)
  else atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble), xx: ddouble )

// The arc-tangent of `x`. Returns the angle in radians.
pub fun atanstd/num/ddouble/atan: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble(xx: ddouble,onestd/num/ddouble/one: ddouble)

// The arc-tangent of a point (`x`,`y`). Returns the angle with respect to the x-axis in radians between -&pi; and &pi;.
pub fun atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( yy: ddouble : ddoublestd/num/ddouble/ddouble: V, xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then
    if yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble /* dd-nan */ else dd-pi2std/num/ddouble/dd-pi2: ddouble.with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble(yy: ddouble)
  elif yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then
    if xx: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else dd-pistd/num/ddouble/dd-pi: ddouble
  elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool yy: ddouble then
    if yy: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then dd-pi4std/num/ddouble/dd-pi4: ddouble else ~std/num/ddouble/(~): (x : ddouble) -> ddoubledd-pi34std/num/ddouble/dd-pi34: ddouble
  elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool ~std/num/ddouble/(~): (x : ddouble) -> ddoubleyy: ddouble then
    if yy: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then dd-pi34std/num/ddouble/dd-pi34: ddouble else ~std/num/ddouble/(~): (x : ddouble) -> ddoubledd-pi4std/num/ddouble/dd-pi4: ddouble
  else
    val rr: ddouble = sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble yy: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble )
    val xrxr: ddouble = xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble rr: ddouble
    val yryr: ddouble = yy: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble rr: ddouble
    val zz: ddouble = atan2std/num/float64/atan2: (x : float64, y : float64) -> float64( yy: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64, xx: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 ).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble // approximation
    val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinzsinz: ddouble,coszcosz: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble)(zz: ddouble)
    if xrxr: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool yryr: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64
     then zz: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ((yryr: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinzsinz: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble coszcosz: ddouble)
     else zz: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ((xrxr: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble coszcosz: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble sinzsinz: ddouble)


// The hyperbolic sine of `x`.
pub fun sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble
  elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5 then
    val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble
    if !std/core/types/bool/(!): (b : bool) -> boolexex: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then exex: ddouble else ( exex: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble
  else
    // small x so use Taylor series to avoid cancellation
    val x2x2: ddouble = xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble
    xx: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble6literal: int
dec  = 6
hex8 = 0x06
bit8 = 0b00000110.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble20literal: int
dec  = 20
hex8 = 0x14
bit8 = 0b00010100.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble42literal: int
dec  = 42
hex8 = 0x2A
bit8 = 0b00101010.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble))))


// The hyperbolic cosine of `x`.
pub fun coshstd/num/ddouble/cosh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble
  elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5 then
    val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble
    if !std/core/types/bool/(!): (b : bool) -> boolexex: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then exex: ddouble else (exex: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble

  else
    val ss: ddouble = xx: ddouble.sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble
    sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble )


// The hyperbolic tangent of `x`.
pub fun tanhstd/num/ddouble/tanh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble
  elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5 then
    val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble
    val iexiex: ddouble = onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble
    if exex: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble
    elif exex: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble
    else ((exex: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble iexiex: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (exex: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble iexiex: ddouble))

  else
    val ss: ddouble = xx: ddouble.sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble
    val cc: ddouble = sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble )
    returnreturn: ddouble (ss: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble cc: ddouble)


// The area hyperbolic sine of `x`.
pub fun asinhstd/num/ddouble/asinh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  lnstd/num/ddouble/ln: (x : ddouble) -> ddouble( xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble))

// The area hyperbolic cosine of `x`.
pub fun acoshstd/num/ddouble/acosh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble else lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble))

// The area hyperbolic tangent of `x`.
pub fun atanhstd/num/ddouble/atanh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  if xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble
   else lnstd/num/ddouble/ln: (x : ddouble) -> ddouble((onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble