/*
* Euler - a numerical lab
*
* platform : neutral
*
* file : scalp.c -- long accumulator
*/
#include <math.h>
#include <float.h>
#include <string.h>
#include "interval.h"
#include "express.h"
#include "output.h"
#include "scalp.h"
/************************ Exact Scalar Product ******************/
#define LR_SIZE 128
#define LR_BASE 64
#define LR_U_MIN -64
#define LR_E_MAX 63
#define LR_DOUBLE_LENGTH 8
#define error(s) error=1, output(s"\n")
typedef struct
{ int E,U,S;
unsigned short M[LR_SIZE];
} Accu;
Accu A1,A2,h,g,h1,g1,hh;
int accumode=s_real;
#define accu(a,n) (((a)->M)[LR_BASE+(n)])
#define iszero(a) ((a)->U==(a)->E && accu((a),(a)->U)==0)
#define MI 65536.0
static double EXP2 [17] =
{1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384.0,
32768.0,65536.0};
static void accutruncate (Accu *a)
{ while (accu(a,a->E)==0 && a->E>a->U) a->E--;
while (accu(a,a->U)==0 && a->E>a->U) a->U++;
if (a->E==a->U && accu(a,a->E)==0)
{ a->E=a->U=0; accu(a,0)=0; a->S=1;
}
}
static void accuload (Accu *a, double x)
{ double y;
int n;
if (x==0) { a->S=1; a->E=a->U=0; accu(a,0)=0; return; }
if (x<0) { a->S=-1; x=-x; }
else a->S=1;
if (x>1)
{ x=frexp(x,&n);
a->E=n/16; n-=a->E*16; x*=EXP2[n];
a->U=a->E;
while (x!=0.0)
{ accu(a,a->U)=(unsigned short)(y=floor(x)); x-=y; a->U--;
x*=MI;
}
a->U++;
if (a->U<=LR_U_MIN || a->U<a->E-LR_DOUBLE_LENGTH) return;
}
else if (x<1)
{ x=frexp(x,&n);
a->E=n/16-1; n-=a->E*16; x*=EXP2[n];
a->U=a->E;
while (x!=0.0)
{ accu(a,a->U)=(unsigned short)(y=floor(x)); x-=y; a->U--;
x*=MI;
}
a->U++;
if (a->U<=LR_U_MIN || a->U<a->E-LR_DOUBLE_LENGTH) return;
}
else
{ accu(a,0)=1; a->E=0; a->U=0;
}
accutruncate(a);
}
#define exp16(n) ldexp(1,(n)*16)
static double accuget (Accu *a)
{ double y;
int j;
y=accu(a,a->E);
for (j=a->E-1; j>=a->U && j>=a->E-LR_DOUBLE_LENGTH; j--)
{ y=y*MI+accu(a,j);
}
return a->S*y*exp16(j+1);
}
static void accucopy (Accu *a, Accu *b)
{ int i;
a->E=b->E; a->U=b->U; a->S=b->S;
for (i=a->U; i<=a->E; i++) accu(a,i)=accu(b,i);
}
static void accuadd (Accu *a, Accu *b);
static void accusub (Accu *a, Accu *b)
{ unsigned long l,l1;
int i,j;
int o=0;
if (iszero(b)) return;
if (iszero(a)) { accucopy(a,b); a->S=-a->S; return; }
if (b->S!=a->S)
{ if (a->S<0) { b->S=-1; accuadd(a,b); b->S=1; return; }
else { b->S=1; accuadd(a,b); b->S=-1; return; }
}
else
{ if (b->E<a->E)
{ greater :
if (b->U<a->U)
{ for (j=b->U; j<a->U; j++) accu(a,j)=0;
a->U=b->U;
}
i=b->U;
while (i<=b->E)
{ l=accu(a,i); l1=(unsigned long)accu(b,i)+o;
if (l<l1) { l+=65536l; o=1; }
else o=0;
accu(a,i)=(unsigned short)(l-l1);
i++;
}
if (o) while (1)
{ if (accu(a,i)>0) { accu(a,i)--; break; }
accu(a,i)=65535u;
i++;
}
}
else if (b->E>a->E)
{ less :
i=a->U;
if (b->U<a->U)
{ for (j=b->U; j<a->U; j++) accu(a,j)=accu(b,j);
a->U=b->U;
}
while (i<=a->E)
{ if (i>=b->U) l=accu(b,i);
else l=0;
l1=(unsigned long)accu(a,i)+o;
if (l<l1) { l+=65536l; o=1; }
else o=0;
accu(a,i)=(unsigned short)(l-l1);
i++;
}
if (o) while (1)
{ if (i>=b->U && accu(b,i)>0)
{ accu(a,i)=accu(b,i)-1; i++; break;
}
accu(a,i)=65535u;
i++;
}
for (; i<=b->E; i++) accu(a,i)=accu(b,i);
a->E=b->E;
a->S=-a->S;
}
else
{ i=a->E;
while (i>=a->U)
{ if (i<b->U) goto greater;
if (accu(a,i)<accu(b,i)) goto less;
if (accu(a,i)>accu(b,i)) goto greater;
i--;
}
if (i>=b->U) goto less;
accu(a,0)=0; a->E=a->U=0; a->S=1;
return;
}
}
accutruncate(a);
}
static void accuadd (Accu *a, Accu *b)
{ unsigned long l;
int i,j;
int o=0;
if (iszero(b)) return;
if (iszero(a)) { accucopy(a,b); return; }
if (b->S==a->S)
{ if (b->U>a->E)
{ for (j=a->E+1; j<b->U; j++) accu(a,j)=0;
a->E=b->U-1;
}
i=b->U;
while (i<=b->E)
{ if (i<a->U)
{ accu(a,i)=accu(b,i);
}
else
{ if (i>a->E)
{ l=(unsigned long)accu(b,i)+o;
a->E=i;
}
else l=((unsigned long)accu(a,i)+accu(b,i))+o;
accu(a,i)=(unsigned short)(l&65535l);
if (l>=65536l) o=1;
else o=0;
}
i++;
}
while (o)
{ if (i>LR_E_MAX)
{ output("Accumulator overflow.\n");
error=1; return;
}
if (i>a->E) { accu(a,i)=1; a->E=i; i++; break; }
if (i<a->U) { accu(a,i)=1; i++; break; }
l=(unsigned long)accu(a,i)+1;
accu(a,i)=(unsigned short)(l&65535l);
if (l>=65536l) o=1;
else o=0;
i++;
}
if (i<a->U) for (j=i; j<a->U; j++) accu(a,j)=0;
if (b->U<a->U) a->U=b->U;
}
else if (b->S<0) { b->S=1; accusub(a,b); b->S=-1; return; }
else if (a->S<0) { b->S=-1; accusub(a,b); b->S=1; return; }
accutruncate(a);
}
static void accumult (Accu *a, Accu *b)
{ unsigned long l,l1,o1;
static Accu hm;
int i,j,k,o,start=1;
if (iszero(a)) return;
if (iszero(b)) { accuload(a,0.0); return; }
for (i=a->U; i<=a->E; i++)
{ for (o1=0,o=0,j=b->U,k=i+j; j<=b->E; j++,k++)
{ if (k>LR_E_MAX || k<LR_U_MIN)
{ output("Accumulator overflow.\n");
error=1; return;
}
l=(unsigned long)accu(a,i)*accu(b,j)+o1;
l1=(start?0:(unsigned long)accu(&hm,k))+(l&65535l)+o;
if (l1>=65536l)
{ accu(&hm,k)=(unsigned short)(l1&65535l); o=1;
}
else
{ accu(&hm,k)=(unsigned short)l1; o=0;
}
o1=(l>>16)&65535l;
}
start=0;
if (k>LR_E_MAX || k<LR_U_MIN)
{ output("Accumulator overflow.\n");
error=1; return;
}
accu(&hm,k)=(unsigned short)(o1+o);
}
if (a->S==b->S) a->S=1;
else a->S=-1;
a->U=b->U+a->U;
a->E=b->E+a->E;
if (accu(&hm,a->E+1)!=0) a->E++;
for (i=a->U; i<=a->E; i++) accu(a,i)=accu(&hm,i);
accutruncate(a);
}
static void accuimult (double a, double b, double a1, double b1)
{ if (a>=0)
{ if (a1>=0)
{ accuload(&h,a); accuload(&g,a1);
accuload(&h1,b); accuload(&g1,b1);
}
else if (b1>=0)
{ accuload(&h,b); accuload(&g,a1);
accuload(&h1,b); accuload(&g1,b1);
}
else
{ accuload(&h,b); accuload(&g,a1);
accuload(&h1,a); accuload(&g1,b1);
}
}
else if (b>=0)
{ if (a1>=0)
{ accuload(&h,a); accuload(&g,b1);
accuload(&h1,b); accuload(&g1,b1);
}
else if (b1>=0)
{ if (a*b1 <= b*a1)
{ accuload(&h,a); accuload(&g,b1);
}
else
{ accuload(&h,b); accuload(&g,a1);
}
if (b*b1 >= a*a1)
{ accuload(&h1,b); accuload(&g1,b1);
}
else
{ accuload(&h1,a); accuload(&g1,a1);
}
} else
{ accuload(&h,b); accuload(&g,a1);
accuload(&h1,a); accuload(&g1,b1);
}
}
else
{ if (a1>=0)
{ accuload(&h,a); accuload(&g,b1);
accuload(&h1,b); accuload(&g1,a1);
}
else if (b1>=0)
{ accuload(&h,a); accuload(&g,b1);
accuload(&h1,a); accuload(&g1,a1);
}
else
{ accuload(&h,b); accuload(&g,b1);
accuload(&h1,a); accuload(&g1,a1);
}
}
accumult(&g,&h); accuadd(&A1,&g); if (error) return;
accumult(&g1,&h1); accuadd(&A2,&g1);
}
void maccuload (header *hd)
{ header *st=hd,*result;
int r,c;
double *m;
hd=getvalue(hd); if (error) return;
if (hd->type==s_real || hd->type==s_matrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&A1,*m++);
accumode=s_real;
while (c>1)
{ accuload(&h,*m++); accuadd(&A1,&h); c--;
if (error) return;
}
result=new_real(accuget(&A1),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_complex || hd->type==s_cmatrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&A1,*m++);
accuload(&A2,*m++);
accumode=s_complex;
while (c>1)
{ accuload(&h,*m++); accuadd(&A1,&h);
if (error) return;
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
c--;
}
result=new_complex(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_interval || hd->type==s_imatrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&A1,*m++);
accuload(&A2,*m++);
accumode=s_interval;
while (c>1)
{ accuload(&h,*m++); accuadd(&A1,&h);
if (error) return;
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
c--;
}
result=new_interval(
round_down(accuget(&A1)),round_up(accuget(&A2)),"");
if (error) return;
moveresult(st,result);
return;
}
else error("Wrong argument for accuload.\n");
}
void maccuadd (header *hd)
{ header *st=hd,*result;
int r,c;
double *m;
hd=getvalue(hd); if (error) return;
if (hd->type==s_real || hd->type==s_matrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
if (accumode==s_real || accumode==s_complex)
{ accuload(&h,*m++); accuadd(&A1,&h);
if (error) return;
}
else
{ accuload(&h,*m++);
accuadd(&A1,&h); accuadd(&A2,&h);
if (error) return;
}
accumode=s_real;
while (c>1)
{ if (accumode==s_real || accumode==s_complex)
{ accuload(&h,*m++); accuadd(&A1,&h);
if (error) return;
}
else
{ accuload(&h,*m++);
accuadd(&A1,&h); accuadd(&A2,&h);
if (error) return;
}
c--;
}
result=new_real(accuget(&A1),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_complex || hd->type==s_cmatrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
if (accumode==s_interval)
{ output("Cannot add complex to an interval accu.\n");
error=1; return;
}
if (accumode==s_real)
{ accuload(&A2,0.0); accumode=s_complex;
}
accuload(&h,*m++); accuadd(&A1,&h);
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
while (c>1)
{ accuload(&h,*m++); accuadd(&A1,&h);
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
c--;
}
result=new_complex(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_interval || hd->type==s_imatrix)
{ getmatrix(hd,&r,&c,&m);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
if (accumode==s_complex)
{ output("Cannot add interval to a complex accu.\n");
error=1; return;
}
if (accumode==s_real)
{ A2=A1; accumode=s_interval;
}
accuload(&h,*m++); accuadd(&A1,&h);
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
while (c>1)
{ accuload(&h,*m++); accuadd(&A1,&h);
accuload(&h,*m++); accuadd(&A2,&h);
if (error) return;
c--;
}
result=new_interval(
round_down(accuget(&A1)),round_up(accuget(&A2)),"");
if (error) return;
moveresult(st,result);
return;
}
else error("Wrong argument for accuload.\n");
}
void maccuload2 (header *hd)
{ header *st=hd,*hd1,*result;
int r,c,r1,c1;
double *m,*m1;
hd1=nextof(hd);
equal_params_2(&hd,&hd1); if (error) return;
if (hd->type==s_real || hd->type==s_matrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1 || r1!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
if (c==0 || c!=c1)
{ output("Wrong vectors in accuload.\n");
error=1; return;
}
accuload(&A1,*m++);
accuload(&h,*m1++);
accumult(&A1,&h);
accumode=s_real;
if (error) return;
while (c>1)
{ accuload(&h,*m++); accuload(&h1,*m1++);
accumult(&h,&h1);
accuadd(&A1,&h); c--;
if (error) return;
}
result=new_real(accuget(&A1),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_complex || hd->type==s_cmatrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1) { output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&h,*m++);
accuload(&h1,*m++);
accuload(&g,*m1++);
accuload(&g1,*m1++);
accucopy(&hh,&h); accumult(&hh,&g); accucopy(&A1,&hh);
accucopy(&hh,&h1); accumult(&hh,&g1); accusub(&A1,&hh);
accucopy(&hh,&h); accumult(&hh,&g1); accucopy(&A2,&hh);
accucopy(&hh,&h1); accumult(&hh,&g); accuadd(&A2,&hh);
if (error) return;
accumode=s_complex;
while (c>1)
{ accuload(&h,*m++);
accuload(&h1,*m++);
accuload(&g,*m1++);
accuload(&g1,*m1++);
accucopy(&hh,&h); accumult(&hh,&g); accuadd(&A1,&hh);
accucopy(&hh,&h1); accumult(&hh,&g1); accusub(&A1,&hh);
accucopy(&hh,&h); accumult(&hh,&g1); accuadd(&A2,&hh);
accucopy(&hh,&h1); accumult(&hh,&g); accuadd(&A2,&hh);
if (error) return;
c--;
}
result=new_complex(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_interval || hd->type==s_imatrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&A1,0);
accuload(&A2,0);
accuimult(*m,*(m+1),*m1,*(m1+1)); m+=2; m1+=2;
if (error) return;
while (c>1)
{ accuimult(*m,*(m+1),*m1,*(m1+1)); m+=2; m1+=2;
if (error) return;
c--;
}
result=new_interval(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else error("Wrong argument for accuload.\n");
}
void maccuadd2 (header *hd)
{ header *st=hd,*hd1,*result;
int r,c,r1,c1;
double *m,*m1;
hd1=nextof(hd);
equal_params_2(&hd,&hd1); if (error) return;
if (hd->type==s_real || hd->type==s_matrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1 || r1!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
if (c==0 || c!=c1)
{ output("Wrong vectors in accuload.\n");
error=1; return;
}
accuload(&h,*m++);
accuload(&h1,*m1++);
accumult(&h,&h1);
accuadd(&A1,&h);
accumode=s_real;
if (error) return;
while (c>1)
{ accuload(&h,*m++); accuload(&h1,*m1++);
accumult(&h,&h1);
accuadd(&A1,&h); c--;
if (error) return;
}
result=new_real(accuget(&A1),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_complex || hd->type==s_cmatrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuload(&h,*m++);
accuload(&h1,*m++);
accuload(&g,*m1++);
accuload(&g1,*m1++);
accucopy(&hh,&h); accumult(&hh,&g); accuadd(&A1,&hh);
accucopy(&hh,&h1); accumult(&hh,&g1); accusub(&A1,&hh);
accucopy(&hh,&h); accumult(&hh,&g1); accuadd(&A2,&hh);
accucopy(&hh,&h1); accumult(&hh,&g); accuadd(&A2,&hh);
accumode=s_complex;
if (error) return;
while (c>1) { accuload(&h,*m++);
accuload(&h1,*m++);
accuload(&g,*m1++);
accuload(&g1,*m1++);
accucopy(&hh,&h); accumult(&hh,&g); accuadd(&A1,&hh);
accucopy(&hh,&h1); accumult(&hh,&g1); accusub(&A1,&hh);
accucopy(&hh,&h); accumult(&hh,&g1); accuadd(&A2,&hh);
accucopy(&hh,&h1); accumult(&hh,&g); accuadd(&A2,&hh);
if (error) return;
c--;
}
result=new_complex(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else if (hd->type==s_interval || hd->type==s_imatrix)
{ getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
if (r!=1)
{ output("Accuload needs one or two vectors.\n");
error=1; return;
}
accuimult(*m,*(m+1),*m1,*(m1+1)); m+=2; m1+=2;
if (error) return;
while (c>1)
{ accuimult(*m,*(m+1),*m1,*(m1+1)); m+=2; m1+=2;
if (error) return;
c--;
}
result=new_interval(accuget(&A1),accuget(&A2),"");
if (error) return;
moveresult(st,result);
return;
}
else error("Wrong argument for accuload.\n");
}
void mresiduum (header *hd)
{ header *st=hd,*hd1,*hd2,*result;
int i,j,k,c,r,c1,r1,c2,r2,scalar=0;
double *m,*m1,*m2,*mr;
hd1=nextof(st); hd2=nextof(hd1);
equal_params_3(&hd,&hd1,&hd2); if (error) return;
getmatrix(hd,&r,&c,&m);
getmatrix(hd1,&r1,&c1,&m1);
getmatrix(hd2,&r2,&c2,&m2);
if (r2==1 && c2==1)
{ scalar=1;
if (c!=r1)
{ error("Wrong matris size for residuum.\n");
return;
}
}
else if (c1!=c2 || r!=r2 || c!=r1)
{ error("Wrong matris size for residuum.\n");
return;
}
if (isreal(hd))
{ result=new_matrix(r,c1,""); if (error) return;
mr=matrixof(result);
for (i=0; i<r; i++)
for (j=0; j<c1; j++)
{ if (scalar) accuload(&A1,-*m2);
else accuload(&A1,-*mat(m2,c2,i,j));
for (k=0; k<c; k++)
{ accuload(&h,*mat(m,c,i,k));
accuload(&h1,*mat(m1,c1,k,j));
accumult(&h,&h1);
accuadd(&A1,&h);
if (error) return;
}
*mat(mr,c1,i,j)=accuget(&A1);
}
moveresult(st,result);
}
else if (iscomplex(hd))
{ result=new_cmatrix(r,c1,""); if (error) return;
mr=matrixof(result);
for (i=0; i<r; i++)
for (j=0; j<c1; j++)
{ if (scalar)
{ accuload(&A1,-*m2);
accuload(&A2,-*(m2+1));
}
else
{ accuload(&A1,-*cmat(m2,c2,i,j));
accuload(&A2,-*(cmat(m2,c2,i,j)+1));
}
for (k=0; k<c; k++)
{ accuload(&h,*cmat(m,c,i,k));
accuload(&h1,*(cmat(m,c,i,k)+1));
accuload(&g,*cmat(m1,c1,k,j));
accuload(&g1,*(cmat(m1,c1,k,j)+1));
accucopy(&hh,&h); accumult(&hh,&g); accuadd(&A1,&hh);
accucopy(&hh,&h1); accumult(&hh,&g1); accusub(&A1,&hh);
accucopy(&hh,&h); accumult(&hh,&g1); accuadd(&A2,&hh);
accucopy(&hh,&h1); accumult(&hh,&g); accuadd(&A2,&hh);
if (error) return;
}
*cmat(mr,c1,i,j)=accuget(&A1);
*(cmat(mr,c1,i,j)+1)=accuget(&A2);
}
moveresult(st,result);
}
else if (isinterval(hd))
{ result=new_imatrix(r,c1,""); if (error) return;
mr=matrixof(result);
for (i=0; i<r; i++)
for (j=0; j<c1; j++)
{ if (scalar)
{ accuload(&A2,-*m2);
accuload(&A1,-*(m2+1));
}
else
{ accuload(&A2,-(*imat(m2,c2,i,j)));
accuload(&A1,-(*(imat(m2,c2,i,j)+1)));
}
for (k=0; k<c; k++)
{ accuimult(*imat(m,c,i,k),*(imat(m,c,i,k)+1),
*imat(m1,c1,k,j),*(imat(m1,c1,k,j)+1));
if (error) return;
}
*imat(mr,c1,i,j)=accuget(&A1);
*(imat(mr,c1,i,j)+1)=accuget(&A2);
}
moveresult(st,result);
}
else error("Wrong argument for residuum!\n");
}
void accuinit ()
{ accuload(&A1,0.0); accuload(&A2,0.0);
accumode=s_real;
}
void maccu1 (header *hd)
{ header *result;
double *mr;
int i;
result=new_matrix(1,A1.E-A1.U+1,"");
mr=matrixof(result);
for (i=A1.E; i>=A1.U; i--) *mr++=A1.S*accu(&A1,i)*exp16(i);
}
void maccu2 (header *hd)
{ header *result;
double *mr;
int i;
result=new_matrix(1,A2.E-A2.U+1,"");
if (accumode==s_real)
{ output("Accu is real.\n");
error=1; return;
}
mr=matrixof(result);
for (i=A2.E; i>=A2.U; i--) *mr++=A2.S*accu(&A2,i)*exp16(i);
}
syntax highlighted by Code2HTML, v. 0.9.1