Yes, I had one - full float arithmetics without integer multiplication-division.
Floating Point Systems - I remember its computers from long ago, hosted by an IBM mainframe. I'd compile code for them with a cross-compiler on the host system. FPS's computers were Vector processor ones, with instructions that do the same thing on the contents of a sequence of memory locations. These were common in the 1970's to 1990's, but they went into decline as more conventional chips because cheap competition for them.
#include<stdarg.h>
#include <strings.h>
void print_arrays(int nrow,char *prspec, ...){
//print arrays of same length to the console
// prpec forma stringas as in printf()
// narg numer of arrays passed to function,
//picnt pdcnt counters for number of int and double array pointers
// pi pd arrays of pointers to arrays
cout<<"print arrays "<<endl;
int i,j,narg = 0,picnt = 0,pdcnt = 0;
int * pi[20];
double *pd[20];
//--------------------------------------
//split format string into substrigs for printf()
char fdel[2] = " "; //split delimiter
char pfmt[10][10];
char * token = strtok(prspec,fdel);
while( token != NULL ) {
strcpy(pfmt[narg++],token);
//printf( " %s\n", token ); //printing each token
token = strtok(NULL, fdel);
}
//-------------------------------------
// make var_args format string
int nfmt = 100;
char *fmt = new char[nfmt];
for(i=0;i<nfmt;i++)fmt[i] = '\0'; //null the array so strcat() starts at fmt[0]
for(i=0;i<narg;i++){
if(strchr(pfmt[i],'d') != NULL)strcat(fmt," d");
if(strchr(pfmt[i],'f') != NULL)strcat(fmt," f");
if(strchr(pfmt[i],'e') != NULL)strcat(fmt," f");
}
//cout<<"fmt "<<fmt<<endl;
//--------------------------------------
//get pointers to arrays
va_list args;
va_start(args, fmt);
while (*fmt != '\0') {
if(*fmt == 'd')pi[picnt++] = va_arg(args, int*);
if(*fmt == 'f')pd[pdcnt++] = va_arg(args, double*);
if(*fmt == 'e')pd[pdcnt++] = va_arg(args, double*);
//cout<<*fmt;
++fmt;
}//while
//cout<<endl;
va_end(args);
//printf("counts %d %d\n",picnt,pdcnt);
//--------------------------------------------
//print arrays
for(i = 0;i<nrow;i++){
picnt = 0;
pdcnt = 0;
for(j=0;j<narg;j++){
if(strchr(pfmt[j],'d') != NULL)printf(pfmt[j],*(pi[picnt++]+i));
if(strchr(pfmt[j],'f') != NULL)printf(pfmt[j],*(pd[pdcnt++]+i));
if(strchr(pfmt[j],'e') != NULL)printf(pfmt[j],*(pd[pdcnt++]+i));
if(j<narg-1)printf("\t"); //no tab after last column
}//for j
printf("\n");
}//for i
}//print_arrays()
int main(){
int n = 4;
int xi2[n] = {10,20,30,40},xi1[n] = {1,2,3,4};
double xd1[n] = {1.1,2.2,3.3,4.1};
double xd2[n] = {10.1,20.2,30.3,40.4};
char fmt[] = "%10d %5.3f %5d %1.2e";
print_arrays(n,fmt,xi1,xd1,xi2,xd2);
return 0;
}//fft class
#include <iostream>
#include "stdlib.h"
#include "stdio.h"
#include <complex>
using namespace std;
typedef complex<double> compvar;
const double _PI= 3.14159265358979323846;
class FFT{
public:
int n = 0;
compvar *tfc;
double *t,*y ,*freqs,*mag,*phase,fs;
int ffti(void);
int fft(void);
void scale(void);
void mag_phase(void);
};//FFT
void FFT::scale(void){
int i;
double dt = 1./fs,df = fs/n,f = 0;
for(i=0;i<n;i++){
t[i] = i * dt;
freqs[i] = f;
f += df;
}
}//fft_scale()
void FFT::mag_phase(void){
int i;
for(i=0;i<n;i++){
if(i==0){
mag[0] = real(tfc[0])/n;
phase[0] = 0;
}
else{
mag[i] = 2*abs(tfc[i])/n;
phase[i] = atan(imag(tfc[i])/real(tfc[i]))*180./_PI;
}
}
}//mag_phase()
int FFT::fft(void){
// n length of transform power of 2
// y real input signal
// tfc compplex transform out
cout<<"FFT Forward"<<endl;
if(n%2 != 0){cout<<"Not Power Of 2 "<<endl;return 1;}
compvar w,temp;
int i,j,k,n1,n2,a;
int m = log(n)/log(2); //exponet of 2^m
for(i=0;i<n;i++)tfc[i] = y[i] + 0.i;
// Bit reversal
j = 0; n2 = n/2;
for (i=1; i < n - 1; i++) {
n1 = n2;
while ( j >= n1 ) {j = j - n1;n1 = n1/2;}
j = j + n1;
if (i < j) {temp = tfc[i];tfc[i] = tfc[j];tfc[j] = temp;}
}//i
// Transform
n1 = 0;n2 = 1;
for (i=0; i < m; i++) {
n1 = n2;n2 += n2;a = 0;
for (j=0; j < n1; j++) {
w = exp(-2.i*_PI*double(a)/double(n)); //smpled sine wave
a += 1 << (m-i-1);
for (k=j; k < n; k=k+n2) {
temp = tfc[k+n1]*w;
tfc[k+n1] = tfc[k] - temp;
tfc[k] = tfc[k] + temp;
}//k
}//j
}//i
return 0;
}//fft()
int FFT::ffti(void){
// n length otransform power of 2
// yr real input signal
// tfc compplex transform out
cout<<"FFT Inverse"<<endl;
if(n%2 != 0){cout<<"Not Power Of 2 "<<endl;return 1;}
double tre,tim;
compvar w,temp;
compvar *tfcc = new compvar[n];
int i,j,k,n1,n2,a;
int m = log(n)/log(2); //exponet of 2^m
for(i=0;i<n;i++)tfcc[i] = tfc[i];
// Bit reversal
j = 0; n2 = n/2;
for (i=1; i < n - 1; i++) {
n1 = n2;
while ( j >= n1 ) {j = j - n1;n1 = n1/2;}
j = j + n1;
if (i < j) {temp = tfc[i];tfc[i] = tfc[j];tfc[j] = temp;}
}//i
// Transform
n1 = 0;n2 = 1;
for (i=0; i < m; i++) {
n1 = n2;n2 += n2;a = 0;
for (j=0; j < n1; j++) {
w = exp(2.i*_PI*double(a)/double(n));
a += 1 << (m-i-1);
for (k=j; k < n; k=k+n2) {
temp = tfc[k+n1]*w;
tfc[k+n1] = tfc[k] - temp;
tfc[k] = tfc[k] + temp;
y[k+n1] = real(tfc[k+n1])/n;
y[k] = real(tfc[k])/n;
}//k
}//j
}//i
for(i=0;i<n;i++)tfc[i] = tfcc[i];
delete [] tfcc;
return 0;
}//ffti()#include <iostream>
#include "stdlib.h"
#include "stdio.h"
#include <complex>
using namespace std;
const double _PI= 3.14159265358979323846;
typedef complex<double> compvar;
int nmax = pow(2,18);
compvar *tfc = new compvar[nmax];
double *yr = new double[nmax];
double *yi = new double[nmax];
double *mag = new double[nmax];
double *ph = new double[nmax];
double *t = new double[nmax];
double *f = new double[nmax];
int save_fft(int n,double *f,double *mag,double *phase,double *t,double *yr,double *yi){
cout<<"Save "<<endl;
int i;
FILE *p1 = fopen("fft.dat","w");
FILE *p2 = fopen("signals.dat","w");
if(!p1 || !p2){cout<<"fail";return 1;}
for(i=0;i<n;i++){
fprintf(p1,"% .20f %.20f %.20f\n",f[i],mag[i],phase[i]);
fprintf(p2,"% .20f %.20f %.20f\n",t[i],yr[i],yi[i]);
}
fclose(p1);fclose(p2);
return 0;
}//save()
class FFT{
public:
int n = 0;
compvar *tfc;
double *t,*y ,*freqs,*mag,*phase,fs;
int ffti(void);
int fft(void);
void scale(void);
void mag_phase(void);
};
void am(int n,double *t,double *y,double fc,double fm,double m,double a){
//amplitude modulation
int i;
for(i=0;i<n;i++)y[i] = (a+m*cos(2*_PI*fm*t[i]))*sin(2*_PI*fc*t[i]);
}
void rect(int n,double *t,double *y,double tlo,double thi,
double vlo,double vhi){
// recngular wave
double dt = t[2]-t[1], tcount = 0;
int i = 0;
while(1){
tcount = 0;
while(tcount <= thi) {
if(i == n) break;
y[i] = vhi;tcount += dt;i += 1;}
tcount = 0;
while(tcount <= tlo){
if(i == n) break;
y[i] = vlo ;tcount += dt;i += 1;}
if(i == n) break;
}
}//rect()
int main(void){
int i,n = pow(2,16);
cout<<n<<endl;
FFT fft;
fft.n = n;fft.fs = 1000;fft.t = t;
fft.y = yr;fft.tfc = tfc;fft.freqs = f;fft.mag = mag;fft.phase = ph;
fft.scale();
rect(n,t,yr,.025,.025,-1,1);
//am(n,t,yr,100,20,.8,4);
//for(i=0;i<n;i++)yr[i] = sin(2*_PI*20*t[i])-2; //sine wave
fft.fft();
fft.y = yi;
fft.ffti();
fft.mag_phase();
save_fft(n,f,mag,ph,t,yr,yi);
//excute Gnuplot scripts
string fnfft = path + "fft.plt";
string fnsig = path + "signals.plt";
system(fnfft.c_str());
system(fnsig.c_str());
return 0;int save(int n,double *x,double *y,string fname){
int i;
ofstream ofile;
ofile.open(fname);
if(!ofile) return 1;
ofile<<setprecision(10);
for(i=0;i<n;i++)ofile<<x[i]<<setw(10)<<"\t"<<y[i]<<endl;
ofile.close();
return 0;
}//save()
class GNUPLOTXY{
public:
string xlo,xhi,ylo,yhi;
string plot_file = "gnp.plt";
string data_file = "signals.dat";
int save(void);
void plot(void);
};
int GNUPLOTXY::save(void){
int i,n = 7;
string s[n];
s[0] = "set term windows background rgb 'white' fontscale 1 \n";
s[1] = "reset\n";
s[2] = "set grid lt 1 lw 1 lc rgb 'black' dashtype solid\n";
s[3] = "set zeroaxis\n";
s[4] = "set xrange["+xlo+":"+xhi+"]\n";
s[5] = "set yrange["+ylo+":"+yhi+"]\n";
s[6] = "plot '" + data_file + "' using 1:2 with lines ls 4 lt -1 lw 3\n";
ofstream ofile;
ofile.open(plot_file);
if(!ofile) return 1;
for(i=0;i<n;i++)ofile<<s[i];
ofile.close();
return 0;
}
void GNUPLOTXY::plot(void){system(plot_file.c_str());}
int main()
{
GNUPLOTXY gp;
int i,n = 200;
double t[n],y[n],dt = 1./(n-1);
for(i=0;i<n;i++){
t[i] = i*dt;
y[i] = sin(2*_PI*t[i]);
}
save(n,t, y,"sig.dat");
gp.plot_file = "gp.plt";
gp.data_file = "sig.dat";
gp.xlo = '0';gp.xhi = "1";
gp.ylo = "-1";gp.yhi = "1";
gp.save();
gp.plot();
return 0;
}int rand_int(unsigned seed){
//uniform distribution ints 0 - 2147483646
static unsigned rand_num = 1;
unsigned m = 0x7fffffff;
unsigned k1 = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
unsigned k2 = 8*ceil(m*atan(1.)/8.) + 5;
static unsigned temp_seed = 1;
if(seed != temp_seed){
temp_seed = seed;
rand_num = seed;
}
rand_num = rand_num*k2 + k1;
return rand_num%m;
}
double rand_dec(unsigned seed){
// unform distribution decimals 0 <= x < 1
static unsigned rand_num = 1;
unsigned m = 0x7fffffff;
unsigned k1 = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
unsigned k2 = 8*ceil(m*atan(1.)/8.) + 5;
static unsigned temp_seed = 1;
if(seed != temp_seed){
temp_seed = seed;
rand_num = seed;
}
rand_num = rand_num*k2 + k1;
return double(rand_num%m)/(m+1);
}
void rand_dec_test(void){
cout<<"rand_dec_test"<<endl;
long long i,n = pow(10,9);
double k = double(n)/10.,s = 0;
double x,min,max,pc[10];
int count[10];
for(i=0;i<10;i++)count[i] = 0;
min =1;
max = 0;
unsigned seed = time(NULL);
srand(seed);
for(i=0;i<n;i++){
//x = rand_dec(seed);
x = (double(rand()%RAND_MAX))/double(RAND_MAX+1);
if(x> 0 && x < .1)count[0]++;
if(x >= .1 && x < .2)count[1]++;
if(x >= .2 && x < .3)count[2]++;
if(x >= .3 && x < .4)count[3]++;
if(x >= .4 && x < .5)count[4]++;
if(x >= .5 && x < .6)count[5]++;
if(x >= .6 && x < .7)count[6]++;
if(x >= .7 && x < .8)count[7]++;
if(x >= .8 && x < .9)count[8]++;
if(x >= .9 && x <1.)count[9]++;
if(x > max)max = x;
if(x < min)min = x;
}
for(i=0;i<10;i++){
pc[i] = 100*(double(count[i]) - k)/k;
s += abs(pc[i]);
}
s = s/10; //average deviation
for(i=0;i<10;i++)printf("%10d \t%+.5f\n",count[i],pc[i]);
printf("\nAverage Deviation %.5f\n",s);
printf(" min %.15f max %.15f\n",min,max);
}//rand_dec_test()
int main(){
cout<<"RUNNING"<<endl;
rand_dec_test();
return 0;
#include <random>
#include <chrono>
using namespace std;
int main(){
int i, n=100,yi[100];
double yd[100];
typedef mt19937_64 default_random_engine;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator (seed);
uniform_real_distribution<double> distr1 (0,1.);
uniform_int_distribution<int> distr2 (-5,5);
for(i=0;i<n;i++) {
yi[i] = distr2(generator) ;
yd[i] = distr1(generator) ;
}
for(i=0;i<n;i++)printf("%+5d \t %+5.15f\n",yi[i],yd[i]);
return 0;
}#include <random>
#include <chrono>
void save(int n,string fname, double *x1,double *x2){
int i;
FILE *p = fopen(fname.c_str(),"w");
for(i=0;i<n;i++)fprintf(p,"%.10f\t %.10f\t \n",x1[i],x2[i]);
fclose(p);
}//save()
int comp_doub_ascend(const void * a, const void * b) {
if (*(double*)a > *(double*)b) return 1;
else if (*(double*)a < *(double*)b) return -1;
else return 0;
}
void cum_dist(int n,double *y,double *cd){
int i;
qsort(y,n,sizeof(double),comp_doub_ascend);
for(i = 0;i < n;i++)cd[i] = 100.* double(i)/double(n-1);
}//cum_dist()
double median(int n,double *y){
//y is sorted ascending
int i ;
double c;
for(i=0;i<n;i++){
c = double(i)/double(n-1);
if(c >= 0.5)break;
}
return y[i];
}//median()
double avg(int n,double *y){
int i;
double s = 0,av;
for(i=0;i<n;i++)s += y[i];
return s/double(n);
}//avg()
double sdev(int n,double *y){
int i;
double s = 0,av;
for(i=0;i<n;i++)s += y[i];
av = s/double(n);
s = 0;
for(i=0;i<n;i++)s += pow(av-y[i],2);
return sqrt(s/double(n));
}//sdev()
double minf(int n,double *y){
int i;
double ymin = y[0];
for(i=0;i<n;i++)if(y[i] < ymin)ymin = y[i];
return ymin;
}
double maxf(int n,double *y){
int i;
double ymax = y[0];
for(i=0;i<n;i++)if(y[i] > ymax)ymax = y[i];
return ymax;
}
void normal(int n,double u,double std,double *y){
int i;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator (seed);
normal_distribution<double> dist(u, std);
for(i=0;i<n;i++)y[i] = dist(generator);
}
void expon(int n,double mean,double *y){
int i;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator (seed);
exponential_distribution<double> dist(1./mean);
for(i=0;i<n;i++)y[i] = dist(generator);
}
int main(){
int n2=10000;
int i, n =100,yi[n];
double *yn = new double[n2];
double *ye = new double[n2];
double *cde = new double[n2];
double *cdn = new double[n2];
double minn,maxn,mine,maxe;
double avgn,stdn,medn,avge,stde,mede;
normal(n2,125,20,yn);
expon(n2,20,ye);
cum_dist(n2,ye,cde);
avgn = avg(n2,yn);
cum_dist(n2,ye,cde);
cum_dist(n2,yn,cdn);
stdn = sdev(n2,yn);
medn = median(n2,yn);
minn = minf(n2,yn);
maxn = maxf(n2,yn);
printf("normal %.5f \t %.5f \t %.5f \t %.5f\n",avgn,stdn,medn,(maxn-minn)/stdn);
avge = avg(n2,ye);
stde = sdev(n2,ye);
mede = median(n2,ye);
printf("exponential %.5f \t %.5f \t %.5f %.5f\n",avge,stde,mede,log(2)*avge);
save(n2,"norm.dat", yn,cdn);
save(n2,"exp.dat", ye,cde);
return 0;
} unsigned m = 0x7fffffff;
unsigned k1 = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
unsigned k2 = 8*ceil(m*atan(1.)/8.) + 5; rand_num = rand_num*k2 + k1;