Hex Artifact Content
Not logged in

Artifact c57a811ec0083746c3b209dd2a46460040689086:


0000: 2f 2a 20 54 68 69 73 20 66 69 6c 65 20 69 73 20  /* This file is 
0010: 64 69 73 74 72 69 62 75 74 65 64 20 75 6e 64 65  distributed unde
0020: 72 20 74 68 65 20 55 6e 69 76 65 72 73 69 74 79  r the University
0030: 20 6f 66 20 49 6c 6c 69 6e 6f 69 73 20 4f 70 65   of Illinois Ope
0040: 6e 20 53 6f 75 72 63 65 0a 20 2a 20 4c 69 63 65  n Source. * Lice
0050: 6e 73 65 2e 20 53 65 65 20 4c 49 43 45 4e 53 45  nse. See LICENSE
0060: 2e 54 58 54 20 66 6f 72 20 64 65 74 61 69 6c 73  .TXT for details
0070: 2e 0a 20 2a 2f 0a 0a 2f 2a 20 6c 6f 6e 67 20 64  .. */../* long d
0080: 6f 75 62 6c 65 20 5f 5f 67 63 63 5f 71 64 69 76  ouble __gcc_qdiv
0090: 28 6c 6f 6e 67 20 64 6f 75 62 6c 65 20 78 2c 20  (long double x, 
00a0: 6c 6f 6e 67 20 64 6f 75 62 6c 65 20 79 29 3b 0a  long double y);.
00b0: 20 2a 20 54 68 69 73 20 66 69 6c 65 20 69 6d 70   * This file imp
00c0: 6c 65 6d 65 6e 74 73 20 74 68 65 20 50 6f 77 65  lements the Powe
00d0: 72 50 43 20 31 32 38 2d 62 69 74 20 64 6f 75 62  rPC 128-bit doub
00e0: 6c 65 2d 64 6f 75 62 6c 65 20 64 69 76 69 73 69  le-double divisi
00f0: 6f 6e 20 6f 70 65 72 61 74 69 6f 6e 2e 0a 20 2a  on operation.. *
0100: 20 54 68 69 73 20 69 6d 70 6c 65 6d 65 6e 74 61   This implementa
0110: 74 69 6f 6e 20 69 73 20 73 68 61 6d 65 6c 65 73  tion is shameles
0120: 73 6c 79 20 63 72 69 62 62 65 64 20 66 72 6f 6d  sly cribbed from
0130: 20 41 70 70 6c 65 27 73 20 44 44 52 54 2c 20 63   Apple's DDRT, c
0140: 69 72 63 61 20 31 39 39 33 28 21 29 0a 20 2a 2f  irca 1993(!). */
0150: 0a 0a 23 69 6e 63 6c 75 64 65 20 22 44 44 2e 68  ..#include "DD.h
0160: 22 0a 0a 6c 6f 6e 67 20 64 6f 75 62 6c 65 20 5f  "..long double _
0170: 5f 67 63 63 5f 71 64 69 76 28 6c 6f 6e 67 20 64  _gcc_qdiv(long d
0180: 6f 75 62 6c 65 20 61 2c 20 6c 6f 6e 67 20 64 6f  ouble a, long do
0190: 75 62 6c 65 20 62 29 0a 7b 09 0a 09 73 74 61 74  uble b).{...stat
01a0: 69 63 20 63 6f 6e 73 74 20 75 69 6e 74 33 32 5f  ic const uint32_
01b0: 74 20 69 6e 66 69 6e 69 74 79 48 69 20 3d 20 55  t infinityHi = U
01c0: 49 4e 54 33 32 5f 43 28 30 78 37 66 66 30 30 30  INT32_C(0x7ff000
01d0: 30 30 29 3b 0a 09 44 44 20 64 73 74 20 3d 20 7b  00);..DD dst = {
01e0: 20 2e 6c 64 20 3d 20 61 20 7d 2c 20 73 72 63 20   .ld = a }, src 
01f0: 3d 20 7b 20 2e 6c 64 20 3d 20 62 20 7d 3b 0a 09  = { .ld = b };..
0200: 0a 09 72 65 67 69 73 74 65 72 20 64 6f 75 62 6c  ..register doubl
0210: 65 20 78 20 3d 20 64 73 74 2e 73 2e 68 69 2c 20  e x = dst.s.hi, 
0220: 78 31 20 3d 20 64 73 74 2e 73 2e 6c 6f 2c 0a 09  x1 = dst.s.lo,..
0230: 09 09 09 09 79 20 3d 20 73 72 63 2e 73 2e 68 69  ....y = src.s.hi
0240: 2c 20 79 31 20 3d 20 73 72 63 2e 73 2e 6c 6f 3b  , y1 = src.s.lo;
0250: 0a 09 0a 20 20 20 20 64 6f 75 62 6c 65 20 79 48  ...    double yH
0260: 69 2c 20 79 4c 6f 2c 20 71 48 69 2c 20 71 4c 6f  i, yLo, qHi, qLo
0270: 3b 0a 20 20 20 20 64 6f 75 62 6c 65 20 79 71 2c  ;.    double yq,
0280: 20 74 6d 70 2c 20 71 3b 0a 09 0a 20 20 20 20 71   tmp, q;...    q
0290: 20 3d 20 78 20 2f 20 79 3b 0a 09 0a 09 2f 2a 20   = x / y;..../* 
02a0: 44 65 74 65 63 74 20 73 70 65 63 69 61 6c 20 63  Detect special c
02b0: 61 73 65 73 20 2a 2f 0a 09 69 66 20 28 71 20 3d  ases */..if (q =
02c0: 3d 20 30 2e 30 29 20 7b 0a 09 09 64 73 74 2e 73  = 0.0) {...dst.s
02d0: 2e 68 69 20 3d 20 71 3b 0a 09 09 64 73 74 2e 73  .hi = q;...dst.s
02e0: 2e 6c 6f 20 3d 20 30 2e 30 3b 0a 09 09 72 65 74  .lo = 0.0;...ret
02f0: 75 72 6e 20 64 73 74 2e 6c 64 3b 0a 09 7d 0a 09  urn dst.ld;..}..
0300: 0a 09 63 6f 6e 73 74 20 64 6f 75 62 6c 65 62 69  ..const doublebi
0310: 74 73 20 71 42 69 74 73 20 3d 20 7b 20 2e 64 20  ts qBits = { .d 
0320: 3d 20 71 20 7d 3b 0a 09 69 66 20 28 28 28 75 69  = q };..if (((ui
0330: 6e 74 33 32 5f 74 29 28 71 42 69 74 73 2e 78 20  nt32_t)(qBits.x 
0340: 3e 3e 20 33 32 29 20 26 20 69 6e 66 69 6e 69 74  >> 32) & infinit
0350: 79 48 69 29 20 3d 3d 20 69 6e 66 69 6e 69 74 79  yHi) == infinity
0360: 48 69 29 20 7b 0a 09 09 64 73 74 2e 73 2e 68 69  Hi) {...dst.s.hi
0370: 20 3d 20 71 3b 0a 09 09 64 73 74 2e 73 2e 6c 6f   = q;...dst.s.lo
0380: 20 3d 20 30 2e 30 3b 0a 09 09 72 65 74 75 72 6e   = 0.0;...return
0390: 20 64 73 74 2e 6c 64 3b 0a 09 7d 0a 09 0a 20 20   dst.ld;..}...  
03a0: 20 20 79 48 69 20 3d 20 68 69 67 68 32 36 62 69    yHi = high26bi
03b0: 74 73 28 79 29 3b 0a 20 20 20 20 71 48 69 20 3d  ts(y);.    qHi =
03c0: 20 68 69 67 68 32 36 62 69 74 73 28 71 29 3b 0a   high26bits(q);.
03d0: 09 0a 20 20 20 20 79 71 20 3d 20 79 20 2a 20 71  ..    yq = y * q
03e0: 3b 0a 20 20 20 20 79 4c 6f 20 3d 20 79 20 2d 20  ;.    yLo = y - 
03f0: 79 48 69 3b 0a 20 20 20 20 71 4c 6f 20 3d 20 71  yHi;.    qLo = q
0400: 20 2d 20 71 48 69 3b 0a 09 0a 20 20 20 20 74 6d   - qHi;...    tm
0410: 70 20 3d 20 4c 4f 57 4f 52 44 45 52 28 79 71 2c  p = LOWORDER(yq,
0420: 20 79 48 69 2c 20 79 4c 6f 2c 20 71 48 69 2c 20   yHi, yLo, qHi, 
0430: 71 4c 6f 29 3b 0a 20 20 20 20 74 6d 70 20 3d 20  qLo);.    tmp = 
0440: 28 78 20 2d 20 79 71 29 20 2d 20 74 6d 70 3b 0a  (x - yq) - tmp;.
0450: 20 20 20 20 74 6d 70 20 3d 20 28 28 74 6d 70 20      tmp = ((tmp 
0460: 2b 20 78 31 29 20 2d 20 79 31 20 2a 20 71 29 20  + x1) - y1 * q) 
0470: 2f 20 79 3b 0a 20 20 20 20 78 20 3d 20 71 20 2b  / y;.    x = q +
0480: 20 74 6d 70 3b 0a 09 0a 20 20 20 20 64 73 74 2e   tmp;...    dst.
0490: 73 2e 6c 6f 20 3d 20 28 71 20 2d 20 78 29 20 2b  s.lo = (q - x) +
04a0: 20 74 6d 70 3b 0a 20 20 20 20 64 73 74 2e 73 2e   tmp;.    dst.s.
04b0: 68 69 20 3d 20 78 3b 0a 09 0a 20 20 20 20 72 65  hi = x;...    re
04c0: 74 75 72 6e 20 64 73 74 2e 6c 64 3b 0a 7d 0a     turn dst.ld;.}.