"Fossies" - the Fresh Open Source Software Archive

Member "asymptote-2.61/runmath.in" (18 Nov 2019, 8900 Bytes) of package /linux/misc/asymptote-2.61.src.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "runmath.in": 2.53_vs_2.54.

    1 /*****
    2  * runmath.in
    3  *
    4  * Runtime functions for math operations.
    5  *
    6  *****/
    7 
    8 pair     => primPair()
    9 realarray* => realArray()
   10 pairarray* => pairArray()
   11 
   12 #include <inttypes.h>
   13 
   14 #include "mathop.h"
   15 #include "path.h"
   16 
   17 #ifdef __CYGWIN__
   18 extern "C" double yn(int, double);
   19 extern "C" double jn(int, double);
   20 extern "C" int __signgam;
   21 #define signgam __signgam
   22 #endif
   23 
   24 using namespace camp;
   25 
   26 typedef array realarray;
   27 typedef array pairarray;
   28 
   29 using types::realArray;
   30 using types::pairArray;
   31 
   32 using run::integeroverflow;
   33 using vm::frame;
   34 
   35 const char *invalidargument="invalid argument";
   36 
   37 extern uint32_t CLZ(uint32_t a);
   38 
   39 inline unsigned intbits() {
   40   static unsigned count=0;
   41   if(count > 0) return count;
   42   while((1ULL << count) < Int_MAX)
   43     ++count;
   44   ++count;
   45   return count;
   46 }
   47 
   48 static const unsigned char BitReverseTable8[256]=
   49 {
   50 #define R2(n)     n,    n+2*64,    n+1*64,    n+3*64
   51 #define R4(n) R2(n),R2(n+2*16),R2(n+1*16),R2(n+3*16)
   52 #define R6(n) R4(n),R4(n+2*4 ),R4(n+1*4 ),R4(n+3*4 )
   53 R6(0),R6(2),R6(1),R6(3)
   54 };
   55 #undef R2
   56 #undef R4
   57 #undef R6
   58 
   59 unsigned long long bitreverse8(unsigned long long a)
   60 {
   61   return
   62     (unsigned long long) BitReverseTable8[a]; 
   63 }
   64 
   65 unsigned long long bitreverse16(unsigned long long a)
   66 {
   67   return
   68     ((unsigned long long) BitReverseTable8[a & 0xff] << 8) | 
   69     ((unsigned long long) BitReverseTable8[(a >> 8)]);
   70 }
   71 
   72 unsigned long long bitreverse24(unsigned long long a)
   73 {
   74   return
   75     ((unsigned long long) BitReverseTable8[a & 0xff] << 16) | 
   76     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 8) | 
   77     ((unsigned long long) BitReverseTable8[(a >> 16)]);
   78 }
   79 
   80 unsigned long long bitreverse32(unsigned long long a)
   81 {
   82   return
   83     ((unsigned long long) BitReverseTable8[a & 0xff] << 24) | 
   84     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 16) | 
   85     ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 8) |
   86     ((unsigned long long) BitReverseTable8[(a >> 24)]);
   87 }
   88 
   89 unsigned long long bitreverse40(unsigned long long a)
   90 {
   91   return
   92     ((unsigned long long) BitReverseTable8[a & 0xff] << 32) | 
   93     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 24) | 
   94     ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 16) |
   95     ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 8) |
   96     ((unsigned long long) BitReverseTable8[(a >> 32)]);
   97 }
   98 
   99 unsigned long long bitreverse48(unsigned long long a)
  100 {
  101   return
  102     ((unsigned long long) BitReverseTable8[a & 0xff] << 40) | 
  103     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 32) | 
  104     ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 24) |
  105     ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 16) |
  106     ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 8) |
  107     ((unsigned long long) BitReverseTable8[(a >> 40)]);
  108 }
  109 
  110 unsigned long long bitreverse56(unsigned long long a)
  111 {
  112   return
  113     ((unsigned long long) BitReverseTable8[a & 0xff] << 48) | 
  114     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 40) | 
  115     ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 32) |
  116     ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 24) |
  117     ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 16) |
  118     ((unsigned long long) BitReverseTable8[(a >> 40) & 0xff] << 8) | 
  119     ((unsigned long long) BitReverseTable8[(a >> 48)]);
  120 }
  121 
  122 unsigned long long bitreverse64(unsigned long long a)
  123 {
  124   return
  125     ((unsigned long long) BitReverseTable8[a & 0xff] << 56) | 
  126     ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 48) | 
  127     ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 40) |
  128     ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 32) |
  129     ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 24) |
  130     ((unsigned long long) BitReverseTable8[(a >> 40) & 0xff] << 16) | 
  131     ((unsigned long long) BitReverseTable8[(a >> 48) & 0xff] << 8) |
  132     ((unsigned long long) BitReverseTable8[(a >> 56)]);
  133 }
  134 
  135 // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
  136 #define T unsignedInt
  137 Int popcount(T a)
  138 {
  139   a=a-((a >> 1) & (T)~(T)0/3);
  140   a=(a & (T)~(T)0/15*3)+((a >> 2) & (T)~(T)0/15*3);
  141   a=(a+(a >> 4)) & (T)~(T)0/255*15;
  142 return (T)(a*((T)~(T)0/255)) >> (sizeof(T)-1)*CHAR_BIT;
  143 }
  144 #undef T
  145 
  146 // Return the factorial of a non-negative integer using a lookup table.
  147 Int factorial(Int n)
  148 {
  149   static Int *table;
  150   static Int size=0;
  151   if(size == 0) {
  152     Int f=1;
  153     size=2;
  154     while(f <= Int_MAX/size)
  155       f *= (size++);
  156     table=new Int[size];
  157     table[0]=f=1;
  158     for(Int i=1; i < size; ++i) {
  159       f *= i;
  160       table[i]=f;
  161     }
  162   }
  163   if(n >= size) integeroverflow(0);
  164   return table[n];
  165 }
  166 
  167 static inline Int Round(double x) 
  168 {
  169   return Int(x+((x >= 0) ? 0.5 : -0.5));
  170 }
  171 
  172 inline Int sgn(double x) 
  173 {
  174   return (x > 0.0 ? 1 : (x < 0.0 ? -1 : 0));
  175 }
  176 
  177 static bool initializeRandom=true;
  178 
  179 void Srand(Int seed)
  180 { 
  181   initializeRandom=false;
  182   const int n=256;
  183   static char state[n];
  184   initstate(intcast(seed),state,n);
  185 }  
  186 
  187 // Autogenerated routines:
  188 
  189 
  190 real ^(real x, Int y)
  191 {
  192   return pow(x,y);
  193 }
  194 
  195 pair ^(pair z, Int y)
  196 {
  197   return pow(z,y);
  198 }
  199 
  200 Int quotient(Int x, Int y)
  201 { 
  202   return quotient<Int>()(x,y);
  203 }  
  204 
  205 Int abs(Int x)
  206 { 
  207   return Abs(x);
  208 }  
  209 
  210 Int sgn(real x)
  211 { 
  212   return sgn(x);
  213 }  
  214 
  215 Int rand()
  216 { 
  217   if(initializeRandom)
  218     Srand(1);
  219   return random();
  220 }  
  221 
  222 void srand(Int seed)
  223 { 
  224   Srand(seed);
  225 }  
  226 
  227 // a random number uniformly distributed in the interval [0,1]
  228 real unitrand()
  229 {                         
  230   return ((real) random())/RANDOM_MAX;
  231 }
  232 
  233 Int ceil(real x)
  234 { 
  235   return Intcast(ceil(x));
  236 }
  237 
  238 Int floor(real x)
  239 { 
  240   return Intcast(floor(x));
  241 }
  242 
  243 Int round(real x)
  244 { 
  245   if(validInt(x)) return Round(x);
  246   integeroverflow(0);
  247 }
  248 
  249 Int Ceil(real x)
  250 { 
  251   return Ceil(x);
  252 }
  253 
  254 Int Floor(real x)
  255 { 
  256   return Floor(x);
  257 }
  258 
  259 Int Round(real x)
  260 { 
  261   return Round(Intcap(x));
  262 }
  263 
  264 real fmod(real x, real y)
  265 {
  266   if (y == 0.0) dividebyzero();
  267   return fmod(x,y);
  268 }
  269 
  270 real atan2(real y, real x)
  271 { 
  272   return atan2(y,x);
  273 }  
  274 
  275 real hypot(real x, real y)
  276 { 
  277   return hypot(x,y);
  278 }  
  279 
  280 real remainder(real x, real y)
  281 { 
  282   return remainder(x,y);
  283 }  
  284 
  285 real Jn(Int n, real x)
  286 {
  287   return jn(n,x);
  288 }
  289 
  290 real Yn(Int n, real x)
  291 {
  292   return yn(n,x);
  293 }
  294 
  295 real erf(real x)
  296 {
  297   return erf(x);
  298 }
  299 
  300 real erfc(real x)
  301 {
  302   return erfc(x);
  303 }
  304 
  305 Int factorial(Int n) {
  306   if(n < 0) error(invalidargument);
  307   return factorial(n);
  308 }
  309 
  310 Int choose(Int n, Int k) {
  311   if(n < 0 || k < 0 || k > n) error(invalidargument);
  312   Int f=1;
  313   Int r=n-k;
  314   for(Int i=n; i > r; --i) {
  315     if(f > Int_MAX/i) integeroverflow(0);
  316     f=(f*i)/(n-i+1);
  317   }
  318   return f;
  319 }
  320 
  321 real gamma(real x)
  322 {
  323 #ifdef HAVE_TGAMMA
  324   return tgamma(x);
  325 #else
  326   real lg = lgamma(x);
  327   return signgam*exp(lg);
  328 #endif
  329 }
  330 
  331 realarray *quadraticroots(real a, real b, real c)
  332 {
  333   quadraticroots q(a,b,c);
  334   array *roots=new array(q.roots);
  335   if(q.roots >= 1) (*roots)[0]=q.t1;
  336   if(q.roots == 2) (*roots)[1]=q.t2;
  337   return roots;
  338 }
  339 
  340 pairarray *quadraticroots(explicit pair a, explicit pair b, explicit pair c)
  341 {
  342   Quadraticroots q(a,b,c);
  343   array *roots=new array(q.roots);
  344   if(q.roots >= 1) (*roots)[0]=q.z1;
  345   if(q.roots == 2) (*roots)[1]=q.z2;
  346   return roots;
  347 }
  348 
  349 realarray *cubicroots(real a, real b, real c, real d)
  350 {
  351   cubicroots q(a,b,c,d);
  352   array *roots=new array(q.roots);
  353   if(q.roots >= 1) (*roots)[0]=q.t1;
  354   if(q.roots >= 2) (*roots)[1]=q.t2;
  355   if(q.roots == 3) (*roots)[2]=q.t3;
  356   return roots;
  357 }
  358 
  359 
  360 // Logical operations
  361 
  362 bool !(bool b)
  363 {
  364   return !b;
  365 }
  366 
  367 bool :boolMemEq(frame *a, frame *b)
  368 {
  369   return a == b;
  370 }
  371 
  372 bool :boolMemNeq(frame *a, frame *b)
  373 {
  374   return a != b;
  375 }
  376 
  377 bool :boolFuncEq(callable *a, callable *b)
  378 {
  379   return a->compare(b);
  380 }
  381 
  382 bool :boolFuncNeq(callable *a, callable *b)
  383 {
  384   return !(a->compare(b));
  385 }
  386 
  387 
  388 // Bit operations
  389 
  390 Int AND(Int a, Int b) 
  391 {
  392   return a & b;
  393 }
  394 
  395 Int OR(Int a, Int b) 
  396 {
  397   return a | b;
  398 }
  399 
  400 Int XOR(Int a, Int b) 
  401 {
  402   return a ^ b;
  403 }
  404 
  405 Int NOT(Int a)
  406 {
  407   return ~a;
  408 }
  409 
  410 Int CLZ(Int a) 
  411 {
  412   if((unsigned long long) a > 0xFFFFFFFF)
  413     return CLZ((uint32_t) ((unsigned long long) a >> 32));
  414   else {
  415     int bits=intbits();
  416     if(a != 0) return bits-32+CLZ((uint32_t) a);
  417     return bits;
  418   }
  419 }
  420 
  421 Int popcount(Int a) 
  422 {
  423   return popcount(a);
  424 }
  425 
  426 Int CTZ(Int a) 
  427 {
  428   return popcount((a&-a)-1);
  429 }
  430 
  431 // bitreverse a within a word of length bits.
  432 Int bitreverse(Int a, Int bits)
  433 {
  434   typedef unsigned long long Bitreverse(unsigned long long a);
  435   static Bitreverse *B[]={bitreverse8,bitreverse16,bitreverse24,bitreverse32,
  436                      bitreverse40,bitreverse48,bitreverse56,bitreverse64};
  437   int maxbits=intbits()-1; // Drop sign bit
  438 #if Int_MAX2 >= 0x7fffffffffffffffLL
  439   --maxbits;               // Drop extra bit for reserved values
  440 #endif  
  441   if(bits <= 0 || bits > maxbits || a < 0 ||
  442      (unsigned long long) a >= (1ULL << bits))
  443     return -1;
  444   unsigned int bytes=(bits+7)/8;
  445   return B[bytes-1]((unsigned long long) a) >> (8*bytes-bits);
  446 }