* as Nikolai proposed, start of translation of

    transcendent helpers functions of Bochs for use in the softfloat unit

rtl/inc/f80constant.inc

@@ -0,0 +1,91 @@
+{ The code is translated from Bochs by Florian Klaempfl and contains 
+  the copyright notice below. The original filename is fpu_constant.h
+
+   All changes to the code are copyrighted by the Free Pascal development team and
+   released under the same license as the orginal code, see below.
+}
+
+{============================================================================
+This source file is an extension to the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2b, written for Bochs = x86 achitecture simulator);
+floating point emulation.
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
+been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
+RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
+AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
+COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
+EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
+INSTITUTE = possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
+OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) the source code for the derivative work includes prominent notice that
+the work is derivative, and = 2) the source code includes prominent notice with
+these four paragraphs for those parts of this code that are retained.
+=============================================================================}
+
+// Pentium CPU uses only 68-bit precision M_PI approximation
+//   BETTER_THAN_PENTIUM
+
+{============================================================================
+ * Written for Bochs = x86 achitecture simulator) by
+ *            Stanislav Shwartsman [sshwarts at sourceforge net]
+ * ==========================================================================}
+
+{ the pascal translation is based on https://github.com/lubomyr/bochs/blob/8e0b9abcd81cd24d4d9c68f7fdef2f53bc180d33/cpu/fpu/fpu_constant.h
+
+  if you update it, please make a note here
+
+}
+
+//////////////////////////////
+// PI, PI/2, PI/4 constants
+//////////////////////////////
+
+const
+  FLOATX80_PI_EXP  = $4000;
+
+// 128-bit PI fraction
+{$ifdef BETTER_THAN_PENTIUM}
+  FLOAT_PI_HI = QWord($c90fdaa22168c234);
+  FLOAT_PI_LO = QWord($c4c6628b80dc1cd1);
+{$else}
+  FLOAT_PI_HI = QWord($c90fdaa22168c234);
+  FLOAT_PI_LO = QWord($C000000000000000);
+{$endif}
+
+{$ifdef dummy}
+  FLOATX80_PI2_EXP  = $3FFF);
+  FLOATX80_PI4_EXP  = $3FFE);
+
+//////////////////////////////
+// 3PI/4 constant
+//////////////////////////////
+
+  FLOATX80_3PI4_EXP = $4000;
+
+// 128-bit 3PI/4 fraction
+{$ifdef BETTER_THAN_PENTIUM}
+  FLOAT_3PI4_HI = QWord($96cbe3f9990e91a7);
+  FLOAT_3PI4_LO = QWord($9394c9e8a0a5159c));
+{$else}
+  FLOAT_3PI4_HI = QWord($96cbe3f9990e91a7);
+  FLOAT_3PI4_LO = QWord($9000000000000000);
+{$endif}
+
+//////////////////////////////
+// 1/LN2 constant
+//////////////////////////////
+
+  FLOAT_LN2INV_EXP  = $3FFF);
+
+// 128-bit 1/LN2 fraction
+{$ifdef BETTER_THAN_PENTIUM}
+  FLOAT_LN2INV_HI = QWord($b8aa3b295c17f0bb);
+  FLOAT_LN2INV_LO = QWord($be87fed0691d3e89);
+{$else}
+  FLOAT_LN2INV_HI = QWord($b8aa3b295c17f0bb);
+  FLOAT_LN2INV_LO = QWord($C000000000000000);
+{$endif}
+{$endif dummy}

rtl/inc/f80poly.inc

