libm

Author: ktye

k+/l.c

@@ -1,17 +1,21 @@
+enum fn {
+#ifdef BLAS
+  nrm2, asum, imax,
+#endif
+};
 void reg(int64_t id, int64_t name, int32_t arity){
  int64_t l=ti(14,(uint32_t)l2(id, Ku(name)));
  SetI32((int32_t)l-12,arity);
  dx(Asn(sc(Ku(name)),l));
 }
 void libs(void){ //encode strings with: https://play.golang.org/p/4ethx6OEVCR
 #ifdef BLAS
- reg(0, 846033518ull,1); //nrm2
- reg(1,1836413793ull,1); //asum
- reg(2,2019650921ull,1); //imax
+ reg(nrm2, 846033518ull,1); //id k-name arity (name is same as id-enum)
+ reg(asum,1836413793ull,1);
+ reg(imax,2019650921ull,1);
 #endif
- //reg(1,...);
- //reg(2,...);
 }
+
 int64_t fzvec(int64_t x){
  int32_t xt=tp(x);
  if(xt<16){x=Enl(x);xt=tp(x);}
@@ -29,18 +33,18 @@ int cblas_izamax(int,double*,int);
 
 #define fp(x) ((double*)(M_+(int32_t)x))
 
-int64_t cnative(int64_t x, int64_t y){
+int64_t cnative(int64_t x, int64_t y){ //function-id, list of args(length arity)
  int32_t xt,xn,i;
  double r;
  switch(x){ //switch registered function id
 #ifdef BLAS
- case 0:x=fzvec(Fst(y));xn=nn(x);if(xt==Zt)xn*=2;
+ case nrm2:x=fzvec(Fst(y));xn=nn(x);if(xt==Zt)xn*=2;
   r=cblas_dnrm2(xn,fp(x),1);
   dx(x);return Kf(r);
- case 1:x=fzvec(Fst(y));xn=nn(x);if(xt==Zt)xn*=2;
+ case asum:x=fzvec(Fst(y));xn=nn(x);if(xt==Zt)xn*=2;
   r=cblas_dasum(xn,fp(x),1);
   dx(x);return Kf(r);
- case 2:x=fzvec(Fst(y));xn=nn(x);
+ case imax:x=fzvec(Fst(y));xn=nn(x);
   i=tp(x)==Ft?cblas_idamax(xn,fp(x),1):cblas_izamax(xn,fp(x),1);
   dx(x);return Ki(i);
 #endif

k+/l.h

@@ -1,3 +1,8 @@
 #define NATIVE 1
 void libs(void);
 int64_t cnative(int64_t,int64_t);
+static uint64_t se(double);
+void sincos(double,double*,double*);
+static void sin1(int32_t,int32_t,int32_t);
+static void cos1(int32_t,int32_t,int32_t);
+#define cosin_(x,r,_) sincos(x,(double*)(M_+8+r),(double*)(M_+r));

k+/m.awk

@@ -0,0 +1,29 @@
+# patch out built-in math library, replace with -lm */
+
+BEGIN                {x=1;    }
+/double hypot_\(/    {x=match($0,/;/);next} #cut both declaration&function
+/void cosin_\(/      {x=match($0,/;/);next}
+/void sin1\(/        {x=match($0,/;/);next}
+/void cos1\(/        {x=match($0,/;/);next}
+/double atan2_\(/    {x=match($0,/;/);next}
+/double atan_\(/     {x=match($0,/;/);next}
+/double satan\(/     {x=match($0,/;/);next}
+/double xatan\(/     {x=match($0,/;/);next}
+/double exp_\(/      {x=match($0,/;/);next}
+/double expmulti\(/  {x=match($0,/;/);next}
+/double ldexp_\(/    {x=match($0,/;/);next}
+/int32_t frexp1\(/   {x=match($0,/;/);next}
+/double frexp2\(/    {x=match($0,/;/);next}
+/int64_t frexp3\(/   {x=match($0,/;/);next}
+/double normalize\(/ {x=match($0,/;/);next}
+/double log_\(/      {x=match($0,/;/);next}
+/double modabsfi\(/  {x=match($0,/;/);next}
+/double pow_\(/      {x=match($0,/;/);next}
+/uint64_t se\(/      {x=match($0,/;/);next}
+/exp_/  { sub(/exp_/  ,"exp"  )}
+/log_/  { sub(/log_/  ,"log"  )}
+/pow_/  { sub(/pow_/  ,"pow"  )}
+/atan2_/{ sub(/atan2_/,"atan2")}
+/hypot_/{gsub(/hypot_/,"hypot")}
+{        if          (x) print}
+/^}/                 {x=1;    }

k+/m.c

@@ -0,0 +1,42 @@
+//reinsert lost functions implemented with libm
+
+/*
+static uint64_t se(double x){
+ int32_t ei;
+ double f;
+ int64_t e;
+ f=x;
+ e=(int64_t)(0ll);
+ if(frexp1(x)){
+  f=frexp2(x);
+  e=frexp3(x);
+ }
+ x=(0.3010299956639812*(double)(e));
+ ei=(int32_t)(F64floor(x));
+ x=(x-(double)(ei));
+ return ucat(cat1(sf((f*pow_(10.,x))),Kc(101)),si(ei));
+}
+static void sin1(int32_t xp, int32_t yp, int32_t rp){
+ cosin_(F64(xp),rp,1);
+}
+static void cos1(int32_t xp, int32_t yp, int32_t rp){
+ cosin_(F64(xp),rp,2);
+}
+*/
+static uint64_t se(double x){
+ int32_t ei;
+ double f;
+ int e;
+ f=frexp(x,&e);
+ x=(0.3010299956639812*(double)(e));
+ ei=(int32_t)(F64floor(x));
+ x=(x-(double)(ei));
+ return ucat(cat1(sf((f*pow(10,x))),Kc(101)),si(ei));
+}
+static void sin1(int32_t xp, int32_t yp, int32_t rp){
+ SetF64(rp,sin(F64(xp)));
+}
+static void cos1(int32_t xp, int32_t yp, int32_t rp){
+ SetF64(rp,cos(F64(xp)));
+}
+

k+/readme

@@ -1,17 +1,36 @@
-# extend k.c with some c-libraries using the native function type
-# blas lapack fftw lp_solve ..
+# k+ extends k.c with some c-libraries using the native function type.
+# blas lapack fftw lp_solve ode specfun ..
+#
+# m.awk also replaces the builtin math functions with libm.
+#
+# example gemm             C := a*op(A)*op(B) + beta*C
+# gemm[a;A;rA;B;rB;op;b;C] C is reused only if refcount 1, op is 012 notrans/trans/herm
+# matrices(A, B) are passed as vectors, e.g. ,/3^?12
+# rA is the number of rows, if negative A is stored in column major format.
+# e.g. A:3^?12 is a list of length 3 containing length-4 vectors.
+# if these represent A matrix with 3 columns, rA is -4.
+# if they are a 3 rows of a matrix, rA is 3 
+
+# BLAS    /usr/include/x86_64-linux-gnu/cblas-netlib.h
+# nrm2[x] FZ l2 vector norm
+# asum[x] FZ l1 vector norm
+# imax[x] FZ index of l1 max element
 
 set -e
 set -x
 if [ "$1" = "get" ];then
  wget -Okc.orig https://github.com/ktye/i/releases/download/latest/k.c 
 fi
 
-sed -e 's/math.h>$/math.h>\n#include"l.h"/' -e 's/kinit()/kinit();libs()/' kc.orig>k.c
-cat l.c                                                                          >>k.c
-
-CC=cc
-DD="-D BLAS"
-
+awk -f m.awk kc.orig |\
+sed -e 's/math.h>$/math.h>\n#include"l.h"/' \
+    -e 's/kinit()/kinit();libs()/' \
+        >k.c
+cat m.c>>k.c
+cat l.c>>k.c
 
-$CC -O3 -march=native -ok+ $DD k.c -lopenblas -lm
+ C=cc
+ O="-O0 -march=native"
+ D="-D BLAS   "
+ L="-lopenblas"
+$C $O $D -ok+ k.c $L -lm