home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qflta.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  5.2 KB  |  429 lines
  1. /*        qflta.c
  2.  * Utilities for extended precision arithmetic, called by qflt.c.
  3.  * These should all be written in machine language for speed.
  4.  *
  5.  * addm( x, y )        add significand of x to that of y
  6.  * shdn1( x )        shift significand of x down 1 bit
  7.  * shdn8( x )        shift significand of x down 8 bits
  8.  * shdn16( x )        shift significand of x down 16 bits
  9.  * shup1( x )        shift significand of x up 1 bit
  10.  * shup8( x )        shift significand of x up 8 bits
  11.  * shup16( x )        shift significand of x up 16 bits
  12.  * divm( a, b )        divide significand of a into b
  13.  * mulm( a, b )        multiply significands, result in b
  14.  * mdnorm( x )        normalize and round off
  15.  *
  16.  * Copyright (c) 1984 - 1988 by Stephen L. Moshier.  All rights reserved.
  17.  */
  18.  
  19. #include "qhead.h"
  20. #define N (NQ-2)
  21.  
  22. /*
  23. ;    Shift mantissa down by 1 bit
  24. */
  25.  
  26. shdn1(x)
  27. register unsigned short *x;
  28. {
  29. register short bits;
  30. int i;
  31.  
  32. x += 2;    /* point to mantissa area */
  33.  
  34. bits = 0;
  35. for( i=0; i<NQ-1; i++ )
  36.     {
  37.     if( *x & 1 )
  38.         bits |= 1;
  39.     *x >>= 1;
  40.     if( bits & 2 )
  41.         *x |= 0x8000;
  42.     bits <<= 1;
  43.     ++x;
  44.     }    
  45. }
  46.  
  47. /*
  48. ;    Shift mantissa up by 1 bit
  49. */
  50. shup1(x)
  51. register unsigned short *x;
  52. {
  53. register short bits;
  54. int i;
  55.  
  56. x += NQ;
  57. bits = 0;
  58.  
  59. for( i=0; i<NQ-1; i++ )
  60.     {
  61.     if( *x & 0x8000 )
  62.         bits |= 1;
  63.     *x <<= 1;
  64.     if( bits & 2 )
  65.         *x |= 1;
  66.     bits <<= 1;
  67.     --x;
  68.     }
  69. }
  70.  
  71.  
  72.  
  73. /*
  74. ;    Shift mantissa down by 8 bits
  75. */
  76.  
  77. shdn8(x)
  78. register unsigned short *x;
  79. {
  80. register unsigned short newbyt, oldbyt;
  81. int i;
  82.  
  83. x += 2;
  84. oldbyt = 0;
  85. for( i=0; i<NQ-1; i++ )
  86.     {
  87.     newbyt = *x << 8;
  88.     *x >>= 8;
  89.     *x |= oldbyt;
  90.     oldbyt = newbyt;
  91.     ++x;
  92.     }
  93. }
  94.  
  95. /*
  96. ;    Shift mantissa up by 8 bits
  97. */
  98.  
  99. shup8(x)
  100. register unsigned short *x;
  101. {
  102. int i;
  103. register unsigned short newbyt, oldbyt;
  104.  
  105. x += NQ;
  106. oldbyt = 0;
  107.  
  108. for( i=0; i<NQ-1; i++ )
  109.     {
  110.     newbyt = *x >> 8;
  111.     *x <<= 8;
  112.     *x |= oldbyt;
  113.     oldbyt = newbyt;
  114.     --x;
  115.     }
  116. }
  117.  
  118.  
  119.  
  120. /*
  121. ;    Shift mantissa up by 16 bits
  122. */
  123.  
  124. shup16(x)
  125. register short *x;
  126. {
  127. int i;
  128. register short *p;
  129.  
  130. p = x+2;
  131. x += 3;
  132.  
  133. for( i=0; i<NQ-2; i++ )
  134.     *p++ = *x++;
  135.  
  136. *p = 0;
  137. }
  138.  
  139. /*
  140. ;    Shift mantissa down by 16 bits
  141. */
  142.  
  143. shdn16(x)
  144. register unsigned short *x;
  145. {
  146. int i;
  147. register unsigned short *p;
  148.  
  149. x += NQ;
  150. p = x+1;
  151.  
  152. for( i=0; i<NQ-2; i++ )
  153.     *(--p) = *(--x);
  154.  
  155. *(--p) = 0;
  156. }
  157.  
  158.  
  159.  
  160. /*
  161. ;    add mantissas
  162. ;    x + y replaces y
  163. */
  164.  
  165. addm( x, y )
  166. unsigned short *x, *y;
  167. {
  168. register unsigned long a;
  169. int i;
  170. unsigned int carry;
  171.  
  172. x += NQ;
  173. y += NQ;
  174. carry = 0;
  175. for( i=0; i<NQ-1; i++ )
  176.     {
  177.     a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
  178.     if( a & 0x10000 )
  179.         carry = 1;
  180.     else
  181.         carry = 0;
  182.     *y = a;
  183.     --x;
  184.     --y;
  185.     }
  186. }
  187.  
  188. /*
  189. ;    subtract mantissas
  190. ;    y - x replaces y
  191. */
  192.  
  193. subm( x, y )
  194. unsigned short *x, *y;
  195. {
  196. register unsigned long a;
  197. int i;
  198. unsigned int carry;
  199.  
  200. x += NQ;
  201. y += NQ;
  202. carry = 0;
  203. for( i=0; i<NQ-1; i++ )
  204.     {
  205.     a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
  206.     if( a & 0x10000 )
  207.         carry = 1;
  208.     else
  209.         carry = 0;
  210.     *y = a;
  211.     --x;
  212.     --y;
  213.     }
  214. }
  215.  
  216.  
  217. divm( a, ac3 )
  218. unsigned short a[], ac3[];
  219. {
  220. unsigned short act[NQ+1], ac1[NQ+1];
  221. int i, ctr, lost;
  222. unsigned short d, *p, *q, *r;
  223.  
  224. qmovz( a, ac1 );
  225. qclear( act );
  226. act[0] = ac3[0];
  227. act[1] = ac3[1];
  228. act[NQ] = 0;
  229. ac3[NQ] = 0;
  230.  
  231. /* test for word-precision denominator */
  232. for( i=4; i<NQ; i++ )
  233.     {
  234.     if( ac1[i] )
  235.         goto longdiv;
  236.     }
  237. /* Do divide with faster compare and subtract */
  238. d = ac1[3];
  239. p = &ac3[3];
  240. q = &ac3[2];
  241. r = &act[NQ];
  242. for( i=0; i<NBITS+2; i++ )
  243.     {
  244.     if( (*q != 0) || (*p >= d) )
  245.         {
  246.         *p -= d;
  247.         *q = 0;
  248.         *r = 0x8000;
  249.         }
  250.     shup1( ac3 );
  251.     shup1( act );
  252.     }
  253. goto divdon;
  254.  
  255. /* Slower compare and subtract required */
  256. longdiv:
  257. for( i=0; i<NBITS+2; i++ )
  258.     {
  259.     if( cmpm( ac3, ac1 ) >= 0 )
  260.         {
  261.         subm( ac1, ac3 );
  262.         ctr = 1;
  263.         }
  264.     else
  265.         ctr = 0;
  266.     shup1( ac3 );
  267.     shup1( act );
  268.     act[NQ-1] |= ctr;
  269.     }
  270.  
  271. divdon:
  272.  
  273. p = &ac3[2];
  274. lost = 0;
  275. for( i=2; i<NQ; i++ )
  276.     {
  277.     if( *p++ != 0 )
  278.         {
  279.         lost = 1;
  280.         break;
  281.         }
  282.     }
  283. shdn1( act );
  284. shdn1( act );
  285. act[1] -= 1;
  286. if( act[1] & 0x8000 )
  287.     act[1] = 0;
  288. mdnorm( act, lost );
  289. qmov( act, ac3 );
  290. }
  291.  
  292.  
  293.  
  294.  
  295. mulm( b, ac3 )
  296. unsigned short b[], ac3[];
  297. {
  298. unsigned short act[NQ+1], ac2[NQ+1];
  299. unsigned short *p, *q;
  300. int ctr, nct, lost;
  301.  
  302. qmov( b, ac2 );
  303. qclear( act );
  304. act[0] = ac3[0];
  305. act[1] = ac3[1];
  306. act[NQ] = 0;
  307. /* strip trailing zero bits */
  308. nct = NBITS;
  309. p = &ac2[NQ-1];
  310. while( *p == 0 )
  311.     {
  312.     shdn16( ac2 );
  313.     nct -= 16;
  314.     }
  315. if( (*p & 0xff) == 0 )
  316.     {
  317.     shdn8( ac2 );
  318.     nct -= 8;
  319.     }
  320. lost = 0;
  321. q = &act[NQ];
  322. for( ctr=0; ctr<nct; ctr++ )
  323.     {
  324.     if( *p & 1 )
  325.         addm(ac3, act);
  326.     shdn1(ac2);
  327.     lost |= *q & 1;
  328.     shdn1(act);
  329.     }
  330. mdnorm( act, lost );
  331. qmov( act, ac3 );
  332. }
  333.  
  334.  
  335.  
  336.  
  337. mulin( a, b )
  338. unsigned short a[], b[];
  339. {
  340. mulm(a,b);
  341. }
  342.  
  343.  
  344. unsigned short rndbit[NQ+1] = {0};
  345. short rndset = 0;
  346.  
  347. mdnorm( x, lost )
  348. unsigned short x[];
  349. int lost;
  350. {
  351. int i;
  352. register short *p;
  353.  
  354. if( rndset == 0 )
  355.     {
  356.     qclear( rndbit );
  357.     rndbit[NQ-1] = 1;
  358.     rndbit[NQ] = 0;
  359.     rndset = 1;
  360.     }
  361. p = (short *)&x[1];
  362. for( i=0; i<3; i++ )
  363.     {
  364.     if( x[2] == 0 )
  365.         break;
  366.     shdn1(x);
  367.     *p += 1;
  368.     if( *p < 0 )
  369.         *p = 32767; 
  370.     }
  371. for( i=0; i<3; i++ )
  372.     {
  373.     if( x[3] & 0x8000 )
  374.         break;
  375.     shup1(x);
  376.     *p -= 1;
  377.     if( *p < 0 )
  378.         *p = 0;
  379.     }
  380. i = x[NQ] & 0xffff;
  381. if( i & 0x8000 )
  382.     {
  383.     if( (i == 0x8000) && (lost == 0) )
  384.         {
  385.         if( (x[NQ-1] & 1) == 0 )
  386.             goto nornd;
  387.         }
  388.     addm( rndbit, x );
  389.     }
  390. if( x[2] )
  391.     {
  392.     shdn1( x );
  393.     *p += 1;
  394.     if( *p < 0 )
  395.         *p = 32767;
  396.     }
  397. nornd:
  398. x[NQ] = 0;
  399. }
  400.  
  401. /*
  402. ;    move a to b
  403. ;
  404. ;    short a[NQ], b[NQ];
  405. ;    qmov( a, b );
  406. */
  407.  
  408. qmov( a, b )
  409. register short *a, *b;
  410. {
  411. int i;
  412.  
  413. for( i=0; i<NQ; i++ )
  414.     *b++ = *a++;
  415. }
  416.  
  417.  
  418. qmovz( a, b )
  419. register short *a, *b;
  420. {
  421. int i;
  422.  
  423. for( i=0; i<NQ; i++ )
  424.     *b++ = *a++;
  425. *b++ = 0;
  426. }
  427.  
  428.  
  429.