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