gcc_qmul.c 1.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
  2. // See https://llvm.org/LICENSE.txt for license information.
  3. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
  4. // long double __gcc_qmul(long double x, long double y);
  5. // This file implements the PowerPC 128-bit double-double multiply operation.
  6. // This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!)
  7. #include "DD.h"
  8. long double __gcc_qmul(long double x, long double y) {
  9. static const uint32_t infinityHi = UINT32_C(0x7ff00000);
  10. DD dst = {.ld = x}, src = {.ld = y};
  11. register double A = dst.s.hi, a = dst.s.lo, B = src.s.hi, b = src.s.lo;
  12. double aHi, aLo, bHi, bLo;
  13. double ab, tmp, tau;
  14. ab = A * B;
  15. // Detect special cases
  16. if (ab == 0.0) {
  17. dst.s.hi = ab;
  18. dst.s.lo = 0.0;
  19. return dst.ld;
  20. }
  21. const doublebits abBits = {.d = ab};
  22. if (((uint32_t)(abBits.x >> 32) & infinityHi) == infinityHi) {
  23. dst.s.hi = ab;
  24. dst.s.lo = 0.0;
  25. return dst.ld;
  26. }
  27. // Generic cases handled here.
  28. aHi = high26bits(A);
  29. bHi = high26bits(B);
  30. aLo = A - aHi;
  31. bLo = B - bHi;
  32. tmp = LOWORDER(ab, aHi, aLo, bHi, bLo);
  33. tmp += (A * b + a * B);
  34. tau = ab + tmp;
  35. dst.s.lo = (ab - tau) + tmp;
  36. dst.s.hi = tau;
  37. return dst.ld;
  38. }