mandelps.c

/* -----------------------------------------------------------------------------
 *  Copyright (C) 2019-2021 daiteq s.r.o.                http://www.daiteq.com
 *
 *  This program is distributed WITHOUT ANY WARRANTY; without even
 *  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 *  PURPOSE.
 *
 * -----------------------------------------------------------------------------
 *  Filename    : mandelps.c
 *  Authors     : Martin Danek
 *  Description : Mandelbrot set with packed single-precision ops
 *  Release     :
 *  Version     : 1.5
 *  Date        : 27.4.2021
 * -----------------------------------------------------------------------------
 */

#include <stdint.h>
#include <leon-myuart.h>
#include <my_printf.h>

#include <sfmath.h>
#include <sysregs.h>
#include <lconf.h>

/* for clearing BSS section */
extern unsigned long __bss_start;  /* start address for the .bss section. defined in linker script */
extern unsigned long _end;         /* end address for the .bss section. defined in linker script */


#define NOIOTIME

#if (defined(PACKED) && (!defined(GCC)))

  typedef float cpx_t __attribute__((ext_vector_type(2)));

#else

  typedef struct cpx_t {
    float x;
    float y;
  } cpx_t;

#endif

#if (!defined(GCC))

#define FZERO 0.F
#define FONE 1.F
#define FTHREE 3.F
#define FFOUR 4.F
#define FMAX 1.e+35F

#else

#define FZERO ((float)0.)
#define FONE ((float)1.)
#define FTHREE ((float)3.)
#define FFOUR ((float)4.)
#define FMAX ((float)1.e+35)

#endif

// typedef union fp64 {
// /*
//   struct {
//     float up;
//     float dn;
//   } stf;
// */
//   cpx_t ps;
//   double d;
//   float f[2];
//   unsigned u[2];
// } fp64;


#if defined(PACKED)

#if defined(GCC)

  #define opcode_FADDPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "faddps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

  #define opcode_FSUBPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fsubps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

  #define opcode_FMULPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fmulps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \


#endif

  #define opcode_FADDRPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "faddrps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

  #define opcode_FSUBRPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fsubrps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

  #define opcode_FMULXPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fmulxps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

  #define opcode_FSUBADDRPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fsubaddrps\t%1, %2, %0 \n" \
      : "=f" (RES)          \
      : "If" (A), "If" (B)  \
    )                       \

  #define opcode_FMULRPS(A, B, RES)   \
  __asm__ __volatile__ (       \
      "fmulrps\t%1, %2, %0 \n" \
      : "=f" (RES)                                                             \
      : "If" (A), "If" (B)                                                     \
    )                       \

#endif


#if defined(PACKED)

#if defined(GCC)

#define dcadd(a,b,z) \
  opcode_FADDPS(*((double*)&a),*((double*)&b),*((double*)&z));

#define dcsub(a,b,z) \
  opcode_FSUBPS(*((double*)&a),*((double*)&b),*((double*)&z));


// 1. C = FMULP(A,B) = (ar*br, ai*bi)
// 2. D = FMULXP(A,B) = (ar*bi, ai*br)
// 3. H = FSUBADDRP(C,D) = (ar*br-ai*bi, ar*bi+ai*br)
#define dcmul(a,b,z) \
  do { \
  cpx_t tc,td; \
  opcode_FMULPS(*((double*)&a),*((double*)&b),*((double*)&tc)); \
  opcode_FMULXPS(*((double*)&a),*((double*)&b),*((double*)&td)); \
  opcode_FSUBADDRPS(*((double*)&tc),*((double*)&td),*((double*)&z)); \
  } while (0)

#define dcabss(a,z) \
  do { \
  cpx_t tc, td, th; \
  th.x=FZERO; th.y=FZERO; \
  opcode_FMULPS(*((double*)&a),*((double*)&a),*((double*)&tc)); \
  opcode_FADDRPS(*((double*)&tc),*((double*)&th),*((double*)&td)); \
  z=td.x; \
  } while (0)


