aboutsummaryrefslogtreecommitdiff
path: root/src/benchmarks/cfrac/pmul.c
blob: e69a366299f03a86215fb624f51bc425bed7cf2b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#include "pdefs.h"
#include "precision.h"
#include <string.h>

#ifdef ASM_16BIT
#include "asm16bit.h"
#endif

/*
 *   Multiply u by v (assumes normalized)
 */
precision pmul(u, v)
   register precision v;			/* register a5 on 68000 */
#ifdef ASM_16BIT
   register precision u;			/* register a4 */
{
#else
   precision u;
{
   digitPtr vPtr;
   register digitPtr uPtr, wPtr, HiDigit;
   register accumulator	  temp;	       /* 0 <= temp < base * base */   /* d7 */
   register digit 	  vdigit;				       /* d6 */
#endif
   register digit	  hi; 			/* 0 <= hi < base */   /* d5 */
   precision w;

   (void) pparm(u);
   (void) pparm(v);
   /*
    * Check for multiply by zero.  Helps prevent wasted storage and -0
    */
   if (peqz(u) || peqz(v)) {
      w = palloc(1);
      if (w == pUndef) return w;

      w->sign	  = false;
      w->value[0] = 0;
   } else {
      if (u->size < v->size) { /* u is biggest number (for inner loop speed) */
	 w = u; u = v; v = w;
      }

      w = palloc(u->size + v->size);
      if (w == pUndef) return w;

      w->sign = (u->sign != v->sign);

#ifndef ASM_16BIT
      uPtr = u->value;
      vPtr = v->value;
      wPtr = w->value + u->size;			 /* this is correct! */
      do {
	 *--wPtr = 0;
      } while (wPtr > w->value);

      vPtr    = v->value;
      HiDigit = u->value + u->size;
      do {
	 uPtr	= u->value;
	 wPtr	= w->value + (vPtr - v->value);
	 hi	= 0;
	 vdigit = *vPtr;
	 do {	  
	    temp    = uMul(vdigit, *uPtr++);	/* 0 <= temp <= (base-1)^2   */
	    temp   += *wPtr;			/* 0 <= temp <= base(base-1) */
	    temp   += hi;			/* 0 <= temp < base * base   */
	    hi	    = divBase(temp);		/* 0 <= hi < base	     */
	    *wPtr++ = modBase(temp); 
	 } while (uPtr < HiDigit);
	 *wPtr++    = hi;
      } while (++vPtr < v->value + v->size);
#else
      hi = memmulw(w->value, u->value, u->size, v->value, v->size);
#endif
      if (hi == 0) {
	 --(w->size);				/* normalize */
      }
   }

   pdestroy(u);
   pdestroy(v);
   return presult(w);
}