summaryrefslogtreecommitdiff
path: root/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
diff options
context:
space:
mode:
Diffstat (limited to 'texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h')
-rw-r--r--texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h788
1 files changed, 678 insertions, 110 deletions
diff --git a/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h b/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
index bfcbbf23..27a5434f 100644
--- a/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
+++ b/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h
@@ -1,76 +1,8 @@
-int length2(double gdata[][2]);
-int length3(double gdata[][3]);
-int length4(double gdata[][4]);
-int length5(double gdata[][5]);
-int length2i(int gdata[][2]);
-void copyd(int from, int upto, double q[], double p[]);
-void copyi(int from, int upto, int q[], int p[]);
-int addptd2(double gdata[][2], double pt[2]);
-int addptd3(double gdata[][3], double pt[3]);
-int addptd4(double gdata[][4], double pt[4]);
-int addptd5(double gdata[][5], double pt[5]);
-int add2(double gdata[][2], double x, double y);
-int add3(double gdata[][3], double x, double y, double z);
-int push2(double x, double y, double ptL[][2]);
-int push3(double x, double y, double z, double ptL[][3]);
-int push4(double x, double y, double z, double w, double ptL[][4]);
-int push5(double x, double y, double z, double w, double u, double ptL[][5]);
-void pull2(int num, double ptL[][2],double pt[2]);
-void pull3(int num, double ptL[][3],double pt[3]);
-void pull4(int num, double ptL[][4],double pt[4]);
-void pull5(int num, double ptL[][5],double pt[5]);
-int appendd2(int level, int from,int upto, double gdata[][2],double mat[][2]);
-int appendd3(int level, int from,int upto, double gdata[][3],double mat[][3]);
-int prependd2(int from,int upto, double gdata[][2],double mat[][2]);
-int appendd1(double numvec[], double num);
-int appendi1(int numvec[],int num);
-int prependi1(int num, int numvec[]);
-int removei2(int nrmv,int mat[][2]);
-int dataindexi1(int out[],int din[][2]);
-int dataindexd1(double out[],int din[][2]);
-int dataindexd2(int level,double out[][2],int din[][2]);
-int dataindexd3(int level, double out[][3],int din[][2]);
-int dataindexd4(int level, double out[][4],int din[][2]);
-void dispmatd1(int from, int upto, double mat[]);
-void dispmatd2(int from, int upto, double mat[][2]);
-void dispmatd3(int from, int upto, double mat[][3]);
-void dispmatd4(int from, int upto, double mat[][4]);
-void dispmatd5(int from, int upto, double mat[][5]);
-void dispmatd6(int from, int upto, double mat[][6]);
-void dispmati2(int from, int upto, int mat[][2]);
-void dispmati1(int from, int upto, int mat[]);
-int readdata3(const char *fname, double out[][3]);
-int readdata2(const char *fname, double out[][2]);
-int output2(const char *var, const char *fname, int level, double out[][2]);
-void simplesort(double number[]);
-int memberi(int a, int list[]);
-double dotprod(int dim, double p[], double q[]);
-double dist(int dim, double p[], double q[]);
-void pointoncurve(double t, double gdata[][2], double eps, double pt[]);
-int partcrv(double a,double b,double pkL[][2], double pL[][2]);
-int partcrv3(double a,double b,double pkL[][3], double pL[][3]);
-int connectseg(double pdata[][2], double eps, double out[][2]);
-void koutenseg(double a[2], double b[2], double c[2], double d[2],
- double eps, double eps2, double out[4]);
-double paramoncurve(double p[2], int n, double ptL[][2]);
-void nearestptpt(double pa[2], double pL[][2], double ans[4]);
-int intersectselfPp(double data1[][2],double eps,double eps1s,double kL[][4]);
-int intersectcrvsPp(double data1[][2], double data2[][2],
- double eps,double eps1s, double kL[][4]);
-int dropnumlistcrv(double qd[][2], double eps1, int out[]);
-int increasept2(double ptL[][2], int nn, double out[][2]);
-int increasept3(double ptL[][3], int nn, double out[][3]);
-void bezierpt(double t, double p0[2], double p3[2],
- double p1[2], double p2[2], double ptout[2]);
-int bezier(double ptL[][2], double ctrL[][4], int num, double out[][2]);
-int ospline(int level, double ptL[][2], int num, double out[][2]);
-void bezierpt3(double t, double p0[3], double p3[3],
- double p1[3], double p2[3], double ptout[3]);
-int bezier3(double ptL[][3], double ctrL[][6], int num, double out[][3]);
-int ospline3(int level, double ptL[][3], int num, double out[][3]);
-
-const double Inf=99999;
-const int Infint=30000;
+//20180308 norm added
+
+int length1(double gdata[]){
+ return floor(gdata[0]+0.5);
+}
int length2(double gdata[][2]){
return floor(gdata[0][0]+0.5);
@@ -88,6 +20,10 @@ int length5(double gdata[][5]){
return floor(gdata[0][0]+0.5);
}
+int length6(double gdata[][6]){
+ return floor(gdata[0][0]+0.5);
+}
+
int length2i(int gdata[][2]){
return gdata[0][0];
}
@@ -140,6 +76,15 @@ int push5(double x, double y, double z, double w, double u, double ptL[][5]){
return nall;
}
+int push6(double x, double y, double z, double w, double u, double v,double ptL[][6]){
+ int nall=ptL[0][0];
+ nall++;
+ ptL[nall][0]=x; ptL[nall][1]=y;ptL[nall][2]=z;
+ ptL[nall][3]=w; ptL[nall][4]=u; ptL[nall][5]=v;
+ ptL[0][0]=nall;
+ return nall;
+}
+
void pull3(int num, double ptL[][3],double pt[3]){
pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2];
}
@@ -158,6 +103,11 @@ void pull5(int num, double ptL[][5],double pt[5]){
pt[2]=ptL[num][2]; pt[3]=ptL[num][3]; pt[4]=ptL[num][4];
}
+void pull6(int num, double ptL[][6],double pt[6]){
+ pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2];
+ pt[3]=ptL[num][3]; pt[4]=ptL[num][4]; pt[5]=ptL[num][5];
+}
+
int add2(double gdata[][2], double x, double y){
int nall;
nall=push2(x, y, gdata);
@@ -194,6 +144,12 @@ int addptd5(double gdata[][5], double pt[5]){
return nall;
}
+int addptd6(double gdata[][6], double pt[6]){
+ int nall;
+ nall=push6(pt[0],pt[1],pt[2],pt[3],pt[4],pt[5],gdata);
+ return nall;
+}
+
int appendd2(int level, int from,int upto, double gdata[][2],double mat[][2]){
int j, np, nall=length2(mat);
np=upto-from+1;
@@ -423,11 +379,20 @@ int dataindexi1(int out[],int din[][2]){
return ctr;
}
-void dispmatd1(int from, int upto, double mat[]){
- int i;
- for(i=from; i<upto+1; i++){
- printf("%d %7.5f\n",i,mat[i]);
- }
+void dispvec(int n,double dt[]){
+ int i;
+ for(i=0; i<n; i++){
+ printf("%f ",dt[i]);
+ }
+ printf("\n");
+}
+
+void dispveci(int n,int dt[]){
+ int i;
+ for(i=0; i<n; i++){
+ printf("%d ",dt[i]);
+ }
+ printf("\n");
}
void dispmatd2(int from, int upto, double mat[][2]){
@@ -468,6 +433,49 @@ void dispmatd6(int from, int upto, double mat[][6]){
}
}
+void dispmatd2all(double mat[][2]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f\n",i,mat[i][0],mat[i][1]);
+ }
+}
+
+void dispmatd3all(double mat[][3]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f\n",i,mat[i][0],mat[i][1],mat[i][2]);
+ }
+}
+
+void dispmatd4all(double mat[][4]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f\n",mat[i][2],mat[i][3]);
+ }
+}
+
+void dispmatd5all(double mat[][5]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][2],mat[i][3],mat[i][4]);
+ }
+}
+
+void dispmatd6all(double mat[][6]){
+ int i;
+ if(mat[0][0]==0){printf("%s\n","no data");return;}
+ for(i=1; i<mat[0][0]+1; i++){
+ printf("%3d %7.5f %7.5f %7.5f ",i,mat[i][0],mat[i][1],mat[i][2]);
+ printf("%7.5f %7.5f %7.5f\n",mat[i][3],mat[i][4],mat[i][5]);
+ }
+}
+
void dispmati2(int from, int upto, int mat[][2]){
int i;
for(i=from; i<upto+1; i++){
@@ -482,19 +490,19 @@ void dispmati1(int from, int upto, int mat[]){
}
}
-int readdata3(const char *fname, double out[][3]){
+void readdata3(const char *fname, double out[][3]){
double data[DsizeL][4];
int i, j, ndata, nall;
FILE *fp;
fp=fopen(fname,"r");
if(fp==NULL){
printf("file not found\n");
- return 0;
+ return;
}
ndata=0;
i=1;
while( ! feof(fp) && i<DsizeL){
- for(j=0;j<=3;j++){
+ for(j=0;j<4;j++){
fscanf(fp,"%lf",&data[i][j]);
}
i++;
@@ -504,9 +512,36 @@ int readdata3(const char *fname, double out[][3]){
fclose(fp);
out[0][0]=0; nall=0;
for(i=1; i<=ndata; i++){
- nall=add3(out, data[i][1], data[i][2], data[i][3]);
+ nall=add3(out, data[i][0], data[i][1], data[i][2]);
}
- return nall;
+ return;
+}
+
+void readdata2(const char *fname, double out[][2]){
+ double data[DsizeL][3];
+ int i, j, ndata, nall;
+ FILE *fp;
+ fp=fopen(fname,"r");
+ if(fp==NULL){
+ printf("file not found\n");
+ return;
+ }
+ ndata=0;
+ i=1;
+ while( ! feof(fp) && i<DsizeL){
+ for(j=0;j<3;j++){
+ fscanf(fp,"%lf",&data[i][j]);
+ }
+ i++;
+ ndata++;
+ }
+ ndata--;
+ fclose(fp);
+ out[0][0]=0; nall=0;
+ for(i=1; i<=ndata; i++){
+ nall=add2(out, data[i][0], data[i][1]);
+ }
+ return;
}
int output2(const char *var, const char *fname, int level, double out[][2]){
@@ -542,33 +577,6 @@ int output2(const char *var, const char *fname, int level, double out[][2]){
return 0;
}
-int readdata2(const char *fname, double out[][2]){
- double data[DsizeL][3];
- int i, j, ndata, nall;
- FILE *fp;
- fp=fopen(fname,"r");
- if(fp==NULL){
- printf("file not found\n");
- return 0;
- }
- ndata=0;
- i=1;
- while( ! feof(fp) && i<DsizeL){
- for(j=0;j<=2;j++){
- fscanf(fp,"%lf",&data[i][j]);
- }
- i++;
- ndata++;
- }
- ndata--;
- fclose(fp);
- out[0][0]=0; nall=0;
- for(i=1; i<=ndata; i++){
- nall=add2(out, data[i][1], data[i][2]);
- }
- return nall;
-}
-
void simplesort(double number[]){
int i,j, total=number[0];
double tmp;
@@ -602,6 +610,12 @@ double dotprod(int dim, double p[], double q[]){
return ans;
}
+double norm(int dim, double p[]){
+ double ans;
+ ans=dotprod(dim,p,p);
+ return sqrt(ans);
+}
+
double dist(int dim, double p[], double q[]){
double tmp;
tmp=pow(q[0]-p[0],2)+pow(q[1]-p[1],2);
@@ -611,6 +625,34 @@ double dist(int dim, double p[], double q[]){
return sqrt(tmp);
}
+void crossprod(int dim,double a[3],double b[3], double out[3]){
+ if(dim==3){
+ out[0]=a[1]*b[2]-a[2]*b[1];
+ out[1]=a[2]*b[0]-a[0]*b[2];
+ out[2]=a[0]*b[1]-a[1]*b[0];
+ }else{
+ out[0]=a[0]*b[1]-a[1]*b[0];
+ out[1]=Inf;
+ }
+}
+
+void reflectpoint(double p1[2],double regard[][2],double p2[2]){
+ double pa[2],pb[2],u,v,a,b,u2,v2;
+ pull2(1,regard,pa);
+ if(length2(regard)==1){
+ p2[0]=2*pa[0]-p1[0]; p2[1]=2*pa[1]-p1[1];
+ }else{
+ pull2(2,regard,pb);
+ u=pb[0]-pa[0]; u2=pow(u,2.0);
+ v=pb[1]-pa[1]; v2=pow(v,2.0);
+ a=pa[0]; b=pa[1];
+ p2[0]=(u2-v2)/(u2+v2)*p1[0]+2*u*v/(u2+v2)*p1[1];
+ p2[0]=p2[0]-2*v*(u*b-v*a)/(u2+v2);
+ p2[1]=2*u*v/(u2+v2)*p1[0]-(u2-v2)/(u2+v2)*p1[1];
+ p2[1]=p2[1]+2*u*(u*b-v*a)/(u2+v2);
+ }
+}
+
void pointoncurve(double t, double gdata[][2], double eps, double pt[]){
double pa[2],pb[2], s;
int n, nall;
@@ -1557,3 +1599,529 @@ int ospline3(int level, double ptL[][3], int num, double out[][3]){
}
return nall;
}
+
+void segplot(double x1,double y1,double x2,double y2,double seg[][2]){
+ seg[0][0]=0;
+ push2(x1,y1,seg);
+ push2(x2,y2,seg);
+}
+
+void addsegplot(double x1,double y1,double seg[][2]){
+ push2(x1,y1,seg);
+}
+
+void circledata(double cx,double cy,double r, int num,double out[][2]){
+ double t,dt,pt[2];
+ int i;
+ dt=2*M_PI/num;
+ out[0][0]=0;
+ for(i=0; i<=num; i++){
+ push2(cx+r*cos(i*dt),cy+r*sin(i*dt),out);
+ }
+}
+
+void intersectline(double p1[2],double v1[2], double p2[2],double v2[2],double out[4]){
+ double tmp,d,dt,ds,t,s,tmp1[2],tmp2[2],tmp3[2],pt[2];
+ tmp=dotprod(2,v1,v2);
+ tmp1[0]=tmp; tmp1[1]=pow(norm(2,v1),2.0);
+ tmp2[0]=-pow(norm(2,v2),2.0);tmp2[1]=-tmp;
+ pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1];
+ tmp3[0]=dotprod(2,pt,v2);tmp3[1]=dotprod(2,pt,v1);
+ crossprod(2,tmp1,tmp2,pt); d=pt[0];
+ crossprod(2,v1,v2,pt); tmp=fabs(pt[0]);
+ if(tmp>Eps){
+ crossprod(2,tmp3,tmp2,pt); dt=pt[0];
+ crossprod(2,tmp1,tmp3,pt); ds=pt[0];
+ t=dt/d;
+ s=ds/d;
+ pt[0]=p1[0]+v1[0]*t; pt[1]=p1[1]+v1[1]*t;
+ out[0]=pt[0];out[1]=pt[1];out[2]=t; out[3]=s;
+ }else{
+ pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1];
+ crossprod(2,pt,v1,tmp1);
+ tmp=tmp1[0]/norm(2,v1);
+ out[0]=fabs(tmp); out[1]=Inf;
+ }
+}
+
+void intersectseg(double seg1[][2],double seg2[][2], double out[5]){
+ double p1[2],q1[2],p2[2],q2[2],pt[2],pt1[2],nv[2],pts[DsizeS][4];
+ double q1p1[2],q2p2[2],p1q2[2],q1p2[2],p1p2[2],q1q2[2];
+ double tmp[2],tmp1[2],tmp2[2],tmp3[2],tmp4[2];
+ double tmpd1,x,y,t,s,t1,s1,distance,tmpd[4],tmp3d[7];
+ int ctr=0, i;
+ pull2(1,seg1,p1); pull2(2,seg1,q1);
+ pull2(1,seg2,p2); pull2(2,seg2,q2);
+ q1p1[0]=q1[0]-p1[0];q1p1[1]=q1[1]-p1[1];
+ q2p2[0]=q2[0]-p2[0];q2p2[1]=q2[1]-p2[1];
+ p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1];
+ q1p2[0]=q1[0]-p2[0];q1p2[1]=q1[1]-p2[1];
+ p1p2[0]=p1[0]-p2[0];p1p2[1]=p1[1]-p2[1];
+ p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1];
+ q1q2[0]=q1[0]-q2[0];q1q2[1]=q1[1]-q2[1];
+ if((norm(2,q1p1)<Eps)||(norm(2,q2p2)<Eps)){
+ out[0]=-1; out[1]=Inf;
+ }else{
+ intersectline(p1,q1p1,p2,q2p2,tmpd);
+ if(tmpd[1]<Inf-Eps){
+ pt[0]=tmpd[0];pt[1]=tmpd[1];t=tmpd[2];s=tmpd[3];
+ if((t*(t-1)<Eps)&&(s*(s-1)<Eps)){
+ out[0]=0;out[1]=pt[0];out[2]=pt[1];out[3]=t;out[4]=s;
+ }else{
+ if(t<0){t=0;}; if(t>1){t=1;};
+ if(s<0){s=0;}; if(s>1){s=1;};
+ tmp3d[ctr]=norm(2,p1p2);ctr++;
+ tmp3d[ctr]=norm(2,p1q2);ctr++;
+ tmp3d[ctr]=norm(2,q1p2);ctr++;
+ tmp3d[ctr]=norm(2,q1q2);ctr++;
+ tmp1[0]=q2p2[1];tmp1[1]=-q2p2[0];
+ intersectline(p1,tmp1,p2,q2p2,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,p1);ctr++;
+ }
+ }
+ intersectline(q1,tmp1,p2,q2p2,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,q1);ctr++;
+ }
+ }
+ tmp1[0]=q1p1[1];tmp1[1]=-q1p1[0];
+ intersectline(p2,tmp1,p1,q1p1,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,p2);ctr++;
+ }
+ }
+ intersectline(q2,tmp1,p1,q1p1,tmpd);
+ if(tmpd[1]!=Inf){
+ pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3];
+ if(s1*(s1-1)<Eps){
+ tmp3d[ctr]=dist(2,pt1,q2);ctr++;
+ }
+ }
+ tmpd1=tmp3d[0];
+ for(i=1;i<ctr;i++){
+ if(tmp3d[i]<tmpd1){tmpd1=tmp3d[i];};
+ }
+ out[0]=tmpd1;out[1]=pt[0];out[2]=pt[1];
+ out[3]=t;out[4]=s;
+ }
+ }else{
+ distance=tmpd[0];
+ nv[0]=q1p1[1]/norm(2,q1p1);nv[1]=-q1p1[0]/norm(2,q1p1);
+ pts[0][0]=0;
+ intersectline(p1,nv,p2,q2p2,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1];
+ push4(x,y,0,s1,pts);
+ }
+ intersectline(q1,nv,p2,q2p2,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1];
+ push4(x,y,1,s1,pts);
+ }
+ intersectline(p2,nv,p1,q1p1,tmpd);
+ s1=tmpd[3];
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1];
+ push4(x,y,s1,0,pts);
+ }
+ intersectline(q2,nv,p1,q1p1,tmpd);
+ if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){
+ x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1];
+ push4(x,y,s1,2,pts);
+ }
+ if(length4(pts)==0){
+ tmpd1=norm(2,p1p2);
+ if(norm(2,p1q2)<tmpd1){tmpd1=norm(2,p1q2);}
+ if(norm(2,q1p2)<tmpd1){tmpd1=norm(2,q1p2);}
+ if(norm(2,q1q2)<tmpd1){tmpd1=norm(2,q1q2);}
+ out[0]=tmpd1;out[1]=Inf;
+ }else{
+ if(distance>Eps1){
+ out[0]=distance;out[1]=Inf;
+ }else{
+ x=0; y=0;
+ for(i=1; i<=length4(pts); i++){
+ x=x+pts[i][0];
+ y=y+pts[i][1];
+ }
+ tmp3[0]=x/length4(pts); tmp3[1]=y/length4(pts);
+ nearestptpt(tmp3,seg1,tmpd);
+ t=tmpd[2];
+ nearestptpt(tmp3,seg2,tmpd);
+ s=tmpd[2];
+ out[0]=distance; out[1]=tmp3[0]; out[2]=tmp3[1];
+ out[3]=t;out[4]=s;
+ }
+ }
+ }
+ }
+}
+
+int osplineseg(double ptlist[5][2],int num,double out[][2]){
+ double p0[2],p1[2],p2[2],p3[2],p2p0[2],p3p1[2],p2p1[2];
+ double tmp,cc,pq[2],pr[2],knots[3][2],ctrlist[2][4];
+ pull2(1,ptlist,p0); pull2(2,ptlist,p1);
+ pull2(3,ptlist,p2); pull2(4,ptlist,p3);
+ p2p0[0]=p2[0]-p0[0];p2p0[1]=p2[1]-p0[1];
+ p3p1[0]=p3[0]-p1[0];p3p1[1]=p3[1]-p1[1];
+ p2p1[0]=p2[0]-p1[0];p2p1[1]=p2[1]-p1[1];
+ tmp=norm(2,p2p0)*norm(2,p3p1);
+ tmp=1+sqrt((1+dotprod(2,p2p0,p3p1)/tmp)/2);
+ cc=4*norm(2,p2p1)/3/(norm(2,p2p0)+norm(2,p3p1))/tmp;
+ pq[0]=p1[0]+cc*p2p0[0]; pq[1]=p1[1]+cc*p2p0[1];
+ pr[0]=p2[0]-cc*p3p1[0];pr[1]=p2[1]-cc*p3p1[1];
+ ctrlist[0][0]=0;
+ push4(pq[0],pq[1],pr[0],pr[1],ctrlist);
+ knots[0][0]=0;
+ addptd2(knots,p1);addptd2(knots,p2);
+ bezier(knots,ctrlist,num,out);
+ return length2(out);
+}
+
+void intersectpartseg(double crv1[][2],double crv2[][2],int ii, int jj,int num,double out[6]){
+ double distance=10*Eps2,seg1[3][2],seg2[3][2],tmpd[2],tmp1d[2],tmp2d[2];
+ double tmp,snang,dst,tmpd5[5],os1[DsizeS][2],os2[DsizeS][2];
+ double p0[2],p1[2],p2[2],p3[2],ptlist[5][2],regard[3][2];
+ double tmp1md[20][2],tmpd3[3],tmp3md3[20][3],tmp2md3[20][3];
+ double tmpd4[4],Eps00=pow(10,-8.0);
+ double sx,sy;
+ int kk,ll,nn;
+ out[0]=Inf;
+ for(kk=1;kk<6;kk++){out[kk]=0;}
+ seg1[0][0]=0;seg2[0][0]=0;
+ appendd2(1,ii,ii+1,crv1,seg1);
+ appendd2(1,jj,jj+1,crv2,seg2);
+ pull2(2,seg1,tmp1d);
+ pull2(1,seg1,tmpd);
+ tmp1d[0]=tmp1d[0]-tmpd[0]; tmp1d[1]=tmp1d[1]-tmpd[1];
+ pull2(2,seg2,tmp2d);
+ pull2(1,seg2,tmpd);
+ tmp2d[0]=tmp2d[0]-tmpd[0]; tmp2d[1]=tmp2d[1]-tmpd[1];
+ crossprod(2,tmp1d,tmp2d,tmpd);
+ snang=fabs(tmpd[0])/(norm(2,tmp1d)*norm(2,tmp2d));
+ intersectseg(seg1,seg2,tmpd5);
+ dst=tmpd5[0];
+ if(dst<Eps){
+ if(dst>-Eps){
+ out[0]=tmpd5[1];out[1]=tmpd5[2];out[2]=ii+tmpd5[3];
+ out[3]=jj+tmpd5[4];out[4]=dst;out[5]=snang;
+ }
+ }else{
+ if(dst<Eps2){
+ os1[0][0]=0;os2[0][0]=0;
+ pull2(1,seg1,tmp1d);
+ pull2(2,seg1,tmp2d);
+ tmpd[0]=tmp2d[0]-tmp1d[0];tmpd[1]=tmp2d[1]-tmp1d[1];
+ if((length2(crv1)==2)||(norm(2,tmpd)>distance-Eps)){
+ appendd2(1,1,2,seg1,os1);
+ }else{
+ pull2(1,seg1,p1);pull2(2,seg1,p2);
+ if(ii==1){
+ pull2(3,crv1,p3);
+ tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p3,regard,p0);
+ }else{
+ if(ii==length2(crv1)-1){
+ pull2(ii-1,crv1,p0);
+ tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p0,regard,p3);
+ }else{
+ pull2(ii-1,crv1,p0);pull2(ii+2,crv1,p3);
+ }
+ }
+ ptlist[0][0]=0;
+ addptd2(ptlist,p0);addptd2(ptlist,p1);
+ addptd2(ptlist,p2);addptd2(ptlist,p3);
+ osplineseg(ptlist,num,os1);
+ }
+ pull2(1,seg2,tmp1d);
+ pull2(2,seg2,tmp2d);
+ if((length2(crv2)==2)||(dist(2,tmp2d,tmp1d)>distance-Eps)){
+ appendd2(1,1,2,seg2,os2);
+ }else{
+ pull2(1,seg2,p1);pull2(2,seg2,p2);
+ if(jj==1){
+ pull2(3,crv2,p3);
+ tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p3,regard,p0);
+ }else{
+ if(jj==length2(crv2)-1){
+ pull2(jj-1,crv2,p0);
+ tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1];
+ tmp1d[0]=(p1[0]+p2[0])/2;
+ tmp1d[1]=(p1[1]+p2[1])/2;
+ tmp2d[0]=tmp1d[0]+tmpd[1];
+ tmp2d[1]=tmp1d[1]-tmpd[0];
+ regard[0][0]=0;
+ addptd2(regard,tmp1d); addptd2(regard,tmp2d);
+ reflectpoint(p0,regard,p3);
+ }else{
+ pull2(jj-1,crv2,p0);pull2(jj+2,crv2,p3);
+ }
+ }
+ ptlist[0][0]=0;
+ addptd2(ptlist,p0);addptd2(ptlist,p1);
+ addptd2(ptlist,p2);addptd2(ptlist,p3);
+ osplineseg(ptlist,num,os2);
+ }
+ tmp2md3[0][0]=0;
+ for(kk=1;kk<length2(os1);kk++){
+ for(ll=1;ll<length2(os2);ll++){
+ seg1[0][0]=0; seg2[0][0]=0;
+ pull2(kk,os1,tmpd);addptd2(seg1,tmpd);
+ pull2(kk+1,os1,tmpd);addptd2(seg1,tmpd);
+ pull2(ll,os2,tmpd);addptd2(seg2,tmpd);
+ pull2(ll+1,os2,tmpd);addptd2(seg2,tmpd);
+ intersectseg(seg1,seg2,tmpd5);
+ if((tmpd5[0]<Eps1)&&(tmpd5[1]<Inf-Eps)){
+ if(tmpd5[0]<dst+Eps00){
+ dst=tmpd5[0];
+ tmp3md3[0][0]=0;
+ for(nn=1;nn<=length3(tmp2md3);nn++){
+ pull3(nn,tmp2md3,tmpd3);
+ if(tmpd3[0]<dst+Eps00){
+ addptd3(tmp3md3,tmpd3);
+ }
+ }
+ tmp2md3[0][0]=0;
+ for(nn=1;nn<=length3(tmp3md3);nn++){
+ pull3(nn,tmp3md3,tmpd3);
+ addptd3(tmp2md3,tmpd3);
+ }
+ push3(dst,tmpd5[1],tmpd5[2],tmp2md3);
+ }
+ }
+ }
+ }
+ if(length3(tmp2md3)>0){
+ tmp1md[0][0]=0;
+ sx=0; sy=0;
+ for(nn=1;nn<=length3(tmp2md3);nn++){
+ pull3(nn,tmp2md3,tmpd3);
+ push2(tmpd3[1],tmpd3[2],tmp1md);
+ sx=sx+tmpd3[1];
+ sy=sy+tmpd3[2];
+ }
+ sx=sx/length3(tmp2md3);
+ sy=sy/length3(tmp2md3);
+ p0[0]=sx; p0[1]=sy;
+ pull2(ii,crv1,p1); pull2(ii+1,crv1,p2);
+ tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1];
+ tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0];
+ intersectline(p0,tmpd,p1,tmp1d,tmpd4);
+ tmp=tmpd4[3];
+ if(tmp<0){tmp=0;}
+ if(tmp>1){tmp=1;}
+ out[0]=sx;out[1]=sy;out[2]=ii+tmp;
+ pull2(jj,crv2,p1); pull2(jj+1,crv2,p2);
+ tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1];
+ tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0];
+ intersectline(p0,tmpd,p1,tmp1d,tmpd4);
+ tmp=tmpd4[3];
+ if(tmp<0){tmp=0;}
+ if(tmp>1){tmp=1;}
+ out[3]=jj+tmp;out[4]=dst;out[5]=snang;
+ }
+ }
+ }
+}
+
+void collectsameseg(double ptdL[][6],double rL[][6]){
+ double Eps00=pow(10,-8.0),tmp1d[2],tmp2d[2],tmp1d6[6],tmp2d6[6];
+ int ii,jj,kk,nr,ctr,numL[40],diffL[40];
+ double s1,e1,s2,e2,dst,tmp1,tmp2,tmp1md6[40][6];
+ rL[0][0]=0;
+ if(length6(ptdL)==0){
+ return;
+ }
+ for(ii=2;ii<=length6(ptdL);ii++){
+ pull6(ii,ptdL,tmp1d6);addptd6(rL,tmp1d6);
+ }
+ tmp1md6[0][0]=0;
+ pull6(1,ptdL,tmp1d6);
+ addptd6(tmp1md6,tmp1d6);
+ kk=tmp1d6[2];
+ if(tmp1d6[2]<kk+Eps00){
+ s1=kk-1-Eps00;e1=s1+2+2*Eps00;
+ }else{
+ s1=kk-Eps00; e1=s1+1+2*Eps00;
+ }
+ kk=tmp1d6[3];
+ if(tmp1d6[3]<kk+Eps00){
+ s2=kk-1-Eps00;e2=s2+2+2*Eps00;
+ }else{
+ s2=kk-Eps00; e2=s2+1+2*Eps00;
+ }
+ ctr=1;numL[0]=0;
+ for(ii=1;ii<=length6(rL);ii++){
+ pull6(ii,rL,tmp1d6);
+ tmp1=tmp1d6[2]; tmp2=tmp1d6[3];
+ if((tmp1>s1)&&(tmp1<e1)&&(tmp2>s2)&&(tmp2<e2)){
+ addptd6(tmp1md6,tmp1d6);
+ numL[ctr]=ii;ctr++;
+ }
+ }
+ dst=100;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmp1d6);
+ if(tmp1d6[4]<dst){
+ dst=tmp1d6[4];
+ }
+ }
+ ptdL[0][0]=0;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmp1d6);
+ if(tmp1d6[4]<dst+Eps00){
+ addptd6(ptdL,tmp1d6);
+ }
+ }
+ for(ii=0;ii<=length6(rL);ii++){
+ diffL[ii]=0;
+ }
+ for(ii=1;ii<=ctr;ii++){
+ kk=numL[ii];
+ diffL[kk]=1;
+ }
+ nr=length6(rL);
+ rL[0][0]=0;
+ for(ii=1;ii<=nr;ii++){
+ if(diffL[ii]==0){
+ pull6(ii,rL,tmp1d6);addptd6(rL,tmp1d6);
+ }
+ }
+}
+
+void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][6]){
+ double Eps00=pow(10,-8.0),Dist=10*Eps2, crv1[DsizeL][2],crv2[DsizeL][2];
+ double tmp1md6[DsizeS][6],tmp2md6[DsizeS][6];
+ double tmpd[2],tmp1d[2],tmp2d[2],tmpd6[6],tmp1d6[6];
+ double sx,sy,dst,ans[4];
+ int ii,jj,kk,nall,self,js,je;
+ crv1[0][0]=0;
+ pull2(1,crv1s,tmpd);addptd2(crv1,tmpd);
+ for(ii=2;ii<=length2(crv1s);ii++){
+ pull2(length2(crv1),crv1,tmp1d);
+ pull2(ii,crv1s,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>Eps){
+ addptd2(crv1,tmp2d);
+ }
+ }
+ crv2[0][0]=0;
+ pull2(1,crv2s,tmpd);
+ addptd2(crv2,tmpd);
+ for(ii=2;ii<=length2(crv2s);ii++){
+ pull2(length2(crv2),crv2,tmp1d);
+ pull2(ii,crv2s,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>Eps){
+ addptd2(crv2,tmp2d);
+ }
+ }
+ if(length2(crv1)!=length2(crv2)){
+ self=0;
+ }else{
+ self=1;
+ for(ii=1;ii<=length2(crv1);ii++){
+ pull2(ii,crv1,tmp1d);pull2(ii,crv2,tmp2d);
+ if(dist(2,tmp1d,tmp2d)>0){
+ self=0;
+ break;
+ }
+ }
+ }
+ tmp1md6[0][0]=0;
+ for(ii=1;ii<=length2(crv1)-1;ii++){
+ if(self==0){
+ js=1; je=length2(crv2)-1;
+ }else{
+ js=ii+2; je=length2(crv2)-1;
+ }
+ for(jj=js;jj<=je;jj++){
+ intersectpartseg(crv1,crv2,ii,jj,num,tmpd6);
+ if(tmpd6[0]<Inf-Eps){
+ tmpd[0]=tmpd6[0];tmpd[1]=tmpd6[1];
+ if(length6(tmp1md6)==0){
+ addptd6(tmp1md6,tmpd6);
+ }else{
+ pull6(length6(tmp1md6),tmp1md6,tmp1d6);
+ tmp1d[0]=tmp1d6[0];tmp1d[1]=tmp1d6[1];
+ if(dist(2,tmp1d,tmpd)>Eps1){
+ addptd6(tmp1md6,tmpd6);
+ }
+ }
+ if(self==1){
+ push6(tmpd6[0],tmpd6[1],tmpd6[3],tmpd6[2],tmpd6[4],tmpd6[5],tmp1md6);
+ }
+ }
+ }
+ }
+ out[0][0]=0;
+ tmp2md6[0][0]=0;
+ while(length6(tmp1md6)>0){
+ collectsameseg(tmp1md6,tmp2md6);
+ dst=100;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmpd6);
+ if(tmpd6[4]<dst){
+ dst=tmpd6[4];
+ }
+ }
+ sx=0; sy=0;nall=0;
+ for(ii=1;ii<=length6(tmp1md6);ii++){
+ pull6(ii,tmp1md6,tmpd6);
+ if(tmpd6[4]<dst+Eps00){
+ nall++;
+ sx=sx+tmpd6[0]; sy=sy+tmpd6[1];
+ }
+ tmp1d[0]=sx/nall; tmp1d[1]=sy/nall;
+ tmp1d6[0]=tmp1d[0];tmp1d6[1]=tmp1d[1];
+ nearestptpt(tmp1d,crv1,ans);
+ tmp1d6[2]=ans[2];
+ nearestptpt(tmp1d,crv2,ans);
+ tmp1d6[3]=ans[2];
+ tmp1d6[4]=dst;
+ tmp1d6[5]=tmpd6[5];
+ addptd6(out,tmp1d6);
+ }
+ tmp1md6[0][0]=0;
+ for(ii=1;ii<=length6(tmp2md6);ii++){
+ pull6(ii,tmp2md6,tmpd6);addptd6(tmp1md6,tmpd6);
+ }
+ }
+}
+
+void intersectcurves(double crv1[][2],double crv2[][2],int num,double out[][2]){
+ double out6[DsizeS][6],tmpd[2],tmpd6[6];
+ int jj;
+ intersectcurvesPp(crv1,crv2,num,out6);
+ out[0][0]=0;
+ for(jj=1;jj<=length6(out6);jj++){
+ pull6(jj,out6,tmpd6);push2(tmpd6[0],tmpd6[1],out);
+ }
+}