Chapter 2
Numbers

Naturals

The naturals (kNatural) are either 32 or 64 bits in size, and store non-negative integers. We start with an error declaration. These constants will be used as error codes throughout this module.

typedef enum
{
  kNaturalErrorOk,
  kNaturalErrorDivZero,
  kNaturalErrorOverflow
} kNaturalErrorKind;

The simple functions that follow are defined.

inline kBool kNaturalIsZero(kNatural x)
{
  kError = kNaturalErrorOk;

  return x == 0;
}

inline kBool kNaturalIsEven(kNatural x)
{
  kError = kNaturalErrorOk;

  return x % 2 == 0;
}

inline kBool kNaturalIsOdd(kNatural x)
{
  kError = kNaturalErrorOk;

  return x % 2 == 1;
}

inline kNatural kNaturalMin(kNatural x, kNatural y)
{
  kError = kNaturalErrorOk;

  return x <= y ? x : y;
}

inline kNatural kNaturalMax(kNatural x, kNatural y)
{
  kError = kNaturalErrorOk;

  return x >= y ? x : y;
}

There are functions that tell whether a natural is zero or even or odd, and a couple of functions that return the minimum and maximum of two naturals.

Integers

The integers (kInteger) are numbers that are either positive, or negative, or zero. Their size is either 32 or 64 bits, depending on the CPU. As above, we have error codes as follows.

typedef enum
{
  kIntegerErrorOk,
  kIntegerErrorDivZero,
  kIntegerErrorOverflow
} kIntegerErrorKind;

We define functions on this type as follows.

inline kBool kIntegerIsZero(kInteger x)
{
  kError = kIntegerErrorOk;

  return x == 0;
}

inline kBool kIntegerIsNegative(kInteger x)
{
  kError = kIntegerErrorOk;

  return x < 0;
}

inline kBool kIntegerIsPositive(kInteger x)
{
  kError = kIntegerErrorOk;

  return x > 0;
}

inline kBool kIntegerIsEven(kInteger x)
{
  kError = kIntegerErrorOk;

  return x % 2 == 0;
}

inline kBool kIntegerIsOdd(kInteger x)
{
  kError = kIntegerErrorOk;

  return x % 2 == 1;
}

inline kInteger kIntegerAbs(kInteger x)
{
  kError = kIntegerErrorOk;

  return x < 0 ? -x : x;
}

inline kInteger kIntegerNeg(kInteger x)
{
  kError = kIntegerErrorOk;

  return -x;
}

inline kInteger kIntegerMin(kInteger x, kInteger y)
{
  kError = kIntegerErrorOk;

  return x <= y ? x : y;
}

inline kInteger kIntegerMax(kInteger x, kInteger y)
{
  kError = kIntegerErrorOk;

  return x >= y ? x : y;
}

There are functions that say whether the integer is zero or positive or negative or even or odd, the function that returns the absolute value of an integer, the one that returns the integer with its sign reversed, and two functions that return the minimum or maximum of two integers.

Reals

Another set of numbers we are interested in are the reals. These are represented by the type kReal, which is either 32 or 64 bits in size. We start with a couple of enums, used to classify our real numbers.

typedef enum
{
  kRealSgnKindNaN,
  kRealNegKindInfinity,
  kRealNegKindNormal,
  kRealNegKindDenormal,
  kRealNegKindZero,
  kRealPosKindZero,
  kRealPosKindDenormal,
  kRealPosKindNormal,
  kRealPosKindInfinity
} kRealSgnKind;

typedef enum
{
  kRealAbsKindNaN,
  kRealAbsKindInfinity,
  kRealAbsKindNormal,
  kRealAbsKindDenormal,
  kRealAbsKindZero
} kRealAbsKind;

Another enum follows, used in rounding operations:

typedef enum
{
  kRealRoundModeNearest,
  kRealRoundModeNegInfinity,
  kRealRoundModePosInfinity,
  kRealRoundModeZero
} kRealRoundModeKind;

We need error codes as follows.