#else

#define dcadd(a,b,z) \
  z=a+b;

#define dcsub(a,b,z) \
  z=a-b;


// 1. C = FMULP(A,B) = (ar*br, ai*bi)
// 2. D = FMULXP(A,B) = (ar*bi, ai*br)
// 3. H = FSUBADDRP(C,D) = (ar*br-ai*bi, ar*bi+ai*br)
#define dcmul(a,b,z) \
  do { \
  cpx_t tc,td; \
  tc=a*b; \
  opcode_FMULXPS(a,b,td); \
  opcode_FSUBADDRPS(tc,td,z); \
  } while (0)

#define dcabss(a,z) \
  do { \
  cpx_t tc, td, th; \
  th.x=FZERO; th.y=FZERO; \
  tc=a*a; \
  opcode_FADDRPS(tc,th,td); \
  z=td.x; \
  } while (0)


#endif

#else

#define dcadd(a,b,z) \
    z.x=a.x+b.x; \
    z.y=a.y+b.y;

#define dcsub(a,b,z) \
    z.x=a.x-b.x; \
    z.y=a.y-b.y;


#define dcmul(a,b,z) \
  do { \
  cpx_t th; \
  th.x=a.x*b.x-a.y*b.y; \
  th.y=a.y*b.x+a.x*b.y; \
  z=th; \
  } while (0)

#define dcabss(a,z) \
  z=a.x*a.x+a.y*a.y;


#endif


#if defined(FCALL)

cpx_t cadd(cpx_t a, cpx_t b) {
  cpx_t c;
#if defined(PACKED)
  c=a+b;
#else
  c.x=a.x+b.x;
  c.y=a.y+b.y;
#endif
  return c;
}



cpx_t csub(cpx_t a, cpx_t b) {
  cpx_t c;
#if defined(PACKED)
  c=a-b;
#else
  c.x=a.x-b.x;
  c.y=a.y-b.y;
#endif
  return c;
}



cpx_t cmul(cpx_t a, cpx_t b) {
// 1. C = FMULP(A,B) = (ar*br, ai*bi)
// 2. D = FMULXP(A,B) = (ar*bi, ai*br)
// 3. H = FSUBADDRP(C,D) = (ar*br-ai*bi, ar*bi+ai*br)
//

  cpx_t x,y,c,d,e,f,g,h,zero;
#if defined(PACKED)
//   zero.d=0.;
  x=a;
  y=b;
  c=x*y;
  opcode_FMULXPS(x,y,d);
  opcode_FSUBADDRPS(c,d,h);
#else
  h.x=a.x*b.x-a.y*b.y;
  h.y=a.y*b.x+a.x*b.y;
#endif
  return h;
}

float cabss(cpx_t a) {
  cpx_t x,c,d;
  const cpx_t zero={FZERO,FZERO};
#if defined(PACKED)
  c=a*a;
  opcode_FADDRPS(c,zero,d);
#else
  d.x=a.x*a.x+a.y*a.y;
#endif
  return d.x;
}

#endif


#define RES_X 80
#define RES_Y 80
#define MAX_ITER 100

#define MAX_COL 80
#define MAX_ROW 80

#define printpacked(A) \
  c=A.x; \
  d=A.y; \
  printf(" (%f,%f)",c,d);


