asymptote 1.72 sources (not integrated into build yet)

git-svn-id: svn://tug.org/texlive/trunk@13110 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 221 insertions, 0 deletions

diff --git a/Build/source/utils/asymptote/simpson.cc b/Build/source/utils/asymptote/simpson.cc
new file mode 100644
index 00000000000..c17f3545c03
--- /dev/null
+++ b/Build/source/utils/asymptote/simpson.cc
@@ -0,0 +1,221 @@
+#include <cmath>
+#include <cassert>
+#include <cfloat>
+
+// Compute a numerical approximation to an integral via adaptive Simpson's Rule
+// This routine ignores underflow.
+
+const int nest=DBL_MANT_DIG;
+
+typedef struct {
+  bool left;                    // left interval?
+  double psum, f1t, f2t, f3t, dat, estr;
+} TABLE;
+
+bool                            // Returns true iff successful.
+simpson(double& integral,       // Approximate value of the integral.
+        double (*f)(double),    // Pointer to function to be integrated.
+        double a, double b,     // Lower, upper limits of integration (a <= b).
+        double acc,             // Desired relative accuracy of integral.
+                                // Try to make |error| <= acc*abs(integral).
+        double dxmax)           // Maximum limit on the width of a subinterval
+// For periodic functions, dxmax should be
+// set to the period or smaller to prevent
+// premature convergence of Simpson's rule. 
+{
+  double diff,area,estl,estr,alpha,da,dx,wt,est,fv[5];
+  TABLE table[nest],*p,*pstop;
+  static const double sixth=1.0/6.0;
+        
+  bool success=true;
+  p=table;
+  pstop=table+nest-1;
+  p->left=true;
+  p->psum=0.0;
+  alpha=a;
+  da=b-a;
+  fv[0]=(*f)(alpha);
+  fv[2]=(*f)(alpha+0.5*da);
+  fv[4]=(*f)(alpha+da);
+  wt=sixth*da;
+  est=wt*(fv[0]+4.0*fv[2]+fv[4]);
+  area=est;
+
+  // Have estimate est of integral on (alpha, alpha+da).
+  // Bisect and compute estimates on left and right half intervals.
+  // integral is the best value for the integral.
+
+  for(;;) {
+    dx=0.5*da;
+    double arg=alpha+0.5*dx;
+    fv[1]=(*f)(arg);
+    fv[3]=(*f)(arg+dx);
+    wt=sixth*dx;
+    estl=wt*(fv[0]+4.0*fv[1]+fv[2]);
+    estr=wt*(fv[2]+4.0*fv[3]+fv[4]);
+    integral=estl+estr;
+    diff=est-integral;
+    area -= diff;
+
+    if(p >= pstop) success=false;
+    if(!success || (fabs(diff) <= acc*fabs(area) && da <= dxmax)) {
+      // Accept approximate integral.
+      // If it was a right interval, add results to finish at this level.
+      // If it was a left interval, process right interval.
+
+      for(;;) {
+        if(p->left == false) { // process right-half interval
+          alpha += da;
+          p->left=true;
+          p->psum=integral;
+          fv[0]=p->f1t;
+          fv[2]=p->f2t;
+          fv[4]=p->f3t;
+          da=p->dat;
+          est=p->estr;
+          break;
+        }
+        integral += p->psum;
+        if(--p <= table) return success;
+      }
+
+    } else {
+      // Raise level and store information for processing right-half interval.
+      ++p;
+      da=dx;
+      est=estl;
+      p->left=false;
+      p->f1t=fv[2];
+      p->f2t=fv[3];
+      p->f3t=fv[4];
+      p->dat=dx;
+      p->estr=estr;
+      fv[4]=fv[2];
+      fv[2]=fv[1];
+    }
+  }
+}
+
+// Use adaptive Simpson integration to determine the upper limit of
+// integration required to make the definite integral of a continuous
+// non-negative function close to a user specified sum.
+// This routine ignores underflow.
+
+bool                            // Returns true iff successful.
+unsimpson(double integral,      // Given value for the integral.
+          double (*f)(double),  // Pointer to function to be integrated.
+          double a, double& b,  // Lower, upper limits of integration (a <= b).
+                                // The value of b provided on entry is used
+                                // as an initial guess; somewhat faster if the
+                                // given value is an underestimation.
+          double acc,           // Desired relative accuracy of b.
+                                // Try to make |integral-area| <= acc*integral.
+          double& area,         // Computed integral of f(x) on [a,b].
+          double dxmax,         // Maximum limit on the width of a subinterval
+                                // For periodic functions, dxmax should be
+                                // set to the period or smaller to prevent
+                                // premature convergence of Simpson's rule. 
+          double dxmin=0)       // Lower limit on sampling width.
+{
+  double diff,estl,estr,alpha,da,dx,wt,est,fv[5];
+  double sum,parea,pdiff,b2;
+  TABLE table[nest],*p,*pstop;
+  static const double sixth=1.0/6.0;
+        
+  p=table;
+  pstop=table+nest-1;
+  p->psum=0.0;
+  alpha=a;
+  parea=0.0;
+  pdiff=0.0;
+  
+  for(;;) {
+    p->left=true;
+    da=b-alpha;
+    fv[0]=(*f)(alpha);
+    fv[2]=(*f)(alpha+0.5*da);
+    fv[4]=(*f)(alpha+da);
+    wt=sixth*da;
+    est=wt*(fv[0]+4.0*fv[2]+fv[4]);
+    area=est;
+
+    // Have estimate est of integral on (alpha, alpha+da).
+    // Bisect and compute estimates on left and right half intervals.
+    // Sum is better value for integral.
+
+    bool cont=true;
+    while(cont) {
+      dx=0.5*da;
+      double arg=alpha+0.5*dx;
+      fv[1]=(*f)(arg);
+      fv[3]=(*f)(arg+dx);
+      wt=sixth*dx;
+      estl=wt*(fv[0]+4.0*fv[1]+fv[2]);
+      estr=wt*(fv[2]+4.0*fv[3]+fv[4]);
+      sum=estl+estr;
+      diff=est-sum;
+
+      assert(sum >= 0.0);
+      area=parea+sum;
+      b2=alpha+da;
+      if(fabs(fabs(integral-area)-fabs(pdiff))+fabs(diff) <= fv[4]*acc*(b2-a)){
+        b=b2;
+        return true;
+      }
+      if(fabs(integral-area) > fabs(pdiff+diff)) {
+        if(integral <= area) {
+          p=table;
+          p->left=true;
+          p->psum=parea;
+        } else {
+          if((fabs(diff) <= fv[4]*acc*da || dx <= dxmin) && da <= dxmax) {
+            // Accept approximate integral sum.
+            // If it was a right interval, add results to finish at this level.
+            // If it was a left interval, process right interval.
+      
+            pdiff += diff;
+            for(;;) {
+              if(p->left == false) { // process right-half interval
+                parea += sum;
+                alpha += da;
+                p->left=true;
+                p->psum=sum;
+                fv[0]=p->f1t;
+                fv[2]=p->f2t;
+                fv[4]=p->f3t;
+                da=p->dat;
+                est=p->estr;
+                break;
+              }
+              sum += p->psum;
+              parea -= p->psum;
+              if(--p <= table) {
+                p=table;
+                p->psum=parea=sum;
+                alpha += da;
+                b += b-a;
+                cont=false;
+                break;
+              }
+            }
+            continue;
+          }
+        }
+      }
+      if(p >= pstop) return false;
+// Raise level and store information for processing right-half interval.
+      ++p;
+      da=dx;
+      est=estl;
+      p->psum=0.0;
+      p->left=false;
+      p->f1t=fv[2];
+      p->f2t=fv[3];
+      p->f3t=fv[4];
+      p->dat=dx;
+      p->estr=estr;
+      fv[4]=fv[2];
+      fv[2]=fv[1];
+    }
+  }
+}