typedef enum
{
  kRealErrorOk,
  kRealErrorInvalid,
  kRealErrorDenormal,
  kRealErrorDivZero,
  kRealErrorOverflow,
  kRealErrorUnderflow,
  kRealErrorPrecision,
  kRealErrorDomain
} kRealErrorKind;

Our final types are used in the classification of reals.

typedef union
{
  kSingle f;
  kULong w;
  struct
  {
    kULong mantissa : 23;
    kULong exponent : 8;
    kULong negative : 1;
  } b;
} kSingleComponents;

typedef union
{
  kDouble f;
  kULLong w;
  struct
  {
    kULLong mantissa : 52;
    kULLong exponent : 11;
    kULLong negative : 1;
  } b;
} kDoubleComponents;

#ifdef kBit32
typedef kSingleComponents kRealComponents;
#endif

#ifdef kBit64
typedef kDoubleComponents kRealComponents;
#endif

Next, we define a few useful constants: π and variants, Neper’s number, and Euler’s gamma constant.

const kReal kRealHalfPi = 3.14159265358979323846264338327950288 / 2.0,
  kRealPi            = 3.14159265358979323846264338327950288,
  kRealThreeHalfsPi  = 3.14159265358979323846264338327950288 * 1.5,
  kRealDoublePi      = 3.14159265358979323846264338327950288 * 2.0,
  kRealE             = 2.71828182845904523536028747135266249,
  kRealGamma         = 0.57721566490153286060651209008240243;

The following are sometimes useful:

const kReal kRealNaN = 0.0 / 0.0,
  kRealNegInfinity = -1.0 / 0.0,
  kRealPosInfinity = +1.0 / 0.0;

Our last definitions are as follows.

#ifdef kBit32
#define kRealEpsilon   0X1P-23F
#define kRealPrecision 9
#define kRealMinExp    -37
#define kRealMaxExp    +38
#define kRealMinNum    0X1P-126F
#define kRealMaxNum    0X1.FFFFFEP127F
#endif

#ifdef kBit64
#define kRealEpsilon   0X1P-52
#define kRealPrecision 17
#define kRealMinExp    -307
#define kRealMaxExp    +308
#define kRealMinNum    0X1P-1022
#define kRealMaxNum    0X1.FFFFFFFFFFFFFP1023
#endif

We need a private variable to store the current rounding mode:

kUShort kRealRoundMode;

And, finally, here are the functions. We start with a private one, the one that deals with errors:

void kRealSetError(void)
{
  #ifdef kCPUIntel
  kUShort sw;

  __asm__
  (
    "fstsw %[sw];"
    : [sw] "=m" (sw)
  );

  if(sw & 0x0001)
    kError = kRealErrorInvalid;
  else if(sw & 0x0002)
    kError = kRealErrorDenormal;
  else if(sw & 0x0004)
    kError = kRealErrorDivZero;
  else if(sw & 0x0008)
    kError = kRealErrorOverflow;
  else if(sw & 0x0010)
    kError = kRealErrorUnderflow;
  else if(sw & 0x0020)
    kError = kRealErrorPrecision;
  else
    kError = kRealErrorOk;
  #endif
}

Another private function is used further below in the calculation of the sine of a real. It was coded by Konstantinos Gkizinos, and is accurate up to a few decimal digits. For the moment it will have to do.

kReal kRealDoSin(kReal x)
{
  kSInt n = 20;
  kReal sum, fa, pow;

  sum = 0.0;
  for(kSInt i = 0; i <= n; i++)
  {
    fa = 1.0;
    pow = 1.0;
    for(kSInt j = 1; j <= 2 * i + 1; j++)
    {
      fa *= j;
      pow *= x;
    }
    sum += ((i % 2 ? -1.0 : 1.0) / fa) * pow;
  }

  return sum;
}

Next are the classification of reals, sign-wise and absolute.