@@ -0,0 +1,107 @@
+{ The code is translated from Bochs by Florian Klaempfl and contains 
+  the copyright notice below. The original filename is poly.cc
+
+   All changes to the code are copyrighted by the Free Pascal development team and
+   released under the same license as the orginal code, see below.
+}
+
+{============================================================================
+This source file is an extension to the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
+floating point emulation.
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
+been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
+RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
+AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
+COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
+EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
+INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
+OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) the source code for the derivative work includes prominent notice that
+the work is derivative, and (2) the source code includes prominent notice with
+these four paragraphs for those parts of this code that are retained.
+=============================================================================}
+
+{============================================================================
+ * Written for Bochs (x86 achitecture simulator) by
+ *            Stanislav Shwartsman [sshwarts at sourceforge net]
+ * ==========================================================================}
+
+{ the pascal translation is based on https://github.com/lubomyr/bochs/blob/8e0b9abcd81cd24d4d9c68f7fdef2f53bc180d33/cpu/fpu/poly.cc
+
+  if you update it, please make a note here
+
+  Translation notes:
+    * the use of status paramteters has been removed, this is handled directly by the status variable
+}
+
+{$define FLOAT128}
+
+//                            2         3         4               n
+// f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+//         0    1         2         3         4               n
+//
+//          --       2k                --        2k+1
+//   p(x) = >  C  * x           q(x) = >  C   * x
+//          --  2k                     --  2k+1
+//
+//   f(x) ~ [ p(x) + x * q(x) ]
+//
+
+function EvalPoly(x: float128; arr: pfloat128; n: Integer): float128;
+var
+  r: float128;
+begin
+    dec(n);
+    r := arr[n];
+
+    repeat
+        r := float128_mul(r, x);
+        dec(n);
+        r := float128_add(r, arr[n]);
+    until not(n > 0);
+
+    Result := r;
+end;
+
+
+//                  2         4         6         8               2n
+// f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+//         0    1         2         3         4               n
+//
+//          --       4k                --        4k+2
+//   p(x) = >  C  * x           q(x) = >  C   * x
+//          --  2k                     --  2k+1
+//
+//                    2
+//   f(x) ~ [ p(x) + x * q(x) ]
+//
+
+function EvenPoly(x: float128; arr: pfloat128; n: Integer): float128;
+begin
+     Result := EvalPoly(float128_mul(x, x), arr, n);
+end;
+
+
+//                        3         5         7         9               2n+1
+// f(x) ~ (C * x) + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+//          0         1         2         3         4               n
+//                        2         4         6         8               2n
+//      = x * [ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+//               0    1         2         3         4               n
+//
+//          --       4k                --        4k+2
+//   p(x) = >  C  * x           q(x) = >  C   * x
+//          --  2k                     --  2k+1
+//
+//                        2
+//   f(x) ~ x * [ p(x) + x * q(x) ]
+//
+
+function OddPoly(x: float128; arr: pfloat128; n: Integer) : float128;
+begin
+     Result := float128_mul(x, EvenPoly(x, arr, n));
+end;

rtl/inc/f80sincos.inc

