1/*
2** emfloat.c
3** Source for emulated floating-point routines.
4** BYTEmark (tm)
5** BYTE's Native Mode Benchmarks
6** Rick Grehan, BYTE Magazine.
7**
8** Created:
9** Last update: 3/95
10**
11** DISCLAIMER
12** The source, executable, and documentation files that comprise
13** the BYTEmark benchmarks are made available on an "as is" basis.
14** This means that we at BYTE Magazine have made every reasonable
15** effort to verify that the there are no errors in the source and
16** executable code.  We cannot, however, guarantee that the programs
17** are error-free.  Consequently, McGraw-HIll and BYTE Magazine make
18** no claims in regard to the fitness of the source code, executable
19** code, and documentation of the BYTEmark.
20**  Furthermore, BYTE Magazine, McGraw-Hill, and all employees
21** of McGraw-Hill cannot be held responsible for any damages resulting
22** from the use of this code or the results obtained from using
23** this code.
24*/
25
26#include "../pub/libvex_basictypes.h"
27
28static HWord (*serviceFn)(HWord,HWord) = 0;
29
30
31/////////////////////////////////////////////////////////////////////
32/////////////////////////////////////////////////////////////////////
33
34static char* my_strcpy ( char* dest, const char* src )
35{
36   char* dest_orig = dest;
37   while (*src) *dest++ = *src++;
38   *dest = 0;
39   return dest_orig;
40}
41
42static void* my_memcpy ( void *dest, const void *src, int sz )
43{
44   const char *s = (const char *)src;
45   char *d = (char *)dest;
46
47   while (sz--)
48      *d++ = *s++;
49
50   return dest;
51}
52
53static void* my_memmove( void *dst, const void *src, unsigned int len )
54{
55    register char *d;
56    register char *s;
57    if ( dst > src ) {
58        d = (char *)dst + len - 1;
59        s = (char *)src + len - 1;
60        while ( len >= 4 ) {
61            *d-- = *s--;
62            *d-- = *s--;
63            *d-- = *s--;
64            *d-- = *s--;
65            len -= 4;
66        }
67        while ( len-- ) {
68            *d-- = *s--;
69        }
70    } else if ( dst < src ) {
71        d = (char *)dst;
72        s = (char *)src;
73        while ( len >= 4 ) {
74            *d++ = *s++;
75            *d++ = *s++;
76            *d++ = *s++;
77            *d++ = *s++;
78            len -= 4;
79        }
80        while ( len-- ) {
81            *d++ = *s++;
82        }
83    }
84    return dst;
85}
86
87/////////////////////////////////////////////////////////////////////
88
89static void vexxx_log_bytes ( char* p, int n )
90{
91   int i;
92   for (i = 0; i < n; i++)
93      (*serviceFn)( 1, (int)p[i] );
94}
95
96/*---------------------------------------------------------*/
97/*--- vexxx_printf                                        ---*/
98/*---------------------------------------------------------*/
99
100/* This should be the only <...> include in the entire VEXXX library.
101   New code for vexxx_util.c should go above this point. */
102#include <stdarg.h>
103
104static HChar vexxx_toupper ( HChar c )
105{
106   if (c >= 'a' && c <= 'z')
107      return toHChar(c + ('A' - 'a'));
108   else
109      return c;
110}
111
112static Int vexxx_strlen ( const HChar* str )
113{
114   Int i = 0;
115   while (str[i] != 0) i++;
116   return i;
117}
118
119Bool vexxx_streq ( const HChar* s1, const HChar* s2 )
120{
121   while (True) {
122      if (*s1 == 0 && *s2 == 0)
123         return True;
124      if (*s1 != *s2)
125         return False;
126      s1++;
127      s2++;
128   }
129}
130
131/* Some flags.  */
132#define VG_MSG_SIGNED    1 /* The value is signed. */
133#define VG_MSG_ZJUSTIFY  2 /* Must justify with '0'. */
134#define VG_MSG_LJUSTIFY  4 /* Must justify on the left. */
135#define VG_MSG_PAREN     8 /* Parenthesize if present (for %y) */
136#define VG_MSG_COMMA    16 /* Add commas to numbers (for %d, %u) */
137
138/* Copy a string into the buffer. */
139static UInt
140myvprintf_str ( void(*send)(HChar), Int flags, Int width, HChar* str,
141                Bool capitalise )
142{
143#  define MAYBE_TOUPPER(ch) toHChar(capitalise ? vexxx_toupper(ch) : (ch))
144   UInt ret = 0;
145   Int i, extra;
146   Int len = vexxx_strlen(str);
147
148   if (width == 0) {
149      ret += len;
150      for (i = 0; i < len; i++)
151         send(MAYBE_TOUPPER(str[i]));
152      return ret;
153   }
154
155   if (len > width) {
156      ret += width;
157      for (i = 0; i < width; i++)
158         send(MAYBE_TOUPPER(str[i]));
159      return ret;
160   }
161
162   extra = width - len;
163   if (flags & VG_MSG_LJUSTIFY) {
164      ret += extra;
165      for (i = 0; i < extra; i++)
166         send(' ');
167   }
168   ret += len;
169   for (i = 0; i < len; i++)
170      send(MAYBE_TOUPPER(str[i]));
171   if (!(flags & VG_MSG_LJUSTIFY)) {
172      ret += extra;
173      for (i = 0; i < extra; i++)
174         send(' ');
175   }
176
177#  undef MAYBE_TOUPPER
178
179   return ret;
180}
181
182/* Write P into the buffer according to these args:
183 *  If SIGN is true, p is a signed.
184 *  BASE is the base.
185 *  If WITH_ZERO is true, '0' must be added.
186 *  WIDTH is the width of the field.
187 */
188static UInt
189myvprintf_int64 ( void(*send)(HChar), Int flags, Int base, Int width, ULong pL)
190{
191   HChar buf[40];
192   Int   ind = 0;
193   Int   i, nc = 0;
194   Bool  neg = False;
195   HChar *digits = "0123456789ABCDEF";
196   UInt  ret = 0;
197   UInt  p = (UInt)pL;
198
199   if (base < 2 || base > 16)
200      return ret;
201
202   if ((flags & VG_MSG_SIGNED) && (Int)p < 0) {
203      p   = - (Int)p;
204      neg = True;
205   }
206
207   if (p == 0)
208      buf[ind++] = '0';
209   else {
210      while (p > 0) {
211         if ((flags & VG_MSG_COMMA) && 10 == base &&
212             0 == (ind-nc) % 3 && 0 != ind)
213         {
214            buf[ind++] = ',';
215            nc++;
216         }
217         buf[ind++] = digits[p % base];
218         p /= base;
219      }
220   }
221
222   if (neg)
223      buf[ind++] = '-';
224
225   if (width > 0 && !(flags & VG_MSG_LJUSTIFY)) {
226      for(; ind < width; ind++) {
227	//vassert(ind < 39);
228         buf[ind] = toHChar((flags & VG_MSG_ZJUSTIFY) ? '0': ' ');
229      }
230   }
231
232   /* Reverse copy to buffer.  */
233   ret += ind;
234   for (i = ind -1; i >= 0; i--) {
235      send(buf[i]);
236   }
237   if (width > 0 && (flags & VG_MSG_LJUSTIFY)) {
238      for(; ind < width; ind++) {
239	 ret++;
240         send(' ');  // Never pad with zeroes on RHS -- changes the value!
241      }
242   }
243   return ret;
244}
245
246
247/* A simple vprintf().  */
248static
249UInt vprintf_wrk ( void(*send)(HChar), const HChar *format, va_list vargs )
250{
251   UInt ret = 0;
252   int i;
253   int flags;
254   int width;
255   Bool is_long;
256
257   /* We assume that vargs has already been initialised by the
258      caller, using va_start, and that the caller will similarly
259      clean up with va_end.
260   */
261
262   for (i = 0; format[i] != 0; i++) {
263      if (format[i] != '%') {
264         send(format[i]);
265	 ret++;
266         continue;
267      }
268      i++;
269      /* A '%' has been found.  Ignore a trailing %. */
270      if (format[i] == 0)
271         break;
272      if (format[i] == '%') {
273         /* `%%' is replaced by `%'. */
274         send('%');
275	 ret++;
276         continue;
277      }
278      flags = 0;
279      is_long = False;
280      width = 0; /* length of the field. */
281      if (format[i] == '(') {
282	 flags |= VG_MSG_PAREN;
283	 i++;
284      }
285      /* If ',' follows '%', commas will be inserted. */
286      if (format[i] == ',') {
287         flags |= VG_MSG_COMMA;
288         i++;
289      }
290      /* If '-' follows '%', justify on the left. */
291      if (format[i] == '-') {
292         flags |= VG_MSG_LJUSTIFY;
293         i++;
294      }
295      /* If '0' follows '%', pads will be inserted. */
296      if (format[i] == '0') {
297         flags |= VG_MSG_ZJUSTIFY;
298         i++;
299      }
300      /* Compute the field length. */
301      while (format[i] >= '0' && format[i] <= '9') {
302         width *= 10;
303         width += format[i++] - '0';
304      }
305      while (format[i] == 'l') {
306         i++;
307         is_long = True;
308      }
309
310      switch (format[i]) {
311         case 'd': /* %d */
312            flags |= VG_MSG_SIGNED;
313            if (is_long)
314               ret += myvprintf_int64(send, flags, 10, width,
315				      (ULong)(va_arg (vargs, Long)));
316            else
317               ret += myvprintf_int64(send, flags, 10, width,
318				      (ULong)(va_arg (vargs, Int)));
319            break;
320         case 'u': /* %u */
321            if (is_long)
322               ret += myvprintf_int64(send, flags, 10, width,
323				      (ULong)(va_arg (vargs, ULong)));
324            else
325               ret += myvprintf_int64(send, flags, 10, width,
326				      (ULong)(va_arg (vargs, UInt)));
327            break;
328         case 'p': /* %p */
329	    ret += 2;
330            send('0');
331            send('x');
332            ret += myvprintf_int64(send, flags, 16, width,
333				   (ULong)((HWord)va_arg (vargs, void *)));
334            break;
335         case 'x': /* %x */
336            if (is_long)
337               ret += myvprintf_int64(send, flags, 16, width,
338				      (ULong)(va_arg (vargs, ULong)));
339            else
340               ret += myvprintf_int64(send, flags, 16, width,
341				      (ULong)(va_arg (vargs, UInt)));
342            break;
343         case 'c': /* %c */
344	    ret++;
345            send(toHChar(va_arg (vargs, int)));
346            break;
347         case 's': case 'S': { /* %s */
348            char *str = va_arg (vargs, char *);
349            if (str == (char*) 0) str = "(null)";
350            ret += myvprintf_str(send, flags, width, str,
351                                 toBool(format[i]=='S'));
352            break;
353	 }
354#        if 0
355	 case 'y': { /* %y - print symbol */
356	    Addr a = va_arg(vargs, Addr);
357
358            HChar *name;
359	    if (VG_(get_fnname_w_offset)(a, &name)) {
360               HChar buf[1 + VG_strlen(name) + 1 + 1];
361	       if (flags & VG_MSG_PAREN) {
362                  VG_(sprintf)(str, "(%s)", name):
363	       } else {
364                  VG_(sprintf)(str, "%s", name):
365               }
366	       ret += myvprintf_str(send, flags, width, buf, 0);
367	    }
368	    break;
369	 }
370#        endif
371         default:
372            break;
373      }
374   }
375   return ret;
376}
377
378
379/* A general replacement for printf().  Note that only low-level
380   debugging info should be sent via here.  The official route is to
381   to use vg_message().  This interface is deprecated.
382*/
383static HChar myprintf_buf[1000];
384static Int   n_myprintf_buf;
385
386static void add_to_myprintf_buf ( HChar c )
387{
388   if (c == '\n' || n_myprintf_buf >= 1000-10 /*paranoia*/ ) {
389      (*vexxx_log_bytes)( myprintf_buf, vexxx_strlen(myprintf_buf) );
390      n_myprintf_buf = 0;
391      myprintf_buf[n_myprintf_buf] = 0;
392   }
393   myprintf_buf[n_myprintf_buf++] = c;
394   myprintf_buf[n_myprintf_buf] = 0;
395}
396
397static UInt vexxx_printf ( const char *format, ... )
398{
399   UInt ret;
400   va_list vargs;
401   va_start(vargs,format);
402
403   n_myprintf_buf = 0;
404   myprintf_buf[n_myprintf_buf] = 0;
405   ret = vprintf_wrk ( add_to_myprintf_buf, format, vargs );
406
407   if (n_myprintf_buf > 0) {
408      (*vexxx_log_bytes)( myprintf_buf, n_myprintf_buf );
409   }
410
411   va_end(vargs);
412
413   return ret;
414}
415
416/*---------------------------------------------------------------*/
417/*--- end                                          vexxx_util.c ---*/
418/*---------------------------------------------------------------*/
419
420
421/////////////////////////////////////////////////////////////////////
422/////////////////////////////////////////////////////////////////////
423
424//#include <stdio.h>
425//#include <string.h>
426//#include <malloc.h>
427
428typedef unsigned char uchar;
429typedef unsigned int uint;
430typedef unsigned short ushort;
431typedef unsigned long ulong;
432typedef int int32;              /* Signed 32 bit integer */
433
434#define INTERNAL_FPF_PRECISION 4
435#define CPUEMFLOATLOOPMAX 500000L
436#define EMFARRAYSIZE 3000L
437
438typedef struct {
439        int adjust;             /* Set adjust code */
440        ulong request_secs;     /* # of seconds requested */
441        ulong arraysize;        /* Size of array */
442        ulong loops;            /* Loops per iterations */
443        double emflops;         /* Results */
444} EmFloatStruct;
445
446
447
448/* Is this a 64 bit architecture? If so, this will define LONG64 */
449/* Uwe F. Mayer 15 November 1997                                 */
450// #include "pointer.h"
451
452#define u8 unsigned char
453#define u16 unsigned short
454#ifdef LONG64
455#define u32 unsigned int
456#else
457#define u32 unsigned long
458#endif
459#define uchar unsigned char
460#define ulong unsigned long
461
462#define MAX_EXP 32767L
463#define MIN_EXP (-32767L)
464
465#define IFPF_IS_ZERO 0
466#define IFPF_IS_SUBNORMAL 1
467#define IFPF_IS_NORMAL 2
468#define IFPF_IS_INFINITY 3
469#define IFPF_IS_NAN 4
470#define IFPF_TYPE_COUNT 5
471
472#define ZERO_ZERO                       0
473#define ZERO_SUBNORMAL                  1
474#define ZERO_NORMAL                     2
475#define ZERO_INFINITY                   3
476#define ZERO_NAN                        4
477
478#define SUBNORMAL_ZERO                  5
479#define SUBNORMAL_SUBNORMAL             6
480#define SUBNORMAL_NORMAL                7
481#define SUBNORMAL_INFINITY              8
482#define SUBNORMAL_NAN                   9
483
484#define NORMAL_ZERO                     10
485#define NORMAL_SUBNORMAL                11
486#define NORMAL_NORMAL                   12
487#define NORMAL_INFINITY                 13
488#define NORMAL_NAN                      14
489
490#define INFINITY_ZERO                   15
491#define INFINITY_SUBNORMAL              16
492#define INFINITY_NORMAL                 17
493#define INFINITY_INFINITY               18
494#define INFINITY_NAN                    19
495
496#define NAN_ZERO                        20
497#define NAN_SUBNORMAL                   21
498#define NAN_NORMAL                      22
499#define NAN_INFINITY                    23
500#define NAN_NAN                         24
501#define OPERAND_ZERO                    0
502#define OPERAND_SUBNORMAL               1
503#define OPERAND_NORMAL                  2
504#define OPERAND_INFINITY                3
505#define OPERAND_NAN                     4
506
507typedef struct
508{
509        u8 type;        /* Indicates, NORMAL, SUBNORMAL, etc. */
510        u8 sign;        /* Mantissa sign */
511        short exp;      /* Signed exponent...no bias */
512        u16 mantissa[INTERNAL_FPF_PRECISION];
513} InternalFPF;
514
515static
516void SetupCPUEmFloatArrays(InternalFPF *abase,
517        InternalFPF *bbase, InternalFPF *cbase, ulong arraysize);
518static
519ulong DoEmFloatIteration(InternalFPF *abase,
520        InternalFPF *bbase, InternalFPF *cbase,
521        ulong arraysize, ulong loops);
522
523static void SetInternalFPFZero(InternalFPF *dest,
524                        uchar sign);
525static void SetInternalFPFInfinity(InternalFPF *dest,
526                        uchar sign);
527static void SetInternalFPFNaN(InternalFPF *dest);
528static int IsMantissaZero(u16 *mant);
529static void Add16Bits(u16 *carry,u16 *a,u16 b,u16 c);
530static void Sub16Bits(u16 *borrow,u16 *a,u16 b,u16 c);
531static void ShiftMantLeft1(u16 *carry,u16 *mantissa);
532static void ShiftMantRight1(u16 *carry,u16 *mantissa);
533static void StickyShiftRightMant(InternalFPF *ptr,int amount);
534static void normalize(InternalFPF *ptr);
535static void denormalize(InternalFPF *ptr,int minimum_exponent);
536static void RoundInternalFPF(InternalFPF *ptr);
537static void choose_nan(InternalFPF *x,InternalFPF *y,InternalFPF *z,
538                int intel_flag);
539static void AddSubInternalFPF(uchar operation,InternalFPF *x,
540                InternalFPF *y,InternalFPF *z);
541static void MultiplyInternalFPF(InternalFPF *x,InternalFPF *y,
542                        InternalFPF *z);
543static void DivideInternalFPF(InternalFPF *x,InternalFPF *y,
544                        InternalFPF *z);
545
546static void Int32ToInternalFPF(int32 mylong,
547                InternalFPF *dest);
548static int InternalFPFToString(char *dest,
549                InternalFPF *src);
550
551static int32 randnum(int32 lngval);
552
553static int32 randwc(int32 num)
554{
555	return(randnum((int32)0)%num);
556}
557
558static int32 randw[2] = { (int32)13 , (int32)117 };
559static int32 randnum(int32 lngval)
560{
561	register int32 interm;
562
563	if (lngval!=(int32)0)
564	{	randw[0]=(int32)13; randw[1]=(int32)117; }
565
566	interm=(randw[0]*(int32)254754+randw[1]*(int32)529562)%(int32)999563;
567	randw[1]=randw[0];
568	randw[0]=interm;
569	return(interm);
570}
571
572
573static
574void SetupCPUEmFloatArrays(InternalFPF *abase,
575                InternalFPF *bbase,
576                InternalFPF *cbase,
577                ulong arraysize)
578{
579ulong i;
580InternalFPF locFPF1,locFPF2;
581
582randnum((int32)13);
583
584for(i=0;i<arraysize;i++)
585{/*       LongToInternalFPF(randwc(50000L),&locFPF1); */
586        Int32ToInternalFPF(randwc((int32)50000),&locFPF1);
587 /*       LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */
588        Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2);
589        DivideInternalFPF(&locFPF1,&locFPF2,abase+i);
590 /*       LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */
591        Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2);
592        DivideInternalFPF(&locFPF1,&locFPF2,bbase+i);
593}
594return;
595}
596
597
598static char* str1 = "loops %d\n";
599static
600ulong DoEmFloatIteration(InternalFPF *abase,
601                InternalFPF *bbase,
602                InternalFPF *cbase,
603                ulong arraysize, ulong loops)
604{
605static uchar jtable[16] = {0,0,0,0,1,1,1,1,2,2,2,2,2,3,3,3};
606ulong i;
607int number_of_loops;
608 loops = 100;
609number_of_loops=loops-1; /* the index of the first loop we run */
610
611vexxx_printf(str1, (int)loops);
612
613/*
614** Each pass through the array performs operations in
615** the followingratios:
616**   4 adds, 4 subtracts, 5 multiplies, 3 divides
617** (adds and subtracts being nearly the same operation)
618*/
619
620{
621        for(i=0;i<arraysize;i++)
622                switch(jtable[i % 16])
623                {
624                        case 0: /* Add */
625                                AddSubInternalFPF(0,abase+i,
626                                  bbase+i,
627                                  cbase+i);
628                                break;
629                        case 1: /* Subtract */
630                                AddSubInternalFPF(1,abase+i,
631                                  bbase+i,
632                                  cbase+i);
633                                break;
634                        case 2: /* Multiply */
635                                MultiplyInternalFPF(abase+i,
636                                  bbase+i,
637                                  cbase+i);
638                                break;
639                        case 3: /* Divide */
640                                DivideInternalFPF(abase+i,
641                                  bbase+i,
642                                  cbase+i);
643                                break;
644                }
645{
646  ulong j[8];   /* we test 8 entries */
647  int k;
648  ulong i;
649  char buffer[1024];
650  if (100==loops) /* the first loop */
651    {
652      j[0]=(ulong)2;
653      j[1]=(ulong)6;
654      j[2]=(ulong)10;
655      j[3]=(ulong)14;
656      j[4]=(ulong)(arraysize-14);
657      j[5]=(ulong)(arraysize-10);
658      j[6]=(ulong)(arraysize-6);
659      j[7]=(ulong)(arraysize-2);
660      for(k=0;k<8;k++){
661	i=j[k];
662	InternalFPFToString(buffer,abase+i);
663	vexxx_printf("%6d: (%s) ",i,buffer);
664	switch(jtable[i % 16])
665	  {
666	  case 0: my_strcpy(buffer,"+"); break;
667	  case 1: my_strcpy(buffer,"-"); break;
668	  case 2: my_strcpy(buffer,"*"); break;
669	  case 3: my_strcpy(buffer,"/"); break;
670	  }
671	vexxx_printf("%s ",buffer);
672	InternalFPFToString(buffer,bbase+i);
673	vexxx_printf("(%s) = ",buffer);
674	InternalFPFToString(buffer,cbase+i);
675	vexxx_printf("%s\n",buffer);
676      }
677return 0;
678    }
679}
680}
681return 0;
682}
683
684/***********************
685** SetInternalFPFZero **
686************************
687** Set an internal floating-point-format number to zero.
688** sign determines the sign of the zero.
689*/
690static void SetInternalFPFZero(InternalFPF *dest,
691                        uchar sign)
692{
693int i;          /* Index */
694
695dest->type=IFPF_IS_ZERO;
696dest->sign=sign;
697dest->exp=MIN_EXP;
698for(i=0;i<INTERNAL_FPF_PRECISION;i++)
699        dest->mantissa[i]=0;
700return;
701}
702
703/***************************
704** SetInternalFPFInfinity **
705****************************
706** Set an internal floating-point-format number to infinity.
707** This can happen if the exponent exceeds MAX_EXP.
708** As above, sign picks the sign of infinity.
709*/
710static void SetInternalFPFInfinity(InternalFPF *dest,
711                        uchar sign)
712{
713int i;          /* Index */
714
715dest->type=IFPF_IS_INFINITY;
716dest->sign=sign;
717dest->exp=MIN_EXP;
718for(i=0;i<INTERNAL_FPF_PRECISION;i++)
719        dest->mantissa[i]=0;
720return;
721}
722
723/**********************
724** SetInternalFPFNaN **
725***********************
726** Set an internal floating-point-format number to Nan
727** (not a number).  Note that we "emulate" an 80x87 as far
728** as the mantissa bits go.
729*/
730static void SetInternalFPFNaN(InternalFPF *dest)
731{
732int i;          /* Index */
733
734dest->type=IFPF_IS_NAN;
735dest->exp=MAX_EXP;
736dest->sign=1;
737dest->mantissa[0]=0x4000;
738for(i=1;i<INTERNAL_FPF_PRECISION;i++)
739        dest->mantissa[i]=0;
740
741return;
742}
743
744/*******************
745** IsMantissaZero **
746********************
747** Pass this routine a pointer to an internal floating point format
748** number's mantissa.  It checks for an all-zero mantissa.
749** Returns 0 if it is NOT all zeros, !=0 otherwise.
750*/
751static int IsMantissaZero(u16 *mant)
752{
753int i;          /* Index */
754int n;          /* Return value */
755
756n=0;
757for(i=0;i<INTERNAL_FPF_PRECISION;i++)
758        n|=mant[i];
759
760return(!n);
761}
762
763/**************
764** Add16Bits **
765***************
766** Add b, c, and carry.  Retult in a.  New carry in carry.
767*/
768static void Add16Bits(u16 *carry,
769                u16 *a,
770                u16 b,
771                u16 c)
772{
773u32 accum;              /* Accumulator */
774
775/*
776** Do the work in the 32-bit accumulator so we can return
777** the carry.
778*/
779accum=(u32)b;
780accum+=(u32)c;
781accum+=(u32)*carry;
782*carry=(u16)((accum & 0x00010000) ? 1 : 0);     /* New carry */
783*a=(u16)(accum & 0xFFFF);       /* Result is lo 16 bits */
784return;
785}
786
787/**************
788** Sub16Bits **
789***************
790** Additive inverse of above.
791*/
792static void Sub16Bits(u16 *borrow,
793                u16 *a,
794                u16 b,
795                u16 c)
796{
797u32 accum;              /* Accumulator */
798
799accum=(u32)b;
800accum-=(u32)c;
801accum-=(u32)*borrow;
802*borrow=(u32)((accum & 0x00010000) ? 1 : 0);    /* New borrow */
803*a=(u16)(accum & 0xFFFF);
804return;
805}
806
807/*******************
808** ShiftMantLeft1 **
809********************
810** Shift a vector of 16-bit numbers left 1 bit.  Also provides
811** a carry bit, which is shifted in at the beginning, and
812** shifted out at the end.
813*/
814static void ShiftMantLeft1(u16 *carry,
815                        u16 *mantissa)
816{
817int i;          /* Index */
818int new_carry;
819u16 accum;      /* Temporary holding placed */
820
821for(i=INTERNAL_FPF_PRECISION-1;i>=0;i--)
822{       accum=mantissa[i];
823        new_carry=accum & 0x8000;       /* Get new carry */
824        accum=accum<<1;                 /* Do the shift */
825        if(*carry)
826                accum|=1;               /* Insert previous carry */
827        *carry=new_carry;
828        mantissa[i]=accum;              /* Return shifted value */
829}
830return;
831}
832
833/********************
834** ShiftMantRight1 **
835*********************
836** Shift a mantissa right by 1 bit.  Provides carry, as
837** above
838*/
839static void ShiftMantRight1(u16 *carry,
840                        u16 *mantissa)
841{
842int i;          /* Index */
843int new_carry;
844u16 accum;
845
846for(i=0;i<INTERNAL_FPF_PRECISION;i++)
847{       accum=mantissa[i];
848        new_carry=accum & 1;            /* Get new carry */
849        accum=accum>>1;
850        if(*carry)
851                accum|=0x8000;
852        *carry=new_carry;
853        mantissa[i]=accum;
854}
855return;
856}
857
858
859/*****************************
860** StickyShiftMantRight **
861******************************
862** This is a shift right of the mantissa with a "sticky bit".
863** I.E., if a carry of 1 is shifted out of the least significant
864** bit, the least significant bit is set to 1.
865*/
866static void StickyShiftRightMant(InternalFPF *ptr,
867                        int amount)
868{
869int i;          /* Index */
870u16 carry;      /* Self-explanatory */
871u16 *mantissa;
872
873mantissa=ptr->mantissa;
874
875if(ptr->type!=IFPF_IS_ZERO)     /* Don't bother shifting a zero */
876{
877        /*
878        ** If the amount of shifting will shift everyting
879        ** out of existence, then just clear the whole mantissa
880        ** and set the lowmost bit to 1.
881        */
882        if(amount>=INTERNAL_FPF_PRECISION * 16)
883        {
884                for(i=0;i<INTERNAL_FPF_PRECISION-1;i++)
885                        mantissa[i]=0;
886                mantissa[INTERNAL_FPF_PRECISION-1]=1;
887        }
888        else
889                for(i=0;i<amount;i++)
890                {
891                        carry=0;
892                        ShiftMantRight1(&carry,mantissa);
893                        if(carry)
894                                mantissa[INTERNAL_FPF_PRECISION-1] |= 1;
895                }
896}
897return;
898}
899
900
901/**************************************************
902**         POST ARITHMETIC PROCESSING            **
903**  (NORMALIZE, ROUND, OVERFLOW, AND UNDERFLOW)  **
904**************************************************/
905
906/**************
907** normalize **
908***************
909** Normalize an internal-representation number.  Normalization
910** discards empty most-significant bits.
911*/
912static void normalize(InternalFPF *ptr)
913{
914u16     carry;
915
916/*
917** As long as there's a highmost 0 bit, shift the significand
918** left 1 bit.  Each time you do this, though, you've
919** gotta decrement the exponent.
920*/
921while ((ptr->mantissa[0] & 0x8000) == 0)
922{
923        carry = 0;
924        ShiftMantLeft1(&carry, ptr->mantissa);
925        ptr->exp--;
926}
927return;
928}
929
930/****************
931** denormalize **
932*****************
933** Denormalize an internal-representation number.  This means
934** shifting it right until its exponent is equivalent to
935** minimum_exponent. (You have to do this often in order
936** to perform additions and subtractions).
937*/
938static void denormalize(InternalFPF *ptr,
939                int minimum_exponent)
940{
941long exponent_difference;
942
943if (IsMantissaZero(ptr->mantissa))
944{
945        vexxx_printf("Error:  zero significand in denormalize\n");
946}
947
948exponent_difference = ptr->exp-minimum_exponent;
949if (exponent_difference < 0)
950{
951        /*
952        ** The number is subnormal
953        */
954        exponent_difference = -exponent_difference;
955        if (exponent_difference >= (INTERNAL_FPF_PRECISION * 16))
956        {
957                /* Underflow */
958                SetInternalFPFZero(ptr, ptr->sign);
959        }
960        else
961        {
962                ptr->exp+=exponent_difference;
963                StickyShiftRightMant(ptr, exponent_difference);
964        }
965}
966return;
967}
968
969
970/*********************
971** RoundInternalFPF **
972**********************
973** Round an internal-representation number.
974** The kind of rounding we do here is simplest...referred to as
975** "chop".  "Extraneous" rightmost bits are simply hacked off.
976*/
977void RoundInternalFPF(InternalFPF *ptr)
978{
979/* int i; */
980
981if (ptr->type == IFPF_IS_NORMAL ||
982        ptr->type == IFPF_IS_SUBNORMAL)
983{
984        denormalize(ptr, MIN_EXP);
985        if (ptr->type != IFPF_IS_ZERO)
986        {
987
988                /* clear the extraneous bits */
989                ptr->mantissa[3] &= 0xfff8;
990/*              for (i=4; i<INTERNAL_FPF_PRECISION; i++)
991                {
992                        ptr->mantissa[i] = 0;
993                }
994*/
995                /*
996                ** Check for overflow
997                */
998/*              Does not do anything as ptr->exp is a short and MAX_EXP=37268
999		if (ptr->exp > MAX_EXP)
1000                {
1001                        SetInternalFPFInfinity(ptr, ptr->sign);
1002                }
1003*/
1004        }
1005}
1006return;
1007}
1008
1009/*******************************************************
1010**  ARITHMETIC OPERATIONS ON INTERNAL REPRESENTATION  **
1011*******************************************************/
1012
1013/***************
1014** choose_nan **
1015****************
1016** Called by routines that are forced to perform math on
1017** a pair of NaN's.  This routine "selects" which NaN is
1018** to be returned.
1019*/
1020static void choose_nan(InternalFPF *x,
1021                InternalFPF *y,
1022                InternalFPF *z,
1023                int intel_flag)
1024{
1025int i;
1026
1027/*
1028** Compare the two mantissas,
1029** return the larger.  Note that we will be emulating
1030** an 80387 in this operation.
1031*/
1032for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1033{
1034        if (x->mantissa[i] > y->mantissa[i])
1035        {
1036                my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1037                return;
1038        }
1039        if (x->mantissa[i] < y->mantissa[i])
1040        {
1041                my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1042                return;
1043        }
1044}
1045
1046/*
1047** They are equal
1048*/
1049if (!intel_flag)
1050        /* if the operation is addition */
1051        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1052else
1053        /* if the operation is multiplication */
1054        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1055return;
1056}
1057
1058
1059/**********************
1060** AddSubInternalFPF **
1061***********************
1062** Adding or subtracting internal-representation numbers.
1063** Internal-representation numbers pointed to by x and y are
1064** added/subtracted and the result returned in z.
1065*/
1066static void AddSubInternalFPF(uchar operation,
1067                InternalFPF *x,
1068                InternalFPF *y,
1069                InternalFPF *z)
1070{
1071int exponent_difference;
1072u16 borrow;
1073u16 carry;
1074int i;
1075InternalFPF locx,locy;  /* Needed since we alter them */
1076
1077/*
1078** Following big switch statement handles the
1079** various combinations of operand types.
1080*/
1081switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1082{
1083case ZERO_ZERO:
1084        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1085        if (x->sign ^ y->sign ^ operation)
1086        {
1087                z->sign = 0; /* positive */
1088        }
1089        break;
1090
1091case NAN_ZERO:
1092case NAN_SUBNORMAL:
1093case NAN_NORMAL:
1094case NAN_INFINITY:
1095case SUBNORMAL_ZERO:
1096case NORMAL_ZERO:
1097case INFINITY_ZERO:
1098case INFINITY_SUBNORMAL:
1099case INFINITY_NORMAL:
1100        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1101        break;
1102
1103
1104case ZERO_NAN:
1105case SUBNORMAL_NAN:
1106case NORMAL_NAN:
1107case INFINITY_NAN:
1108        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1109        break;
1110
1111case ZERO_SUBNORMAL:
1112case ZERO_NORMAL:
1113case ZERO_INFINITY:
1114case SUBNORMAL_INFINITY:
1115case NORMAL_INFINITY:
1116        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1117        z->sign ^= operation;
1118        break;
1119
1120case SUBNORMAL_SUBNORMAL:
1121case SUBNORMAL_NORMAL:
1122case NORMAL_SUBNORMAL:
1123case NORMAL_NORMAL:
1124        /*
1125        ** Copy x and y to locals, since we may have
1126        ** to alter them.
1127        */
1128        my_memmove((void *)&locx,(void *)x,sizeof(InternalFPF));
1129        my_memmove((void *)&locy,(void *)y,sizeof(InternalFPF));
1130
1131        /* compute sum/difference */
1132        exponent_difference = locx.exp-locy.exp;
1133        if (exponent_difference == 0)
1134        {
1135                /*
1136                ** locx.exp == locy.exp
1137                ** so, no shifting required
1138                */
1139                if (locx.type == IFPF_IS_SUBNORMAL ||
1140                  locy.type == IFPF_IS_SUBNORMAL)
1141                        z->type = IFPF_IS_SUBNORMAL;
1142                else
1143                        z->type = IFPF_IS_NORMAL;
1144
1145                /*
1146                ** Assume that locx.mantissa > locy.mantissa
1147                */
1148                z->sign = locx.sign;
1149                z->exp= locx.exp;
1150        }
1151        else
1152                if (exponent_difference > 0)
1153                {
1154                        /*
1155                        ** locx.exp > locy.exp
1156                        */
1157                        StickyShiftRightMant(&locy,
1158                                 exponent_difference);
1159                        z->type = locx.type;
1160                        z->sign = locx.sign;
1161                        z->exp = locx.exp;
1162                }
1163                else    /* if (exponent_difference < 0) */
1164                {
1165                        /*
1166                        ** locx.exp < locy.exp
1167                        */
1168                        StickyShiftRightMant(&locx,
1169                                -exponent_difference);
1170                        z->type = locy.type;
1171                        z->sign = locy.sign ^ operation;
1172                        z->exp = locy.exp;
1173                }
1174
1175                if (locx.sign ^ locy.sign ^ operation)
1176                {
1177                        /*
1178                        ** Signs are different, subtract mantissas
1179                        */
1180                        borrow = 0;
1181                        for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1182                                Sub16Bits(&borrow,
1183                                        &z->mantissa[i],
1184                                        locx.mantissa[i],
1185                                        locy.mantissa[i]);
1186
1187                        if (borrow)
1188                        {
1189                                /* The y->mantissa was larger than the
1190                                ** x->mantissa leaving a negative
1191                                ** result.  Change the result back to
1192                                ** an unsigned number and flip the
1193                                ** sign flag.
1194                                */
1195                                z->sign = locy.sign ^ operation;
1196                                borrow = 0;
1197                                for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1198                                {
1199                                        Sub16Bits(&borrow,
1200                                                &z->mantissa[i],
1201                                                0,
1202                                                z->mantissa[i]);
1203                                }
1204                        }
1205                        else
1206                        {
1207                                /* The assumption made above
1208                                ** (i.e. x->mantissa >= y->mantissa)
1209                                ** was correct.  Therefore, do nothing.
1210                                ** z->sign = x->sign;
1211                                */
1212                        }
1213
1214                        if (IsMantissaZero(z->mantissa))
1215                        {
1216                                z->type = IFPF_IS_ZERO;
1217                                z->sign = 0; /* positive */
1218                        }
1219                        else
1220                                if (locx.type == IFPF_IS_NORMAL ||
1221                                         locy.type == IFPF_IS_NORMAL)
1222                                {
1223                                        normalize(z);
1224                                }
1225                }
1226                else
1227                {
1228                        /* signs are the same, add mantissas */
1229                        carry = 0;
1230                        for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1231                        {
1232                                Add16Bits(&carry,
1233                                        &z->mantissa[i],
1234                                        locx.mantissa[i],
1235                                        locy.mantissa[i]);
1236                        }
1237
1238                        if (carry)
1239                        {
1240                                z->exp++;
1241                                carry=0;
1242                                ShiftMantRight1(&carry,z->mantissa);
1243                                z->mantissa[0] |= 0x8000;
1244                                z->type = IFPF_IS_NORMAL;
1245                        }
1246                        else
1247                                if (z->mantissa[0] & 0x8000)
1248                                        z->type = IFPF_IS_NORMAL;
1249        }
1250        break;
1251
1252case INFINITY_INFINITY:
1253        SetInternalFPFNaN(z);
1254        break;
1255
1256case NAN_NAN:
1257        choose_nan(x, y, z, 1);
1258        break;
1259}
1260
1261/*
1262** All the math is done; time to round.
1263*/
1264RoundInternalFPF(z);
1265return;
1266}
1267
1268
1269/************************
1270** MultiplyInternalFPF **
1271*************************
1272** Two internal-representation numbers x and y are multiplied; the
1273** result is returned in z.
1274*/
1275static void MultiplyInternalFPF(InternalFPF *x,
1276                        InternalFPF *y,
1277                        InternalFPF *z)
1278{
1279int i;
1280int j;
1281u16 carry;
1282u16 extra_bits[INTERNAL_FPF_PRECISION];
1283InternalFPF locy;       /* Needed since this will be altered */
1284/*
1285** As in the preceding function, this large switch
1286** statement selects among the many combinations
1287** of operands.
1288*/
1289switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1290{
1291case INFINITY_SUBNORMAL:
1292case INFINITY_NORMAL:
1293case INFINITY_INFINITY:
1294case ZERO_ZERO:
1295case ZERO_SUBNORMAL:
1296case ZERO_NORMAL:
1297        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1298        z->sign ^= y->sign;
1299        break;
1300
1301case SUBNORMAL_INFINITY:
1302case NORMAL_INFINITY:
1303case SUBNORMAL_ZERO:
1304case NORMAL_ZERO:
1305        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1306        z->sign ^= x->sign;
1307        break;
1308
1309case ZERO_INFINITY:
1310case INFINITY_ZERO:
1311        SetInternalFPFNaN(z);
1312        break;
1313
1314case NAN_ZERO:
1315case NAN_SUBNORMAL:
1316case NAN_NORMAL:
1317case NAN_INFINITY:
1318        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1319        break;
1320
1321case ZERO_NAN:
1322case SUBNORMAL_NAN:
1323case NORMAL_NAN:
1324case INFINITY_NAN:
1325        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1326        break;
1327
1328
1329case SUBNORMAL_SUBNORMAL:
1330case SUBNORMAL_NORMAL:
1331case NORMAL_SUBNORMAL:
1332case NORMAL_NORMAL:
1333        /*
1334        ** Make a local copy of the y number, since we will be
1335        ** altering it in the process of multiplying.
1336        */
1337        my_memmove((void *)&locy,(void *)y,sizeof(InternalFPF));
1338
1339        /*
1340        ** Check for unnormal zero arguments
1341        */
1342        if (IsMantissaZero(x->mantissa) || IsMantissaZero(y->mantissa))
1343                SetInternalFPFInfinity(z, 0);
1344
1345        /*
1346        ** Initialize the result
1347        */
1348        if (x->type == IFPF_IS_SUBNORMAL ||
1349            y->type == IFPF_IS_SUBNORMAL)
1350                z->type = IFPF_IS_SUBNORMAL;
1351        else
1352                z->type = IFPF_IS_NORMAL;
1353
1354        z->sign = x->sign ^ y->sign;
1355        z->exp = x->exp + y->exp ;
1356        for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1357        {
1358                z->mantissa[i] = 0;
1359                extra_bits[i] = 0;
1360        }
1361
1362        for (i=0; i<(INTERNAL_FPF_PRECISION*16); i++)
1363        {
1364                /*
1365                ** Get rightmost bit of the multiplier
1366                */
1367                carry = 0;
1368                ShiftMantRight1(&carry, locy.mantissa);
1369                if (carry)
1370                {
1371                        /*
1372                        ** Add the multiplicand to the product
1373                        */
1374                        carry = 0;
1375                        for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--)
1376                                Add16Bits(&carry,
1377                                        &z->mantissa[j],
1378                                        z->mantissa[j],
1379                                        x->mantissa[j]);
1380                }
1381                else
1382                {
1383                        carry = 0;
1384                }
1385
1386                /*
1387                ** Shift the product right.  Overflow bits get
1388                ** shifted into extra_bits.  We'll use it later
1389                ** to help with the "sticky" bit.
1390                */
1391                ShiftMantRight1(&carry, z->mantissa);
1392                ShiftMantRight1(&carry, extra_bits);
1393        }
1394
1395        /*
1396        ** Normalize
1397        ** Note that we use a "special" normalization routine
1398        ** because we need to use the extra bits. (These are
1399        ** bits that may have been shifted off the bottom that
1400        ** we want to reclaim...if we can.
1401        */
1402        while ((z->mantissa[0] & 0x8000) == 0)
1403        {
1404                carry = 0;
1405                ShiftMantLeft1(&carry, extra_bits);
1406                ShiftMantLeft1(&carry, z->mantissa);
1407                z->exp--;
1408        }
1409
1410        /*
1411        ** Set the sticky bit if any bits set in extra bits.
1412        */
1413        if (IsMantissaZero(extra_bits))
1414        {
1415                z->mantissa[INTERNAL_FPF_PRECISION-1] |= 1;
1416        }
1417        break;
1418
1419case NAN_NAN:
1420        choose_nan(x, y, z, 0);
1421        break;
1422}
1423
1424/*
1425** All math done...do rounding.
1426*/
1427RoundInternalFPF(z);
1428return;
1429}
1430
1431
1432/**********************
1433** DivideInternalFPF **
1434***********************
1435** Divide internal FPF number x by y.  Return result in z.
1436*/
1437static void DivideInternalFPF(InternalFPF *x,
1438                        InternalFPF *y,
1439                        InternalFPF *z)
1440{
1441int i;
1442int j;
1443u16 carry;
1444u16 extra_bits[INTERNAL_FPF_PRECISION];
1445InternalFPF locx;       /* Local for x number */
1446
1447/*
1448** As with preceding function, the following switch
1449** statement selects among the various possible
1450** operands.
1451*/
1452switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1453{
1454case ZERO_ZERO:
1455case INFINITY_INFINITY:
1456        SetInternalFPFNaN(z);
1457        break;
1458
1459case ZERO_SUBNORMAL:
1460case ZERO_NORMAL:
1461        if (IsMantissaZero(y->mantissa))
1462        {
1463                SetInternalFPFNaN(z);
1464                break;
1465        }
1466
1467case ZERO_INFINITY:
1468case SUBNORMAL_INFINITY:
1469case NORMAL_INFINITY:
1470        SetInternalFPFZero(z, x->sign ^ y->sign);
1471        break;
1472
1473case SUBNORMAL_ZERO:
1474case NORMAL_ZERO:
1475        if (IsMantissaZero(x->mantissa))
1476        {
1477                SetInternalFPFNaN(z);
1478                break;
1479        }
1480
1481case INFINITY_ZERO:
1482case INFINITY_SUBNORMAL:
1483case INFINITY_NORMAL:
1484        SetInternalFPFInfinity(z, 0);
1485        z->sign = x->sign ^ y->sign;
1486        break;
1487
1488case NAN_ZERO:
1489case NAN_SUBNORMAL:
1490case NAN_NORMAL:
1491case NAN_INFINITY:
1492        my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1493        break;
1494
1495case ZERO_NAN:
1496case SUBNORMAL_NAN:
1497case NORMAL_NAN:
1498case INFINITY_NAN:
1499        my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1500        break;
1501
1502case SUBNORMAL_SUBNORMAL:
1503case NORMAL_SUBNORMAL:
1504case SUBNORMAL_NORMAL:
1505case NORMAL_NORMAL:
1506        /*
1507        ** Make local copy of x number, since we'll be
1508        ** altering it in the process of dividing.
1509        */
1510        my_memmove((void *)&locx,(void *)x,sizeof(InternalFPF));
1511
1512        /*
1513        ** Check for unnormal zero arguments
1514        */
1515        if (IsMantissaZero(locx.mantissa))
1516        {
1517                if (IsMantissaZero(y->mantissa))
1518                        SetInternalFPFNaN(z);
1519                else
1520                        SetInternalFPFZero(z, 0);
1521                break;
1522        }
1523        if (IsMantissaZero(y->mantissa))
1524        {
1525                SetInternalFPFInfinity(z, 0);
1526                break;
1527        }
1528
1529        /*
1530        ** Initialize the result
1531        */
1532        z->type = x->type;
1533        z->sign = x->sign ^ y->sign;
1534        z->exp = x->exp - y->exp +
1535                        ((INTERNAL_FPF_PRECISION * 16 * 2));
1536        for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1537        {
1538                z->mantissa[i] = 0;
1539                extra_bits[i] = 0;
1540        }
1541
1542        while ((z->mantissa[0] & 0x8000) == 0)
1543        {
1544                carry = 0;
1545                ShiftMantLeft1(&carry, locx.mantissa);
1546                ShiftMantLeft1(&carry, extra_bits);
1547
1548                /*
1549                ** Time to subtract yet?
1550                */
1551                if (carry == 0)
1552                        for (j=0; j<INTERNAL_FPF_PRECISION; j++)
1553                        {
1554                                if (y->mantissa[j] > extra_bits[j])
1555                                {
1556                                        carry = 0;
1557                                        goto no_subtract;
1558                                }
1559                                if (y->mantissa[j] < extra_bits[j])
1560                                        break;
1561                        }
1562                /*
1563                ** Divisor (y) <= dividend (x), subtract
1564                */
1565                carry = 0;
1566                for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--)
1567                        Sub16Bits(&carry,
1568                                &extra_bits[j],
1569                                extra_bits[j],
1570                                y->mantissa[j]);
1571                carry = 1;      /* 1 shifted into quotient */
1572        no_subtract:
1573                ShiftMantLeft1(&carry, z->mantissa);
1574                z->exp--;
1575        }
1576        break;
1577
1578case NAN_NAN:
1579        choose_nan(x, y, z, 0);
1580        break;
1581}
1582
1583/*
1584** Math complete...do rounding
1585*/
1586RoundInternalFPF(z);
1587}
1588
1589/**********************
1590** LongToInternalFPF **
1591** Int32ToInternalFPF **
1592***********************
1593** Convert a signed (long) 32-bit integer into an internal FPF number.
1594*/
1595/* static void LongToInternalFPF(long mylong, */
1596static void Int32ToInternalFPF(int32 mylong,
1597                InternalFPF *dest)
1598{
1599int i;          /* Index */
1600u16 myword;     /* Used to hold converted stuff */
1601/*
1602** Save the sign and get the absolute value.  This will help us
1603** with 64-bit machines, since we use only the lower 32
1604** bits just in case. (No longer necessary after we use int32.)
1605*/
1606/* if(mylong<0L) */
1607if(mylong<(int32)0)
1608{       dest->sign=1;
1609        mylong=(int32)0-mylong;
1610}
1611else
1612        dest->sign=0;
1613/*
1614** Prepare the destination floating point number
1615*/
1616dest->type=IFPF_IS_NORMAL;
1617for(i=0;i<INTERNAL_FPF_PRECISION;i++)
1618        dest->mantissa[i]=0;
1619
1620/*
1621** See if we've got a zero.  If so, make the resultant FP
1622** number a true zero and go home.
1623*/
1624if(mylong==0)
1625{       dest->type=IFPF_IS_ZERO;
1626        dest->exp=0;
1627        return;
1628}
1629
1630/*
1631** Not a true zero.  Set the exponent to 32 (internal FPFs have
1632** no bias) and load the low and high words into their proper
1633** locations in the mantissa.  Then normalize.  The action of
1634** normalizing slides the mantissa bits into place and sets
1635** up the exponent properly.
1636*/
1637dest->exp=32;
1638myword=(u16)((mylong >> 16) & 0xFFFFL);
1639dest->mantissa[0]=myword;
1640myword=(u16)(mylong & 0xFFFFL);
1641dest->mantissa[1]=myword;
1642normalize(dest);
1643return;
1644}
1645
1646#if 1
1647/************************
1648** InternalFPFToString **
1649*************************
1650** FOR DEBUG PURPOSES
1651** This routine converts an internal floating point representation
1652** number to a string.  Used in debugging the package.
1653** Returns length of converted number.
1654** NOTE: dest must point to a buffer big enough to hold the
1655**  result.  Also, this routine does append a null (an effect
1656**  of using the sprintf() function).  It also returns
1657**  a length count.
1658** NOTE: This routine returns 5 significant digits.  Thats
1659**  about all I feel safe with, given the method of
1660**  conversion.  It should be more than enough for programmers
1661**  to determine whether the package is properly ported.
1662*/
1663static int InternalFPFToString(char *dest,
1664                InternalFPF *src)
1665{
1666InternalFPF locFPFNum;          /* Local for src (will be altered) */
1667InternalFPF IFPF10;             /* Floating-point 10 */
1668InternalFPF IFPFComp;           /* For doing comparisons */
1669int msign;                      /* Holding for mantissa sign */
1670int expcount;                   /* Exponent counter */
1671int ccount;                     /* Character counter */
1672int i,j,k;                      /* Index */
1673u16 carryaccum;                 /* Carry accumulator */
1674u16 mycarry;                    /* Local for carry */
1675
1676/*
1677** Check first for the simple things...Nan, Infinity, Zero.
1678** If found, copy the proper string in and go home.
1679*/
1680switch(src->type)
1681{
1682        case IFPF_IS_NAN:
1683                my_memcpy(dest,"NaN",3);
1684                return(3);
1685
1686        case IFPF_IS_INFINITY:
1687                if(src->sign==0)
1688                        my_memcpy(dest,"+Inf",4);
1689                else
1690                        my_memcpy(dest,"-Inf",4);
1691                return(4);
1692
1693        case IFPF_IS_ZERO:
1694                if(src->sign==0)
1695                        my_memcpy(dest,"+0",2);
1696                else
1697                        my_memcpy(dest,"-0",2);
1698                return(2);
1699}
1700
1701/*
1702** Move the internal number into our local holding area, since
1703** we'll be altering it to print it out.
1704*/
1705my_memcpy((void *)&locFPFNum,(void *)src,sizeof(InternalFPF));
1706
1707/*
1708** Set up a floating-point 10...which we'll use a lot in a minute.
1709*/
1710/* LongToInternalFPF(10L,&IFPF10); */
1711Int32ToInternalFPF((int32)10,&IFPF10);
1712
1713/*
1714** Save the mantissa sign and make it positive.
1715*/
1716msign=src->sign;
1717
1718/* src->sign=0 */ /* bug, fixed Nov. 13, 1997 */
1719(&locFPFNum)->sign=0;
1720
1721expcount=0;             /* Init exponent counter */
1722
1723/*
1724** See if the number is less than 10.  If so, multiply
1725** the number repeatedly by 10 until it's not.   For each
1726** multiplication, decrement a counter so we can keep track
1727** of the exponent.
1728*/
1729
1730while(1)
1731{       AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp);
1732        if(IFPFComp.sign==0) break;
1733        MultiplyInternalFPF(&locFPFNum,&IFPF10,&IFPFComp);
1734        expcount--;
1735        my_memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF));
1736}
1737/*
1738** Do the reverse of the above.  As long as the number is
1739** greater than or equal to 10, divide it by 10.  Increment the
1740** exponent counter for each multiplication.
1741*/
1742
1743while(1)
1744{
1745        AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp);
1746        if(IFPFComp.sign!=0) break;
1747        DivideInternalFPF(&locFPFNum,&IFPF10,&IFPFComp);
1748        expcount++;
1749        my_memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF));
1750}
1751
1752/*
1753** About time to start storing things.  First, store the
1754** mantissa sign.
1755*/
1756ccount=1;               /* Init character counter */
1757if(msign==0)
1758        *dest++='+';
1759else
1760        *dest++='-';
1761
1762/*
1763** At this point we know that the number is in the range
1764** 10 > n >=1.  We need to "strip digits" out of the
1765** mantissa.  We do this by treating the mantissa as
1766** an integer and multiplying by 10. (Not a floating-point
1767** 10, but an integer 10.  Since this is debug code and we
1768** could care less about speed, we'll do it the stupid
1769** way and simply add the number to itself 10 times.
1770** Anything that makes it to the left of the implied binary point
1771** gets stripped off and emitted.  We'll do this for
1772** 5 significant digits (which should be enough to
1773** verify things).
1774*/
1775/*
1776** Re-position radix point
1777*/
1778carryaccum=0;
1779while(locFPFNum.exp>0)
1780{
1781        mycarry=0;
1782        ShiftMantLeft1(&mycarry,locFPFNum.mantissa);
1783        carryaccum=(carryaccum<<1);
1784        if(mycarry) carryaccum++;
1785        locFPFNum.exp--;
1786}
1787
1788while(locFPFNum.exp<0)
1789{
1790        mycarry=0;
1791        ShiftMantRight1(&mycarry,locFPFNum.mantissa);
1792        locFPFNum.exp++;
1793}
1794
1795for(i=0;i<6;i++)
1796        if(i==1)
1797        {       /* Emit decimal point */
1798                *dest++='.';
1799                ccount++;
1800        }
1801        else
1802        {       /* Emit a digit */
1803                *dest++=('0'+carryaccum);
1804                ccount++;
1805
1806                carryaccum=0;
1807                my_memcpy((void *)&IFPF10,
1808                        (void *)&locFPFNum,
1809                        sizeof(InternalFPF));
1810
1811                /* Do multiply via repeated adds */
1812                for(j=0;j<9;j++)
1813                {
1814                        mycarry=0;
1815                        for(k=(INTERNAL_FPF_PRECISION-1);k>=0;k--)
1816                                Add16Bits(&mycarry,&(IFPFComp.mantissa[k]),
1817                                        locFPFNum.mantissa[k],
1818                                        IFPF10.mantissa[k]);
1819                        carryaccum+=mycarry ? 1 : 0;
1820                        my_memcpy((void *)&locFPFNum,
1821                                (void *)&IFPFComp,
1822                                sizeof(InternalFPF));
1823                }
1824        }
1825
1826/*
1827** Now move the 'E', the exponent sign, and the exponent
1828** into the string.
1829*/
1830*dest++='E';
1831
1832/* sprint is supposed to return an integer, but it caused problems on SunOS
1833 * with the native cc. Hence we force it.
1834 * Uwe F. Mayer
1835 */
1836if (expcount < 0) {
1837     *dest++ = '-';
1838     expcount =- expcount;
1839}
1840else *dest++ = ' ';
1841
1842*dest++ = (char)(expcount + '0');
1843*dest++ = 0;
1844
1845ccount += 3;
1846/*
1847** All done, go home.
1848*/
1849return(ccount);
1850
1851}
1852
1853#endif
1854
1855
1856
1857////////////////////////////////////////////////////////////////////////
1858static
1859void* AllocateMemory ( unsigned long n, int* p )
1860{
1861  *p = 0;
1862  void* r = (void*) (*serviceFn)(2,n);
1863  return r;
1864}
1865static
1866void FreeMemory ( void* p, int* zz )
1867{
1868  *zz = 0;
1869  // free(p);
1870}
1871
1872
1873
1874/**************
1875** DoEmFloat **
1876***************
1877** Perform the floating-point emulation routines portion of the
1878** CPU benchmark.  Returns the operations per second.
1879*/
1880static
1881void DoEmFloat(void)
1882{
1883EmFloatStruct *locemfloatstruct;        /* Local structure */
1884InternalFPF *abase;             /* Base of A array */
1885InternalFPF *bbase;             /* Base of B array */
1886InternalFPF *cbase;             /* Base of C array */
1887ulong tickcount;                /* # of ticks */
1888char *errorcontext;             /* Error context string pointer */
1889int systemerror;                /* For holding error code */
1890ulong loops;                    /* # of loops */
1891
1892/*
1893** Link to global structure
1894*/
1895EmFloatStruct global_emfloatstruct;
1896 global_emfloatstruct.adjust = 0;
1897 global_emfloatstruct.request_secs = 0;
1898 global_emfloatstruct.arraysize = 100;
1899 global_emfloatstruct.loops = 1;
1900 global_emfloatstruct.emflops = 0.0;
1901locemfloatstruct=&global_emfloatstruct;
1902
1903/*
1904** Set the error context
1905*/
1906errorcontext="CPU:Floating Emulation";
1907
1908
1909abase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1910		&systemerror);
1911
1912bbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1913		&systemerror);
1914
1915cbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1916		&systemerror);
1917
1918/*
1919** Set up the arrays
1920*/
1921SetupCPUEmFloatArrays(abase,bbase,cbase,locemfloatstruct->arraysize);
1922
1923 loops=100;
1924	       tickcount=DoEmFloatIteration(abase,bbase,cbase,
1925			locemfloatstruct->arraysize,
1926			loops);
1927
1928FreeMemory((void *)abase,&systemerror);
1929FreeMemory((void *)bbase,&systemerror);
1930FreeMemory((void *)cbase,&systemerror);
1931
1932return;
1933}
1934
1935//////////////////
1936void entry ( HWord(*f)(HWord,HWord) )
1937{
1938  serviceFn = f;
1939  vexxx_printf("starting\n");
1940  DoEmFloat();
1941  (*serviceFn)(0,0);
1942}
1943