kRealSgnKind kRealSgnClassify(kReal x)
{
  kRealComponents r;

  r.f = x;

  kError = kRealErrorOk;

  #ifdef kBit32
  kULong negative = r.b.negative, exponent = r.b.exponent, mantissa = r.b.mantissa;

  if((exponent == 0) && (mantissa == 0))
  {
    if(negative)
      return kRealNegKindZero;
    else
      return kRealPosKindZero;
  }
  else if(negative == 0)
  {
    if(exponent == 0)
      return kRealPosKindDenormal;
    else if((exponent == 255) && (mantissa == 0))
      return kRealPosKindInfinity;
    else if(exponent == 255)
      return kRealSgnKindNaN;
    else
      return kRealPosKindNormal;
  }
  else
  {
    if(exponent == 0)
      return kRealNegKindDenormal;
    else if((exponent == 255) && (mantissa == 0))
      return kRealNegKindInfinity;
    else if(exponent == 255)
      return kRealSgnKindNaN;
    else
      return kRealNegKindNormal;
  }
  #endif

  #ifdef kBit64
  kULLong negative = r.b.negative, exponent = r.b.exponent, mantissa = r.b.mantissa;

  if((exponent == 0) && (mantissa == 0))
  {
    if(negative)
      return kRealNegKindZero;
    else
      return kRealPosKindZero;
  }
  else if(negative == 0)
  {
    if(exponent == 0)
      return kRealPosKindDenormal;
    else if((exponent == 2047) && (mantissa == 0))
      return kRealPosKindInfinity;
    else if(exponent == 2047)
      return kRealSgnKindNaN;
    else
      return kRealPosKindNormal;
  }
  else
  {
    if(exponent == 0)
      return kRealNegKindDenormal;
    else if((exponent == 2047) && (mantissa == 0))
      return kRealNegKindInfinity;
    else if(exponent == 2047)
      return kRealSgnKindNaN;
    else
      return kRealNegKindNormal;
  }
  #endif
}

kRealAbsKind kRealAbsClassify(kReal x)
{
  kRealComponents r;

  r.f = x;

  kError = kRealErrorOk;

  #ifdef kBit32
  kULong exponent = r.b.exponent, mantissa = r.b.mantissa;

  if((exponent == 0) && (mantissa == 0))
    return kRealAbsKindZero;
  else
  {
    if(exponent == 0)
      return kRealAbsKindDenormal;
    else if((exponent == 255) && (mantissa == 0))
      return kRealAbsKindInfinity;
    else if(exponent == 255)
      return kRealAbsKindNaN;
    else
      return kRealAbsKindNormal;
  }
  #endif

  #ifdef kBit64
  kULLong exponent = r.b.exponent, mantissa = r.b.mantissa;

  if((exponent == 0) && (mantissa == 0))
    return kRealAbsKindZero;
  else
  {
    if(exponent == 0)
      return kRealAbsKindDenormal;
    else if((exponent == 2047) && (mantissa == 0))
      return kRealAbsKindInfinity;
    else if(exponent == 2047)
      return kRealAbsKindNaN;
    else
      return kRealAbsKindNormal;
  }
  #endif
}

The functions that control the rounding method are next.

void kRealSetRoundMode(kRealRoundModeKind mode)
{
  kError = kRealErrorOk;

  #ifdef kCPUIntel
  kRealRoundMode = mode << 10;
  #endif
}

kRealRoundModeKind kRealGetRoundMode(void)
{
  kError = kRealErrorOk;

  #ifdef kCPUIntel
  return kRealRoundMode >> 10;
  #endif
}

Next are the absolute value and the sign-reversing functions.

kReal kRealAbs(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fabs;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealNeg(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fchs;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

The following functions return the fractional part, the integer part, the reciprocal, and the remainder of a real.

kReal kRealFracPart(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x0e7f, %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fldt %[aux];\
    frndint;\
    fsubr %%st(1), %%st;\
    fabs;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw)
  );

  #endif

  kRealSetError();

  return aux;
}