@@ -0,0 +1,512 @@
+{ The code is translated from Bochs by Florian Klaempfl and contains 
+  the copyright notice below. The original filename is fsincos.cc
+
+   All changes to the code are copyrighted by the Free Pascal development team and
+   released under the same license as the orginal code, see below.
+}
+
+{============================================================================
+This source file is an extension to the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
+floating point emulation.
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
+been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
+RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
+AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
+COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
+EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
+INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
+OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) the source code for the derivative work includes prominent notice that
+the work is derivative, and (2) the source code includes prominent notice with
+these four paragraphs for those parts of this code that are retained.
+=============================================================================}
+
+{============================================================================
+ * Written for Bochs (x86 achitecture simulator) by
+ *            Stanislav Shwartsman [sshwarts at sourceforge net]
+ * ==========================================================================}
+
+{ the pascal translation is based on https://github.com/lubomyr/bochs/blob/8e0b9abcd81cd24d4d9c68f7fdef2f53bc180d33/cpu/fpu/fsincos.cc
+
+  if you update it, please make a note here
+
+}
+
+{$define FLOAT128}
+
+{$define USE_estimateDiv128To64}
+
+{$I f80constant.inc}
+
+const
+  floatx80_one : floatx80 = 
+     (
+         low: QWord($8000000000000000);
+         high: $3fff;
+     );
+
+{ reduce trigonometric function argument using 128-bit precision
+   M_PI approximation }
+function argument_reduction_kernel(aSig0: Bit64u; Exp: Integer; var zSig0, zSig1: Bit64u) : Bit64u;
+var
+    aSig1, term0, term1, term2, q : Bit64u;
+begin
+    aSig1 := 0;
+
+    shortShift128Left(aSig1, aSig0, Exp, aSig1, aSig0);
+    q := estimateDiv128To64(aSig1, aSig0, FLOAT_PI_HI);
+    mul128By64To192(FLOAT_PI_HI, FLOAT_PI_LO, q, term0, term1, term2);
+    sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
+    while (Bit64s(zSig1) < 0) do
+    begin
+        dec(q);
+        add192(zSig1, zSig0, term2, 0, FLOAT_PI_HI, FLOAT_PI_LO, zSig1, zSig0, term2);
+    end;
+    zSig1 := term2;
+    Result := q;
+end;
+
+function reduce_trig_arg(expDiff: Integer; var zSign: Integer; var aSig0, aSig1: Bit64u) : Integer;
+var
+  term0, term1, q : Bit64u;
+  eq, lt : Integer;
+begin
+    q := 0;
+
+    if (expDiff < 0) then
+    begin
+        shift128Right(aSig0, 0, 1, aSig0, aSig1);
+        expDiff := 0;
+    end;
+    if (expDiff > 0) then
+    begin
+        q := argument_reduction_kernel(aSig0, expDiff, aSig0, aSig1);
+    end
+    else 
+    begin
+        if (FLOAT_PI_HI <= aSig0) then
+        begin
+            dec(aSig0, FLOAT_PI_HI);
+            q := 1;
+        end;
+    end;
+
+    shift128Right(FLOAT_PI_HI, FLOAT_PI_LO, 1, &term0, &term1);
+    if (not(lt128(aSig0, aSig1, term0, term1))<>0) then
+    begin
+        lt := lt128(term0, term1, aSig0, aSig1);
+        eq := eq128(aSig0, aSig1, term0, term1);
+
+        if (((eq<>0) and ((q and 1)<>0)) or (lt<>0)) then
+        begin
+            zSign := not(zSign);
+            inc(q);
+        end;
+        if (lt<>0) then
+          sub128(FLOAT_PI_HI, FLOAT_PI_LO, aSig0, aSig1, &aSig0, &aSig1);
+    end;
+
+    Result:= Integer(q and 3);
+end;
+
+
+const
+  sin_arr: array[0..10] of float128 = (    
+    (low: $3fff000000000000; high: $0000000000000000), {  1 }
+    (low: $bffc555555555555; high: $5555555555555555), {  3 }
+    (low: $3ff8111111111111; high: $1111111111111111), {  5 }
+    (low: $bff2a01a01a01a01; high: $a01a01a01a01a01a), {  7 }
+    (low: $3fec71de3a556c73; high: $38faac1c88e50017), {  9 }
+    (low: $bfe5ae64567f544e; high: $38fe747e4b837dc7), { 11 }
+    (low: $3fde6124613a86d0; high: $97ca38331d23af68), { 13 }
+    (low: $bfd6ae7f3e733b81; high: $f11d8656b0ee8cb0), { 15 }
+    (low: $3fce952c77030ad4; high: $a6b2605197771b00), { 17 }
+    (low: $bfc62f49b4681415; high: $724ca1ec3b7b9675), { 19 }
+    (low: $3fbd71b8ef6dcf57; high: $18bef146fcee6e45)  { 21 }
+  );
+
+
+  cos_arr: array[0..10] of float128 = (
+    (low: $3fff000000000000; high: $0000000000000000), {  0 }
+    (low: $bffe000000000000; high: $0000000000000000), {  2 }
+    (low: $3ffa555555555555; high: $5555555555555555), {  4 }
+    (low: $bff56c16c16c16c1; high: $6c16c16c16c16c17), {  6 }
+    (low: $3fefa01a01a01a01; high: $a01a01a01a01a01a), {  8 }
+    (low: $bfe927e4fb7789f5; high: $c72ef016d3ea6679), { 10 }
+    (low: $3fe21eed8eff8d89; high: $7b544da987acfe85), { 12 }
+    (low: $bfda93974a8c07c9; high: $d20badf145dfa3e5), { 14 }
+    (low: $3fd2ae7f3e733b81; high: $f11d8656b0ee8cb0), { 16 }
+    (low: $bfca6827863b97d9; high: $77bb004886a2c2ab), { 18 }
+    (low: $3fc1e542ba402022; high: $507a9cad2bf8f0bb)  { 20 }
+  );
+
+
+{$I f80poly.inc}
+
+{ 0 <= x <= pi/4 }
+function poly_sin(x: float128): float128;inline;
+begin
+    //                 3     5     7     9     11     13     15
+    //                x     x     x     x     x      x      x
+    // sin (x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- =
+    //                3!    5!    7!    9!    11!    13!    15!
+    //
+    //                 2     4     6     8     10     12     14
+    //                x     x     x     x     x      x      x
+    //   := x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- ] =
+    //                3!    5!    7!    9!    11!    13!    15!
+    //
+    //           3                          3
+    //          --       4k                --        4k+2
+    //   p(x) := >  C  * x   > 0     q(x) := >  C   * x     < 0
+    //          --  2k                     --  2k+1
+    //          k=0                        k=0
+    //
+    //                          2
+    //   sin(x) ~ x * [ p(x) + x * q(x) ]
+    //
+
+    Result := OddPoly(x, sin_arr, length(sin_arr));
+end;
+
+
+{ 0 <= x <= pi/4 }
+function poly_cos(x: float128): float128;inline;
+begin
+    //                 2     4     6     8     10     12     14
+    //                x     x     x     x     x      x      x
+    // cos (x) ~ 1 - --- + --- - --- + --- - ---- + ---- - ----
+    //                2!    4!    6!    8!    10!    12!    14!
+    //
+    //           3                          3
+    //          --       4k                --        4k+2
+    //   p(x) := >  C  * x   > 0     q(x) := >  C   * x     < 0
+    //          --  2k                     --  2k+1
+    //          k=0                        k=0
+    //
+    //                      2
+    //   cos(x) ~ [ p(x) + x * q(x) ]
+    //
+
+    Result := EvenPoly(x, cos_arr, length(cos_arr));
+end;
+
+
+procedure sincos_invalid(sin_a: pfloatx80; cos_a: pfloatx80; a: floatx80);
+begin
+    if assigned(sin_a) then
+      sin_a^ := a;
+    if assigned(cos_a) then
+      cos_a^ := a;
+end;
+
+
+procedure sincos_tiny_argument(sin_a, cos_a: pfloatx80; a: floatx80);
+begin
+    if assigned(sin_a) then
+      sin_a^ := a;
+    if assigned(cos_a) then
+      cos_a^ := floatx80_one;
+end;
+
+
+function floatx80_chs(const a: floatx80): floatx80;
+begin
+    Result := a;
+    Result.high := Result.high xor $8000;
+end;
+
+
+function sincos_approximation(neg: Integer; r: float128; quotient: Bit64u): floatx80;
+begin
+    if (quotient and $1)<>0 then
+    begin
+        r := poly_cos(r);
+        neg := 0;
+    end
+    else  
+    begin
+        r := poly_sin(r);
+    end;
+
+    result := float128_to_floatx80(r);
+    if (quotient and $2)<>0 then
+        neg := neg xor 1;
+
+    if (neg<>0) then
+        floatx80_chs(result);
+end;
+
+// =================================================
+// FSINCOS               Compute sin(x) and cos(x)
+// =================================================
+
+//
+// Uses the following identities:
+// ----------------------------------------------------------
+//
+//  sin(-x) := -sin(x)
+//  cos(-x) :=  cos(x)
+//
+//  sin(x+y) := sin(x)*cos(y)+cos(x)*sin(y)
+//  cos(x+y) := sin(x)*sin(y)+cos(x)*cos(y)
+//
+//  sin(x+ pi/2)  :=  cos(x)
+//  sin(x+ pi)    := -sin(x)
+//  sin(x+3pi/2)  := -cos(x)
+//  sin(x+2pi)    :=  sin(x)
+//
+
+function fsincos(a: floatx80;sin_a, cos_a: pfloatx80): Integer;
+var
+    aSig0, aSig1: Bit64u;
+    aExp, zExp, expDiff: Bit32s;
+    aSign, zSign, q: Integer;
+    r: float128;
+label
+    invalid;    
+begin
+    aSig1 := 0;
+    q := 0;
+
+    // handle unsupported extended double-precision floating encodings
+    if (floatx80_is_unsupported(a)<>0) then
+    begin
+        goto invalid;
+    end;
+
+    aSig0 := extractFloatx80Frac(a);
+    aExp := extractFloatx80Exp(a);
+    aSign := extractFloatx80Sign(a);
+
+    { invalid argument }
+    if (aExp = $7FFF) then
+    begin
+        if (Bit64u(aSig0 shl 1)<>0) then
+        begin
+            sincos_invalid(sin_a, cos_a, propagateFloatx80NaN(a, a)); // not sure if this is correctly translated (FK)
+            Result := 0;
+            Exit;
+        end;
+
+    invalid:
+        float_raise(float_flag_invalid);
+        sincos_invalid(sin_a, cos_a, floatx80_default_nan);
+        Result := 0;
+        Exit;
+    end;
+
+    if (aExp = 0) then
+    begin
+        if (aSig0 = 0) then
+        begin
+            sincos_tiny_argument(sin_a, cos_a, a);
+            Result := 0;
+            Exit;
+        end;
+
+        float_raise(float_flag_denormal);
+
+        { handle pseudo denormals }
+        if (not(aSig0 and QWord($8000000000000000))<>0) then
+        begin
+            float_raise(float_flag_inexact);
+            if assigned(sin_a) then
+                float_raise(float_flag_underflow);
+            sincos_tiny_argument(sin_a, cos_a, a);
+            Result := 0;
+            Exit;
+        end;
+
+        normalizeFloatx80Subnormal(aSig0, aExp, aSig0);
+    end;
+
+    zSign := aSign;
+    zExp := FLOATX80_EXP_BIAS;
+    expDiff := aExp - zExp;
+
+    { argument is out-of-range }
+    if (expDiff >= 63) then    
+        Exit(-1);
+
+    float_raise(float_flag_inexact);
+
+    if (expDiff < -1) then
+    begin    // doesn't require reduction
+        if (expDiff <= -68) then
+        begin
+            a := packFloatx80(aSign, aExp, aSig0);
+            sincos_tiny_argument(sin_a, cos_a, a);
+            Result := 0;
+            Exit;
+        end;
+        zExp := aExp;
+    end
+    else begin
+        q := reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+    end;
+
+    { **************************** }
+    { argument reduction completed }
+    { **************************** }
+
+    { using float128 for approximation }
+    r := normalizeRoundAndPackFloat128(0, zExp-$10, aSig0, aSig1);
+
+    if (aSign<>0) then
+      q := -q;
+    if assigned(sin_a) then
+      sin_a^ := sincos_approximation(zSign, r,   q);
+    if assigned(cos_a) then
+      cos_a^ := sincos_approximation(zSign, r, q+1);
+
+    Result := 0;
+end;
+
+function fsin(var a: floatx80): Integer;
+begin
+    Result := fsincos(a, @a, nil);
+end;
+
+function fcos(var a: floatx80): Integer;
+begin
+    Result := fsincos(a, nil, @a);
+end;
+
+
+// =================================================
+// FPTAN                 Compute tan(x)
+// =================================================
+
+//
+// Uses the following identities:
+//
+// 1. ----------------------------------------------------------
+//
+//  sin(-x) := -sin(x)
+//  cos(-x) :=  cos(x)
+//
+//  sin(x+y) := sin(x)*cos(y)+cos(x)*sin(y)
+//  cos(x+y) := sin(x)*sin(y)+cos(x)*cos(y)
+//
+//  sin(x+ pi/2)  :=  cos(x)
+//  sin(x+ pi)    := -sin(x)
+//  sin(x+3pi/2)  := -cos(x)
+//  sin(x+2pi)    :=  sin(x)
+//
+// 2. ----------------------------------------------------------
+//
+//           sin(x)
+//  tan(x) := ------
+//           cos(x)
+//
+
+function ftan(var a: floatx80): Integer;
+var
+    aSig0, aSig1: Bit64u;
+    aExp, zExp, expDiff: Bit32s;
+    aSign, zSign, q: Integer;    
+    r, cos_r, sin_r: float128;
+label
+    invalid;
+begin
+    aSig1 := 0;
+    q := 0;
+
+    // handle unsupported extended double-precision floating encodings
+    if (floatx80_is_unsupported(a)<>0) then
+    begin
+        goto invalid;
+    end;
+
+    aSig0 := extractFloatx80Frac(a);
+    aExp := extractFloatx80Exp(a);
+    aSign := extractFloatx80Sign(a);
+
+    { invalid argument }
+    if (aExp = $7FFF) then
+    begin
+        if (Bit64u(aSig0 shl 1)<>0) then
+        begin
+            a := propagateFloatx80NaN(a, a); // not sure if this is correctly translated (FK)
+            Result := 0;
+            Exit;
+        end;
+
+    invalid:
+        float_raise(float_flag_invalid);
+        a := floatx80_default_nan;
+        Result := 0;
+        Exit;
+    end;
+
+    if (aExp = 0) then
+    begin
+        if (aSig0 = 0) then
+            Exit(0);
+        float_raise(float_flag_denormal);
+        { handle pseudo denormals }
+        if (not(aSig0 and QWord($8000000000000000))<>0) then
+        begin
+            float_raise([float_flag_inexact,float_flag_underflow]);
+            Result := 0;
+            Exit;
+        end;
+        normalizeFloatx80Subnormal(aSig0, aExp, aSig0);
+    end;
+
+    zSign := aSign;
+    zExp := FLOATX80_EXP_BIAS;
+    expDiff := aExp - zExp;
+
+    { argument is out-of-range }
+    if (expDiff >= 63) then
+        Exit(-1);
+
+    float_raise(float_flag_inexact);
+
+    if (expDiff < -1) then
+    begin    // doesn't require reduction
+        if (expDiff <= -68) then
+        begin
+            a := packFloatx80(aSign, aExp, aSig0);
+            Result := 0;
+            exit;
+        end;
+        zExp := aExp;
+    end
+    else 
+    begin
+        q := reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+    end;
+
+    { **************************** }
+    { argument reduction completed }
+    { **************************** }
+
+    { using float128 for approximation }
+    r := normalizeRoundAndPackFloat128(0, zExp-$10, aSig0, aSig1);
+
+    sin_r := poly_sin(r);
+    cos_r := poly_cos(r);
+
+    if (q and $1)<>0 then
+    begin
+        r := float128_div(cos_r, sin_r);
+        zSign :=  zSign xor 1;
+    end
+    else
+    begin
+        r := float128_div(sin_r, cos_r);
+    end;
+
+    a := float128_to_floatx80(r);
+    if (zSign<>0) then
+        floatx80_chs(a);
+
+    Result := 0;
+end;