int main(void) {

    float a,b,c,d;
    cpx_t pa,pb,pc,point,next;

    int i, j;
    unsigned k,l,m,hit;


  /* reset BSS section to zero */
  for(unsigned long *pMem = &__bss_start; pMem < &_end;) {
    *(pMem++) = 0;
  }

  init_uart();

  init_softfloat();

#ifndef NOHWCONF
  do {
    lconf_type lconf;
    lconf_fpu_type lconf_fpu;
    lconf_swar_type lconf_swar;
    get_lconf(&lconf);
    get_lconf_fpu(&lconf_fpu);
    get_lconf_swar(&lconf_swar);
    report_lconf(&lconf);
    if ((lconf.fpu==FPU_DAIFPU)||(lconf.fpu==FPU_MEIKO))
      report_lconf_fpu(&lconf_fpu);
    report_lconf_swar(&lconf_swar);
    sysregs_init();
  } while(0);
#endif /* NOHWCONF */

#if defined (PACKED)
    m_print("MANDELREF - SINGLE PRECISION, PACKED\n");
#else
    m_print("MANDELREF - SINGLE PRECISION, NORMAL\n");
#endif

// #define INIT_TEST
#if defined (INIT_TEST)
    a = FTHREE;
    b = FFOUR;

    pa.x = a;
    pa.y = b;
    printf("A=(%f,%f)\n",pa.x,pa.y);
    tx= pa;
    printf("A=(0x%08X, 0x%08X)\n",tx.x,tx.y);

    pb.x = FZERO;
    pb.y = FONE;
    printf("B=(%f,%f)\n",pb.x,pb.y);

#if defined(FCALL)
    pc=cadd(pa,pb);
#else
    dcadd(pa,pb,pc);
#endif

    printf("pa+pb (%f,%f)\n",pc.x,pc.y);
#if defined(FCALL)
    pc=csub(pa,pb);
#else
    dcsub(pa,pb,pc);
#endif
    printf("pa-pb (%f,%f)\n",pc.x,pc.y);
#if defined(FCALL)
    pc=cmul(pa,pb);
#else
    dcmul(pa,pb,pc);
#endif
    printf("paxpb (%f,%f)\n",pc.x,pc.y);
#endif

#ifdef NOIOTIME
    sysregs_init();
#endif

    for (j=-RES_Y;j<RES_Y;j++) {
      for (i=-2*RES_X;i<RES_X;i++) {
        point.x=((float)i)/(((float)RES_X));point.y=((float)j)/(((float)RES_X));
//         point.x=i/(RES_X*1.);point.y=j/(RES_X*1.);
        k=0;
        next.x=FZERO;
        next.y=FZERO;
#if defined(FCALL)
        while ((k<MAX_ITER)&&(cabss(next)<FMAX)) {
#else
        dcabss(next,a);
        while ((k<MAX_ITER)&&(a<FMAX)) {
#endif
#if defined(FCALL)
          next=cadd(cmul(next,next),point);
#else
          dcmul(next,next,pa);
          dcadd(pa,point,next);
          dcabss(next,a);
#endif
          k++;
        }
        if (k==MAX_ITER) {
#ifdef NOIOTIME
          sysregs_stop();
#endif
          m_putchar('*');
#ifdef NOIOTIME
          sysregs_start();
#endif
        }
        else {
          if (k<10) {
#ifdef NOIOTIME
            sysregs_stop();
#endif
            m_putchar('0'+k);
#ifdef NOIOTIME
            sysregs_start();
#endif
          }
          else {
#ifdef NOIOTIME
            sysregs_stop();
#endif
            m_putchar('A'+((k>35)?35:k)-10);
#ifdef NOIOTIME
            sysregs_start();
#endif
          }
        }
      }
      NL;
    }

    NL;

#ifndef NOHWCONF
    do {
      unsigned ticks_lo, ticks_hi, insns_lo, insns_hi;
      sysregs_stop();
      sysregs_read(&ticks_lo, &ticks_hi, &insns_lo, &insns_hi);
      if (ticks_hi)
        printf("TEST EXECUTED IN (%u * 2^32 + %u) TICKS AND (%u * 2^32 + %u) INSTRUCTIONS\n\n", ticks_hi, ticks_lo, insns_hi, insns_lo);
      else
        printf("TEST EXECUTED IN %u TICKS AND %u INSTRUCTIONS\n\n", ticks_lo, insns_lo);
    } while(0);
#endif /* NOHWCONF */


    printf("\n\nEND OF TEST\n\n");

    __asm__("ta 1; nop; nop;");


    return 0;
}