kReal kRealIntgPart(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x0e7f, %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    frndint;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealReciprocal(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fld1;\
    fdiv %%st, %%st(1);\
    fstpt %[aux];\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealRemainder(kReal x, kReal y)
{
  kExtended auxx = x, auxy = y;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[auxx];\
    fldt %[auxy];\
    fprem1;\
    fstpt %[auxx];"
    : [auxx] "=mr" (auxx), [auxy] "=mr" (auxy), [cw] "=m" (cw) : "m" (auxx), "m" (auxy), "m" (cw),
      [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return auxx;
}

There follow the functions to calculate a real raised to the power of another real and to calculate the square root of a real.

kReal kRealPow(kReal base, kReal x)
{
  kExtended auxbase = base, aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fldt %[auxbase];\
    fyl2x;\
    fld %%st;\
    frndint;\
    fst %%st(2);\
    fsubrp;\
    f2xm1;\
    fld1;\
    faddp;\
    fscale;\
    fstp %%st(1);\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), [auxbase] "m" (auxbase), "m" (cw),
      [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealSqrt(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fsqrt;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

The transcendental functions follow. There is the exponential function, the natural logarithm function, and the function that calculates the logarithm to the base of a real number.

kReal kRealExp(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldl2e;\
    fldt %[aux];\
    fmul;\
    fld %%st;\
    frndint;\
    fst %%st(2);\
    fsubrp;\
    f2xm1;\
    fld1;\
    faddp;\
    fscale;\
    fstp %%st(1);\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealLn(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldl2e;\
    fld1;\
    fdiv %%st(1), %%st;\
    fldt %[aux];\
    fyl2x;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealLog(kReal base, kReal x)
{
  kExtended auxbase = base, aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fld1;\
    fldt %[auxbase];\
    fyl2x;\
    fld1;\
    fdiv %%st(1), %%st;\
    fldt %[aux];\
    fyl2x;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), [auxbase] "m" (auxbase), "m" (cw),
      [kRealRoundMode] "r" (kRealRoundMode)
  );
  #endif

  kRealSetError();

  return aux;
}

Then, there are two functions to convert between radians and degrees.

kReal kRealRad2Deg(kReal x)
{
  kError = kRealErrorOk;

  return 180.0 * x / kRealPi;
}

kReal kRealDeg2Rad(kReal x)
{
  kError = kRealErrorOk;

  return kRealPi * x / 180.0;
}

We then have the trigonometric functions.

kReal kRealSin(kReal x)
{
  kRealAbsKind kind = kRealAbsClassify(x);

  if((kind == kRealAbsKindNaN) || (kind == kRealAbsKindInfinity))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else if(x < 0.0)
    return -kRealSin(-x);
  else if(x == 0.0)
  {
    kError = kRealErrorOk;

    return 0.0;
  }
  else if(x < kRealHalfPi)
  {
    kError = kRealErrorOk;

    return kRealDoSin(x);
  }
  else if(x == kRealHalfPi)
  {
    kError = kRealErrorOk;

    return 1.0;
  }
  else if(x < kRealPi)
    return kRealSin(kRealPi - x);
  else if(x == kRealPi)
  {
    kError = kRealErrorOk;

    return 0.0;
  }
  else if(x < kRealThreeHalfsPi)
    return -kRealSin(x - kRealPi);
  else if(x == kRealThreeHalfsPi)
  {
    kError = kRealErrorOk;

    return -1.0;
  }
  else if(x < kRealDoublePi)
    return -kRealSin(kRealDoublePi - x);
  else if(x == kRealDoublePi)
  {
    kError = kRealErrorOk;

    return 0.0;
  }
  else
    return kRealSin(kRealAbs(kRealRemainder(x, kRealDoublePi)));
}

kReal kRealCos(kReal x)
{
  kRealAbsKind kind = kRealAbsClassify(x);

  if((kind == kRealAbsKindNaN) || (kind == kRealAbsKindInfinity))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else if(x >= 0.0)
    return -kRealSin(x - kRealHalfPi);
  else
    return kRealSin(x + kRealHalfPi);
}

kReal kRealTan(kReal x)
{
  kReal vCos = kRealCos(x);
  kRealAbsKind vKind = kRealAbsClassify(vCos);

  if((vKind == kRealAbsKindNaN) || (vKind == kRealAbsKindInfinity) || (vKind == kRealAbsKindZero))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else
    return kRealSin(x) / vCos;
}

kReal kRealCot(kReal x)
{
  kReal vSin = kRealSin(x);
  kRealAbsKind vKind = kRealAbsClassify(vSin);

  if((vKind == kRealAbsKindNaN) || (vKind == kRealAbsKindInfinity) || (vKind == kRealAbsKindZero))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else
    return kRealCos(x) / vSin;
}

kReal kRealSec(kReal x)
{
  kReal vCos = kRealCos(x);
  kRealAbsKind vKind = kRealAbsClassify(vCos);

  if((vKind == kRealAbsKindNaN) || (vKind == kRealAbsKindInfinity) || (vKind == kRealAbsKindZero))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else
    return 1.0 / vCos;
}

kReal kRealCsc(kReal x)
{
  kReal vSin = kRealSin(x);
  kRealAbsKind vKind = kRealAbsClassify(vSin);

  if((vKind == kRealAbsKindNaN) || (vKind == kRealAbsKindInfinity) || (vKind == kRealAbsKindZero))
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
  else
    return 1.0 / vSin;
}

The inverse trigonometric functions are as follows.

kReal kRealASin(kReal x)
{
  kError = kRealErrorOk;

  if(x == -1.0)
    return -kRealPi / 2.0;
  else if(x == 1.0)
    return kRealPi / 2.0;
  else if((x > -1.0) && (x < 1.0))
    return kRealATan(x / kRealSqrt(1.0 - x * x));
  else
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
}

kReal kRealACos(kReal x)
{
  if(x == 0.0)
  {
    kError = kRealErrorOk;

    return kRealPi / 2.0;
  }
  else if((x > 0.0) && (x <= 1.0))
    return kRealATan(kRealSqrt(1.0 - x * x) / x);
  else if((x >= -1.0) && (x < 0.0))
    return kRealPi + kRealATan(kRealSqrt(1.0 - x * x) / x);
  else
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
}

kReal kRealATan(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel

  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x033f, %[cw];\
    orw %[kRealRoundMode], %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    fld1;\
    fpatan;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw), [kRealRoundMode] "r" (kRealRoundMode)
  );

  #endif

  kRealSetError();

  return aux;
}

kReal kRealACot(kReal x)
{
  if(x == 0.0)
  {
    kError = kRealErrorOk;

    return kRealPi / 2.0;
  }
  else
    return kRealATan(1.0 / x);
}

kReal kRealASec(kReal x)
{
  if((x <= -1.0) || (x >= 1.0))
    return kRealATan(kRealSqrt(x * x - 1.0));
  else
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
}

kReal kRealACsc(kReal x)
{
  kError = kRealErrorOk;

  if(x == -1.0)
    return -kRealPi / 2.0;
  else if(x == 1.0)
    return kRealPi / 2.0;
  else if((x < -1.0) || (x > 1.0))
    return kRealATan(1.0 / kRealSqrt(x * x - 1.0));
  else
  {
    kError = kRealErrorDomain;

    return kRealNaN;
  }
}

The functions that round a real are next. They are the floor, the ceiling, and the round to nearest functions.

kReal kRealFloor(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x063f, %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    frndint;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealCeil(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x0a3f, %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    frndint;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw)
  );
  #endif

  kRealSetError();

  return aux;
}

kReal kRealRound(kReal x)
{
  kExtended aux = x;

  #ifdef kCPUIntel
  kUShort cw;

  __asm__
  (
    "finit;\
    fstcw %[cw];\
    andw $0xf0c0, %[cw];\
    orw $0x023f, %[cw];\
    fldcw %[cw];\
    fldt %[aux];\
    frndint;\
    fstpt %[aux];"
    : [aux] "=mr" (aux), [cw] "=m" (cw) : "m" (aux), "m" (cw)
  );
  #endif

  kRealSetError();

  return aux;
}

Finally, the functions that return the minimum and the maximum of two reals close this section.

inline kReal kRealMin(kReal x, kReal y)
{
  kError = kRealErrorOk;

  return x <= y ? x : y;
}

inline kReal kRealMax(kReal x, kReal y)
{
  kError = kRealErrorOk;

  return x >= y ? x : y;
}

Copyright © 2023 Rui Cuco. All rights reserved.