rtl/inc/softfpu.pp

@@ -71,6 +71,11 @@ these four paragraphs for those parts of this code that are retained.
 
 { $define FPC_SOFTFLOAT_FLOATX80}
 { $define FPC_SOFTFLOAT_FLOAT128}
+{ $define FPC_SOFTFLOAT_FLOATX80_FUNCS}
+
+{$ifdef FPC_SOFTFLOAT_FLOATX80_TRIG}
+{$define FPC_SOFTFLOAT_FLOAT128}
+{$endif FPC_SOFTFLOAT_FLOATX80_TRIG}
 
 { the softfpu unit can be also embedded directly into the system unit }
 
@@ -112,7 +117,9 @@ TYPE
   bits16 = word;
   sbits16 = smallint;
   sbits32 = longint;
+  Bit32s = sbits32;
   bits32 = longword;
+  Bit32u = bits32;
 {$ifndef fpc}
   qword = int64;
 {$endif}
@@ -120,7 +127,9 @@ TYPE
   uint64 = qword;
   }
   bits64 = qword;
+  Bit64u = bits64;
   sbits64 = int64;
+  Bit64s = sbits64;
 
 {$ifdef ENDIAN_LITTLE}
 {$ifndef FPC_SYSTEM_HAS_float64}
@@ -181,6 +190,9 @@ TYPE
   end;
 {$endif}
 
