home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / cephes22 / qfloat / qfltb.c < prev    next >
C/C++ Source or Header  |  1992-11-17  |  8KB  |  605 lines
  1. /* 
  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, b )
  218. unsigned short a[], b[];
  219. {
  220. unsigned short sqr[NQ+2], prod[NQ+2], quot[NQ+2];
  221. int i, prec, k;
  222. unsigned short d, qu, *p, *q, *r;
  223. unsigned long u;
  224.  
  225. /* Test if denominator has only 16 bits of significance. */
  226. p = &a[4];
  227. i = NQ-4;
  228. do
  229.     {
  230.     if( *p++ != 0 )
  231.         goto longdiv;
  232.     }
  233. while( --i );
  234.  
  235. /* Do single precision divides if so. */
  236. qmov( b, prod );
  237. prod[NQ] = 0;
  238. prod[NQ+1] = 0;
  239. shdn1( prod );
  240. shdn1( prod );
  241. d = a[3];
  242. u = ((unsigned long )prod[3] << 16) | prod[4];
  243. for( i=3; i<NQ; i++ )
  244.     {
  245.     qu = u / d;
  246.     prod[i] = qu;
  247.     u = ((u - (unsigned long )d * qu) << 16) | prod[i+2];
  248.     }
  249. prod[NQ] = u / d;
  250. goto divdon;
  251.  
  252.  
  253. longdiv:
  254.  
  255. /* Slower procedure is required */
  256. qclear(quot);
  257. quot[NQ] = 0;
  258. qclear(prod);
  259. qclear(sqr);
  260. quot[3] = ((unsigned long )0x40000000) / (unsigned long )a[3];
  261. prec = 2;
  262. k = 1;
  263. while( prec < NQ-2 )
  264.     {
  265.     k = 2 * k;
  266.     if( k > NQ-2 )
  267.         prec = NQ - 2;
  268.     else
  269.         prec = k;
  270.     squarev( quot, sqr, prec );
  271.     mulv( a, sqr, prod, prec );
  272.     subm( prod, quot );
  273.     shup1( quot );
  274.     }
  275. mulv( quot, b, prod, NQ-2 );
  276. prod[0] = b[0];
  277. prod[1] = b[1];
  278.  
  279. divdon:
  280.  
  281. mdnorm( prod );
  282. qmov( prod, b );
  283. }
  284.  
  285. /*
  286. prtemp( s, z )
  287. char *s;
  288. unsigned short z[];
  289. {
  290. int i;
  291.  
  292. printf( "%s ", s );
  293. for( i=0; i<8; i++ )
  294.     printf( "%04x ", z[i+2] );
  295. printf( "\n" );
  296. }
  297. */
  298.  
  299.  
  300. /* Variable precision multiply of significands.
  301.  * c must not be in the same location as either a or b.
  302.  */
  303. static int mulv( a, b, c, prec )
  304. unsigned short a[], b[], c[];
  305. int prec;
  306. {
  307. register unsigned short *p, *q, *r;
  308. register unsigned long u, lp;
  309. unsigned short carry;
  310. int k, i;
  311.  
  312. k = prec+2;
  313. p = &c[2];
  314. do
  315.     *p++ = 0;
  316. while( --k );
  317.  
  318. r = &c[prec+3];
  319. for( k=prec+2; k>=3; k-- )
  320. {
  321. q = &b[3];
  322. p = &a[k];
  323. for( i=k; i>=3; i-- )
  324.     {
  325.     if( (*p == 0) || (*q == 0) )
  326.         {
  327.         --p;
  328.         ++q;
  329.         continue;
  330.         }
  331.     lp = (unsigned long )(*p--) * (unsigned long )(*q++);
  332.     u = (unsigned long )(*r) + (lp & 0xffff);
  333.     if( u & 0x10000 )
  334.         carry = 1;
  335.     else
  336.         carry = 0;
  337.     *r-- = u;
  338.     u = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  339.     *r-- = u;
  340.     if( u & 0x10000 )
  341.         *r += 1;
  342.     r += 2;
  343.     }
  344. --r;
  345. }
  346. }
  347.  
  348.  
  349.  
  350. /* Variable precision square.
  351.  * b must be in a different location from a.
  352.  */
  353. static squarev( a, b, prec )
  354. unsigned short a[], b[];
  355. int prec;
  356. {
  357. register unsigned short *p, *q, *r;
  358. register unsigned long u, lp;
  359. unsigned short carry;
  360. int k, i;
  361.  
  362. k = prec+2;
  363. p = &b[2];
  364. do
  365.     *p++ = 0;
  366. while( --k );
  367.  
  368. r = &b[prec+3];
  369. for( k=prec+2; k>=3; k-- )
  370. {
  371. q = &a[3];
  372. p = &a[k];
  373. while( p >= q )    
  374.     {
  375.     if( (*p == 0) || (*q == 0) )
  376.         {
  377.         --p;
  378.         ++q;
  379.         continue;
  380.         }
  381. /*    printf( "%d %d %d\n", p - &a[3], q - &a[3], r - &b[3] );*/
  382.     lp = (unsigned long )(*p) * (unsigned long )(*q);
  383.     if( p != q )
  384.         {
  385.         if( lp & 0x80000000 )
  386.             *(r-2) += 1;
  387.         lp <<= 1;
  388.         }
  389.     --p;
  390.     ++q;
  391.     u = (unsigned long )(*r) + (lp & 0xffff);
  392.     if( u & 0x10000 )
  393.         carry = 1;
  394.     else
  395.         carry = 0;
  396.     *r-- = u;
  397.     u = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  398.     *r-- = u;
  399.     if( u & 0x10000 )
  400.         *r += 1;
  401.     r += 2;
  402.     }
  403. --r;
  404. }
  405. shup1(b);
  406. }
  407.  
  408.  
  409.  
  410.  
  411.  
  412. mulm( b, ac3 )
  413. unsigned short b[], ac3[];
  414. {
  415. register unsigned short *p, *q;
  416. unsigned short act[NQ+2];
  417. unsigned short carry, *r;
  418. unsigned long lp, a, y;
  419. int i, k, m, o;
  420.  
  421. qclear( act );
  422. act[0] = ac3[0];
  423. act[1] = ac3[1];
  424. act[NQ] = 0;
  425. act[NQ+1] = 0;
  426. r = &act[NQ+1];
  427. for( k=NQ; k>=3; k-- )
  428. {
  429. if( k == NQ )
  430.     {
  431.     m = NQ-1;
  432.     o = 4;
  433.     }
  434. else
  435.     {
  436.     m = k;
  437.     o = 3;
  438.     }
  439. q = &ac3[o];
  440. p = &b[m];
  441.  
  442. for( i=m; i>=o; i-- )
  443.     {
  444.     if( (*p == 0) || (*q == 0) )
  445.         {
  446.         --p;
  447.         ++q;
  448.         continue;
  449.         }
  450.     lp = (unsigned long )(*p--) * (unsigned long )(*q++);
  451.     a = (unsigned long )(*r) + (lp & 0xffff);
  452.     if( a & 0x10000 )
  453.         carry = 1;
  454.     else
  455.         carry = 0;
  456.     *r-- = a;
  457.     a = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  458.     *r-- = a;
  459.     if( a & 0x10000 )
  460.         *r += 1;
  461.     r += 2;
  462.     }
  463. --r;
  464. }
  465. mdnorm( act );
  466. qmov( act, ac3 );
  467. }
  468.  
  469.  
  470.  
  471.  
  472. mulin( b, ac3 )
  473. unsigned short b[], ac3[];
  474. {
  475. register unsigned short *p, *r;
  476. unsigned short act[NQ+1];
  477. unsigned short carry, y;
  478. unsigned long lp, a;
  479. int i;
  480.  
  481. qclear( act );
  482. act[0] = ac3[0];
  483. act[1] = ac3[1];
  484. act[NQ] = 0;
  485. r = &act[NQ];
  486. y = b[3];
  487. p = &ac3[NQ-1];
  488. for( i=NQ-1; i>=3; i-- )
  489.     {
  490.     if( *p == 0 )
  491.         {
  492.         --p;
  493.         --r;
  494.         continue;
  495.         }
  496.     lp = (unsigned long )(*p--) * y;
  497.     a = (unsigned long )(*r) + (lp & 0xffff);
  498.     if( a & 0x10000 )
  499.         carry = 1;
  500.     else
  501.         carry = 0;
  502.     *r-- = a;
  503.     a = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  504.     *r-- = a;
  505.     if( a & 0x10000 )
  506.         *r += 1;
  507.     ++r;
  508.     }
  509. mdnorm( act );
  510. qmov( act, ac3 );
  511. }
  512.  
  513.  
  514. unsigned short rndbit[NQ+1] = {0};
  515. short rndset = 0;
  516.  
  517. mdnorm( x )
  518. unsigned short x[];
  519. {
  520. int i;
  521. register short *p;
  522.  
  523. if( rndset == 0 )
  524.     {
  525.     qclear( rndbit );
  526.     rndbit[NQ-1] = 1;
  527.     rndbit[NQ] = 0;
  528.     rndset = 1;
  529.     }
  530. p = (short *)&x[1];
  531. for( i=0; i<3; i++ )
  532.     {
  533.     if( x[2] == 0 )
  534.         break;
  535.     shdn1(x);
  536.     *p += 1;
  537.     if( *p < 0 )
  538.         *p = 32767; 
  539.     }
  540. for( i=0; i<3; i++ )
  541.     {
  542.     if( x[3] & 0x8000 )
  543.         break;
  544.     shup1(x);
  545.     *p -= 1;
  546.     if( *p < 0 )
  547.         *p = 0;
  548.     }
  549.  
  550. if( x[NQ] & 0x8000 )
  551.     {
  552. /*
  553.     if(  ((x[NQ] & 0x8000) == 0x8000) && ((x[NQ-1] & 1) == 0)  )
  554.         goto nornd;
  555. */
  556.     addm( rndbit, x );
  557.     }
  558. if( x[2] )
  559.     {
  560.     shdn1( x );
  561.     *p += 1;
  562.     if( *p < 0 )
  563.         *p = 32767;
  564.     }
  565. nornd:
  566. x[NQ] = 0;
  567. }
  568.  
  569. /*
  570. ;    move a to b
  571. ;
  572. ;    short a[NQ], b[NQ];
  573. ;    qmov( a, b );
  574. */
  575.  
  576. qmov( a, b )
  577. register short *a, *b;
  578. {
  579. register int i;
  580.  
  581. i = NQ;
  582. do
  583.     {
  584.     *b++ = *a++;
  585.     }
  586. while( --i );
  587. }
  588.  
  589.  
  590. qmovz( a, b )
  591. register short *a, *b;
  592. {
  593. register int i;
  594.  
  595. i = NQ;
  596. do
  597.     {
  598.     *b++ = *a++;
  599.     }
  600. while( --i );
  601. *b++ = 0;
  602. }
  603.  
  604.  
  605.