+pfloat128 = ^float128;
+pfloatx80 = ^floatx80;
+
 {$define FPC_SYSTEM_HAS_float64}
 
 
@@ -6724,6 +6736,13 @@ const
     floatx80_default_nan_high = $FFFF;
     floatx80_default_nan_low  = bits64( $C000000000000000 );
 
+    FLOATX80_EXP_BIAS = $3FFF;
+
+    floatx80_default_nan : floatx80 = 
+      (low:  floatx80_default_nan_low;
+       high: floatx80_default_nan_high);
+
+
 {*----------------------------------------------------------------------------
 | Returns 1 if the extended double-precision floating-point value `a' is a
 | signaling NaN; otherwise returns 0.
@@ -6740,6 +6759,18 @@ begin
        and ( a.low = aLow ) );
 end;
 
+
+{----------------------------------------------------------------------------
+| Returns 1 if the extended double-precision floating-point value `a' is an
+| unsupported; otherwise returns 0.
+*----------------------------------------------------------------------------}
+
+function floatx80_is_unsupported(a: floatx80): Integer;
+begin
+    Result := ord(((a.high and $7FFF)<>0) and not((a.low and QWord($8000000000000000))<>0));
+end;
+
+
 {*----------------------------------------------------------------------------
 | Returns the result of converting the extended double-precision floating-
 | point NaN `a' to the canonical NaN format.  If `a' is a signaling NaN, the
@@ -9376,6 +9407,13 @@ end;
 
 {$endif FPC_SOFTFLOAT_FLOAT128}
 
+
+{$ifdef FPC_SOFTFLOAT_FLOATX80_FUNCS}
+
+{$I f80sincos.inc}
+
+{$endif FPC_SOFTFLOAT_FLOATX80_FUNCS}
+
 {$endif not(defined(fpc_softfpu_interface))}
 
 {$if not(defined(fpc_softfpu_interface)) and not(defined(fpc_softfpu_implementation))}

rtl/linux/Makefile

@@ -179,6 +179,8 @@ else
 ARCH=$(CPU_TARGET)
 endif
 endif
+ifeq ($(FULL_TARGET),aarch64-embedded)
+endif
 ifeq ($(FULL_TARGET),arm-embedded)
 ifeq ($(SUBARCH),)
 $(error When compiling for arm-embedded, a sub-architecture (e.g. SUBARCH=armv4t or SUBARCH=armv7m) must be defined)
@@ -4564,7 +4566,7 @@ ctypes$(PPUEXT) :  $(INC)/ctypes.pp unixtype$(PPUEXT) $(SYSTEMUNIT)$(PPUEXT)
 	$(COMPILER) $(INC)/ctypes.pp
 fpcylix$(PPUEXT) : fpcylix.pp cthreads$(PPUEXT) cwstring$(PPUEXT) dynlibs$(PPUEXT) objpas$(PPUEXT) $(SYSTEMUNIT)$(PPUEXT)
 	$(COMPILER) fpcylix.pp
-softfpu$(PPUEXT) : $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT)
+softfpu$(PPUEXT) : $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT) $(INC)/f80sincos.inc $(INC)/f80constant.inc $(INC)/f80poly.inc
 	$(COMPILER) $<
 sfpux80$(PPUEXT) : $(INC)/sfpux80.pp $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT)
 	$(COMPILER) $<

rtl/linux/Makefile.fpc

@@ -504,7 +504,7 @@ ctypes$(PPUEXT) :  $(INC)/ctypes.pp unixtype$(PPUEXT) $(SYSTEMUNIT)$(PPUEXT)
 fpcylix$(PPUEXT) : fpcylix.pp cthreads$(PPUEXT) cwstring$(PPUEXT) dynlibs$(PPUEXT) objpas$(PPUEXT) $(SYSTEMUNIT)$(PPUEXT)
 	$(COMPILER) fpcylix.pp
 
-softfpu$(PPUEXT) : $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT)
+softfpu$(PPUEXT) : $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT) $(INC)/f80sincos.inc $(INC)/f80constant.inc $(INC)/f80poly.inc
 	$(COMPILER) $<
 
 sfpux80$(PPUEXT) : $(INC)/sfpux80.pp $(INC)/softfpu.pp $(SYSTEMUNIT)$(PPUEXT) objpas$(PPUEXT)