math.js

js analyticphysics.com collected 2026-10-02 05:03:33 UTC 150,022 bytes, 6,211 lines download raw bytes

1
2var defaultDecimals = 18;
3
4var pi = Math.PI;
5
6var eulerGamma = .5772156649015329;
7
8var constants = {
9
10  decimals: 50, precisionScale: 10n**50n,
11
12  e: 271828182845904523536028747135266249775724709369995n,
13
14  eulerGamma: 57721566490153286060651209008240243104215933593992n,
15
16  ln10: 230258509299404568401799145468436420760110148862877n,
17
18  pi: 314159265358979323846264338327950288419716939937510n
19
20};
21
22function getConstant( name ) {
23
24  return constants[name] * precisionScale / constants.precisionScale;
25
26}
27
28
29// oeis.org/A000367
30bernoulli2nN = [ 1n, 1n, -1n, 1n, -1n, 5n, -691n, 7n, -3617n, 43867n, -174611n, 854513n, -236364091n, 8553103n, -23749461029n, 8615841276005n, -7709321041217n, 2577687858367n, -26315271553053477373n, 2929993913841559n, -261082718496449122051n, 1520097643918070802691n, -27833269579301024235023n, 596451111593912163277961n, -5609403368997817686249127547n, 495057205241079648212477525n, -801165718135489957347924991853n, 29149963634884862421418123812691n, -2479392929313226753685415739663229n, 84483613348880041862046775994036021n, -1215233140483755572040304994079820246041491n, 12300585434086858541953039857403386151n, -106783830147866529886385444979142647942017n, 1472600022126335654051619428551932342241899101n, -78773130858718728141909149208474606244347001n, 1505381347333367003803076567377857208511438160235n, -5827954961669944110438277244641067365282488301844260429n, 34152417289221168014330073731472635186688307783087n, -24655088825935372707687196040585199904365267828865801n, 414846365575400828295179035549542073492199375372400483487n, -4603784299479457646935574969019046849794257872751288919656867n, 1677014149185145836823154509786269900207736027570253414881613n, -2024576195935290360231131160111731009989917391198090877281083932477n, 660714619417678653573847847426261496277830686653388931761996983n, -1311426488674017507995511424019311843345750275572028644296919890574047n, 1179057279021082799884123351249215083775254949669647116231545215727922535n, -1295585948207537527989427828538576749659341483719435143023316326829946247n, 1220813806579744469607301679413201203958508415202696621436215105284649447n, -211600449597266513097597728109824233673043954389060234150638733420050668349987259n, 67908260672905495624051117546403605607342195728504487509073961249992947058239n, -94598037819122125295227433069493721872702841533066936133385696204311395415197247711n ];
31
32// oeis.org/A002445
33bernoulli2nD = [ 1n, 6n, 30n, 42n, 30n, 66n, 2730n, 6n, 510n, 798n, 330n, 138n, 2730n, 6n, 870n, 14322n, 510n, 6n, 1919190n, 6n, 13530n, 1806n, 690n, 282n, 46410n, 66n, 1590n, 798n, 870n, 354n, 56786730n, 6n, 510n, 64722n, 30n, 4686n, 140100870n, 6n, 30n, 3318n, 230010n, 498n, 3404310n, 6n, 61410n, 272118n, 1410n, 6n, 4501770n, 6n, 33330n ];
34
35
36var factorialCache = [ 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800 ];
37
38
39function complex( x, y=0 ) {
40
41  if ( y === 0 && isArbitrary(x) ) y = 0n;
42
43  return { re: x, im: y };
44
45}
46
47var C = complex;
48
49function isComplex( x ) { return typeof x === 'object' && 're' in x; }
50
51
52var decimals = [], precisionScale;
53var arb1, arb2, onePi, twoPi, halfPi, ln10;
54
55function setPrecisionScale( n ) {
56
57  if ( n === 'reset' ) decimals.shift();
58  else decimals.unshift( n );
59
60  precisionScale = 10n**BigInt( decimals[0] );
61
62  // set some commonly used constants
63  arb1 = arbitrary(1);
64  arb2 = 2n * arb1;
65  onePi = getConstant( 'pi' );
66  twoPi = 2n * onePi;
67  halfPi = onePi / 2n;
68  ln10 = getConstant( 'ln10' );
69
70}
71
72function resetPrecisionScale() { setPrecisionScale( 'reset' ); }
73
74setPrecisionScale( defaultDecimals );
75
76function arbitrary( x ) {
77
78  if ( isComplex(x) ) return { re: arbitrary( x.re ), im: arbitrary( x.im ) };
79
80  if ( isArbitrary(x) ) return Number(x) / 10**decimals[0];
81
82  // BigInt from exponential form includes wrong digits
83  // manual construction from string more accurate
84
85  var parts = x.toExponential().split( 'e' );
86  var mantissa = parts[0].replace( '.', '' );
87  var digits = mantissa.length - ( mantissa[0] === '-' ? 2 : 1 )
88  var padding = +parts[1] + decimals[0] - digits;
89
90  if ( padding < 0 ) return BigInt( Math.round( x * 10**decimals[0] ) );
91
92  return BigInt( mantissa + '0'.repeat(padding) );
93
94}
95
96var A = arbitrary;
97
98function isArbitrary( x ) { return typeof x === 'bigint' || typeof x.re === 'bigint'; }
99
100
101function isZero( x ) {
102
103  if ( isComplex(x) ) return x.re === 0 && x.im === 0;
104  return x === 0;
105
106}
107
108function isUnity( x ) {
109
110  if ( isComplex(x) ) return x.re === 1 && x.im === 0;
111  return x === 1;
112
113}
114
115function isReal( x ) {
116
117  if ( isComplex(x) ) return x.im === 0;
118  return true;
119
120}
121
122function isInteger( x ) {
123
124  if ( isComplex(x) ) return Number.isInteger(x.re) && x.im === 0;
125  return Number.isInteger(x);
126
127}
128
129function isPositiveInteger( x ) {
130
131  if ( isComplex(x) ) return Number.isInteger(x.re) && x.re > 0 && x.im === 0;
132  return Number.isInteger(x) && x > 0;

133
134}
135
136function isPositiveIntegerOrZero( x ) {
137
138  if ( isComplex(x) ) return Number.isInteger(x.re) && x.re >= 0 && x.im === 0;
139  return Number.isInteger(x) && x >= 0;
140
141}
142
143function isNegativeInteger( x ) {
144
145  if ( isComplex(x) ) return Number.isInteger(x.re) && x.re < 0 && x.im === 0;
146  return Number.isInteger(x) && x < 0;
147
148}
149
150function isNegativeIntegerOrZero( x ) {
151
152  if ( isComplex(x) ) return Number.isInteger(x.re) && x.re <= 0 && x.im === 0;
153  return Number.isInteger(x) && x <= 0;
154
155}
156
157function isEqualTo( x, y ) {
158
159  if ( isComplex(x) || isComplex(y) ) {
160    if ( !isComplex(x) ) x = complex(x);
161    if ( !isComplex(y) ) y = complex(y);
162    return x.re === y.re && x.im === y.im;
163  }
164
165  return x === y;
166
167}
168
169
170function re( x ) {
171
172  if ( isComplex(x) ) return x.re;
173  return x;
174
175}
176
177var real = re;
178
179function im( x ) {
180
181  if ( isComplex(x) ) return x.im;
182  return 0;
183
184}
185
186var imag = im;
187
188function abs( x ) {
189
190  if ( isComplex(x) ) {
191
192    if ( isArbitrary(x) ) {
193
194      if ( x.im === 0n ) return abs(x.re);
195
196      if ( x.re === 0n ) return abs(x.im);
197
198      if ( abs(x.re) < abs(x.im) )
199
200        return mul( abs(x.im), sqrt( arb1 + div( div( mul(x.re,x.re), x.im ), x.im ) ) );
201
202      else
203
204         return mul( abs(x.re), sqrt( arb1 + div( div( mul(x.im,x.im), x.re ), x.re ) ) );
205
206    }
207
208    if ( x.im === 0 ) return Math.abs(x.re);
209
210    if ( x.re === 0 ) return Math.abs(x.im);
211
212    if ( Math.abs(x.re) < Math.abs(x.im) )
213
214      return Math.abs(x.im) * Math.sqrt( 1 + ( x.re / x.im )**2 );
215
216    else
217
218      return Math.abs(x.re) * Math.sqrt( 1 + ( x.im / x.re )**2 );
219
220  }
221
222  if ( isArbitrary(x) )
223    if ( x < 0n ) return -x;
224    else return x;
225
226  return Math.abs(x);
227
228}
229
230function arg( x ) {
231
232  if ( isArbitrary(x) ) return ln(x).im;
233
234  // adding zero prevents unexpected behavior for -0
235
236  if ( isComplex(x) ) return Math.atan2( x.im + 0, x.re + 0 );
237
238  return Math.atan2( 0, x + 0 );
239
240}
241
242
243// JavaScript does not support operator overloading
244
245function add( x, y ) {
246
247  if ( arguments.length > 2 ) {
248
249    var z = add( x, y );
250    for ( var i = 2 ; i < arguments.length ; i++ ) z = add( z, arguments[i] );
251    return z; 
252
253  }
254
255  if ( isComplex(x) || isComplex(y) ) {
256
257    if ( !isComplex(x) ) x = complex(x);
258    if ( !isComplex(y) ) y = complex(y);
259
260    return { re: x.re + y.re, im: x.im + y.im };
261
262  }
263
264  return x + y;
265
266}
267
268function sub( x, y ) {
269
270  if ( isComplex(x) || isComplex(y) ) {
271
272    if ( !isComplex(x) ) x = complex(x);
273    if ( !isComplex(y) ) y = complex(y);
274
275    return { re: x.re - y.re, im: x.im - y.im };
276
277  }
278
279  return x - y;
280
281}
282
283function mul( x, y ) {
284
285  if ( arguments.length > 2 ) {
286
287    var z = mul( x, y );
288    for ( var i = 2 ; i < arguments.length ; i++ ) z = mul( z, arguments[i] );
289    return z; 
290
291  }
292
293  if ( isComplex(x) || isComplex(y) ) {
294
295    if ( !isComplex(x) ) x = complex(x);
296    if ( !isComplex(y) ) y = complex(y);
297
298    if ( isArbitrary(x) )
299
300      return { re: ( x.re * y.re - x.im * y.im ) / precisionScale,
301               im: ( x.im * y.re + x.re * y.im ) / precisionScale };
302
303    return { re: x.re * y.re - x.im * y.im,
304             im: x.im * y.re + x.re * y.im };
305
306  }
307
308  if ( isArbitrary(x) ) return x * y / precisionScale;
309
310  return x * y;
311
312}
313
314function neg( x ) {
315
316  if ( isComplex(x) ) return { re: -x.re, im: -x.im };
317
318  return -x;
319
320}
321
322function div( x, y ) {
323
324  if ( isComplex(x) || isComplex(y) ) {
325
326    if ( !isComplex(x) ) x = complex(x);
327    if ( !isComplex(y) ) y = complex(y);
328
329    if ( y.re === 0 && y.im === 0 || y.re === 0n && y.im === 0n ) // operator precedence
330      throw Error( 'Division by zero' );
331
332    if ( isArbitrary(x) ) {
333
334      var N = { re: x.re * y.re + x.im * y.im,
335                im: x.im * y.re - x.re * y.im };
336      var D = y.re * y.re + y.im * y.im;
337
338      return { re: precisionScale * N.re / D,
339               im: precisionScale * N.im / D };
340
341    }
342
343    if ( Math.abs(y.re) < Math.abs(y.im) ) {
344
345      var f = y.re / y.im;
346      return { re: ( x.re * f + x.im ) / ( y.re * f + y.im ),
347               im: ( x.im * f - x.re ) / ( y.re * f + y.im ) };
348
349    } else {
350
351      var f = y.im / y.re;
352      return { re: ( x.re + x.im * f ) / ( y.re + y.im * f ),
353               im: ( x.im - x.re * f ) / ( y.re + y.im * f ) };
354
355    }
356
357  }
358
359  if ( y === 0 || y === 0n ) throw Error( 'Division by zero' );
360
361  if ( isArbitrary(x) ) return precisionScale * x / y;
362
363  return x / y;
364
365}
366
367function inv( x ) {
368
369  if ( isArbitrary(x) ) return div( arb1, x );
370
371  return div( 1, x );
372
373}
374
375function pow( x, y ) {
376
377  if ( isArbitrary(x) || isArbitrary(y) ) {
378
379    if ( !isArbitrary(x) ) x = arbitrary(x);
380    if ( !isArbitrary(y) ) y = arbitrary(y);
381
382    return exp( mul( y, ln(x) ) );
383
384  }
385
386  if ( isComplex(x) || isComplex(y) ) {
387
388    if ( !isComplex(x) ) x = complex(x);
389    if ( !isComplex(y) ) y = complex(y);
390
391    if ( x.re === 0 && x.im === 0 && y.re > 0 )
392      return complex(0);
393    if ( x.re === 0 && x.im === 0 && y.re === 0 && y.im === 0 )
394      return complex(1);
395    if ( x.re === 0 && x.im === 0 && y.re < 0 )
396      throw Error( 'Power singularity' );
397
398    return exp( mul( y, log(x) ) );
399
400  }
401
402  if ( x === 0 && y < 0 ) throw Error( 'Power singularity' );
403
404  if ( x < 0 && !Number.isInteger(y) ) return pow( complex(x), y );
405
406  return x**y;
407
408}
409
410function root( x, y ) { return pow( x, div( 1, y ) ); }
411
412function surd( x, n ) {
413
414  if ( isComplex(x) || isComplex(n) ) throw Error( 'Surd requires real inputs' );
415
416  if ( !isInteger(n) ) throw Error( 'Second parameter of surd must be integer' );
417
418  if ( n & 1 ) {
419    var sign = Math.sign(x); // zero at origin anyway
420    return sign * root( sign*x, n );
421  }
422
423  if ( x < 0 ) throw Error( 'First parameter of surd must be positive for even integers' );
424

425  return root( x, n );
426
427}
428
429function sqrt( x ) {
430
431  if ( isComplex(x) ) {
432
433    if ( isArbitrary(x) ) {
434
435      if ( x.im === 0n )
436        if ( x.re < 0n ) return { re: 0n, im: sqrt(-x.re) };
437        else return { re: sqrt(x.re), im: 0n };
438
439      // need evaluation independent of natural logarithm
440
441      var c = abs(x);
442      var sign = x.im < 0n ? -1n : 1n;
443
444      return { re: sqrt( div( c + x.re, arb2 ) ), im: sign * sqrt( div( c - x.re, arb2 ) ) }
445
446    }
447
448    if ( x.im === 0 )
449      if ( x.re < 0 ) return { re: 0, im: Math.sqrt(-x.re) };
450      else return { re: Math.sqrt(x.re), im: 0 };
451
452    // expression above suffers from rounding errors when not arbitrary,
453    //   especially affecting elliptic integral of third kind
454
455    return exp( mul( .5, log(x) ) );
456
457  }
458
459  if ( isArbitrary(x) ) {
460
461    if ( x === 0n ) return 0n;
462
463    if ( x < 0n ) throw Error( 'Cannot evaluate real square root of ' + x );
464
465    // Brent, Modern Computer Arithmetic, SqrtInt algorithm
466
467    var u = x, s, t;
468
469    while ( u !== s ) {
470      s = u;
471      t = s + div( x, s );
472      u = div( t, arb2 );
473    }
474
475    return s;
476
477  }
478
479  if ( x < 0 ) return { re: 0, im: Math.sqrt(-x) };
480
481  return Math.sqrt(x);
482
483}
484
485function fibonacci( n, real=false ) {
486
487  var sqrt5 = sqrt(5), r = ( 1 + sqrt5 )/2;
488
489  if ( real ) return div( sub( pow(r,n), mul( pow(1/r,n), cos(mul(n,pi)) ) ), sqrt5 );
490
491  return div( sub( pow(r,n), pow(-1/r,n) ), sqrt5 );
492
493}
494
495
496function complexAverage( f, x, offset=1e-5 ) {
497
498  return div( add( f(add(x,offset)), f(sub(x,offset)) ), 2 );
499
500}
501
502
503function complexFromString( s, returnAsString=false ) {
504
505  var lead = '', real, imag;
506
507  if ( s[0] === '+' || s[0] === '-' ) {
508    lead = s[0];
509    s = s.slice(1);
510  }
511
512  if ( s.includes('+') || s.includes('-') ) {
513    if ( s.includes('+') ) {
514      real = lead + s.slice( 0, s.indexOf('+') );
515      imag = s.slice( s.indexOf('+') + 1, s.length - 1 );
516    } else {
517      real = lead + s.slice( 0, s.indexOf('-') );
518      imag = s.slice( s.indexOf('-'), s.length - 1 );
519    }
520  } else {
521    if ( s.includes('i') ) {
522      real = '0';
523      imag = lead + s.slice( 0, s.length - 1 );
524    } else {
525      real = lead + s;
526      imag = '0';
527    }
528  }
529
530  if ( imag === '' || imag === '-' ) imag += '1';
531
532  if ( returnAsString ) return `{ re: ${real}, im: ${imag} }`;
533
534  return { re: +real, im: +imag };
535
536}
537
538function output() {
539
540  for ( var i = 0 ; i < arguments.length ; i++ ) console.log( arguments[i] );
541
542}
543
544
545function besselJ( n, x ) {
546
547  if ( isComplex(n) || isComplex(x) ) {
548
549    if ( isNegativeInteger(n) ) return mul( pow(-1,n), besselJ( neg(n), x ) );
550
551    var product = div( pow( div(x,2), n ), gamma( add(n,1) ) );
552    return mul( product, hypergeometric0F1( add(n,1), mul(-.25, pow(x,2) ) ) );
553
554  } 
555
556  if ( isNegativeInteger(n) ) return (-1)**n * besselJ( -n, x );
557
558  if ( !Number.isInteger(n) && x < 0 ) return besselJ( n, complex(x) );
559
560  return (x/2)**n * hypergeometric0F1( n+1, -.25*x**2 ) / gamma(n+1);
561
562}
563
564function besselJZero( n, m, derivative=false ) {
565
566  if ( isComplex(n) || isComplex(m) )
567    throw Error( 'Complex arguments not supported in Bessel zero' );
568
569  if ( !isPositiveInteger(m) ) throw Error( 'Unsupported index for Bessel zero' );
570
571  // approximations from dlmf.nist.gov/10.21.vi
572  var delta = pi/4;
573
574  if ( derivative ) {
575
576    if ( n === 0 && m === 1 ) return 0;
577
578    var b = ( m + Math.abs(n)/2 - 3/4 ) * pi;
579    var e = b - ( 4*n**2 + 3 ) / ( 8*b );
580
581    // Zeros not simple to enumerate for negative order
582    // Odd double zero for small range of n to be investigated
583    if ( n < 0 ) throw Error( 'Negative order for Bessel derivative zero not supported' );
584
585    return findRoot( x => diff( x => besselJ(n,x), x ), [ e-delta < 0 ? 0 : e-delta, e+delta ] );
586
587  } else {
588
589    var a = ( m + Math.abs(n)/2 - 1/4 ) * pi;
590    var e = a - ( 4*n**2 - 1 ) / ( 8*a );
591
592    if ( n < 0 ) e += fractionalPart(n) * pi;
593
594    return findRoot( x => besselJ(n,x), [ m === 1 ? 1e-10 : e-delta, e+delta ] );
595
596  }
597
598}
599
600function besselY( n, x ) {
601
602  if ( isComplex(n) || isComplex(x) ) {
603
604    if ( isZero(x) ) {
605      if ( !isComplex(n) ) n = complex(n);
606      if ( fractionalPart(n.re) === -.5 && n.im === 0 ) return complex(0);
607    }
608
609    // dlmf.nist.gov/10.2.4
610    if ( isInteger(n) )
611      return div( add( diff( n => besselJ(n,x), n ),
612                       mul( pow(-1,n), diff( n => besselJ(n,x), neg(n) ) ) ), pi );
613
614    var s = sub( mul( besselJ(n,x), cos( mul(n,pi) ) ), besselJ( neg(n), x ) );
615    return div( s, sin( mul(n,pi) ) );
616
617  }
618
619  if ( x < 0 ) return besselY( n, complex(x) );
620
621  if ( x === 0 )
622    if ( fractionalPart(n) === -.5 ) return 0;
623
624  // dlmf.nist.gov/10.2.4
625  if ( Number.isInteger(n) )
626    return ( diff( n => besselJ(n,x), n ) + (-1)**n * diff( n => besselJ(n,x), -n ) ) / pi;
627
628  return ( besselJ(n,x) * cos(n*pi) - besselJ(-n,x) ) / sin(n*pi);
629
630}
631
632function besselYZero( n, m, derivative=false ) {
633
634  if ( isComplex(n) || isComplex(m) )
635    throw Error( 'Complex arguments not supported in Bessel zero' );
636
637  if ( !isPositiveInteger(m) ) throw Error( 'Unsupported index for Bessel zero' );
638
639  // approximations from dlmf.nist.gov/10.21.vi
640  var delta = pi/4;
641
642  if ( derivative ) {
643
644    var b = ( m + Math.abs(n)/2 - 1/4 ) * pi;
645    var e = b - ( 4*n**2 + 3 ) / ( 8*b );
646
647    // Zeros not simple to enumerate for negative order
648    // Odd double zero for small range of n to be investigated
649    if ( n < 0 ) throw Error( 'Negative order for Bessel derivative zero not supported' );
650
651    return findRoot( x => diff( x => besselY(n,x), x ), [ e-delta, e+delta ] );
652
653  } else {
654
655    var a = ( m + Math.abs(n)/2 - 3/4 ) * pi;
656    var e = a - ( 4*n**2 - 1 ) / ( 8*a );
657
658    if ( n < 0 )
659      if ( fractionalPart(n) > -.5 ) e += fractionalPart(n) * pi;
660      else e += ( 1 + fractionalPart(n) ) * pi;
661
662    return findRoot( x => besselY(n,x), [ m === 1 && fractionalPart(n) > -.5 ? 1e-10 : e-delta, e+delta ] );
663
664  }
665
666}
667
668function besselI( n, x ) {
669
670  if ( isComplex(n) || isComplex(x) ) {
671
672    if ( isNegativeInteger(n) ) return besselI( neg(n), x );
673

674    var product = div( pow( div(x,2), n ), gamma( add(n,1) ) );
675    return mul( product, hypergeometric0F1( add(n,1), mul(.25, pow(x,2) ) ) );
676
677  }
678
679  if ( isNegativeInteger(n) ) return besselI( -n, x );
680
681  if ( !Number.isInteger(n) && x < 0 ) return besselI( n, complex(x) );
682
683  return (x/2)**n * hypergeometric0F1( n+1, .25*x**2 ) / gamma(n+1);
684
685}
686
687function besselK( n, x ) {
688
689  var useAsymptotic = 10;
690
691  if ( isComplex(n) || isComplex(x) ) {
692
693    // asymptotic form as per Johansson arxiv.org/abs/1606.06977
694    if ( abs(x) > useAsymptotic ) {
695
696      var t1 = mul( sqrt( div( pi/2, x ) ), exp( neg(x) ) );
697      var t2 = hypergeometric2F0( add(n,.5), sub(.5,n), div(-.5,x) );
698
699      return mul( t1, t2 );
700
701    }
702
703    // dlmf.nist.gov/10.27.5
704    if ( isInteger(n) )
705      return mul( pow(-1,add(n,1)), .5,
706                  add( diff( n => besselI(n,x), n ), diff( n => besselI(n,x), neg(n) ) ) );
707
708    var product = div( pi/2, sin( mul(n,pi) ) );
709    return mul( product, sub( besselI( neg(n), x ), besselI(n,x) ) );
710
711  }
712
713  if ( x > useAsymptotic )
714    return sqrt(pi/2/x) * exp(-x) * hypergeometric2F0( n+.5, .5-n, -.5/x );
715
716  if ( x < 0 ) return besselK( n, complex(x) );
717
718  // dlmf.nist.gov/10.27.5
719  if ( Number.isInteger(n) )
720    return (-1)**(n+1)/2 * ( diff( n => besselI(n,x), n ) + diff( n => besselI(n,x), -n ) );
721
722  return pi/2 * ( besselI(-n,x) - besselI(n,x) ) / sin(n*pi);
723
724}
725
726function hankel1( n, x, tolerance=1e-10 ) {
727
728  var useAsymptotic = 11;
729
730  // dlmf.nist.gov/10.17
731  if ( x.im > 0 && abs(x) > useAsymptotic ) {
732
733    var w = add( x, mul(n,-pi/2), -pi/4 );
734
735    var s = complex(1);
736    var plusI = complex(0,1);
737    var p = complex(1);
738    var i = 1;
739
740    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
741      p = mul( p, plusI, sub( mul(4,n,n), (2*i-1)**2 ), 1/i, 1/8, inv(x) );
742      s = add( s, p );
743      i++;
744    }
745
746    return mul( div( Math.sqrt(2/pi), sqrt(x) ), exp( mul(plusI,w) ), s );
747
748  }
749
750  // not duplicate special case
751  if ( isInteger(n) )
752    return add( besselJ(n,x), mul( complex(0,1), besselY(n,x) ) );
753
754  // one less call to besselJ
755  var s = sub( mul( besselJ(n,x), exp( mul(n,complex(0,-pi)) ) ), besselJ( neg(n), x ) );
756  return div( mul( complex(0,1), s ), sin( mul(n,pi) ) );
757
758}
759
760function hankel2( n, x, tolerance=1e-10 ) {
761
762  var useAsymptotic = 11;
763
764  // dlmf.nist.gov/10.17
765  if ( x.im < 0 && abs(x) > useAsymptotic ) {
766
767    var w = add( x, mul(n,-pi/2), -pi/4 );
768
769    var s = complex(1);
770    var minusI = complex(0,-1);
771    var p = complex(1);
772    var i = 1;
773
774    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
775      p = mul( p, minusI, sub( mul(4,n,n), (2*i-1)**2 ), 1/i, 1/8, inv(x) );
776      s = add( s, p );
777      i++;
778    }
779
780    return mul( div( Math.sqrt(2/pi), sqrt(x) ), exp( mul(minusI,w) ), s );
781
782  }
783
784  // not duplicate special case
785  if ( isInteger(n) )
786    return sub( besselJ(n,x), mul( complex(0,1), besselY(n,x) ) );
787
788  // one less call to besselJ
789  var s = sub( mul( besselJ(n,x), exp( mul(n,complex(0,pi)) ) ), besselJ( neg(n), x ) );
790  return div( mul( complex(0,-1), s ), sin( mul(n,pi) ) );
791
792}
793
794
795// dlmf.nist.gov/9.2.ii and dlmf.nist.gov/9.6.i
796
797function airyAi( x ) {
798
799  if ( isComplex(x) ) {
800
801    if ( isZero(x) ) return complex( 1 / 3**(2/3) / gamma(2/3) );
802
803    if ( x.re < 0 ) {
804
805      var z = mul( 2/3, pow( neg(x), 3/2 ) );
806      return mul( 1/3, sqrt(neg(x)), add( besselJ( 1/3, z ), besselJ( -1/3, z ) ) );
807
808    }
809
810    var z = mul( 2/3, pow( x, 3/2 ) );
811    return mul( 1/pi, sqrt( div( x, 3 ) ), besselK( 1/3, z ) );
812
813  }
814
815  if ( x === 0 ) return 1 / 3**(2/3) / gamma(2/3);
816
817  if ( x < 0 ) {
818
819    var z = 2/3 * (-x)**(3/2);
820    return sqrt(-x) / 3 * ( besselJ( 1/3, z ) + besselJ( -1/3, z ) );
821
822  }
823
824  var z = 2/3 * x**(3/2);
825  return 1/pi * sqrt(x/3) * besselK( 1/3, z );
826
827}
828
829function airyAiPrime( x ) {
830
831  if ( isComplex(x) ) {
832
833    if ( isZero(x) ) return complex( -1 / 3**(1/3) / gamma(1/3) );
834
835    if ( x.re < 0 ) {
836
837      var z = mul( 2/3, pow( neg(x), 3/2 ) );
838      return mul( 1/3, x, sub( besselJ( -2/3, z ), besselJ( 2/3, z ) ) );
839
840    }
841
842    var z = mul( 2/3, pow( x, 3/2 ) );
843    return mul( -1/pi/sqrt(3), x, besselK( 2/3, z ) );
844
845  }
846
847  if ( x === 0 ) return -1 / 3**(1/3) / gamma(1/3);
848
849  if ( x < 0 ) {
850
851    var z = 2/3 * (-x)**(3/2);
852    return x/3 * ( besselJ( -2/3, z ) - besselJ( 2/3, z ) );
853
854  }
855
856  var z = 2/3 * x**(3/2);
857  return -1/pi/sqrt(3) * x * besselK( 2/3, z );
858
859}
860
861function airyBi( x ) {
862
863  if ( isComplex(x) ) {
864
865    if ( isZero(x) ) return complex( 1 / 3**(1/6) / gamma(2/3) );
866
867    if ( x.re < 0 ) {
868
869      var z = mul( 2/3, pow( neg(x), 3/2 ) );
870      return mul( sqrt( div(neg(x),3) ), sub( besselJ( -1/3, z ), besselJ( 1/3, z ) ) );
871
872    }
873
874    var z = mul( 2/3, pow( x, 3/2 ) );
875    return mul( sqrt( div( x, 3 ) ), add( besselI( 1/3, z ), besselI( -1/3, z ) ) );
876
877  }
878
879  if ( x === 0 ) return 1 / 3**(1/6) / gamma(2/3);
880
881  if ( x < 0 ) {
882
883    var z = 2/3 * (-x)**(3/2);
884    return sqrt(-x/3) * ( besselJ( -1/3, z ) - besselJ( 1/3, z ) );
885
886  }
887
888  var z = 2/3 * x**(3/2);

889  return sqrt(x/3) * ( besselI( 1/3, z ) + besselI( -1/3, z ) );
890
891}
892
893function airyBiPrime( x ) {
894
895  if ( isComplex(x) ) {
896
897    if ( isZero(x) ) return complex( 3**(1/6) / gamma(1/3) );
898
899    if ( x.re < 0 ) {
900
901      var z = mul( 2/3, pow( neg(x), 3/2 ) );
902      return mul( 1/sqrt(3), neg(x), add( besselJ( 2/3, z ), besselJ( -2/3, z ) ) );
903
904    }
905
906    var z = mul( 2/3, pow( x, 3/2 ) );
907    return mul( 1/sqrt(3), x, add( besselI( 2/3, z ), besselI( -2/3, z ) ) );
908
909  }
910
911  if ( x === 0 ) return 3**(1/6) / gamma(1/3);
912
913  if ( x < 0 ) {
914
915    var z = 2/3 * (-x)**(3/2);
916    return -x/sqrt(3) * ( besselJ( 2/3, z ) + besselJ( -2/3, z ) );
917
918  }
919
920  var z = 2/3 * x**(3/2);
921  return x/sqrt(3) * ( besselI( 2/3, z ) + besselI( -2/3, z ) );
922
923}
924
925
926function sphericalBesselJ( n, x ) {
927
928  return mul( div( Math.sqrt(pi/2), sqrt(x) ), besselJ( add( n, .5 ), x ) );
929
930}
931
932function sphericalBesselY( n, x ) {
933
934  return mul( div( Math.sqrt(pi/2), sqrt(x) ), besselY( add( n, .5 ), x ) );
935
936}
937
938function sphericalHankel1( n, x ) {
939
940  return mul( div( Math.sqrt(pi/2), sqrt(x) ), hankel1( add( n, .5 ), x ) );
941
942}
943
944function sphericalHankel2( n, x ) {
945
946  return mul( div( Math.sqrt(pi/2), sqrt(x) ), hankel2( add( n, .5 ), x ) );
947
948}
949
950
951function struveH( n, x ) {
952
953  return mul( pow( x, add(n,1) ), inv( mul( pow(2,n), sqrt(pi), gamma( add(n,3/2) ) ) ),
954              hypergeometric1F2( 1, 3/2, add(n,3/2), mul( -1/4, pow(x,2) ) ) );
955
956}
957
958function struveL( n, x ) {
959
960  // one sign different from struveH
961
962  return mul( pow( x, add(n,1) ), inv( mul( pow(2,n), sqrt(pi), gamma( add(n,3/2) ) ) ),
963              hypergeometric1F2( 1, 3/2, add(n,3/2), mul( 1/4, pow(x,2) ) ) );
964
965}
966
967
968function jacobiTheta( n, x, q, tolerance=1e-10 ) {
969
970  if ( abs(q) >= 1 ) throw Error( 'Unsupported elliptic nome' );
971
972  if ( ![1,2,3,4].includes(n) ) throw Error( 'Undefined Jacobi theta index' );
973
974  if ( isComplex(x) || isComplex(q) ) {
975
976    if ( !isComplex(x) ) x = complex(x);
977
978    var piTau = div( log(q), complex(0,1) );
979
980    // dlmf.nist.gov/20.2 to reduce overflow
981    if ( Math.abs(x.im) > Math.abs(piTau.im) || Math.abs(x.re) > Math.PI ) {
982
983      // use floor for consistency with fundamentalParallelogram
984      var pt = Math.floor( x.im / piTau.im );
985      x = sub( x, mul( pt, piTau ) );
986
987      var p = Math.floor( x.re / Math.PI );
988      x = sub( x, p * Math.PI );
989
990      var qFactor = pow( q, -pt*pt );
991      var eFactor = exp( mul( -2 * pt, x, complex(0,1) ) );
992
993      // factors can become huge, so chop spurious parts first
994      switch( n ) {
995
996        case 1:
997
998          return mul( (-1)**(p+pt), qFactor, eFactor, chop( jacobiTheta( n, x, q ), tolerance ) );
999
1000        case 2:
1001
1002          return mul( (-1)**p, qFactor, eFactor, chop( jacobiTheta( n, x, q ), tolerance ) );
1003
1004        case 3:
1005
1006          return mul( qFactor, eFactor, chop( jacobiTheta( n, x, q ), tolerance ) );
1007
1008        case 4:
1009
1010          return mul( (-1)**pt, qFactor, eFactor, chop( jacobiTheta( n, x, q ), tolerance ) );
1011
1012      }
1013
1014    }
1015
1016    switch( n ) {
1017
1018      case 1:
1019
1020        var s = complex(0);
1021        var p = complex(1);
1022        var i = 0;
1023
1024        while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
1025          p = mul( (-1)**i, pow( q, i*i+i ), sin( mul(2*i+1,x) ) );
1026          s = add( s, p );
1027          i++;
1028        }
1029
1030        return mul( 2, pow( q, 1/4 ), s );
1031
1032      case 2:
1033
1034        var s = complex(0);
1035        var p = complex(1);
1036        var i = 0;
1037
1038        while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
1039          p = mul( pow( q, i*i+i ), cos( mul(2*i+1,x) ) );
1040          s = add( s, p );
1041          i++;
1042        }
1043
1044        return mul( 2, pow( q, 1/4 ), s );
1045
1046      case 3:
1047
1048        var s = complex(0);
1049        var p = complex(1);
1050        var i = 1;
1051
1052        while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
1053          p = mul( pow( q, i*i ), cos( mul(2*i,x) ) );
1054          s = add( s, p );
1055          i++;
1056        }
1057
1058        return add( 1, mul(2,s) );
1059
1060      case 4:
1061
1062        var s = complex(0);
1063        var p = complex(1);
1064        var i = 1;
1065
1066        while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
1067          p = mul( pow( neg(q), i*i ), cos( mul(2*i,x) ) );
1068          s = add( s, p );
1069          i++;
1070        }
1071
1072        return add( 1, mul(2,s) );
1073
1074      }
1075
1076  } else {
1077
1078    switch( n ) {
1079
1080      case 1:
1081
1082        if ( q < 0 ) return jacobiTheta( n, x, complex(q) );
1083
1084        var s = 0;
1085        var p = 1;
1086        var i = 0;
1087
1088        while ( Math.abs(p) > tolerance ) {
1089          p = (-1)**i * q**(i*i+i) * sin( (2*i+1) * x );
1090          s += p;
1091          i++;
1092        }
1093
1094        return 2 * q**(1/4) * s;
1095
1096      case 2:
1097
1098        if ( q < 0 ) return jacobiTheta( n, x, complex(q) );
1099
1100        var s = 0;
1101        var p = 1;
1102        var i = 0;
1103
1104        while ( Math.abs(p) > tolerance ) {
1105          p = q**(i*i+i) * cos( (2*i+1) * x );
1106          s += p;
1107          i++;
1108        }
1109
1110        return 2 * q**(1/4) * s;
1111
1112      case 3:
1113
1114        var s = 0;
1115        var p = 1;
1116        var i = 1;
1117
1118        while ( Math.abs(p) > tolerance ) {
1119          p = q**(i*i) * cos( 2*i * x );
1120          s += p;
1121          i++;
1122        }
1123
1124        return 1 + 2 * s;
1125
1126      case 4:
1127
1128        var s = 0;
1129        var p = 1;
1130        var i = 1;
1131
1132        while ( Math.abs(p) > tolerance ) {
1133          p = (-q)**(i*i) * cos( 2*i * x );
1134          s += p;
1135          i++;
1136        }
1137
1138        return 1 + 2 * s;
1139
1140    }
1141
1142  }
1143
1144}
1145
1146
1147function ellipticNome( m ) {
1148
1149  if ( typeof m === 'undefined' ) throw Error( 'Elliptic parameter is undefined' );
1150
1151  if ( isComplex(m) ) return exp( div( mul( -pi, ellipticK( sub(1,m) ) ), ellipticK(m) ) );
1152
1153  if ( m > 1 ) return ellipticNome( complex(m) );
1154
1155  if ( m < 0 ) return -exp( -pi * ellipticK( 1/(1-m) ) / ellipticK( m/(m-1) ) );
1156
1157  return exp( -pi * ellipticK(1-m) / ellipticK(m) );
1158
1159}
1160
1161function am( x, m ) {
1162
1163  if ( m > 1 || isComplex(x) || isComplex(m) ) {
1164
1165    if ( !isComplex(x) ) x = complex(x);
1166    if ( !isComplex(m) ) m = complex(m);
1167
1168    if ( m.im === 0 && m.re <= 1 ) {
1169
1170      var K = ellipticK( m.re );
1171      var n = Math.round( x.re / 2 / K );
1172      x = sub( x, 2 * n * K );
1173
1174      if ( m.re < 0 ) {
1175
1176        var Kp = ellipticK( 1 - m.re );
1177        var p = Math.round( x.im / 2 / Kp.re );
1178
1179        // bitwise test for odd integer
1180        if ( p & 1 ) return sub( n * pi, arcsin( sn(x,m) ) );
1181
1182      }
1183
1184      return add( arcsin( sn(x,m) ), n * pi );
1185
1186    }
1187
1188    return arcsin( sn(x,m) );
1189
1190  } else {
1191
1192    var K = ellipticK(m);
1193    var n = Math.round( x / 2 / K );
1194    x = x - 2 * n * K;
1195
1196    return Math.asin( sn(x,m) ) + n * pi;
1197
1198  }
1199
1200}
1201
1202// Jacobi theta functions grow very quickly in imaginary direction
1203// evaluating elliptic functions without reduction fails accordingly
1204
1205function fundamentalParallelogram( x, p1, p2 ) {
1206
1207  // x = m p1 + n p2, solve for m, n
1208
1209  var det = p1.re * p2.im - p1.im * p2.re;
1210
1211  var m = ( x.re * p2.im - x.im * p2.re ) / det;
1212  var n = ( x.im * p1.re - x.re * p1.im ) / det;
1213
1214  return add( x, mul( -Math.floor(m), p1 ), mul( -Math.floor(n), p2 ) );
1215
1216}
1217
1218function sn( x, m ) {
1219
1220  if ( m > 1 || isComplex(x) || isComplex(m) ) {
1221
1222    if ( !isComplex(m) ) m = complex(m); // ensure K complex
1223
1224    // dlmf.nist.gov/22.17
1225    if ( abs(m) > 1 ) return mul( inv(sqrt(m)), sn( mul(sqrt(m),x), inv(m) ) ); 
1226
1227    // periods 4K, 2iK'
1228    var K = ellipticK(m), p1 = mul( 4, K );
1229    var p2 = mul( complex(0,2), ellipticK( sub(1,m) ) );
1230
1231    x = fundamentalParallelogram( x, p1, p2 );
1232
1233    var q = ellipticNome(m);
1234    var t = div( x, mul( 2/pi, K ) );  // theta(3,0,q)**2 = 2K(m)/pi
1235
1236    return mul( div( jacobiTheta(3,0,q), jacobiTheta(2,0,q) ),
1237                div( jacobiTheta(1,t,q), jacobiTheta(4,t,q) ) );
1238
1239  }
1240
1241  // dlmf.nist.gov/22.5.ii
1242  if ( m === 0 ) return sin(x);

1243  if ( m === 1 ) return tanh(x);
1244
1245  var q = ellipticNome(m);
1246  var t = x / jacobiTheta(3,0,q)**2;
1247
1248  if ( m < 0 )
1249    return jacobiTheta(3,0,q) / jacobiTheta(4,t,q)
1250           * div( jacobiTheta(1,t,q), jacobiTheta(2,0,q) ).re;
1251
1252  return jacobiTheta(3,0,q) / jacobiTheta(2,0,q)
1253         * jacobiTheta(1,t,q) / jacobiTheta(4,t,q);
1254
1255}
1256
1257function cn( x, m ) {
1258
1259  if ( m > 1 || isComplex(x) || isComplex(m) ) {
1260
1261    if ( !isComplex(m) ) m = complex(m); // ensure K complex
1262
1263    // dlmf.nist.gov/22.17
1264    if ( abs(m) > 1 ) return dn( mul(sqrt(m),x), inv(m) ); 
1265
1266    // periods 4K, 2K + 2iK'
1267    var K = ellipticK(m), p1 = mul( 4, K );
1268    var p2 = add( div(p1,2), mul( complex(0,2), ellipticK( sub(1,m) ) ) );
1269
1270    x = fundamentalParallelogram( x, p1, p2 );
1271
1272    var q = ellipticNome(m);
1273    var t = div( x, mul( 2/pi, K ) );  // theta(3,0,q)**2 = 2K(m)/pi
1274
1275    return mul( div( jacobiTheta(4,0,q), jacobiTheta(2,0,q) ),
1276                div( jacobiTheta(2,t,q), jacobiTheta(4,t,q) ) );
1277
1278  }
1279
1280  // dlmf.nist.gov/22.5.ii
1281  if ( m === 0 ) return cos(x);
1282  if ( m === 1 ) return sech(x);
1283
1284  var q = ellipticNome(m);
1285  var t = x / jacobiTheta(3,0,q)**2;
1286
1287  if ( m < 0 )
1288    return jacobiTheta(4,0,q) / jacobiTheta(4,t,q)
1289           * div( jacobiTheta(2,t,q), jacobiTheta(2,0,q) ).re;
1290
1291  return jacobiTheta(4,0,q) / jacobiTheta(2,0,q)
1292         * jacobiTheta(2,t,q) / jacobiTheta(4,t,q);
1293
1294}
1295
1296function dn( x, m ) {
1297
1298  if ( m > 1 || isComplex(x) || isComplex(m) ) {
1299
1300    if ( !isComplex(m) ) m = complex(m); // ensure K complex
1301
1302    // dlmf.nist.gov/22.17
1303    if ( abs(m) > 1 ) return cn( mul(sqrt(m),x), inv(m) ); 
1304
1305    // periods 2K, 4iK'
1306    var K = ellipticK(m), p1 = mul( 2, K );
1307    var p2 = mul( complex(0,4), ellipticK( sub(1,m) ) );
1308
1309    x = fundamentalParallelogram( x, p1, p2 );
1310
1311    var q = ellipticNome(m);
1312    var t = div( x, mul( 2/pi, K ) );  // theta(3,0,q)**2 = 2K(m)/pi
1313
1314    return mul( div( jacobiTheta(4,0,q), jacobiTheta(3,0,q) ),
1315                div( jacobiTheta(3,t,q), jacobiTheta(4,t,q) ) );
1316
1317  }
1318
1319  // dlmf.nist.gov/22.5.ii
1320  if ( m === 0 ) return 1;
1321  if ( m === 1 ) return sech(x);
1322
1323  var q = ellipticNome(m);
1324  var t = x / jacobiTheta(3,0,q)**2;
1325
1326  return jacobiTheta(4,0,q) / jacobiTheta(3,0,q)
1327         * jacobiTheta(3,t,q) / jacobiTheta(4,t,q);
1328
1329}
1330
1331function ns( x, m ) { return inv( sn(x,m) ); }
1332
1333function nc( x, m ) { return inv( cn(x,m) ); }
1334
1335function nd( x, m ) { return inv( dn(x,m) ); }
1336
1337function sc( x, m ) { return div( sn(x,m), cn(x,m) ); }
1338
1339function cs( x, m ) { return div( cn(x,m), sn(x,m) ); }
1340
1341function sd( x, m ) { return div( sn(x,m), dn(x,m) ); }
1342
1343function ds( x, m ) { return div( dn(x,m), sn(x,m) ); }
1344
1345function cd( x, m ) { return div( cn(x,m), dn(x,m) ); }
1346
1347function dc( x, m ) { return div( dn(x,m), cn(x,m) ); }
1348
1349function lemniscateSin( x ) { return sn(x,-1); }
1350
1351function lemniscateCos( x ) { return dn(x,2); }
1352
1353
1354function weierstrassRoots( g2, g3 ) {
1355
1356  function cubicTrigSolution( p, q, n ) {
1357
1358    // p, q both negative in defining cubic
1359
1360    if ( isZero(p) ) return mul( root(q,3), exp( complex(0,2*pi*n/3) ) )
1361
1362    return mul( 2/sqrt(3), sqrt(p),
1363                cos( sub( div( arccos( mul( 3*sqrt(3)/2, q, pow(p,-3/2) ) ), 3 ),
1364                          2*pi*n/3 ) ) );
1365  }
1366
1367  g2 = div( g2, 4 );
1368  g3 = div( g3, 4 );
1369
1370  var e1 = cubicTrigSolution( g2, g3, 0 );
1371  var e2 = cubicTrigSolution( g2, g3, 1 );
1372  var e3 = cubicTrigSolution( g2, g3, 2 );
1373
1374  return [ e1, e2, e3 ];
1375
1376}
1377
1378function weierstrassHalfPeriods( g2, g3 ) {
1379
1380  // Davis, Intro to Nonlinear Diff. & Integral Eqs., pp.157-8
1381  // consistent with periods of Jacobi sine in weierstrassP
1382  // not consistent with Mathematica
1383
1384  var [ e1, e2, e3 ] = weierstrassRoots( g2, g3 );
1385
1386  var lambda = sqrt( sub(e1,e3) );
1387  var m = div( sub(e2,e3), sub(e1,e3) );
1388
1389  var w1 = div( ellipticK(m), lambda );
1390  var w2 = div( mul( complex(0,1), ellipticK( sub(1,m) ) ), lambda );
1391
1392  return [ w1, w2 ];
1393
1394}
1395
1396function weierstrassInvariants( w1, w2 ) {
1397
1398  if ( !isComplex(w1) ) w1 = complex(w1);
1399  if ( !isComplex(w2) ) w2 = complex(w2);
1400

1401  // order half periods by complex slope
1402  if ( w2.im/w2.re < w1.im/w1.re ) [ w1, w2 ] = [ w2, w1 ];
1403
1404  var ratio =  div( w2, w1 ), conjugate;
1405
1406  if ( ratio.im < 0 ) {
1407    ratio.im = -ratio.im;
1408    conjugate = true;
1409  }
1410
1411  var q = exp( mul( complex(0,1), pi, ratio ) );
1412
1413  // en.wikipedia.org/wiki/Weierstrass's_elliptic_functions
1414  // modified for input of half periods
1415
1416  var a = jacobiTheta( 2, 0, q );
1417  var b = jacobiTheta( 3, 0, q );
1418
1419  var g2 = mul( 4/3*pi**4, pow( mul(2,w1), -4 ),
1420                add( pow(a,8), mul( -1, pow(a,4), pow(b,4) ), pow(b,8) ) );
1421
1422  var g3 = mul( 8/27*pi**6, pow( mul(2,w1), -6 ),
1423                add( pow(a,12), mul( -3/2, pow(a,8), pow(b,4) ),
1424                                mul( -3/2, pow(a,4), pow(b,8) ), pow(b,12) ) );
1425
1426  if ( conjugate ) {
1427    g2.im = -g2.im;
1428    g3.im = -g3.im;
1429  }
1430
1431  return [ g2, g3 ];
1432
1433}
1434
1435
1436function weierstrassP( x, g2, g3 ) {
1437
1438  if ( !isComplex(x) ) x = complex(x);
1439
1440  if ( isZero(g2) && isZero(g3) ) return pow( x, -2 );
1441
1442  var [ e1, e2, e3 ] = weierstrassRoots( g2, g3 );
1443
1444  // Whittaker & Watson, Section 22.351
1445
1446  var m = div( sub(e2,e3), sub(e1,e3) );
1447
1448  return add( e3, mul( sub(e1,e3), pow( sn( mul( x, sqrt(sub(e1,e3)) ), m ), -2 ) ) );
1449
1450}
1451
1452function weierstrassPPrime( x, g2, g3 ) {
1453
1454  if ( !isComplex(x) ) x = complex(x);
1455
1456  if ( isZero(g2) && isZero(g3) ) return mul( -2, pow( x, -3 ) );
1457
1458  var [ e1, e2, e3 ] = weierstrassRoots( g2, g3 );
1459
1460  // Whittaker & Watson, Section 22.351
1461
1462  var m = div( sub(e2,e3), sub(e1,e3) );
1463
1464  var argument = mul( x, sqrt(sub(e1,e3)) );
1465
1466  return mul( -2, pow( sub(e1,e3), 3/2 ), cn( argument, m ), dn( argument, m ),
1467              pow( sn( argument, m ), -3 ) );
1468
1469}
1470
1471function inverseWeierstrassP( x, g2, g3 ) {
1472
1473  if ( !isComplex(x) ) x = complex(x);
1474
1475  var [ e1, e2, e3 ] = weierstrassRoots( g2, g3 );
1476
1477  // Johansson arxiv.org/pdf/1806.06725.pdf p.17
1478  // sign of imaginary part on real axis differs from Mathematica
1479
1480  return carlsonRF( sub(x,e1), sub(x,e2), sub(x,e3) );
1481
1482}
1483
1484
1485function kleinJ( x ) {
1486
1487  // from mpmath / elliptic.py
1488
1489  var q = exp( mul( complex(0,pi), x ) );
1490  var t2 = chop( jacobiTheta(2,0,q) );
1491  var t3 = chop( jacobiTheta(3,0,q) );
1492  var t4 = chop( jacobiTheta(4,0,q) );
1493  var P = pow( add( pow(t2,8), pow(t3,8), pow(t4,8) ), 3 );
1494  var Q = mul( 54, pow( mul(t2,t3,t4), 8 ) );
1495
1496  return div( P, Q );
1497
1498}
1499
1500
1501// Carlson symmetric integrals
1502
1503function carlsonRC( x, y ) {
1504
1505  if ( x < 0 || y < 0 || isComplex(x) || isComplex(y) ) {
1506
1507    if ( !isComplex(x) ) x = complex(x);
1508    if ( !isComplex(y) ) y = complex(y);
1509
1510    if ( x.re === y.re && x.im === y.im ) return inv( sqrt(x) );
1511
1512    // return value by continuity
1513
1514    return div( arccos( div( sqrt(x), sqrt(y) ) ), mul( sqrt(y), sqrt( sub( 1, div(x,y) ) ) ) );
1515
1516  }
1517
1518  if ( x === y ) return 1 / Math.sqrt(x);
1519
1520  if ( x < y ) return Math.acos( Math.sqrt(x/y) ) / Math.sqrt(y-x);
1521
1522  return Math.acosh( Math.sqrt(x/y) ) / Math.sqrt(x-y);
1523
1524}
1525
1526function carlsonRD( x, y, z ) {
1527
1528  return carlsonRJ( x, y, z, z );
1529
1530}
1531
1532function carlsonRF( x, y, z, tolerance=1e-10 ) {
1533
1534  if ( isComplex(x) || isComplex(y) || isComplex(z) ) {
1535
1536    var xm = x;
1537    var ym = y;
1538    var zm = z;
1539
1540    var Am = A0 = div( add( x, y, z ), 3 );
1541    var Q = Math.pow( 3*tolerance, -1/6 )
1542            * Math.max( abs( sub(A0,x) ), abs( sub(A0,y) ), abs( sub(A0,z) ) );
1543    var g = .25;
1544    var pow4 = 1;
1545
1546    while ( true ) {
1547      var xs = sqrt(xm);
1548      var ys = sqrt(ym);
1549      var zs = sqrt(zm);
1550      var lm = add( mul(xs,ys), mul(xs,zs), mul(ys,zs) );
1551      var Am1 = mul( add(Am,lm), g );
1552      xm = mul( add(xm,lm), g );
1553      ym = mul( add(ym,lm), g );
1554      zm = mul( add(zm,lm), g );
1555      if ( pow4 * Q < abs(Am) ) break;
1556      Am = Am1;
1557      pow4 *= g;
1558    }
1559
1560    var t = div( pow4, Am );
1561    var X = mul( sub(A0,x), t );
1562    var Y = mul( sub(A0,y), t );
1563    var Z = neg( add(X,Y) );
1564    var E2 = sub( mul(X,Y), mul(Z,Z) );
1565    var E3 = mul(X,Y,Z);
1566
1567    return mul( pow( Am, -.5 ),
1568             add( 9240, mul(-924,E2), mul(385,E2,E2), mul(660,E3), mul(-630,E2,E3) ), 1/9240 );
1569
1570  } else {
1571
1572    if ( y === z ) return carlsonRC( x, y );
1573    if ( x === z ) return carlsonRC( y, x );
1574    if ( x === y ) return carlsonRC( z, x );
1575
1576    if ( x < 0 || y < 0 || z < 0 ) return carlsonRF( complex(x), y, z );
1577
1578    // adapted from mpmath / elliptic.py
1579
1580    var xm = x;
1581    var ym = y;
1582    var zm = z;
1583
1584    var Am = A0 = (x + y + z) / 3;
1585    var Q = Math.pow( 3*tolerance, -1/6 )
1586            * Math.max( Math.abs(A0-x), Math.abs(A0-y), Math.abs(A0-z) );
1587    var g = .25;
1588    var pow4 = 1;
1589
1590    while ( true ) {
1591      var xs = Math.sqrt(xm);
1592      var ys = Math.sqrt(ym);
1593      var zs = Math.sqrt(zm);
1594      var lm = xs*ys + xs*zs + ys*zs;
1595      var Am1 = (Am + lm) * g;
1596      xm = (xm + lm) * g;
1597      ym = (ym + lm) * g;
1598      zm = (zm + lm) * g;
1599      if ( pow4 * Q < Math.abs(Am) ) break;
1600      Am = Am1;
1601      pow4 *= g;
1602    }
1603
1604    var t = pow4 / Am;
1605    var X = (A0-x) * t;
1606    var Y = (A0-y) * t;
1607    var Z = -X-Y;
1608    var E2 = X*Y - Z**2;
1609    var E3 = X*Y*Z;
1610
1611    return Math.pow( Am, -.5 )
1612           * ( 9240 - 924*E2 + 385*E2**2 + 660*E3 - 630*E2*E3 ) / 9240;
1613
1614  }
1615
1616}
1617
1618function carlsonRG( x, y, z ) {
1619
1620  var t1 = mul( z, carlsonRF(x,y,z) );
1621  var t2 = mul( -1/3, sub(x,z), sub(y,z), carlsonRD(x,y,z) );
1622  var t3 = sqrt( mul( x, y, inv(z) ) );
1623
1624  return mul( .5, add( t1, t2, t3 ) );
1625
1626}
1627
1628function carlsonRJ( x, y, z, p, tolerance=1e-10 ) {
1629
1630  if ( isComplex(x) || isComplex(y) || isComplex(z) || isComplex(p) ) {
1631
1632    var xm = x;
1633    var ym = y;
1634    var zm = z;
1635    var pm = p;
1636
1637    var A0 = Am = div( add( x, y, z, mul(2,p) ), 5 );
1638    var delta = mul( sub(p,x), sub(p,y), sub(p,z) );
1639    var Q = Math.pow( .25*tolerance, -1/6 )
1640            * Math.max( abs( sub(A0,x) ), abs( sub(A0,y) ), abs( sub(A0,z) ), abs( sub(A0,p) ) );
1641    var g = .25;
1642    var pow4 = 1;
1643    var S = complex(0);
1644
1645    while ( true ) {
1646      var sx = sqrt(xm);
1647      var sy = sqrt(ym);
1648      var sz = sqrt(zm);
1649      var sp = sqrt(pm);
1650      var lm = add( mul(sx,sy), mul(sx,sz), mul(sy,sz) );
1651      var Am1 = mul( add(Am,lm), g );
1652      xm = mul( add(xm,lm), g );
1653      ym = mul( add(ym,lm), g );
1654      zm = mul( add(zm,lm), g );
1655      pm = mul( add(pm,lm), g );
1656      var dm = mul( add(sp,sx), add(sp,sy), add(sp,sz) );
1657      var em = mul( delta, pow4**3, inv(dm), inv(dm) );
1658      if ( pow4 * Q < abs(Am) ) break;
1659      var T = mul( carlsonRC( 1, add(1,em) ), pow4, inv(dm) );
1660      S = add( S, T );
1661      pow4 *= g;
1662      Am = Am1;
1663    }
1664
1665    var t = div( pow4, Am );
1666    var X = mul( sub(A0,x), t );
1667    var Y = mul( sub(A0,y), t );
1668    var Z = mul( sub(A0,z), t );
1669    var P = div( add(X,Y,Z), -2 );
1670    var E2 = add( mul(X,Y), mul(X,Z), mul(Y,Z), mul(-3,P,P) );
1671    var E3 = add( mul(X,Y,Z), mul(2,E2,P), mul(4,P,P,P) );
1672    var E4 = mul( add( mul(2,X,Y,Z), mul(E2,P), mul(3,P,P,P) ), P );
1673    var E5 = mul(X,Y,Z,P,P);
1674    P = add( 24024, mul(-5148,E2), mul(2457,E2,E2), mul(4004,E3), mul(-4158,E2,E3), mul(-3276,E4), mul(2772,E5) );
1675    var v1 = mul( pow4, pow( Am, -1.5 ), P, 1/24024 );
1676    var v2 = mul(6,S);
1677
1678    return add( v1, v2 );
1679
1680  } else {
1681
1682    if ( x < 0 || y < 0 || z < 0 || p < 0 ) return carlsonRJ( complex(x), y, z, p );
1683
1684    // adapted from mpmath / elliptic.py
1685
1686    var xm = x;
1687    var ym = y;
1688    var zm = z;
1689    var pm = p;
1690

1691    var A0 = Am = (x + y + z + 2*p) / 5;
1692    var delta = (p-x) * (p-y) * (p-z);
1693    var Q = Math.pow( .25*tolerance, -1/6 )
1694            * Math.max( Math.abs(A0-x), Math.abs(A0-y), Math.abs(A0-z), Math.abs(A0-p) );
1695    var g = .25;
1696    var pow4 = 1;
1697    var S = 0;
1698
1699    while ( true ) {
1700      var sx = Math.sqrt(xm);
1701      var sy = Math.sqrt(ym);
1702      var sz = Math.sqrt(zm);
1703      var sp = Math.sqrt(pm);
1704      var lm = sx*sy + sx*sz + sy*sz;
1705      var Am1 = (Am + lm) * g;
1706      xm = (xm + lm) * g;
1707      ym = (ym + lm) * g;
1708      zm = (zm + lm) * g;
1709      pm = (pm + lm) * g;
1710      var dm = (sp+sx) * (sp+sy) * (sp+sz);
1711      var em = delta * pow4**3 / dm**2;
1712      if ( pow4 * Q < Math.abs(Am) ) break;
1713      var T = carlsonRC( 1, 1 + em ) * pow4 / dm;
1714      S += T;
1715      pow4 *= g;
1716      Am = Am1;
1717    }
1718
1719    var t = pow4 / Am;
1720    var X = (A0-x) * t;
1721    var Y = (A0-y) * t;
1722    var Z = (A0-z) * t;
1723    var P = (-X-Y-Z) / 2;
1724    var E2 = X*Y + X*Z + Y*Z - 3*P**2;
1725    var E3 = X*Y*Z + 2*E2*P + 4*P**3;
1726    var E4 = ( 2*X*Y*Z + E2*P + 3*P**3 ) * P;
1727    var E5 = X*Y*Z*P**2;
1728    P = 24024 - 5148*E2 + 2457*E2**2 + 4004*E3 - 4158*E2*E3 - 3276*E4 + 2772*E5;
1729    var v1 = pow4 * Math.pow( Am, -1.5 ) * P / 24024;
1730    var v2 = 6*S;
1731
1732    return v1 + v2;
1733
1734  }
1735
1736}
1737
1738
1739// elliptic integrals
1740
1741function ellipticF( x, m ) {
1742
1743  if ( arguments.length === 1 ) {
1744    m = x;
1745    x = pi/2;
1746  }
1747
1748  if ( isComplex(x) || isComplex(m) ) {
1749
1750    if ( !isComplex(x) ) x = complex(x);
1751
1752    var period = complex(0);
1753    if ( Math.abs(x.re) > pi/2 ) {
1754      var p = Math.round( x.re / pi );
1755      x.re = x.re - p * pi;
1756      period = mul( 2 * p, ellipticK( m ) );
1757    }
1758
1759    return add( mul( sin(x), carlsonRF( mul(cos(x),cos(x)), sub( 1, mul(m,sin(x),sin(x)) ), 1 ) ), period );
1760
1761  } else {
1762
1763    if ( m > 1 && Math.abs(x) > Math.asin( 1 / Math.sqrt(m) ) ) return ellipticF( complex(x), m );
1764
1765    var period = 0;
1766    if ( Math.abs(x) > pi/2 ) {
1767      var p = Math.round( x / pi );
1768      x = x - p * pi;
1769      period = 2 * p * ellipticK( m );
1770    }
1771
1772    return sin(x) * carlsonRF( cos(x)**2, 1 - m * sin(x)**2, 1 ) + period;
1773
1774  }
1775
1776}
1777
1778function ellipticK( m ) {
1779
1780  if ( typeof m === 'undefined' ) throw Error( 'Elliptic parameter is undefined' );
1781
1782  return ellipticF( m );
1783
1784}
1785
1786function ellipticE( x, m ) {
1787
1788  if ( arguments.length === 1 ) {
1789    m = x;
1790    x = pi/2;
1791  }
1792
1793  if ( isComplex(x) || isComplex(m) ) {
1794
1795    if ( !isComplex(x) ) x = complex(x);
1796
1797    var period = complex(0);
1798    if ( Math.abs(x.re) > pi/2 ) {
1799      var p = Math.round( x.re / pi );
1800      x.re = x.re - p * pi;
1801      period = mul( 2 * p,  ellipticE( m ) );
1802    }
1803
1804    return add( mul( sin(x), carlsonRF( mul(cos(x),cos(x)), sub( 1, mul(m,sin(x),sin(x)) ), 1 ) ),
1805                mul( -1/3, m, pow(sin(x),3), carlsonRD( mul(cos(x),cos(x)), sub( 1, mul(m,sin(x),sin(x)) ), 1 ) ),
1806                period );
1807
1808  } else {
1809
1810    if ( m > 1 && Math.abs(x) > Math.asin( 1 / Math.sqrt(m) ) ) return ellipticE( complex(x), m );
1811
1812    var period = 0;
1813    if ( Math.abs(x) > pi/2 ) {
1814      var p = Math.round( x / pi );
1815      x = x - p * pi;
1816      period = 2 * p * ellipticE( m );
1817    }
1818
1819    return sin(x) * carlsonRF( cos(x)**2, 1 - m * sin(x)**2, 1 )
1820           - m / 3 * sin(x)**3 * carlsonRD( cos(x)**2, 1 - m * sin(x)**2, 1 )
1821           + period;
1822
1823  }
1824
1825}
1826
1827function ellipticPi( n, x, m ) {
1828
1829  if ( arguments.length === 2 ) {
1830    m = x;
1831    x = pi/2;
1832  }
1833
1834  // x outside period and abs(n)>1 agrees with mpmath, differs from Mathematica
1835
1836  if ( isComplex(n) || isComplex(x) || isComplex(m) ) {
1837
1838    if ( !isComplex(x) ) x = complex(x);
1839
1840    var period = complex(0);
1841    if ( Math.abs(x.re) > pi/2 ) {
1842      var p = Math.round( x.re / pi );
1843      x.re = x.re - p * pi;
1844      period = mul( 2 * p, ellipticPi( n, m ) );
1845    }
1846
1847    return add( mul( sin(x), carlsonRF( mul(cos(x),cos(x)), sub( 1, mul(m,sin(x),sin(x)) ), 1 ) ),
1848                mul( 1/3, n, pow(sin(x),3),
1849                  carlsonRJ( mul(cos(x),cos(x)), sub( 1, mul(m,sin(x),sin(x)) ), 1,
1850                                 sub( 1, mul(n,sin(x),sin(x)) ) ) ),
1851                period );
1852
1853  } else {
1854
1855    if ( n > 1 && Math.abs(x) > Math.asin( 1 / Math.sqrt(n) ) ) return ellipticPi( n, complex(x), m );
1856
1857    if ( m > 1 && Math.abs(x) > Math.asin( 1 / Math.sqrt(m) ) ) return ellipticPi( n, complex(x), m );
1858
1859    var period = 0;
1860    if ( Math.abs(x) > pi/2 ) {
1861      var p = Math.round( x / pi );
1862      x = x - p * pi;
1863      period = 2 * p * ellipticPi( n, m );
1864    }
1865
1866    return sin(x) * carlsonRF( cos(x)**2, 1 - m * sin(x)**2, 1 )
1867           + n / 3 * sin(x)**3
1868             * carlsonRJ( cos(x)**2, 1 - m * sin(x)**2, 1, 1 - n * sin(x)**2 )
1869           + period;
1870
1871  }
1872
1873}
1874
1875
1876function jacobiZeta( x, m ) {
1877
1878  // using definition matching elliptic integrals
1879  // alternate definition replaces x with am(x,m)
1880
1881  return sub( ellipticE( x, m ), mul( ellipticF(x,m), ellipticE(m), inv( ellipticK(m) ) ) );
1882
1883}
1884
1885
1886function factorial( n ) {
1887
1888  if ( isComplex(n) ) {
1889
1890    if ( isPositiveIntegerOrZero(n) ) return complex( factorial(n.re) );
1891
1892    return gamma( add(n,1) );
1893
1894  }
1895
1896  if ( isPositiveIntegerOrZero(n) ) {
1897
1898    if ( factorialCache[n] ) return factorialCache[n];
1899
1900    var last = factorialCache.length - 1;
1901    var result = factorialCache[last];
1902
1903    for ( var i = last + 1 ; i <= n ; i++ ) {
1904      result *= i;
1905      factorialCache[i] = result;
1906    }
1907
1908    return result;
1909
1910  }
1911
1912  return gamma( n+1 );
1913
1914}
1915
1916function factorial2( n ) {
1917
1918  if ( isComplex(n) ) {
1919
1920    if ( isPositiveIntegerOrZero(n) ) return complex( factorial2(n.re) );
1921
1922    var p1 = pow( 2, div(n,2) );
1923    var p2 = pow( pi/2, div( sub( cos(mul(pi,n)), 1 ), 4 ) );
1924    var p3 = gamma( add( div(n,2) , 1 ) );
1925
1926    return mul( p1, p2, p3 );
1927
1928  }
1929
1930  if ( n === 0 ) return 1;
1931
1932  if ( isPositiveInteger(n) ) {
1933
1934    // bitwise test for odd integer, upward recursion for possible caching
1935    var result = n & 1 ? 1 : 2;
1936    for ( var i = result + 2 ; i <= n ; i += 2 ) result *= i;
1937    return result;
1938
1939  }
1940
1941  // github.com/sympy/sympy/pull/10362#issuecomment-170344298
1942  
1943  return 2**(n/2) * (pi/2)**((cos(pi*n)-1)/4) * gamma( n/2+1 );
1944
1945}
1946
1947function binomial( n, m ) {
1948
1949  if ( isComplex(n) || isComplex(m) )
1950    return div( factorial(n), mul( factorial( sub(n,m) ), factorial(m) ) );
1951
1952  if ( Number.isInteger(m) && m < 0 && n >= 0 ) return 0;
1953
1954  if ( Number.isInteger(n) && Number.isInteger(m) && n >= 0 && m > n ) return 0;
1955
1956  return factorial(n) / factorial(n-m) / factorial(m);
1957
1958}
1959
1960function multinomial( n ) {
1961
1962  if ( arguments.length === 1 ) return isComplex(n) ? complex(1) : 1;
1963
1964  var s = 0, p = 1, nNeg = 0;
1965
1966  for ( var i = 0 ; i < arguments.length ; i++ ) {
1967
1968    var ni = arguments[i];
1969    s = add( s, ni );
1970
1971    // catch negative integers for cancellation in ratio
1972    // reduce all to factorial(-1) divided by factor
1973    // extra overhead with pow() keeps complex in complex out
1974
1975    if ( isNegativeInteger( ni ) ) {
1976      var k = sub( -1, ni );
1977      p = div( p, mul( chop(pow(-1,k)), factorial(k) ) );
1978      nNeg++;
1979    } else
1980      // standard factorial
1981      p = mul( p, factorial( arguments[i] ) );
1982
1983  }
1984
1985  if ( isNegativeInteger(s) ) {
1986    nNeg--;
1987    if ( nNeg > 0 ) return isComplex(s) ? complex(0) : 0;
1988    var k = sub( -1, s );
1989    return inv( mul( p, chop(pow(-1,k)), factorial(k) ) );
1990  }
1991
1992  if ( nNeg > 0 ) return isComplex(s) ? complex(0) : 0;
1993
1994  return div( factorial(s), p );
1995
1996}
1997
1998function pochhammer( x, n ) {
1999
2000  var one = isArbitrary(x) ? arb1 : 1;
2001
2002  if ( isZero(n) ) {
2003    if ( isComplex(x) || isComplex(n) ) return complex(one);
2004    return one;
2005  }
2006
2007  if ( isComplex(n) ) {
2008
2009    if ( isPositiveInteger(n) ) {
2010      var result = pochhammer( x, n.re );
2011      if ( isComplex(result) ) return result;
2012      return complex(result);
2013    }
2014
2015    if ( isArbitrary(x) && !isArbitrary(n) ) n = arbitrary(n);
2016
2017    return div( gamma( add(x,n) ), gamma(x) );
2018
2019  }
2020
2021  if ( isPositiveInteger(n) ) {
2022
2023    var result = x, current = one, count = 1;
2024
2025    while ( count < n ) {
2026      result = mul( result, add( x, current ) );
2027      current = add( current, one );
2028      count++;
2029    }
2030
2031    return result;
2032
2033  }
2034
2035  if ( isArbitrary(x) && !isArbitrary(n) ) n = arbitrary(n);
2036
2037  return div( gamma( add(x,n) ), gamma(x) );
2038
2039}
2040
2041function subfactorial( n ) {
2042
2043  var result = div( gamma( add(n,1), -1 ), exp(1) );
2044
2045  if ( isPositiveInteger(n) ) return round(result);
2046
2047  return result;
2048
2049}
2050
2051
2052function logGamma( x ) {
2053
2054  if ( isArbitrary(x) ) {
2055
2056    var useAsymptotic = 10n * arb1;
2057
2058    if ( isNegativeIntegerOrZero( arbitrary(x) ) ) throw Error( 'Gamma function pole' );
2059
2060    if ( abs(x) < useAsymptotic ) {
2061
2062      // slightly faster to get smaller integer near transition circle
2063      var xRe = isComplex(x) ? x.re : x;
2064      var n = Math.ceil( arbitrary(useAsymptotic) - arbitrary(xRe) );
2065

2066      // ln(pochhammer(x,n))
2067      var lnPochhammer = ln(x), current = arb1, count = 1;
2068
2069      while ( count < n ) {
2070        lnPochhammer = add( lnPochhammer, ln( add( x, current ) ) );
2071        current = add( current, arb1 );
2072        count++;
2073      }
2074
2075      return sub( logGamma(add(x,arbitrary(n))), lnPochhammer );
2076
2077    }
2078
2079    if ( x < 0n ) x = complex(x);
2080
2081    // reflection formula with non-Hare correction to imaginary part
2082    if ( x.re < 0n ) {
2083
2084      // expand sine as exponentials for more accurate result
2085      var t, e, p, k = arb1;
2086
2087      if ( abs(x.im) === 0 ) {
2088
2089        t = neg( ln( sin( mul( x, onePi ) ) ) );
2090
2091      } else if ( x.im < 0n ) {
2092
2093        e = exp( mul( x, complex(0n,-twoPi) ) );
2094        p = complex( e.re, e.im );
2095        t = complex( e.re, e.im );
2096        while ( p.re !== 0n || p.im !== 0n ) {
2097          k += arb1;
2098          p = mul( p, e );
2099          t = add( t, div( p, k ) );
2100        }
2101        t = add( t, mul( x, complex(0n,-onePi) ), ln(arb2) );
2102
2103      } else {
2104
2105        e = exp( mul( x, complex(0n,twoPi) ) );
2106        p = complex( e.re, e.im );
2107        t = complex( e.re, e.im );
2108        while ( p.re !== 0n || p.im !== 0n ) {
2109          k += arb1;
2110          p = mul( p, e );
2111          t = add( t, div( p, k ) );
2112        }
2113        t = add( t, mul( x, complex(0n,onePi) ), complex(0n,-onePi), ln(arb2) );
2114
2115      }
2116
2117      return add( t, ln(onePi), neg( logGamma( sub(arb1,x) ) ), complex(0n,halfPi) );
2118
2119    }
2120

2121    // Johansson arxiv.org/abs/2109.08392
2122
2123    var k = 1, K = arb1, y, p = 1n, s = 0n;
2124
2125    if ( isComplex(x) ) {
2126      function compare() { return p.re !== 0n || p.im !== 0n; };
2127      y = complex( x.re, x.im );
2128    } else {
2129      function compare() { return p !== 0n; };
2130      y = x;
2131    }
2132
2133    while ( compare() ) {
2134
2135      if ( k === bernoulli2nN.length ) {
2136        console.log( 'Not enough Bernoulli numbers for logGamma' );
2137        break;
2138      }
2139
2140      p = div( div( bernoulli2nN[k], bernoulli2nD[k] ),
2141               mul( K+arb1, K, y ) );
2142      s = add( s, p );
2143
2144      y = mul( y, x, x );
2145      K += arb2;
2146      k++;
2147
2148    }
2149
2150    return add( mul( sub(x,arb1/2n), ln(x) ), neg(x), div(ln(twoPi),arb2), s );
2151
2152  }
2153
2154  // log of gamma less likely to overflow than gamma
2155  // Lanczos approximation as evaluated by Paul Godfrey
2156
2157  var c = [ 57.1562356658629235, -59.5979603554754912, 14.1360979747417471,
2158            -.491913816097620199, .339946499848118887e-4, .465236289270485756e-4,
2159            -.983744753048795646e-4, .158088703224912494e-3, -.210264441724104883e-3,
2160            .217439618115212643e-3, -.164318106536763890e-3, .844182239838527433e-4,
2161            -.261908384015814087e-4, .368991826595316234e-5 ];
2162
2163  if ( isNegativeIntegerOrZero(x) ) throw Error( 'Gamma function pole' );
2164
2165  if ( isComplex(x) ) {
2166
2167    // reflection formula with modified Hare correction to imaginary part
2168    if ( x.re < 0 ) {
2169
2170      var logRatio = log( div( pi, sin( mul(pi,x) ) ) );
2171      // rounding errors can lead to wrong side of branch point
2172      if ( isNegativeIntegerOrZero( ( x.re + .5 ) / 2 ) )
2173        logRatio.im = pi * ( x.im > 0 ? 1 : -1 ); // avoid Math.sign(0) = 0
2174
2175      var t = sub( logRatio, logGamma( sub(1,x) ) );
2176      var s = x.im < 0 ? -1 : 1;
2177      var d = x.im === 0 ? 1/4 : 0;
2178      var k = Math.ceil( x.re/2 - 3/4 + d );
2179
2180      return add( t, complex( 0, 2*s*k*pi ) );
2181
2182    }
2183
2184    var t = add( x, 5.24218750000000000 );
2185    t = sub( mul( add( x, .5 ), log(t)), t );
2186    var s = .999999999999997092;
2187    for ( var j = 0 ; j < 14 ; j++ ) s = add( s, div( c[j], add( x, j+1 ) ) );
2188    var u = add( t, log( mul( 2.5066282746310005, div( s, x ) ) ) );
2189
2190    // adjustment to keep imaginary part on same sheet
2191    if ( s.re < 0 ) {
2192      if( x.im < 0 && div(s,x).im < 0 ) u = add( u, complex(0,2*pi) );
2193      if( x.im > 0 && div(s,x).im > 0 ) u = add( u, complex(0,-2*pi) );
2194    }
2195
2196    return u;
2197
2198  } else {
2199
2200    if ( x < 0 ) return logGamma( complex(x) ); 
2201
2202    var t = x + 5.24218750000000000;
2203    t = ( x + .5 ) * log(t) - t;
2204    var s = .999999999999997092;
2205    for ( var j = 0 ; j < 14 ; j++ ) s += c[j] / (x+j+1);
2206    return t + log( 2.5066282746310005 * s / x );
2207
2208  }
2209
2210}
2211
2212function gamma( x, y, z ) {
2213
2214  if ( arguments.length === 2 ) {
2215
2216    if ( isZero(x) ) {
2217
2218      if ( isZero(y) ) throw Error( 'Gamma function pole' );
2219
2220      var result = neg( expIntegralEi( neg(y), true ) );
2221
2222      // complex on negative real axis
2223      if ( y < 0 || y.re < 0 && y.im === 0 ) result = sub( result, complex(0,pi) );
2224
2225      return result;
2226
2227    }
2228
2229    // dlmf.nist.gov/8.4.15
2230    if ( isNegativeInteger(x) ) {
2231
2232      var n = isComplex(x) ? -x.re : -x;
2233      var s = inv(y), p = s; // mul returns new object
2234
2235      for ( var k = 1 ; k < n ; k++ ) {
2236        p = mul( p, -k, inv(y) );
2237        s = add( s, p );
2238      }
2239
2240      var t = mul( exp(neg(y)), s );
2241
2242      // dlmf.nist.gov/8.4.4
2243      var result = sub( neg( expIntegralEi( neg(y), true, 1e-14 ) ), t )
2244
2245      // complex on negative real axis
2246      if ( y < 0 || y.re < 0 && y.im === 0 ) result = sub( result, complex(0,pi) );
2247
2248      result = mul( (-1)**n/factorial(n), result );
2249
2250      if ( isComplex(x) && !isComplex(result) ) return complex(result); // complex in, complex out
2251
2252      return result;
2253
2254    }
2255
2256    return mul( exp(neg(y)), hypergeometricU( sub(1,x), sub(1,x), y ) );
2257
2258  }
2259
2260  if ( arguments.length === 3 ) {
2261
2262    if ( !isZero(y) ) return sub( gamma(x,0,z), gamma(x,0,y) );
2263
2264    return mul( pow(z,x), inv(x), hypergeometric1F1( x, add(x,1), neg(z) ) );
2265
2266  }
2267
2268  if ( isPositiveInteger(x) ) return factorial( sub(x,1) );
2269
2270  // logGamma complex on negative axis
2271  if ( !isComplex(x) && x < 0 ) return exp( logGamma( complex(x) ) ).re;
2272
2273  return exp( logGamma(x) );
2274
2275}
2276
2277function gammaRegularized( x, y, z ) {
2278
2279  if ( arguments.length === 3 ) return div( gamma(x,y,z), gamma(x) );
2280
2281  return div( gamma(x,y), gamma(x) );
2282
2283}
2284
2285function beta( x, y, z, w ) {
2286
2287  if ( arguments.length === 4 )
2288
2289    return sub( beta(y,z,w), beta(x,z,w) );
2290
2291  if ( arguments.length === 3 )
2292
2293    return mul( pow(x,y), inv(y), hypergeometric2F1( y, sub(1,z), add(y,1), x ) );
2294
2295  return div( mul( gamma(x), gamma(y) ), gamma( add(x,y) ) ); 
2296
2297}
2298
2299function betaRegularized( x, y, z, w ) {
2300
2301  if ( arguments.length === 4 )
2302
2303    return div( beta(x,y,z,w), beta(z,w) );
2304
2305  return div( beta(x,y,z), beta(y,z) );
2306
2307}
2308
2309function polygamma( n, x ) {
2310
2311  if ( arguments.length === 1 ) return digamma(x);
2312
2313  if ( !isPositiveInteger(n) ) throw Error( 'Unsupported polygamma index' );
2314
2315  return mul( (-1)**(n+1) * factorial(n), hurwitzZeta( n+1, x ) );
2316
2317}
2318
2319
2320function digamma( x ) {
2321
2322  return diff( x => logGamma(x), x );
2323
2324}
2325
2326
2327function erf( x ) {
2328
2329  var useAsymptotic = 5;
2330
2331  var absArg = Math.abs( arg(x) );
2332
2333  if ( abs(x) > useAsymptotic && ( absArg < pi/4 || absArg > 3*pi/4 ) )
2334
2335    return sub( 1, erfc(x) );
2336
2337  return mul( 2/sqrt(pi), x, hypergeometric1F1( .5, 1.5, neg(mul(x,x)) ) );
2338
2339}
2340
2341function erfc( x ) {
2342
2343  var useAsymptotic = 5;
2344
2345  var absArg = Math.abs( arg(x) );
2346
2347  if ( abs(x) > useAsymptotic && ( absArg < pi/4 || absArg > 3*pi/4 ) ) {
2348
2349    // as per dlmf.nist.gov/7.12.1 this could be an independent sum for minor improvement
2350    // these numbers are tiny and need to stay in this function even though
2351    //   there is some code duplication with erf
2352
2353    var t = mul( 1/sqrt(pi), exp( neg(mul(x,x)) ), inv(x),
2354                 hypergeometric2F0( .5, 1, neg(inv(mul(x,x))) ) );
2355
2356    if ( x.re < 0 || x < 0 ) return add( 2, t );
2357
2358    return t;
2359
2360  }
2361
2362  return sub( 1, erf(x) );
2363
2364}
2365
2366function erfi( x ) {
2367
2368  return mul( complex(0,-1), erf( mul( complex(0,1), x ) ) );
2369
2370}
2371
2372
2373function fresnelS( x ) {
2374
2375  // can also be evaluated with hypergeometric1F2
2376
2377  var m1 = hypergeometric1F1( .5, 1.5, mul( complex(0,pi/2), pow(x,2) ) );
2378  var m2 = hypergeometric1F1( .5, 1.5, mul( complex(0,-pi/2), pow(x,2) ) );
2379
2380  var result = mul( x, sub( m1, m2 ), complex(0,-.5) );
2381
2382  if ( isComplex(x) ) return result;
2383
2384  return result.re;
2385
2386}
2387
2388function fresnelC( x ) {
2389
2390  // can also be evaluated with hypergeometric1F2
2391
2392  var m1 = hypergeometric1F1( .5, 1.5, mul( complex(0,pi/2), pow(x,2) ) );
2393  var m2 = hypergeometric1F1( .5, 1.5, mul( complex(0,-pi/2), pow(x,2) ) );
2394
2395  var result = mul( x, add( m1, m2 ), .5 );
2396
2397  if ( isComplex(x) ) return result;
2398
2399  return result.re;
2400
2401}
2402
2403

2404function expIntegralEi( x, adjustImForGamma=false, tolerance=1e-10 ) {
2405
2406  // direct series hangs at -1 and +i/-i
2407  // complex average until investigate further
2408
2409  if ( isUnity(neg(x)) ) return complexAverage( x => expIntegralEi(x,adjustImForGamma,tolerance), x );
2410
2411  var ix = mul( x, complex(0,1) );
2412  if ( isUnity(ix) || isUnity(neg(ix)) )
2413    return complexAverage( x => expIntegralEi(x,adjustImForGamma,tolerance), x, complex(0,1e-5) );
2414
2415  var useAsymptotic = 26;
2416
2417  if ( isComplex(x) ) {
2418
2419    if ( abs(x) > useAsymptotic ) {
2420
2421      var s = complex(1);
2422      var p = complex(1), pLast = p; // mul returns new object
2423      var i = 1;
2424
2425      while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
2426
2427        p = mul( p, i, inv(x) );
2428        if ( abs(p) > abs(pLast) ) break;
2429
2430        s = add( s, p );
2431        i++;
2432        pLast = p;
2433
2434      }
2435
2436      if ( adjustImForGamma )
2437
2438        return mul( s, exp(x), inv(x) );
2439
2440      else {
2441
2442        // combination of logarithms adds/subtracts complex(0,pi)
2443        var sign = x.im > 0 ? 1 : x.im < 0 ? -1 : 0;
2444
2445        return add( mul( s, exp(x), inv(x) ), complex(0,sign*pi) );
2446
2447      }
2448
2449    }
2450
2451    // determined from pattern on test page
2452    var distanceScale = abs( sub(x,useAsymptotic) ) / useAsymptotic;
2453    var useArbitrary = distanceScale > 1;
2454
2455    // arbitrary precision series is unstable around the origin but not needed
2456
2457    if ( useArbitrary && abs(x) > 1.5 ) {
2458
2459      // use only decimals needed
2460      var n = 17 + Math.round( 10 * ( distanceScale - 1 ) );
2461      setPrecisionScale( n );
2462
2463      var y = arbitrary( x );
2464
2465      var s = arbitrary( complex(0) );
2466      var p = arbitrary( complex(1) );
2467      var i = arb1;
2468
2469      while ( div(p.re,i) !== 0n || div(p.im,i) !== 0n ) {
2470        p = div( mul(p,y), i );
2471        s = add( s, div(p,i) );
2472        i = add( i, arb1 );
2473      }
2474
2475      s = add( s, getConstant( 'eulerGamma' ), ln(y) );
2476
2477      if ( adjustImForGamma ) {
2478        var sign = x.im > 0n ? -1n : x.im < 0n ? 1n : 0n;
2479        s = add( s, complex(0n,sign*onePi) );
2480      }
2481
2482      s = arbitrary( s );
2483
2484      resetPrecisionScale();
2485
2486      // adjust phase on positive imaginary axis according to test
2487      if ( x.im > 0 && x.re === 0 ) s.im += pi;
2488
2489    } else {
2490
2491      var s = complex(0);
2492      var p = complex(1);
2493      var i = 1;
2494
2495      while ( Math.abs(p.re/i) > tolerance || Math.abs(p.im/i) > tolerance ) {
2496        p = mul( p, x, 1/i );
2497        s = add( s, div(p,i) );
2498        i++;
2499      }
2500
2501      s = add( s, eulerGamma, log(x) );
2502
2503      if ( adjustImForGamma ) {
2504        var sign = x.im > 0 ? -1 : x.im < 0 ? 1 : 0;
2505        s = add( s, complex(0,sign*pi) );
2506      }
2507
2508    }
2509
2510    // real on negative real axis, set phase explicitly rather than log combo
2511    if ( x.re < 0 && x.im === 0 ) s.im = 0;
2512
2513    return s;
2514
2515  } else {
2516
2517    if ( x < 0 ) return expIntegralEi( complex(x) ).re;
2518
2519    if ( Math.abs(x) > useAsymptotic ) {
2520
2521      var s = 1;
2522      var p = 1, pLast = p;
2523      var i = 1;
2524
2525      while ( Math.abs(p) > tolerance ) {
2526
2527        p *= i / x;
2528        if ( Math.abs(p) > Math.abs(pLast) ) break;
2529
2530        s += p;
2531        i++;
2532        pLast = p;
2533
2534      }
2535
2536      return s * Math.exp(x) / x;
2537
2538    }
2539
2540    var s = 0;
2541    var p = 1;
2542    var i = 1;
2543
2544    while ( Math.abs(p/i) > tolerance ) {
2545      p *= x / i;
2546      s += p / i;
2547      i++;
2548    }
2549
2550    return s + eulerGamma + Math.log(x);
2551
2552  }
2553
2554}
2555
2556function logIntegral( x ) {
2557
2558  return expIntegralEi( log(x) );
2559
2560}
2561
2562function sinIntegral( x ) {
2563
2564  if ( isZero(x) ) return isComplex(x) ? complex(0) : 0;
2565
2566  var ix = mul( complex(0,1), x );
2567
2568  var result = mul( complex(0,.5), add( gamma(0,neg(ix)), neg(gamma(0,ix)),
2569                                        log(neg(ix)), neg(log(ix)) ) );
2570
2571  if ( isComplex(x) ) return result;
2572
2573  return result.re;
2574
2575}
2576
2577function cosIntegral( x ) {
2578
2579  // complex for negative real argument
2580
2581  var ix = mul( complex(0,1), x );
2582

2583  var result = sub( log(x), mul( .5, add( gamma(0,neg(ix)), gamma(0,ix),
2584                                     log(neg(ix)), log(ix) ) ) );
2585
2586  // real for positive real argument
2587  if ( x > 0 ) return result.re;
2588  if ( x.re > 0 && x.im === 0 ) result.im = chop( result.im );
2589
2590  return result;
2591
2592}
2593
2594function sinhIntegral( x ) {
2595
2596  if ( isZero(x) ) return isComplex(x) ? complex(0) : 0;
2597
2598  var result = mul( .5, add( gamma(0,x), neg(gamma(0,neg(x))), log(x), neg(log(neg(x))) ) );
2599
2600  if ( isComplex(x) ) return result;
2601
2602  return result.re;
2603
2604}
2605
2606function coshIntegral( x ) {
2607
2608  // complex for negative real argument
2609
2610  return mul( -.5, add( gamma(0,x), gamma(0,neg(x)), neg(log(x)), log(neg(x)) ) );
2611
2612}
2613
2614function expIntegralE( n, x ) {
2615
2616  if ( isZero(n) ) return div( exp(neg(x)), x );
2617
2618  if ( isZero(x) && ( n > 1 || n.re > 1 ) ) return inv( sub(n,1) );
2619
2620  var p = pow( x, sub(n,1) );
2621
2622  // real on negative real axis for integer powers
2623  if ( isInteger(n) && x.re < 0 ) p.im = chop( p.im );
2624
2625  return mul( p, gamma( sub(1,n), x ) );
2626
2627}
2628
2629
2630function hypergeometric0F1( a, x, tolerance=1e-10 ) {
2631
2632  var useAsymptotic = 100;
2633
2634  if ( isNegativeIntegerOrZero(a) ) throw Error( 'Hypergeometric function pole' );
2635
2636  if ( isComplex(a) || isComplex(x) ) {
2637
2638    // asymptotic form as per Johansson arxiv.org/abs/1606.06977
2639    if ( abs(x) > useAsymptotic ) {
2640
2641      if ( isNegativeIntegerOrZero( sub(a,.5) ) ) {
2642        var result = complexAverage( a => hypergeometric0F1(a,x), a );
2643        if ( isReal(x) ) result.im = 0;
2644        return result;
2645      }
2646
2647      // transform variables for convenience
2648      var b = sub( mul(2,a), 1 );
2649      a = sub( a, .5 );
2650      x = mul( 4, sqrt(x) );
2651
2652      // copied from hypergeometric1F1
2653      var t1 = mul( gamma(b), pow( neg(x), neg(a) ), inv( gamma(sub(b,a)) ) );
2654      t1 = mul( t1, hypergeometric2F0( a, add(a,neg(b),1), div(-1,x) ) );
2655
2656      var t2 = mul( gamma(b), pow( x, sub(a,b) ), exp(x), inv( gamma(a) ) );
2657      t2 = mul( t2, hypergeometric2F0( sub(b,a), sub(1,a), div(1,x) ) );
2658
2659      var result = mul( exp( div(x,-2) ), add( t1, t2 ) );
2660
2661      // real for real arguments
2662      if ( isReal(a) && isReal(x) ) result.im = 0;
2663
2664      return result;
2665
2666    }
2667
2668    var s = complex(1);
2669    var p = complex(1);
2670    var i = 1;
2671
2672    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
2673      p = mul( p, x, inv(a), 1/i );
2674      s = add( s, p );
2675      a = add( a, 1 );
2676      i++;
2677    }
2678
2679    return s;
2680
2681  } else {
2682
2683    // asymptotic form is complex
2684    if ( Math.abs(x) > useAsymptotic ) return hypergeometric0F1( a, complex(x) ).re;
2685
2686    var s = 1;
2687    var p = 1;
2688    var i = 1;
2689
2690    while ( Math.abs(p) > tolerance ) {
2691      p *= x / a / i;
2692      s += p;
2693      a++;
2694      i++;
2695    }
2696
2697    return s;
2698
2699  }
2700
2701}
2702
2703
2704function hypergeometric1F1( a, b, x, tolerance=1e-10 ) {
2705
2706  if ( isEqualTo(a,b) ) {
2707
2708    if ( ( isComplex(a) || isComplex(b) ) && !isComplex(x) ) x = complex(x);
2709
2710    if ( isNegativeInteger(a) ) {
2711      var n = isComplex(a) ? -a.re : -a;
2712      return summation( k => div( pow(x,k), factorial(k) ), [0,n] );
2713    }
2714
2715    return exp(x);
2716
2717  }
2718
2719  if ( isArbitrary(x) ) {
2720
2721    if ( !isArbitrary(a) ) a = arbitrary(a);
2722    if ( !isArbitrary(b) ) b = arbitrary(b);
2723
2724    if ( isNegativeIntegerOrZero(arbitrary(b)) ) throw Error( 'Hypergeometric function pole' );
2725
2726    // asymptotic forms not yet implemented
2727
2728    if ( isComplex(a) || isComplex(b) || isComplex(x) ) {
2729
2730      if ( !isComplex(x) ) x = complex(x);
2731
2732      // Kummer transformation
2733      if ( x.re < 0n ) return mul( exp(x), hypergeometric1F1( sub(b,a), b, neg(x) ) );
2734
2735      var s = complex(arb1);
2736      var p = complex(arb1);
2737      var i = arb1;
2738
2739      while ( abs(p.re) !== 0n || abs(p.im) !== 0n ) {
2740        p = mul( p, x, a, inv(b), inv(i) );
2741        s = add( s, p );
2742        a = add( a, arb1 );
2743        b = add( b, arb1 );
2744        i += arb1;
2745      }
2746
2747      return s;
2748
2749    } else {
2750
2751      // Kummer transformation
2752      if ( x < 0n ) return mul( exp(x), hypergeometric1F1( b-a, b, -x ) );
2753
2754      var s = arb1;
2755      var p = arb1;
2756      var i = arb1;
2757
2758      while ( abs(p) !== 0n ) {
2759        p = mul( p, x, a, inv(b), inv(i) );
2760        s += p;
2761        a += arb1;
2762        b += arb1;
2763        i += arb1;
2764      }
2765
2766      return s;
2767
2768    }
2769
2770  }
2771
2772  var useAsymptotic = 30;
2773
2774  if ( isNegativeIntegerOrZero(b) ) throw Error( 'Hypergeometric function pole' );
2775
2776  if ( isComplex(a) || isComplex(b) || isComplex(x) ) {
2777
2778    if ( isZero(a) ) return complex(1);
2779
2780    // terminating hypergeometric series holds for entire complex plane
2781    if ( isNegativeInteger(a) ) return hypergeometricSeries( [a], [b], x, true );
2782
2783    if ( !isComplex(x) ) x = complex(x);
2784
2785    // Kummer transformation
2786    if ( x.re < 0 ) return mul( exp(x), hypergeometric1F1( sub(b,a), b, neg(x) ) );
2787
2788    // asymptotic form as per Johansson arxiv.org/abs/1606.06977
2789    if ( abs(x) > useAsymptotic ) {
2790
2791      if ( isNegativeInteger(sub(b,a)) ) // zero/negative a already handled, b is not a
2792        return complexAverage( a =>

2792 hypergeometric1F1(a,b,x), a );
2793
2794      var t1 = mul( gamma(b), pow( neg(x), neg(a) ), inv( gamma(sub(b,a)) ),
2795                    hypergeometric2F0( a, add(a,neg(b),1), div(-1,x) ) );
2796
2797      var t2 = mul( gamma(b), pow( x, sub(a,b) ), exp(x), inv( gamma(a) ),
2798                    hypergeometric2F0( sub(b,a), sub(1,a), div(1,x) ) );
2799
2800      var result = add( t1, t2 );
2801
2802      // real for real arguments
2803      if ( isReal(a) && isReal(b) && isReal(x) ) result.im = 0;
2804
2805      return result;
2806
2807    }
2808
2809    // function to describe wedges along imaginary axis inside
2810    //   direct region as per test pages
2811    function wedge( x ) {
2812      var slope = Math.sqrt(3) - 1;
2813      if ( x.re > 0 ) {
2814        var test = slope*x.re + useAsymptotic/2;
2815        return x.im > test || x.im < -test;
2816      } else {
2817        var test = -slope*x.re + useAsymptotic/2;
2818        return x.im > test || x.im < -test;
2819      }
2820    }
2821
2822    if ( wedge(x) ) {
2823      setPrecisionScale( 22 );
2824      var result = arbitrary( hypergeometric1F1( a, b, arbitrary(x) ) );
2825      resetPrecisionScale();
2826      return result;
2827    }
2828
2829    var s = complex(1);
2830    var p = complex(1);
2831    var i = 1;
2832
2833    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
2834      p = mul( p, x, a, inv(b), 1/i );
2835      s = add( s, p );
2836      a = add( a, 1 );
2837      b = add( b, 1 );
2838      i++;
2839    }
2840
2841    return s;
2842
2843  } else {
2844
2845    if ( a === 0 ) return 1;
2846
2847    // terminating series holds for entire axis
2848    if ( isNegativeInteger(a) ) return hypergeometricSeries( [a], [b], x );
2849
2850    // Kummer transformation
2851    if ( x < 0 ) return exp(x) * hypergeometric1F1( b-a, b, -x );
2852
2853    // asymptotic form is complex
2854    if ( x > useAsymptotic ) return hypergeometric1F1( a, b, complex(x) ).re;
2855
2856    var s = 1;
2857    var p = 1;
2858    var i = 1;
2859
2860    while ( Math.abs(p) > tolerance ) {
2861      p *= x * a / b / i;
2862      s += p;
2863      a++;
2864      b++;
2865      i++;
2866    }
2867
2868    return s;
2869
2870  }
2871
2872}
2873
2874
2875function hypergeometricU( a, b, x ) {
2876
2877  if ( isEqualTo( add(a,1), b ) ) return pow( x, neg(a) );
2878
2879  if ( isZero(a) )
2880    if ( isComplex(a) || isComplex(b) || isComplex(x) ) return complex(1);
2881    else return 1;
2882
2883  var useAsymptotic = 25;
2884
2885  // asymptotic form as per Johansson arxiv.org/abs/1606.06977
2886  if ( abs(x) > useAsymptotic ) {
2887
2888    return mul( pow( x, neg(a) ), hypergeometric2F0( a, add(a,neg(b),1), neg(inv(x)) ) );
2889
2890  }
2891
2892  // special case through independent expIntegralEi for speed
2893  if ( isUnity(a) && isUnity(b) ) return mul( exp(x), gamma(0,x) );
2894
2895  if ( isUnity(b) ) return complexAverage( b => hypergeometricU(a,b,x), b );
2896
2897  // could be combined with previous special case
2898  if ( isPositiveInteger(a) && isEqualTo(a,b) ) return mul( exp(x), gamma(sub(1,a),x) );
2899
2900  // dlmf.nist.gov/13.2.7
2901  if ( isNegativeInteger(a) ) {
2902    var m = isComplex(a) ? -a.re : -a;
2903    return mul( (-1)**m, pochhammer(b,m), hypergeometric1F1(a,b,x) );
2904  }
2905
2906  var t1 = mul( gamma(sub(b,1)), inv( gamma(a) ), pow( x, sub(1,b) ),
2907                hypergeometric1F1( add(a,neg(b),1), sub(2,b), x ) );
2908
2909  var t2 = mul( gamma(sub(1,b)), inv( gamma(add(a,neg(b),1)) ), hypergeometric1F1( a, b, x ) );
2910
2911  var max = Math.max( abs(t1), abs(t2) );
2912  if ( max < 100 ) return add( t1, t2 );
2913
2914  var extraDecimals = +max.toExponential().split( 'e' )[1];
2915  setPrecisionScale( defaultDecimals + extraDecimals );
2916
2917  a = arbitrary(a);
2918  b = arbitrary(b);
2919  x = arbitrary(x);
2920
2921  t1 = mul( gamma(sub(b,arb1)), inv( gamma(a) ), pow( x, sub(arb1,b) ),
2922            hypergeometric1F1( add(a,neg(b),arb1), sub(arb2,b), x ) );
2923
2924  t2 = mul( gamma(sub(arb1,b)), inv( gamma(add(a,neg(b),arb1)) ), hypergeometric1F1( a, b, x ) );
2925
2926  var result = arbitrary( add( t1, t2 ) );
2927
2928  resetPrecisionScale();
2929
2930  return result;
2931
2932}
2933
2934function whittakerM( k, m, x ) {
2935
2936  return mul( exp( mul(-.5,x) ), pow( x, add(m,.5) ),
2937              hypergeometric1F1( add(m,neg(k),.5), add(mul(2,m),1), x ) );
2938
2939}
2940
2941function whittakerW( k, m, x ) {
2942
2943  return mul( exp( mul(-.5,x) ), pow( x, add(m,.5) ),
2944              hypergeometricU( add(m,neg(k),.5), add(mul(2,m),1), x ) );
2945
2946}
2947
2948
2949function hypergeometric2F0( a, b, x, tolerance=1e-10 ) {
2950
2951  var terms = 50;
2952
2953  if ( isComplex(a) || isComplex(b) || isComplex(x) ) {
2954
2955    var s = complex(1);
2956    var p = complex(1), pLast = p; // mul returns new object
2957    var converging = false; // first few terms can be larger than unity
2958    var i = 1;
2959
2960    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
2961
2962      p = mul( p, x, a, b, 1/i );
2963
2964      if ( abs(p) > abs(pLast) && converging ) break; // prevent runaway sum
2965      if ( abs(p) < abs(pLast) ) converging = true;
2966      if ( i > terms ) throw Error( 'Not converging after ' + terms + ' terms' );
2967
2968      s = add( s, p );
2969      a = add( a, 1 );
2970      b = add( b, 1 );
2971      i++;
2972      pLast = p;
2973
2974    }
2975
2976    return s;
2977
2978  } else {
2979
2980    var s = 1;
2981    var p = 1, pLast = p;
2982    var converging = false; // first few terms can be larger than unity
2983    var i = 1;
2984
2985    while ( Math.abs(p) > tolerance ) {
2986
2987      p *= x * a * b / i;
2988
2989      if ( Math.abs(p) > Math.abs(pLast) && converging ) break; // prevent runaway sum
2990      if ( Math.abs(p) < Math.abs(pLast) ) converging = true;
2991      if ( i > terms ) throw Error( 'Not converging after ' + terms + ' terms' );
2992
2993      s += p;
2994      a++;
2995      b++;
2996      i++;
2997      pLast = p;
2998
2999    }
3000
3001    return s;
3002
3003  }
3004
3005}
3006
3007
3008function hypergeometric2F1( a, b, c, x, tolerance=1e-10 ) {
3009
3010  if ( isEqualTo(a,c) ) return pow( sub(1,x), neg(b) );
3011
3012  if ( isEqualTo(b,c) ) return pow( sub(1,x), neg(a) );
3013
3014  if ( !isUnity(neg(x)) && abs(add(x,1)) < tolerance )
3015    return hypergeometric2F1( a, b, c, isComplex(x) ? complex(-1) : -1 );
3016
3017  if ( isNegativeIntegerOrZero(c) ) throw Error( 'Hypergeometric function pole' );
3018
3019  if ( isComplex(a) || isComplex(b) || isComplex(c) || isComplex(x) ) {
3020
3021    if ( isZero(a) || isZero(b) ) return complex(1);
3022
3023    // terminating hypergeometric series holds for entire complex plane
3024    if ( isNegativeInteger(a) || isNegativeInteger(b) )
3025      return hypergeometricSeries( [a,b], [c], x, true );
3026
3027    if ( isUnity(x) )
3028      if ( re(sub(c,add(a,b))) > 0 )
3029        return div( mul( gamma(c), gamma(sub(c,add(a,b))) ),
3030                    mul( gamma(sub(c,a)), gamma(sub(c,b)) ) );
3031      else throw Error( 'Divergent Gauss hypergeometric function' );
3032
3033    // choose smallest absolute value of transformed argument
3034    // transformations from Abramowitz & Stegun p.559
3035    // fewer operations compared to dlmf.nist.gov/15.8
3036
3037    var absArray = [ abs(x), abs(div(x,sub(x,1))), abs(sub(1,x)),
3038                     abs(inv(x)), abs(inv(sub(1,x))), abs(sub(1,inv(x))) ];
3039
3040    var index = absArray.indexOf( Math.min.apply( null, absArray ) );
3041
3042    switch( index ) {
3043
3044      case 0:
3045
3046        break;
3047
3048      case 1:
3049
3050        return mul( pow( sub(1,x), neg(a) ), hypergeometric2F1( a, sub(c,b), c, div(x,sub(x,1)) ) );
3051
3052      case 2:
3053
3054        // zero/negative a or b already handled, c is not a or b
3055        // remaining integer values hit one term or other
3056
3057        if ( isNegativeInteger(sub(c,a)) || isInteger(sub(c,add(a,b))) )
3058          return complexAverage( a =>

3058 hypergeometric2F1(a,b,c,x), a );
3059
3060        if ( isNegativeInteger(sub(c,b)) )
3061          return complexAverage( b => hypergeometric2F1(a,b,c,x), b );
3062
3063        var t1 = mul( gamma(c), gamma( sub(c,add(a,b)) ), 
3064                      inv( gamma(sub(c,a)) ), inv( gamma(sub(c,b)) ),
3065                      hypergeometric2F1( a, b, add(a,b,neg(c),1), sub(1,x) ) );
3066
3067        var t2 = mul( pow( sub(1,x), sub(c,add(a,b)) ),
3068                      gamma(c), gamma( sub(add(a,b),c) ), inv( gamma(a) ), inv( gamma(b) ),
3069                      hypergeometric2F1( sub(c,a), sub(c,b), add(c,neg(a),neg(b),1), sub(1,x) ) );
3070
3071        return add( t1, t2 );
3072
3073      case 3:
3074
3075        if ( isNegativeInteger(sub(c,a)) || isInteger(sub(b,a)) )
3076          return complexAverage( a => hypergeometric2F1(a,b,c,x), a );
3077
3078        if ( isNegativeInteger(sub(c,b)) )
3079          return complexAverage( b => hypergeometric2F1(a,b,c,x), b );
3080
3081        var t1 = mul( gamma(c), gamma(sub(b,a)), inv( gamma(b) ),
3082                      inv( gamma(sub(c,a)) ), pow( neg(x), neg(a) ),
3083                      hypergeometric2F1( a, add(1,neg(c),a), add(1,neg(b),a), inv(x) ) );
3084
3085        var t2 = mul( gamma(c), gamma(sub(a,b)), inv( gamma(a) ),
3086                      inv( gamma(sub(c,b)) ), pow( neg(x), neg(b) ),
3087                      hypergeometric2F1( b, add(1,neg(c),b), add(1,neg(a),b), inv(x) ) );
3088
3089        return add( t1, t2 );
3090
3091      case 4:
3092
3093        if ( isNegativeInteger(sub(c,a)) || isInteger(sub(b,a)) )
3094          return complexAverage( a => hypergeometric2F1(a,b,c,x), a );
3095
3096        if ( isNegativeInteger(sub(c,b)) )
3097          return complexAverage( b => hypergeometric2F1(a,b,c,x), b );
3098
3099        var t1 = mul( pow( sub(1,x), neg(a) ), gamma(c), gamma(sub(b,a)),
3100                      inv( gamma(b) ), inv( gamma(sub(c,a)) ),
3101                      hypergeometric2F1( a, sub(c,b), add(a,neg(b),1), inv(sub(1,x)) ) );
3102
3103        var t2 = mul( pow( sub(1,x), neg(b) ), gamma(c), gamma(sub(a,b)),
3104                      inv( gamma(a) ), inv( gamma(sub(c,b)) ),
3105                      hypergeometric2F1( b, sub(c,a), add(b,neg(a),1), inv(sub(1,x)) ) );
3106
3107        return add( t1, t2 );
3108
3109      case 5:
3110
3111        if ( isNegativeInteger(sub(c,a)) || isInteger(sub(c,add(a,b))) )
3112          return complexAverage( a => hypergeometric2F1(a,b,c,x), a );
3113
3114        if ( isNegativeInteger(sub(c,b)) )
3115          return complexAverage( b => hypergeometric2F1(a,b,c,x), b );
3116
3117        var t1 = mul( gamma(c), gamma( sub(c,add(a,b)) ), inv( gamma(sub(c,a)) ),
3118                      inv( gamma(sub(c,b)) ), pow( x, neg(a) ),
3119                      hypergeometric2F1( a, add(a,neg(c),1), add(a,b,neg(c),1), sub(1,inv(x)) ) );
3120
3121        var t2 = mul( gamma(c), gamma( sub(add(a,b),c) ), inv( gamma(a) ), inv( gamma(b) ),
3122                      pow( sub(1,x), sub(c,add(a,b)) ), pow( x, sub(a,c) ),
3123                      hypergeometric2F1( sub(c,a), sub(1,a), add(c,neg(a),neg(b),1), sub(1,inv(x)) ) );
3124
3125        return add( t1, t2 );
3126
3127    }
3128
3129    var s = complex(1);
3130    var p = complex(1);
3131    var i = 1;
3132
3133    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
3134      p = mul( p, x, a, b, inv(c), 1/i );
3135      s = add( s, p );
3136      a = add( a, 1 );
3137      b = add( b, 1 );
3138      c = add( c, 1 );
3139      i++;
3140    }
3141
3142    return s;
3143
3144  } else {
3145
3146    if ( a === 0 || b === 0 ) return 1;
3147
3148    // terminating series holds for entire axis
3149    if ( isNegativeInteger(a) || isNegativeInteger(b) )
3150      return hypergeometricSeries( [a,b], [c], x );
3151
3152    // transformation from Abramowitz & Stegun p.559
3153    if ( x < -1 ) {
3154
3155      // zero/negative a or b already handled, c is not a or b
3156
3157      if ( isNegativeInteger(c-a) )
3158        return complexAverage( a => hypergeometric2F1(a,b,c,x), a );
3159
3160      if ( isNegativeInteger(c-b) )
3161        return complexAverage( b => hypergeometric2F1(a,b,c,x), b );
3162
3163      var t1 = gamma(c) * gamma(b-a) / gamma(b) / gamma(c-a)
3164                 * (-x)**(-a) * hypergeometric2F1( a, 1-c+a, 1-b+a, 1/x );
3165      var t2 = gamma(c) * gamma(a-b) / gamma(a) / gamma(c-b)
3166                 * (-x)**(-b) * hypergeometric2F1( b, 1-c+b, 1-a+b, 1/x );
3167
3168      return t1 + t2;
3169
3170    }
3171
3172    // faster than direct sum for general case
3173    if ( x === -1 ) return hypergeometric2F1( a, b, c, complex(x) ).re;
3174
3175    if ( x === 1 )
3176      if ( c - a - b > 0 ) return gamma(c) * gamma(c-a-b) / gamma(c-a) / gamma(c-b);
3177      else throw Error( 'Divergent Gauss hypergeometric function' );
3178
3179    if ( x > 1 ) return hypergeometric2F1( a, b, c, complex(x) );
3180
3181    var s = 1;
3182    var p = 1;
3183    var i = 1;
3184
3185    while ( Math.abs(p) > tolerance ) {
3186      p *= x * a * b / c / i;
3187      s += p;
3188      a++;
3189      b++;
3190      c++;
3191      i++;
3192    }
3193
3194    return s;
3195
3196  }
3197
3198}
3199
3200
3201function hypergeometric1F2( a, b, c, x, tolerance=1e-10 ) {
3202
3203  var useAsymptotic = 150;
3204
3205  if ( isComplex(a) || isComplex(b) || isComplex(c) || isComplex(x) ) {
3206
3207    // functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F2/06/02/03/0002/
3208
3209    if ( abs(x) > useAsymptotic ) {
3210
3211      var aNegbNegc = add( a, neg(b), neg(c) );
3212
3213      var p = div( add( aNegbNegc, 1/2 ), 2 );
3214
3215      var ck = [ 1 ];
3216
3217      ck[1] = add( mul( add(mul(3,a),b,c,-2), aNegbNegc, 1/2 ), mul(2,b,c), -3/8 );
3218
3219      var toSquare = add( mul(b,c), mul( aNegbNegc, add(mul(3,a),b,c,-2), 1/4 ), -3/16 );
3220
3221      ck[2] = add( mul( 2, toSquare, toSquare ),
3222                   neg( mul( sub(mul(2,a),3), b, c ) ),
3223                   mul( aNegbNegc, add( mul(-8,mul(a,a)), mul(11,a), b, c, -2 ), 1/4 ),
3224                   -3/16 );
3225
3226      var plusI = complex(0,1), minusI = complex(0,-1);
3227
3228      function w( k ) { return mul( 1/2**k, ck[k], pow(neg(x),-k/2) ); }
3229
3230      var w1 = w(1), w2 = w(2);
3231
3232      var s1 = add( 1, mul(minusI,w1), neg(w2) );
3233      var s2 = add( 1, mul(plusI,w1), neg(w2) );
3234

3235      var k = 2, wk = w(k), powPlusI = powMinusI = -1;
3236
3237      // wk not always complex
3238      while ( abs(wk) > tolerance ) {
3239
3240        k++;
3241        powPlusI = mul( powPlusI, plusI );
3242        powMinusI = mul( powMinusI, minusI );
3243
3244        var t1 = mul( add( 3*k**2, mul( 2*k, add(mul(-3,a),b,c,-2) ),
3245                               mul(3,mul(a,a)), neg(mul(sub(b,c),sub(b,c))),
3246                               neg(mul(2,a,add(b,c,-2))), 1/4 ),
3247                      ck[k-1] );
3248
3249        var t2 = mul( add(k,neg(a),b,neg(c),-1/2), add(k,neg(a),neg(b),c,-1/2),
3250                      add(k,neg(a),b,c,-5/2), ck[k-2] );
3251
3252        ck[k] = div( sub( t1, t2 ), 2*k );
3253
3254        wk = w(k);
3255
3256        s1 = add( s1, mul( powMinusI, wk ) );
3257        s2 = add( s2, mul( powPlusI, wk ) );
3258
3259      }
3260
3261      var exponent = add( mul(pi,p), mul(2,sqrt(neg(x))) );
3262
3263      var u1 = exp( mul( plusI, exponent ) );
3264      var u2 = exp( mul( minusI, exponent ) );
3265
3266      var s = add( mul(u1,s1), mul(u2,s2) );
3267
3268      var t1 = mul( 1/(2*sqrt(pi)), inv(gamma(a)), pow(neg(x),p), s );
3269
3270      var t2 = mul( inv(gamma(sub(b,a))), inv(gamma(sub(c,a))), pow(neg(x),neg(a)),
3271                    hypergeometricSeries( [ a, add(a,neg(b),1), add(a,neg(c),1) ], [], inv(x), true ) );
3272
3273      return mul( gamma(b), gamma(c), add( t1, t2 ) );
3274
3275    }
3276
3277    return hypergeometricSeries( [a], [b,c], x, true );
3278
3279  } else {
3280
3281    // asymptotic form is complex
3282    if ( Math.abs(x) > useAsymptotic ) return hypergeometric1F2( a, b, c, complex(x) ).re;
3283
3284    return hypergeometricSeries( [a], [b,c], x );
3285
3286  }
3287
3288}
3289
3290
3291// convenience function for less-used hypergeometrics
3292// accessing array slower than local variables
3293// for loops faster than forEach or reduce
3294
3295function hypergeometricSeries( A, B, x, complexArguments=false, tolerance=1e-10 ) {
3296
3297  if ( complexArguments ) {
3298
3299    var s = complex(1);
3300    var p = complex(1);
3301    var i = 1;
3302
3303    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
3304
3305      for ( var j = 0 ; j < A.length ; j++ ) {
3306        p = mul( p, A[j] );
3307        A[j] = add( A[j], 1 );
3308      }
3309
3310      for ( var j = 0 ; j < B.length ; j++ ) {
3311        p = div( p, B[j] );
3312        B[j] = add( B[j], 1 );
3313      }
3314
3315      p = mul( p, x, 1/i );
3316      s = add( s, p );
3317      i++;
3318
3319    }
3320
3321    return s;
3322
3323  } else {
3324
3325    var s = 1;
3326    var p = 1;
3327    var i = 1;
3328
3329    while ( Math.abs(p) > tolerance ) {
3330
3331      for ( var j = 0 ; j < A.length ; j++ ) {
3332        p *= A[j];
3333        A[j]++;
3334      }
3335
3336      for ( var j = 0 ; j < B.length ; j++ ) {
3337        p /= B[j];
3338        B[j]++;
3339      }
3340
3341      p *= x / i;
3342      s += p;
3343      i++;
3344
3345    }
3346
3347    return s;
3348
3349  }
3350
3351}
3352
3353
3354function hypergeometricPFQ( A, B, x ) {
3355
3356  // dlmf.nist.gov/16.11 for general transformations
3357
3358  // for B.length > A.length terms can get very large
3359  // roundoff errors even for formally convergent series
3360
3361  if ( abs(x) > 1 ) throw Error( 'General hypergeometric argument currently restricted' );
3362

3363  // check for complex parameters
3364  var cp = false;
3365  A.forEach( a => cp = cp || isComplex(a) );
3366  B.forEach( b => cp = cp || isComplex(b) );
3367
3368  if ( cp || isComplex(x) ) return hypergeometricSeries( A, B, x, true );
3369
3370  return hypergeometricSeries( A, B, x );
3371
3372}
3373
3374
3375function exp( x ) {
3376
3377  if ( isArbitrary(x) ) {
3378
3379    if ( isComplex(x) ) {
3380
3381      var expXre = exp(x.re);
3382
3383      return { re: mul( expXre, cos(x.im) ),
3384               im: mul( expXre, sin(x.im) ) };
3385
3386    }
3387
3388    var m = div( x, ln10 ) / precisionScale;
3389
3390    x -= m * ln10;
3391
3392    // direct sum faster than function inversion
3393    var s = arb1;
3394    var p = arb1;
3395    var i = arb1;
3396
3397    while ( p !== 0n ) {
3398      p = div( mul( p, x ), i );
3399      s += p;
3400      i += arb1;
3401    }
3402
3403    if ( m > 0n ) s *= 10n**m;
3404    if ( m < 0n ) s /= 10n**-m; // value approaches zero for fixed decimals
3405
3406    // could also return as mantissa/exponent
3407    return s;
3408
3409  }
3410
3411  if ( isComplex(x) )
3412
3413    return { re: Math.exp(x.re) * Math.cos(x.im),
3414             im: Math.exp(x.re) * Math.sin(x.im) };
3415
3416  return Math.exp(x);
3417
3418}
3419
3420function logisticSigmoid( x ) { return inv( add( 1, exp(neg(x)) ) ); }
3421
3422
3423function log( x, base ) {
3424
3425  if ( isComplex(x) || isComplex(base) ) {
3426
3427    if ( base === undefined ) return { re: Math.log( abs(x) ), im: arg(x) };
3428
3429    return div( log(x), log(base) );
3430
3431  }
3432
3433  if ( x < 0 ) return log( complex(x), base );
3434
3435  if ( base === undefined ) return Math.log(x);
3436
3437  return Math.log(x) / Math.log(base);
3438
3439}
3440
3441function ln( x ) {
3442
3443  // Brent, Modern Computer Arithmetic, second AGM algorithm
3444
3445  function arbitraryAGM( x, y ) {
3446
3447    var t, u;
3448
3449    if ( isComplex(x) ) {
3450
3451      var maxIter = 20, i = 0;
3452
3453      while( x.re !== y.re || x.im !== y.im ) {
3454        t = x, u = y;
3455        x = div( add( t, u ), arb2 );
3456        y = sqrt( mul( t, u ) );
3457        i++;
3458        if ( i > maxIter ) break; // convergence on complex plane not assured...
3459      }
3460
3461    } else {
3462
3463      while( x !== y ) {
3464        t = x, u = y;
3465        x = div( t + u, arb2 );
3466        y = sqrt( mul( t, u ) );
3467      }
3468
3469    }
3470
3471    return x;
3472
3473  }
3474
3475  function arbitraryTheta2( x ) {
3476
3477    var p = arb1;
3478    var s = arb1;
3479    var i = 1;
3480
3481    if ( isComplex(x) ) {
3482
3483      while ( p.re !== 0n || p.im !== 0n ) {
3484        for ( var j = 0 ; j < 8*i ; j++ ) p = mul( p, x );
3485        s = add( s, p );
3486        i++;
3487      }
3488
3489    } else {
3490
3491      while ( p !== 0n ) {
3492        for ( var j = 0 ; j < 8*i ; j++ ) p = mul( p, x );
3493        s = s + p;
3494        i++;
3495      }
3496
3497    }
3498
3499    return mul( arb2, x, s );
3500
3501  }
3502
3503  function arbitraryTheta3( x ) {
3504
3505    var p = arb1;
3506    var s = 0n;
3507    var i = 1;
3508
3509    if ( isComplex(x) ) {
3510
3511      while ( p.re !== 0n || p.im !== 0n ) {
3512        for ( var j = 0 ; j < 4*(2*i-1) ; j++ ) p = mul( p, x );
3513        s = add( s, p );
3514        i++;
3515      }
3516
3517    } else {
3518
3519      while ( p !== 0n ) {
3520        for ( var j = 0 ; j < 4*(2*i-1) ; j++ ) p = mul( p, x );
3521        s = s + p;
3522        i++;
3523      }
3524
3525    }
3526
3527    return add( arb1, mul( arb2, s ) );
3528
3529  }
3530
3531  function arbitraryLn( x) {
3532
3533    var arb4 = 4n*arb1, arb10 = 10n*arb1;
3534
3535    // convergence near unit circle problematic
3536    // scale argument radially and subtract scaling
3537    // any number will work, convergence faster further out
3538
3539    if ( abs(x) < arb10 ) return sub( arbitraryLn(mul(arb10,x)), ln10 );
3540
3541    x = div( arb1, x );
3542
3543    var t2 = arbitraryTheta2(x);
3544    var t3 = arbitraryTheta3(x);
3545
3546    return div( onePi, mul( arb4, arbitraryAGM( mul(t2,t2), mul(t3,t3) ) ) );
3547
3548  }
3549
3550  if ( isArbitrary(x) ) {
3551
3552    if ( isComplex(x) ) {
3553
3554      if ( x.re === arb1 && x.im === 0n ) return { re: 0n, im: 0n };
3555
3556      if ( x.re === 0n && x.im === 0n )
3557        throw Error( 'Arbitrary natural logarithm singularity' );
3558
3559      // directly set phase on complex axis
3560      if ( x.re === 0n )
3561        if ( x.im > 0n ) return { re: arbitraryLn(x.im), im: halfPi };
3562        else return { re: arbitraryLn(-x.im), im: -halfPi };
3563
3564      var result = arbitraryLn(x);
3565
3566      // adjust imaginary part on left half-plane
3567      if ( x.re < 0n ) {
3568        if ( result.im > 0n ) result.im -= onePi;
3569        else result.im += onePi;
3570      }
3571
3572      return result;
3573
3574    }
3575
3576    if ( x === arb1 ) return 0n;
3577    if ( x === 0n ) throw Error( 'Arbitrary natural logarithm singularity' );
3578
3579    if ( x < 0n ) return { re: arbitraryLn( -x ), im: onePi };
3580
3581    return arbitraryLn(x);
3582
3583  }
3584
3585  return log(x);
3586
3587}
3588
3589
3590function lambertW( k, x, tolerance=1e-10 ) {
3591
3592  if ( arguments.length === 1 ) {
3593    x = k;
3594    k = 0;
3595  }
3596
3597  // restrict to real integers for convenience
3598  if ( !Number.isInteger(k) ) throw Error( 'Unsupported Lambert W index' );
3599
3600  if ( isZero(x) )
3601    if ( k === 0 ) return x;
3602    else throw Error( 'Lambert W pole' );
3603
3604  var expMinusOne = Math.exp(-1);
3605
3606  if ( ( k === 0 || k === -1 ) && abs( add( x, expMinusOne ) ) < tolerance*tolerance )
3607    if ( isComplex(x) ) return complex(-1);
3608    else return -1;
3609
3610  // handle real cases separately, complex otherwise
3611
3612  if ( !isComplex(x) ) {
3613    if ( k === 0 && x > -expMinusOne )
3614      return findRoot( w => w * Math.exp(w) - x, [-1,1000], { tolerance: tolerance } );
3615    if ( k === -1 && x > -expMinusOne && x < 0 )
3616      return findRoot( w => w * Math.exp(w) - x, [-1000,-1], { tolerance: tolerance } );
3617    x = complex(x);
3618  }
3619
3620  // inversion by complex root finding with custom Newton's method
3621  var maxIter = 100;
3622
3623  if ( k === 0 && abs(x) < 1.25 ) {
3624    // unstable along branch cut, add small complex part to avoid error
3625    if ( x.im === 0 &&  x.re < -expMinusOne ) x.im = tolerance;
3626    // based on test page: unstable region jumps between sheets

3627    var w = x.re < .5 ? complex( 0, .5*Math.sign(x.im) ) : complex(0);
3628  } else {
3629    var L = add( log(x), complex(0,2*pi*k) );
3630    var w = sub( L, log(L) );
3631  }
3632
3633  for ( var i = 0; i < maxIter ; i++ ) {
3634     var delta = div( sub( mul(w,exp(w)), x ), mul( exp(w), add(w,1) ) );
3635     w = sub( w, delta );
3636     if ( abs(delta) < tolerance ) return w;
3637  }
3638
3639  throw Error( 'No Lambert W root found at ' + JSON.stringify(x) );
3640
3641}
3642
3643function inverseLambertW( x ) { return mul( x, exp(x) ); }
3644
3645function wrightOmega( x ) {
3646
3647  var unwinding = Math.ceil( ( im(x)/pi - 1 ) / 2 );
3648
3649  return lambertW( unwinding, exp(x) );
3650
3651}
3652
3653
3654function chop( x, tolerance=1e-10 ) {
3655
3656  if ( Array.isArray(x) ) {
3657    var v = vector( x.length );
3658    for ( var i = 0 ; i < x.length ; i++ ) v[i] = chop( x[i] );
3659    return v;
3660  }
3661
3662  if ( isComplex(x) ) return complex( chop(x.re), chop(x.im) );
3663
3664  if ( Math.abs(x) < tolerance ) x = 0;
3665
3666  return x;
3667
3668}
3669
3670function round( x, y ) {
3671
3672  if ( arguments.length === 2 ) return mul( y, round( div(x,y) ) );
3673
3674  if ( isComplex(x) ) return complex( Math.round(x.re), Math.round(x.im) );
3675
3676  return Math.round(x);
3677
3678}
3679
3680function ceiling( x ) {
3681
3682  if ( isComplex(x) ) return complex( Math.ceil(x.re), Math.ceil(x.im) );
3683
3684  return Math.ceil(x);
3685
3686}
3687
3688function floor( x ) {
3689
3690  if ( isComplex(x) ) return complex( Math.floor(x.re), Math.floor(x.im) );
3691
3692  return Math.floor(x);
3693
3694}
3695
3696function sign( x ) {
3697
3698  if ( isZero(x) ) return x;
3699
3700  return div( x, abs(x) );
3701
3702}
3703
3704function integerPart( x ) {
3705
3706  if ( isComplex(x) ) return complex( Math.trunc(x.re), Math.trunc(x.im) );
3707
3708  return Math.trunc(x);
3709
3710}
3711
3712function fractionalPart( x ) { return sub( x, integerPart(x) ); }
3713
3714
3715function random( x, y ) {
3716
3717  if ( arguments.length === 0 ) return Math.random();
3718
3719  if ( arguments.length === 1 )
3720    if ( isComplex(x) )
3721      return { re: x.re * Math.random(), im: x.im * Math.random() };
3722    else return x * Math.random();
3723
3724  if ( isComplex(x) || isComplex(y) ) {
3725    if ( !isComplex(x) ) x = complex(x);
3726    if ( !isComplex(y) ) y = complex(y);
3727    return { re: random( x.re, y.re ), im: random( x.im, y.im ) };
3728  }
3729
3730  return Math.abs(x-y) * Math.random() + ( x > y ? y : x );
3731
3732}
3733
3734
3735function kronecker( i, j ) {
3736
3737  if ( arguments.length === 2 ) {
3738
3739    if ( isComplex(i) || isComplex(j) ) {
3740
3741      if ( !isComplex(i) ) i = complex(i);
3742      if ( !isComplex(j) ) j = complex(j);
3743
3744      return kronecker( i.re, j.re) * kronecker( i.im, j.im );
3745
3746    }
3747
3748    return i === j ? 1 : 0;
3749
3750  }
3751
3752  var result = kronecker( i, j );
3753
3754  for ( var k = 2 ; k < arguments.length ; k++ )
3755    result *= kronecker( i, arguments[k] );
3756
3757  return result;
3758
3759}
3760
3761
3762function piecewise() {
3763
3764  var pieces = arguments;
3765
3766  return function( x ) {
3767
3768    for ( var i = 0 ; i < pieces.length ; i++ ) {
3769      var domain = pieces[i][1];
3770      if ( x >= domain[0] && x <= domain[1] ) return pieces[i][0](x);
3771    }
3772
3773    return 0;
3774
3775  }
3776
3777}
3778
3779
3780function hermite( n, x ) {
3781
3782  function coefficients( n ) {
3783
3784    if ( isComplex(n) ) var returnComplex = true, n = n.re;
3785
3786    var minus2 = [ 1 ];
3787    var minus1 = [ 2, 0 ];
3788    var t, current;
3789

3790    if ( n === 0 ) return minus2;
3791    if ( n === 1 ) return minus1;
3792
3793    for ( var i = 2 ; i <= n ; i++ ) {
3794      current = [];
3795      t = minus1.slice();
3796      t.push( 0 );
3797      minus2.unshift( 0, 0 );
3798      for ( var k = 0 ; k < t.length ; k++ )
3799        current.push( 2*t[k] - 2*(i-1)*minus2[k] );
3800      minus2 = minus1;
3801      minus1 = current;
3802    }
3803
3804    if ( returnComplex ) {
3805      current.forEach( (e,i,a) => a[i] = complex(e) );
3806      return current;
3807    }
3808
3809    return current;
3810
3811  }
3812
3813  if ( isComplex(n) || isComplex(x) ) {
3814
3815    if ( isPositiveIntegerOrZero(n) ) return polynomial( x, coefficients(n) );
3816
3817    var a = div( n, -2 );
3818    var b = div( sub(1,n), 2 );
3819
3820    var s = sub( div( hypergeometric1F1( a, .5, pow(x,2) ), gamma( b ) ),
3821                 mul( 2, x, div( hypergeometric1F1( b, 1.5, pow(x,2) ), gamma( a ) ) ) );
3822
3823    return mul( pow(2,n), sqrt(pi), s );
3824
3825  }
3826
3827  if ( isPositiveIntegerOrZero(n) ) return polynomial( x, coefficients(n) );
3828
3829  var s = hypergeometric1F1( -n/2, .5, x**2 ) / gamma( (1-n)/2 )
3830          - 2 * x * hypergeometric1F1( (1-n)/2, 1.5, x**2 ) / gamma( -n/2 );
3831
3832  return 2**n * sqrt(pi) * s;
3833
3834}
3835
3836
3837function laguerre( n, a, x ) {
3838
3839  // explict recursion unnecessary: hypergeometric series handles integers
3840
3841  if ( arguments.length === 2 ) {
3842    x = a;
3843    return hypergeometric1F1( neg(n), 1, x );
3844  }
3845
3846  return mul( binomial( add(n,a), n ), hypergeometric1F1( neg(n), add(a,1), x ) ); 
3847
3848}
3849
3850
3851function chebyshevT( n, x ) {
3852
3853  return cos( mul( n, arccos(x) ) );
3854
3855}
3856
3857function chebyshevU( n, x ) {
3858
3859  return div( sin( mul( add(n,1), arccos(x) ) ), sin( arccos(x) ) );
3860
3861}
3862
3863
3864function legendreP( l, m, x, renormalized=false ) {
3865
3866  if ( arguments.length === 2 ) {
3867    x = m;
3868    return hypergeometric2F1( neg(l), add(l,1), 1, div(sub(1,x),2) );
3869  }
3870
3871  if ( isZero(m) ) return legendreP( l, x );
3872
3873  if ( Number.isInteger(l) && Number.isInteger(m) && Math.abs(x) <= 1 ) {
3874
3875    var mm = Math.abs(m);
3876    if ( mm > l ) throw Error( 'Invalid spherical harmonic indices' );
3877
3878    if ( !renormalized ) {
3879      var norm = 1;
3880      for ( var i = l-m+1 ; i <= l+m ; i++ ) norm *= i;
3881      norm = Math.sqrt( 4 * pi * norm / (2*l+1) );
3882    }
3883
3884    var legendre1 = (-1)**mm * Math.sqrt( (2*mm+1) / 4 / pi / factorial(2*mm) )
3885                    * factorial2( 2*mm-1 ) * ( 1 - x*x )**(mm/2);
3886
3887    if ( mm === l ) 
3888      if ( renormalized ) return legendre1;
3889      else return norm * legendre1;
3890
3891    var ll = mm + 1;
3892    var factor1 = Math.sqrt( 2*mm+3 );
3893    var legendre2 = factor1 * x * legendre1;
3894
3895    if ( ll === l )
3896      if ( renormalized ) return legendre2;
3897      else return norm * legendre2;
3898
3899    while ( ll < l ) {
3900      ll++
3901      var factor2 = Math.sqrt( ( 4*ll*ll - 1 ) / ( ll*ll - mm*mm ) );
3902      var legendre3 = factor2 * ( x*legendre2 - legendre1/factor1 );
3903      legendre1 = legendre2;
3904      legendre2 = legendre3;
3905      factor1 = factor2;
3906    }
3907
3908    if ( renormalized ) return legendre3;
3909    else return norm * legendre3;
3910
3911  }
3912
3913  // dlmf.nist.gov/14.3.5
3914  if ( isPositiveInteger(m) )
3915    return mul( pow(-1,m), inv( gamma( add(m,1) ) ),
3916                gamma( add(l,m,1) ), inv( gamma( add(l,neg(m),1) ) ),
3917                pow( add(1,x), div(m,-2) ), pow( sub(1,x), div(m,2) ),
3918                hypergeometric2F1( neg(l), add(l,1), add(m,1), div(sub(1,x),2) ) );
3919
3920  return mul( inv( gamma( sub(1,m) ) ),
3921              pow( add(1,x), div(m,2) ), pow( sub(1,x), div(m,-2) ),
3922              hypergeometric2F1( neg(l), add(l,1), sub(1,m), div(sub(1,x),2) ) );
3923
3924}
3925
3926function sphericalHarmonic( l, m, theta, phi, condonShortley=true ) {
3927
3928  var norm = sqrt( mul( div( add(mul(2,l),1), 4*pi ),
3929                        div( factorial(sub(l,m)), factorial(add(l,m)) ) ) );
3930
3931  var cs = condonShortley ? pow(-1,m) : 1;
3932
3933  return mul( cs, norm, legendreP( l, m, cos(theta) ), exp( complex(0,m*phi) ) );
3934
3935}
3936
3937function legendreQ( l, m, x ) {
3938
3939  if ( arguments.length === 2 ) {
3940    x = m;
3941    m = 0;
3942  }
3943
3944  function difference( t ) {
3945
3946    var t1 = mul( cos( mul(pi,t) ), legendreP(l,t,x) );
3947
3948    var t2 = mul( gamma( add(l,t,1) ), inv( gamma( add(l,neg(t),1) ) ), legendreP(l,neg(t),x) );
3949
3950    return sub( t1, t2 );
3951
3952  }
3953
3954  // l'Hopital's rule decent for small m, more accurate might be
3955  // functions.wolfram.com/HypergeometricFunctions/LegendreQ2General/26/01/02/0005/
3956
3957  if ( isInteger(m) ) {
3958
3959    // legendreP is pure real outside unit circle for even integers
3960   if ( abs(x) > 1 && !isComplex(m) && !( m & 1 ) ) m = complex(m);
3961
3962    return mul( .5, pow(-1,m), diff( t => difference(t), m ) );
3963
3964  }
3965
3966  return mul( pi/2, inv( sin(mul(pi,m)) ), difference(m) );
3967
3968}
3969
3970
3971function jacobiP( n, a, b, x ) {
3972
3973  return mul( binomial( add(n,a), n ),
3974              hypergeometric2F1( neg(n), add(n,a,b,1), add(a,1), div(sub(1,x),2) ) );
3975
3976}
3977
3978
3979// This file contains proprietary functions defined on analyticphysics.com
3980// Before each is the title of a presentation describing the function
3981
3982
3983// A Generalized Lambert Function of Two Arguments
3984
3985function doubleLambert( n, x, y, tolerance=1e-10 ) {
3986
3987  if ( arguments.length === 2 ) {
3988    y = x;
3989    x = n;
3990    n = 0;
3991  }
3992
3993  if ( isZero(y) ) return lambertW(n,x);
3994  if ( isZero(x) ) return neg( lambertW( -n, neg(y) ) );
3995
3996  function asymptotic( n, x, y ) {
3997
3998    return add( log(sqrt(x)), neg(log(sqrt(neg(y)))), complex(0,n*pi) );
3999
4000  }
4001
4002  function test( n, x, y ) {
4003
4004    return div( asymptotic(n,x,y), add( 1, mul(2*(-1)**n,sqrt(x),sqrt(neg(y))) ) );
4005
4006  }
4007
4008  function start( n, x, y, tolerance=1e-3 ) {
4009
4010    var a = asymptotic(n,x,y), w1, w2;
4011
4012    var testValue = .9;
4013
4014    if ( abs(test(n,x,y)) < testValue ) return a;
4015
4016    if ( n === 0 ) console.log( 'Using Lambert W on principal branch' );
4017
4018    if ( n & 1 ) {
4019
4020      if ( n > 0 ) {
4021        w1 = lambertW( (n+1)/2, x, tolerance );
4022        w2 = neg( lambertW( -(n+1)/2, neg(y), tolerance ) );
4023      } else {
4024        w1 = lambertW( (n-1)/2, x, tolerance );
4025        w2 = neg( lambertW( -(n-1)/2, neg(y), tolerance ) );
4026      }
4027
4028    } else {
4029
4030      w1 = lambertW( n/2, x, tolerance );
4031      w2 = neg( lambertW( -n/2, neg(y), tolerance ) );
4032
4033    }
4034
4035    if ( Math.abs( w1.im - a.im ) < Math.abs( w2.im - a.im ) )
4036      return w1;
4037    else
4038      return w2;
4039
4040  }
4041
4042  var maxIter = 100;
4043  var root = start(n,x,y);
4044
4045  for ( var i = 0; i < maxIter ; i++ ) {
4046    var N = add( mul(x,exp(neg(root))), mul(y,exp(root)), neg(root) );
4047    var D = add( mul(-1,x,exp(neg(root))), mul(y,exp(root)), -1 );
4048    var delta = div( N, D );
4049    root = sub( root, delta );
4050    if ( abs(delta) < tolerance ) return root;
4051  }
4052
4053  throw Error( 'No double Lambert root found for x = ' + JSON.stringify(x)
4054               + ' and y = ' + JSON.stringify(y) );
4055
4056}
4057
4058
4059// complex circular functions
4060
4061function sin( x ) {
4062
4063  if ( isArbitrary(x) ) {
4064
4065    if ( isComplex(x) )
4066
4067      return { re: mul( sin(x.re), cosh(x.im) ),
4068               im: mul( cos(x.re), sinh(x.im) ) };
4069
4070    x = x % twoPi;
4071
4072    // reduce to [-pi/2,pi/2] with successive reductions
4073    if ( x > halfPi ) return sin( onePi - x );
4074    if ( x < -halfPi ) return sin( -onePi - x );
4075
4076    var s = x;
4077    var p = x;
4078    var i = arb2;
4079
4080    while ( p !== 0n ) {
4081      p = div( mul( p, -arb1, x, x ), mul( i, i + arb1 ) );
4082      s += p;
4083      i += arb2;
4084    }
4085
4086    return s;
4087
4088  }
4089
4090  if ( isComplex(x) )
4091
4092    return { re: Math.sin(x.re) * Math.cosh(x.im),
4093             im: Math.cos(x.re) * Math.sinh(x.im) };
4094
4095  return Math.sin(x);
4096
4097}
4098
4099function cos( x ) {
4100
4101  if ( isArbitrary(x) ) {
4102
4103    if ( isComplex(x) )
4104
4105      return { re: mul( cos(x.re), cosh(x.im) ),
4106               im: -mul( sin(x.re), sinh(x.im) ) };
4107
4108    x = x % twoPi;
4109
4110    // reduce to [-pi/2,pi/2] with successive reductions
4111    if ( x > halfPi ) return -cos( onePi - x );
4112    if ( x < -halfPi ) return -cos( -onePi - x );
4113
4114    var s = arb1;
4115    var p = arb1;
4116    var i = arb1;
4117
4118    while ( p !== 0n ) {
4119      p = div( mul( p, -arb1, x, x ), mul( i, i + arb1 ) );
4120      s += p;
4121      i += arb2;
4122    }
4123
4124    return s;
4125
4126  }
4127
4128  if ( isComplex(x) )
4129

4130    return { re: Math.cos(x.re) * Math.cosh(x.im),
4131             im: -Math.sin(x.re) * Math.sinh(x.im) };
4132
4133  return Math.cos(x);
4134
4135}
4136
4137function tan( x ) {
4138
4139  if ( isComplex(x) ) return div( sin(x), cos(x) );
4140
4141  return Math.tan(x);
4142
4143 }
4144
4145function cot( x ) {
4146
4147  if ( isComplex(x) ) return div( cos(x), sin(x) );
4148
4149  return 1 / Math.tan(x);
4150
4151}
4152
4153function sec( x ) {
4154
4155  if ( isComplex(x) ) return div( 1, cos(x) );
4156
4157  return 1 / Math.cos(x);
4158
4159}
4160
4161function csc( x ) {
4162
4163  if ( isComplex(x) ) return div( 1, sin(x) );
4164
4165  return 1 / Math.sin(x);
4166
4167}
4168
4169
4170// inverse circular functions
4171
4172function arcsin( x ) {
4173
4174  if ( isComplex(x) ) {
4175
4176    var s = sqrt( sub( 1, mul( x, x ) ) );
4177    s = add( mul( complex(0,1), x ), s ); 
4178    return mul( complex(0,-1), log( s ) );
4179
4180  }
4181
4182  if ( Math.abs(x) <= 1 ) return Math.asin(x);
4183
4184  return arcsin( complex(x) );
4185
4186}
4187
4188function arccos( x ) {
4189
4190  if ( isComplex(x) ) {
4191
4192    return sub( pi/2, arcsin(x) );
4193
4194  }
4195
4196  if ( Math.abs(x) <= 1 ) return Math.acos(x);
4197
4198  return arccos( complex(x) );
4199
4200}
4201
4202function arctan( x ) {
4203
4204  if ( isComplex(x) ) {
4205
4206    var s = sub( log( sub( 1, mul( complex(0,1), x ) ) ),
4207                 log( add( 1, mul( complex(0,1), x ) ) ) );
4208    return mul( complex(0,.5), s );
4209
4210  }
4211
4212  return Math.atan(x);
4213
4214}
4215
4216function arccot( x ) {
4217
4218  if ( isComplex(x) ) return arctan( div( 1, x ) );
4219
4220  return Math.atan( 1/x );
4221
4222}
4223
4224function arcsec( x ) {
4225
4226  if ( isComplex(x) ) return arccos( div( 1, x ) );
4227
4228  if ( Math.abs(x) >= 1 ) return Math.acos( 1/x );
4229
4230  return arcsec( complex(x) );
4231
4232}
4233
4234function arccsc( x ) {
4235
4236  if ( isComplex(x) ) return arcsin( div( 1, x ) );
4237
4238  if ( Math.abs(x) >= 1 ) return Math.asin( 1/x );
4239
4240  return arccsc( complex(x) );
4241
4242}
4243
4244
4245// complex hyperbolic functions
4246
4247function sinh( x ) {
4248
4249  if ( isArbitrary(x) ) return div( sub( exp(x), exp(neg(x)) ), arb2 );
4250
4251  if ( isComplex(x) )
4252
4253    return { re: Math.sinh(x.re) * Math.cos(x.im),
4254             im: Math.cosh(x.re) * Math.sin(x.im) };
4255
4256  return Math.sinh(x);
4257
4258}
4259
4260function cosh( x ) {
4261
4262  if ( isArbitrary(x) ) return div( add( exp(x), exp(neg(x)) ), arb2 );
4263
4264  if ( isComplex(x) )
4265

4266    return { re: Math.cosh(x.re) * Math.cos(x.im),
4267             im: Math.sinh(x.re) * Math.sin(x.im) };
4268
4269  return Math.cosh(x);
4270
4271}
4272
4273function tanh( x ) {
4274
4275  if ( isComplex(x) ) return div( sinh(x), cosh(x) );
4276
4277  return Math.tanh(x);
4278
4279}
4280
4281function coth( x ) {
4282
4283  if ( isComplex(x) ) return div( cosh(x), sinh(x) );
4284
4285  return 1 / Math.tanh(x);
4286
4287}
4288
4289function sech( x ) {
4290
4291  if ( isComplex(x) ) return div( 1, cosh(x) );
4292
4293  return 1 / Math.cosh(x);
4294
4295}
4296
4297function csch( x ) {
4298
4299  if ( isComplex(x) ) return div( 1, sinh(x) );
4300
4301  return 1 / Math.sinh(x);
4302
4303}
4304
4305
4306// inverse hyperbolic functions
4307
4308function arcsinh( x ) {
4309
4310  if ( isComplex(x) ) {
4311
4312    var s = sqrt( add( mul( x, x ), 1 ) );
4313    s = add( x, s );
4314    return log( s );
4315
4316  }
4317
4318  return Math.asinh(x);
4319
4320}
4321
4322function arccosh( x ) {
4323
4324  if ( isComplex(x) ) {
4325
4326    var s = mul( sqrt( add( x, 1 ) ), sqrt( sub( x, 1 ) ) );
4327    s = add( x, s ); 
4328    return log( s );
4329
4330  }
4331
4332  if ( x >= 1 ) return Math.acosh(x);
4333
4334  return arccosh( complex(x) );
4335
4336}
4337
4338function arctanh( x ) {
4339
4340  if ( isComplex(x) ) {
4341
4342    var s = sub( log( add( 1, x ) ), log( sub( 1, x ) ) );
4343    return mul( .5, s );
4344
4345  }
4346
4347  if ( Math.abs(x) <= 1 ) return Math.atanh(x);
4348
4349  return arctanh( complex(x) );
4350
4351}
4352
4353function arccoth( x ) {
4354
4355  if ( isComplex(x) ) {
4356
4357    if ( x.re === 0 && x.im === 0 ) throw Error( 'Indeterminate arccoth value' );
4358
4359    return arctanh( div( 1, x ) );
4360
4361  }
4362
4363  if ( Math.abs(x) > 1 ) return Math.atanh( 1/x );
4364
4365  return arccoth( complex(x) );
4366
4367}
4368
4369function arcsech( x ) {
4370
4371  if ( isComplex(x) ) {
4372
4373    if ( x.re === 0 && x.im === 0 ) throw Error( 'Indeterminate arcsech value' );
4374
4375    // adjust for branch cut along negative axis
4376    if ( x.im === 0 ) x.im = -Number.MIN_VALUE;
4377
4378    return arccosh( div( 1, x ) );
4379
4380  }
4381
4382  if ( x > 0 && x < 1 ) return Math.acosh( 1/x );
4383
4384  return arcsech( complex(x) );
4385
4386}
4387
4388function arccsch( x ) {
4389
4390  if ( isComplex(x) ) {
4391
4392    return arcsinh( div( 1, x ) );
4393
4394  }
4395
4396  return Math.asinh( 1/x );
4397
4398}
4399
4400
4401// miscellaneous
4402
4403function sinc( x ) {
4404
4405  if ( isComplex(x) ) {
4406
4407    if ( x.re === 0 && x.im === 0 ) return complex(1);
4408
4409    return div( sin(x), x );
4410
4411  }
4412
4413  if ( x === 0 ) return 1;
4414
4415  return Math.sin(x) / x;
4416
4417}
4418
4419function haversine( x ) { return div( sub( 1, cos(x) ), 2 ); }
4420
4421function inverseHaversine( x ) { return arccos( sub( 1, mul(2,x) ) ); }
4422
4423
4424// analyticphysics.com / The Complex Gudermannian Function
4425
4426function gudermannian( x ) { return mul( 2, arctan( tanh( div(x,2) ) ) ); }
4427
4428function inverseGudermannian( x ) { return mul( 2, arctanh( tan( div(x,2) ) ) ); }
4429
4430
4431function zeta( x, tolerance=1e-10 ) {
4432
4433  if ( isUnity(x) ) throw Error( 'Riemann zeta pole' );
4434
4435  // functional equation dlmf.nist.gov/25.4.2 - connects both half planes
4436  if ( x < 0 || x.re < 0 )
4437    return mul( pow(2,x), pow(pi,sub(x,1)), sin( mul(pi/2,x) ), gamma( sub(1,x) ), zeta( sub(1,x) ) );
4438
4439  // direct summation more accurate in right-hand plane
4440  var directSummation = 5;
4441
4442  if ( x > directSummation || x.re > directSummation ) {
4443
4444    if ( isComplex(x) ) {
4445
4446      var s = complex(1);
4447      var p = complex(1);
4448      var i = 2;
4449
4450      while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance*tolerance ) {
4451        p = pow( i, neg(x) );
4452        s = add( s, p );
4453        i++;
4454      }
4455
4456      return s;
4457
4458    } else {
4459
4460      var s = 1;
4461      var p = 1;
4462      var i = 2;
4463
4464      while ( Math.abs(p) > tolerance ) {
4465        p = 1 / i**x;
4466        s += p;
4467        i++;
4468      }
4469
4470      return s;
4471
4472    }
4473
4474  }
4475
4476  // Borwein, Efficient Algorithm
4477
4478  var n = 14; // from error bound for tolerance
4479
4480  if ( isComplex(x) && x.im !== 0 )
4481    n = Math.max( n, Math.ceil( log( 2 / abs(gamma(x)) / tolerance ) / log( 3 + sqrt(8) ) ) );
4482
4483  var d = [ 1 ], prod = n;
4484  for ( var i = 1 ; i <= n ; i++ ) {
4485    d.push( d[i-1] + n * prod / factorial( 2*i ) * 4**i );
4486    // manually evaluate factorial( n+i-1 ) / factorial( n-i ) to avoid overflow
4487    prod *= (n+i)*(n-i);
4488  }
4489
4490  if ( isComplex(x) ) {
4491
4492    var s = summation( k => div( (-1)**k * ( d[k] - d[n] ), pow( k+1, x ) ), [0,n-1] );
4493
4494    return div( div( s, -d[n] ), sub( 1, pow( 2, sub(1,x) ) ) );
4495
4496  } else {
4497
4498    var s = summation( k => (-1)**k * ( d[k] - d[n] ) / (k+1)**x, [0,n-1] );
4499
4500    return -s / d[n] / ( 1 - 2**(1-x) );
4501
4502  }
4503
4504}
4505
4506function dirichletEta( x ) { return mul( zeta(x), sub( 1, pow( 2, sub(1,x) ) ) ); }
4507
4508function riemannXi( x ) {
4509
4510  if ( isZero(x) || isUnity(x) ) return isComplex(x) ? complex(.5) : .5;
4511
4512  var half = div( x, 2 );
4513
4514  return mul( half, sub(x,1), pow( pi, neg(half) ), gamma(half), zeta(x) );
4515
4516}
4517
4518
4519function bernoulli( n, x ) {
4520
4521  if ( arguments.length === 2 ) {
4522
4523    if ( isZero(n) ) return isComplex(n) || isComplex(x) ? complex(1) : 1;
4524
4525    // avoid Hurwitz zeta parameter poles
4526    if (  isNegativeIntegerOrZero(x) ) return complexAverage( x

4526=> bernoulli(n,x), x );
4527
4528    return mul( neg(n), hurwitzZeta( sub(1,n), x ) );
4529
4530  }
4531
4532  if ( Number.isInteger(n) && n >= 0 ) {
4533
4534    if ( n === 0 ) return 1;
4535
4536    if ( n === 1 ) return -.5;
4537
4538    if ( n & 1 ) return 0;
4539
4540    var m = n/2;
4541    if ( m <= bernoulli2nN.length )
4542      return arbitrary( div( bernoulli2nN[m], bernoulli2nD[m] ) );
4543
4544    return -n * zeta(1-n);
4545
4546  }
4547
4548  if ( isPositiveIntegerOrZero(n) ) return complex( bernoulli( n.re ) );
4549
4550  // generalized Bernoulli number
4551  console.log( 'Returning generalized Bernoulli number' );
4552  return bernoulli( n, 0 );
4553
4554}
4555
4556function harmonic( n, r ) {
4557
4558  if ( arguments.length === 2 ) {
4559
4560    if ( !Number.isInteger(n) || isComplex(r) )
4561      return sub( zeta(r), hurwitzZeta( r, add(n,1) ) );
4562
4563    return summation( i => 1/i**r, [1,n] );
4564
4565  }
4566
4567  if ( !Number.isInteger(n) ) return add( digamma(add(n,1)), eulerGamma );
4568
4569  if ( n > 1e3 ) return log(n) + eulerGamma + 1/2/n - 1/12/n**2;
4570
4571  return summation( i => 1/i, [1,n] );
4572
4573}
4574
4575
4576function hurwitzZeta( x, a, tolerance=1e-10 ) {
4577
4578  if ( isUnity(x) ) throw Error( 'Hurwitz zeta pole' );
4579
4580  if ( isNegativeIntegerOrZero(a) ) throw Error( 'Hurwitz zeta parameter pole' );
4581
4582  if ( isComplex(x) || isComplex(a) ) {
4583
4584    if ( !isComplex(x) ) x = complex(x);
4585    if ( !isComplex(a) ) a = complex(a);
4586
4587    // direct summation more accurate than dlmf.nist.gov/25.11.4 for positive a
4588
4589    if ( a.re < 0 ) {
4590      var m = -Math.floor(a.re);
4591      return add( hurwitzZeta(x,add(a,m)), summation( i => pow( add(a,i), neg(x) ), [0,m-1] ) );
4592    }
4593

4594    // Johansson arxiv.org/abs/1309.2877
4595
4596    var n = 15; // recommendation of Vepstas, Efficient Algorithm, p.12
4597
4598    // Euler-Maclaurin has differences of large values in left-hand plane
4599    var useArbitrary = x.re < 0;
4600
4601    if ( useArbitrary ) {
4602
4603      setPrecisionScale( 20 - Math.round(x.re) );
4604
4605      x = arbitrary(x), a = arbitrary(a);
4606      var arbN = arbitrary(n), arb3 = arbitrary(3);
4607
4608      // at first sight including the factor (a+N)**x in this sum
4609      //   appears to be unnecessary calculational overhead
4610      // turns out to be more accurate by keeping it on a comparable
4611      //   level of numerical precision as the other two terms
4612
4613      var S = 0n;
4614      for ( var i = 0 ; i < n ; i++ )
4615        S = add( S, pow( div( add(a,arbN), add(a,arbitrary(i)) ), x ) );
4616
4617      var I = div( add(a,arbN), sub(x,arb1) );
4618
4619      var p = div( mul( arb1/2n, x ), add(a,arbN) );
4620      var b = div( bernoulli2nN[1], bernoulli2nD[1] );
4621      var t = mul( b, p );
4622      var i = arb2;
4623      var j = 2;
4624
4625      while ( p.re !== 0n || p.im !== 0n ) {
4626        if ( j === bernoulli2nN.length ) break;
4627        p = div( mul( p, add( x, 2n*i - arb2 ), add( x, 2n*i - arb3 ) ),
4628                 mul( 2n*i, 2n*i - arb1, add(a,arbN), add(a,arbN) ) );
4629        b = div( bernoulli2nN[j], bernoulli2nD[j] );
4630        t = add( t, mul( b, p ) );
4631        i += arb1;
4632        j++;
4633      }
4634
4635      var T = add( arb1/2n, t );
4636
4637      var result = arbitrary( mul( add( S, I, T ), pow( add(a,arbN), mul(-arb1,x) ) ) );
4638
4639      resetPrecisionScale();
4640
4641      return result;
4642
4643    }
4644
4645    var S = summation( i => pow( add(a,i), neg(x) ), [0,n-1] );
4646
4647    var I = div( pow( add(a,n), sub(1,x) ), sub(x,1) );
4648
4649    var p = mul( .5, x, inv(add(a,n)) );
4650    var t = mul( bernoulli(2), p );
4651    var i = 2;
4652
4653    // converges rather quickly
4654    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
4655      p = div( mul( p, add( x, 2*i - 2 ), add( x, 2*i - 3 ) ),
4656               mul( 2*i * (2*i-1), add(a,n), add(a,n) ) );
4657      t = add( t, mul( bernoulli(2*i), p ) );
4658      i++;
4659    }
4660
4661    var T = div( add( .5, t ), pow( add(a,n), x ) );
4662
4663    return add( S, I, T );
4664
4665  } else {
4666
4667    if ( a < 0 ) return hurwitzZeta( x, complex(a) );
4668
4669    // direct summation more accurate than dlmf.nist.gov/25.11.4
4670
4671    // Euler-Maclaurin has differences of large values in left-hand plane
4672    // switch to different summation: dlmf.nist.gov/25.11.9
4673
4674    var switchForms = -5;
4675
4676    if ( x < switchForms ) {
4677
4678      x = 1 - x;
4679      var t = Math.cos( pi*x/2 - 2*pi*a );
4680      var s = t;
4681      var i = 1;
4682
4683      while ( Math.abs(t) > tolerance ) {
4684        i++;
4685        t = Math.cos( pi*x/2 - 2*i*pi*a ) / i**x;
4686        s += t;
4687      }
4688
4689      return 2 * gamma(x) / (2*pi)**x * s;
4690
4691    }
4692

4693    // Johansson arxiv.org/abs/1309.2877
4694
4695    var n = 15; // recommendation of Vepstas, Efficient Algorithm, p.12
4696
4697    var S = summation( i => 1 / (a+i)**x, [0,n-1] );
4698
4699    var I = (a+n)**(1-x) / (x-1);
4700
4701    var p = x / 2 / (a+n);
4702    var t = bernoulli(2) * p;
4703    var i = 2;
4704
4705    // converges rather quickly
4706    while ( Math.abs(p) > tolerance ) {
4707      p *= ( x + 2*i - 2 ) * ( x + 2*i - 3 ) / ( 2*i * (2*i-1) * (a+n)**2 );
4708      t += bernoulli(2*i) * p;
4709      i++;
4710    }
4711
4712    var T = ( .5 + t ) / (a+n)**x;
4713
4714    return S + I + T;
4715
4716  }
4717
4718}
4719
4720
4721function polylog( n, x, tolerance=1e-10 ) {
4722
4723  if ( isEqualTo(x,1) ) return zeta(n);
4724
4725  if ( isEqualTo(x,-1) ) return neg( dirichletEta(n) );
4726
4727  if ( isEqualTo(n,1) ) return neg( log( sub(1,x) ) );
4728
4729  if ( isEqualTo(n,0) ) return div( x, sub(1,x) );
4730
4731  if ( isEqualTo(n,-1) ) return div( x, mul( sub(1,x), sub(1,x) ) );
4732
4733  if ( abs(x) > 1 ) {
4734
4735    var twoPiI = complex(0,2*pi);
4736
4737    if ( isPositiveInteger(n) ) {
4738
4739      // Crandall, Note on Fast Polylogarithm Computation
4740
4741      var t1 = mul( (-1)**n, polylog( n, inv(x) ) );
4742
4743      var t2 = mul( div( pow(twoPiI,n), factorial(n) ), bernoulli( n, div(log(x),twoPiI) ) );
4744
4745      var y = isComplex(x) ? x : complex(x); // just for test
4746      var t3 = y.im < 0 || ( y.im === 0 && y.re >= 1 ) ?
4747               mul( twoPiI, div( pow(log(x),n-1), factorial(n-1) ) ) : 0;      
4748
4749      var result = neg( add( t1, t2, t3 ) );
4750
4751      // real on negative real axis
4752      if ( !isComplex(x) && x < 0 ) return result.re;
4753
4754      return result;
4755
4756    }
4757
4758    var v = sub(1,n);
4759    var I = complex(0,1);
4760    var L = div( log(neg(x)), twoPiI );
4761
4762    var z1 = mul( pow(I,v), hurwitzZeta( v, add(.5,L) ) );
4763    var z2 = mul( pow(I,neg(v)), hurwitzZeta( v, sub(.5,L) ) );
4764
4765    return mul( gamma(v), pow(2*pi,neg(v)), add(z1,z2) )
4766
4767  }
4768
4769  if ( isComplex(n) || isComplex(x) ) {
4770
4771    var s = x;
4772    var p = complex(1);
4773    var i = 2;
4774
4775    while ( Math.abs(p.re) > tolerance || Math.abs(p.im) > tolerance ) {
4776      p = div( pow(x,i), pow(i,n) );
4777      s = add( s, p );
4778      i++;
4779    }
4780
4781    return s;
4782
4783  } else {
4784
4785    var s = x;
4786    var p = 1;
4787    var i = 2;
4788
4789    while ( Math.abs(p) > tolerance ) {
4790      p = x**i / i**n;
4791      s += p;
4792      i++;
4793    }
4794
4795    return s;
4796
4797  }
4798
4799}
4800
4801
4802function ode( f, y, [x0,x1], step=.001, method='runge-kutta' ) {
4803
4804  if ( x1 < x0 ) {
4805    var compare = x => x >= x1;
4806    step *= -1;
4807  } else
4808    var compare = x => x <= x1;
4809
4810  // vectorizing first-order real equation works because +[1] = 1
4811  // for complex case +[C(1)] = NaN, so explicit array references
4812  //    are necessary in the input function
4813
4814  if ( !Array.isArray(y) ) {
4815    var g = f;
4816    f = (x,y) => [ g(x,y) ];
4817    y = [ y ];
4818  }
4819
4820  // preparation for complex system
4821  if ( isComplex(x0) || isComplex(x1) || y.some( e => isComplex(e) )
4822         || f(x0,y).some( e => isComplex(e) ) ) {
4823
4824    if ( !isComplex(x0) ) x0 = complex(x0);
4825
4826    y.forEach( (e,i,a) => { if ( !isComplex(e) ) a[i] = complex(e); } );
4827
4828    if ( f(x0,y).every( e => !isComplex(e) ) )
4829      throw Error( 'All functions must handle complex math' );
4830
4831    var d = sub(x1,x0), absD = abs(d);
4832    step = mul( step, div( d, absD ) );
4833    var steps = Math.trunc( absD / abs(step) ), currentStep = 0;
4834
4835  }
4836
4837  var points = [ [x0].concat(y) ];
4838  var size = y.length;
4839
4840  switch( method ) {
4841
4842    case 'euler':
4843
4844      if ( isComplex(x0) ) {
4845
4846        for ( var x = add(x0,step) ; currentStep < steps ; x = add(x,step) ) {
4847
4848          var k = f(x,y);
4849
4850          for ( var i = 0 ; i < size ; i++ ) y[i] = add( y[i], mul( k[i], step ) );
4851
4852          points.push( [x].concat(y) );
4853
4854          currentStep++;
4855
4856        }
4857
4858        return points;
4859
4860      } else {
4861
4862        for ( var x = x0+step ; compare(x) ; x += step ) {
4863
4864          var k = f(x,y);
4865
4866          for ( var i = 0 ; i < size ; i++ ) y[i] += k[i] * step;
4867
4868          points.push( [x].concat(y) );
4869
4870        }
4871
4872        return points;
4873
4874      }
4875
4876    case 'runge-kutta':
4877
4878      if ( isComplex(x0) ) {
4879
4880        var halfStep = div( step, 2 );
4881
4882        for ( var x = add(x0,step) ; currentStep < steps ; x = add(x,step) ) {
4883
4884          var y1 = [], y2 = [], y3 = [];
4885
4886          var k1 = f(x,y);
4887          for ( var i = 0 ; i < size ; i++ ) y1.push( add( y[i], mul( k1[i], halfStep ) ) );
4888          var k2 = f( add( x, halfStep ), y1 );
4889          for ( var i = 0 ; i < size ; i++ ) y2.push( add( y[i], mul( k2[i], halfStep ) ) );
4890          var k3 = f( add( x, halfStep ), y2 );
4891          for ( var i = 0 ; i < size ; i++ ) y3.push( add( y[i], mul( k3[i], step ) ) );
4892          var k4 = f( add( x, step ), y3 );
4893
4894          for ( var i = 0 ; i < size ; i++ )
4895            y[i] = add( y[i], mul( add( k1[i], mul(2,k2[i]), mul(2,k3[i]), k4[i] ), step, 1/6 ) );
4896
4897          points.push( [x].concat(y) );
4898
4899          currentStep++;
4900
4901        }
4902
4903        return points;
4904
4905      } else {
4906
4907        for ( var x = x0+step ; compare(x) ; x += step ) {
4908
4909          var y1 = [], y2 = [], y3 = [];
4910
4911          var k1 = f(x,y);
4912          for ( var i = 0 ; i < size ; i++ ) y1.push( y[i] + k1[i]*step/2 );
4913          var k2 = f( x+step/2, y1 );
4914          for ( var i = 0 ; i < size ; i++ ) y2.push( y[i] + k2[i]*step/2 );
4915          var k3 = f( x+step/2, y2 );
4916          for ( var i = 0 ; i < size ; i++ ) y3.push( y[i] + k3[i]*step );
4917          var k4 = f( x+step, y3 );
4918
4919          for ( var i = 0 ; i < size ; i++ )
4920            y[i] += ( k1[i] + 2*k2[i] + 2*k3[i] + k4[i] ) * step / 6;
4921
4922          points.push( [x].concat(y) );
4923
4924        }
4925
4926        return points;
4927
4928      }
4929
4930    default:
4931
4932      throw Error( 'Unsupported differential equation solver method' );
4933
4934  }
4935
4936}
4937
4938
4939function diff( f, x, n=1, options={} ) {
4940
4941  if ( !Number.isInteger(n) || n < 0 )
4942    throw Error( 'Only positive integer orders for differentiation' );
4943
4944  if ( isComplex(x) || isComplex(f(x)) ) {
4945
4946    if ( !isComplex(f(x)) ) throw Error( 'Function must handle complex math' );
4947
4948    var absX = abs(x);
4949    var normed = absX === 0 ? complex(1) : div( x, absX );
4950
4951    var real = diff( t => f( mul(normed,t) ).re, absX, n, options );
4952    var imag = diff( t => f( mul(normed,t) ).im, absX, n, options );
4953
4954    return div( complex( real, imag ), pow(normed,n) );
4955
4956  }
4957
4958  var method = 'method' in options ? options.method : 'ridders';
4959  var epsilon = 'epsilon' in options ? options.epsilon : 1e-5;
4960

4961  // central differences have h**2 error but division
4962  //   by h**n increases roundoff error
4963  // step sizes chosen as epsilon**(1/(n+2)) to minimize error
4964
4965  function difference() {
4966
4967    var s = 0;
4968    for ( var i = 0 ; i <= n ; i++ )
4969      s += (-1)**i * binomial(n,i) * f( x + (n-2*i)*h );
4970
4971    return s / (2*h)**n
4972
4973  }
4974
4975  switch( method ) {
4976
4977    case 'naive':
4978
4979      // only accurate for first couple derivatives
4980      epsilon /= 1000;
4981      var h = epsilon**(1/(n+2));
4982      return difference();
4983
4984    case 'ridders':
4985
4986      var h = epsilon**(1/(n+2));
4987      var error = Number.MAX_VALUE;
4988      var maxIter = 10;
4989      var result;
4990
4991      var d = [];
4992      for ( var i = 0 ; i < maxIter ; i++ ) d.push( [] );
4993
4994      // Richardson extrapolation as per C. Ridders
4995      d[0][0] = difference();
4996
4997      for ( var i = 1 ; i < maxIter ; i++ ) {
4998
4999        h /= 2;
5000        d[0][i] = difference();
5001
5002        for ( var j = 1 ; j <= i ; j++ ) {
5003
5004          d[j][i] = ( 4**j * d[j-1][i] - d[j-1][i-1] ) / ( 4**j - 1 );
5005
5006          var delta = Math.max( Math.abs( d[j][i] - d[j-1][i] ),
5007                                Math.abs( d[j][i] - d[j-1][i-1] ) );
5008          if ( delta <= error ) {
5009            error = delta;
5010            result = d[j][i];
5011          }
5012
5013        }
5014
5015        if ( Math.abs( d[i][i] - d[i-1][i-1] ) > error ) break;
5016
5017      }
5018
5019      return result;
5020
5021    default:
5022
5023      throw Error( 'Unsupported differentiation method' );
5024
5025  }
5026
5027}
5028
5029var D = diff;
5030
5031
5032function taylorSeries( f, x0, terms=5 ) {
5033
5034  var c = [ f(x0) ];
5035  for ( var i = 1 ; i < terms ; i++ ) c.push( diff( f, x0, i ) );
5036
5037  return function( x ) {
5038
5039    var s = 0;
5040    for ( var i = 0 ; i < c.length ; i++ )
5041      s = add( s, mul( c[i], pow( sub(x,x0), i ), 1/factorial(i) ) );
5042    return s;
5043
5044  }
5045
5046}
5047
5048
5049function gradient( f, point ) {
5050
5051  if ( f.length !== point.length ) throw Error( 'Gradient point length differs from function' );
5052
5053  var result = [];
5054
5055  for ( var i = 0 ; i < point.length ; i++ ) {
5056
5057    var a = point.slice();
5058
5059    result.push( diff( x => { a[i] = x; return f.apply( null, a ); }, a[i] ) );
5060
5061  }
5062
5063  return result;
5064
5065}
5066
5067function findExtremum( f, point, options={} ) {
5068
5069  if ( !Array.isArray(point) ) point = [ point ];
5070
5071  var sign = options.findMaximum ? 1 : -1;
5072  var tolerance = 'tolerance' in options ? options.tolerance : 1e-10;
5073
5074  var maxIter = 1e4;
5075  var gamma = .01 * sign;
5076  var grad, step, test;
5077
5078  for ( var i = 0 ; i < maxIter ; i++ ) {
5079
5080    grad = gradient( f, point );
5081    test = true;
5082
5083    for ( var j = 0 ; j < point.length ; j++ ) {
5084      step = gamma * grad[j];
5085      point[j] += step;
5086      test = test && step < tolerance;
5087    }
5088
5089    if ( test )
5090     if ( point.length === 1 ) return point[0];
5091     else return point;
5092
5093  }
5094
5095  throw Error( 'No extremum found for tolerance ' + tolerance );
5096
5097}
5098
5099
5100function integrate( f, [a,b], options={} ) {
5101
5102  if ( isComplex(a) || isComplex(b) || isComplex(f(a)) ) {
5103
5104    if ( !isComplex(a) ) a = complex(a);
5105    if ( !isComplex(b) ) b = complex(b);
5106
5107    if ( a.re === b.re && a.im === b.im ) return complex(0);
5108
5109    if ( !isComplex(f(a)) ) throw Error( 'Function must handle complex math' );
5110
5111    function lerp( t ) { return add( mul( sub(b,a), t ), a ); }
5112
5113    var real = integrate( t => f( lerp(t) ).re, [0,1], options );
5114    var imag = integrate( t => f( lerp(t) ).im, [0,1], options );
5115
5116    return mul( sub(b,a), complex( real, imag ) );
5117
5118  }
5119
5120  var method = 'method' in options ? options.method : 'adaptive-simpson';
5121  var tolerance = 'tolerance' in options ? options.tolerance : 1e-10;
5122
5123  if ( a === b ) return 0;
5124
5125  if ( options.avoidEndpoints )
5126    if ( a < b ) { a += tolerance; b -= tolerance; }
5127    else { a -= tolerance; b += tolerance; }
5128
5129  function nextEulerIteration() {
5130
5131      h /= 2;
5132      var x = a + h;
5133      while ( x < b ) {
5134        // only add new function evaluations
5135        s += f(x);
5136        x += 2*h;
5137      }
5138
5139  }
5140
5141  switch( method ) {
5142
5143    case 'euler-maclaurin':
5144
5145      // Euler-Maclaurin summation formula
5146
5147      var maxIter = 50;
5148
5149      var h = ( b - a ) / 2;
5150      var s = ( f(a) + f(b) ) / 2 + f( (a+b)/2 );
5151      var result = h * s;
5152      var previous = result;
5153
5154      for ( var i = 0 ; i < maxIter ; i++ ) {
5155
5156        nextEulerIteration();
5157        result = h * s;
5158        if ( Math.abs( result - previous ) < tolerance * Math.abs(previous) ) return result;
5159        previous = result;
5160
5161      }
5162
5163      throw Error( 'Maximum interations reached' );
5164
5165    case 'romberg':
5166
5167      var error = Number.MAX_VALUE;
5168      var maxIter = 30;
5169
5170      var h = ( b - a ) / 2;
5171      var s = ( f(a) + f(b) ) / 2 + f( (a+b)/2 );
5172      var result = h * s;
5173
5174      var d = [];
5175      for ( var i = 0 ; i < maxIter ; i++ ) d.push( [] );
5176
5177      // Richardson extrapolation of Euler-Maclaurin trapezoids
5178      d[0][0] = result;
5179
5180      for ( var i = 1 ; i < maxIter ; i++ ) {
5181
5182        nextEulerIteration();
5183        d[0][i] = h * s;
5184
5185        for ( var j = 1 ; j <= i ; j++ ) {
5186
5187          d[j][i] = ( 4**j * d[j-1][i] - d[j-1][i-1] ) / ( 4**j - 1 );
5188
5189          var delta = Math.max( Math.abs( d[j][i] - d[j-1][i] ),
5190                                Math.abs( d[j][i] - d[j-1][i-1] ) );
5191          if ( delta <= error ) {
5192            error = delta;
5193            result = d[j][i];
5194          }
5195
5196        }
5197
5198        if ( Math.abs( d[i][i] - d[i-1][i-1] ) > error ) break;
5199
5200      }
5201
5202      return result;
5203
5204    case 'adaptive-simpson':
5205
5206      // algorithm by Charles Collins
5207
5208      var maxIter = 50;
5209
5210      function adaptiveSimpson( a, b, fa, fm, fb, s, tolerance, depth ) {
5211
5212        var h = b - a;
5213        var f1 = f( a + h/4 );
5214        var f2 = f( b - h/4 )
5215
5216        if ( isNaN(f1) || isNaN(f2) ) throw Error( 'NaN encountered in integration' );
5217
5218        var s1 = ( fa + 4*f1 + fm ) * h / 12;
5219        var s2 = ( fm + 4*f2 + fb ) * h / 12;
5220        var ss = s1 + s2;
5221        var error = ( ss - s ) / 15;
5222
5223        if ( Math.abs(error) < tolerance  || depth > maxIter ) return ss + error;
5224        else {
5225          var m = a + h/2;
5226          return adaptiveSimpson( a, m, fa, f1, fm, s1, tolerance/2, depth+1 )
5227                 + adaptiveSimpson( m, b, fm, f2, fb, s2, tolerance/2, depth+1 );
5228        }
5229
5230      }
5231
5232      var fa = f(a);
5233      var fm = f( (a+b)/2 );
5234      var fb = f(a);
5235      var s = ( fa + 4*fm + fb ) * (b-a) / 6;
5236      var depth = 0;
5237
5238      return adaptiveSimpson( a, b, fa, fm, fb, s, tolerance, depth );
5239
5240    case 'tanh-sinh':
5241
5242      // based on Borwein & Bailey, Experimentation in Mathematics
5243
5244      var m = 10;
5245      var h = 1 / 2**m;
5246      var x = [], w = [];
5247
5248      for ( var k = 0 ; k <= 20 * 2**m ; k++ ) {
5249        var t = k * h;
5250        x[k] = Math.tanh( Math.PI/2 * Math.sinh(t) );
5251        w[k] = Math.PI/2 * Math.cosh(t) / Math.cosh( Math.PI/2 * Math.sinh(t) )**2;
5252        if ( Math.abs(1-x[k]) < tolerance ) break;
5253      }
5254
5255      var nt = k;

5256      var sum = 0;
5257
5258      // rescale [a,b] to [-1,1]
5259      var len = ( b - a ) / 2;
5260      var mid = ( b + a ) / 2;
5261
5262      for ( var k = 1 ; k <= m ; k++ ) {
5263        for ( var i = 0 ; i < nt ; i += 2**(m-k) ) {
5264          if ( i % 2**(m-k+1) !== 0 || k === 1 ) {
5265            if ( i === 0 ) sum += w[0] * f( mid );
5266            else sum += w[i] * ( f( mid - len*x[i] ) + f( mid + len*x[i] ) );
5267          }
5268        }
5269      }
5270
5271      return len * h * sum;
5272
5273    case 'gaussian':
5274
5275      // based on Borwein & Bailey, Experimentation in Mathematics
5276
5277      var m = 10;
5278      var x = [], w = [];
5279
5280      var n = 3 * 2**m;
5281
5282      for ( var j = 1 ; j <= n/2 ; j++ ) {
5283
5284        var r = Math.cos( Math.PI * (j-.25) / (n+.5) );
5285
5286        while ( true ) {
5287
5288          var t1 = 1, t2 = 0;
5289
5290          for ( var j1 = 1 ; j1 <= n ; j1++ ) {
5291            t3 = t2;
5292            t2 = t1;
5293            t1 = ( (2*j1-1) * r * t2 - (j1-1) * t3 ) / j1;
5294          }
5295
5296          var t4 = n * ( r*t1 - t2 ) / ( r**2 - 1 );
5297          var delta = t1 / t4;
5298          r -= delta;
5299
5300          if ( Math.abs( delta ) < tolerance ) break;
5301
5302        }
5303
5304        x[j] = r;
5305        w[j] = 2 / ( 1 - r**2 ) / t4**2
5306
5307      }
5308
5309      // rescale [a,b] to [-1,1]
5310      var len = ( b - a ) / 2;
5311      var mid = ( b + a ) / 2;
5312
5313      var sum = 0;
5314
5315      for ( var j = 1 ; j <= n/2 ; j++ ) 
5316        sum += w[j] * ( f( mid - len*x[j] ) + f( mid + len*x[j] ) );
5317
5318      return len * sum;
5319
5320    default:
5321
5322      throw Error( 'Unsupported integration method' );
5323
5324  }
5325
5326}
5327
5328
5329function discreteIntegral( values, step ) {
5330
5331  // composite trapezoidal rule
5332
5333  var s = ( values[0] + values[ values.length - 1 ] ) / 2;
5334
5335  for ( var i = 1 ; i < values.length - 1 ; i++ ) s += values[i];
5336
5337  return s * step;
5338
5339}
5340
5341
5342function summation( f, [a,b,step=1] ) {
5343
5344  if ( isComplex( f(a) ) ) {
5345
5346    var s = complex(0);
5347
5348    for ( var i = a ; i <= b ; i += step ) s = add( s, f(i) );
5349
5350    return s;
5351
5352  } else {
5353
5354    var s = 0;
5355
5356    for ( var i = a ; i <= b ; i += step ) s += f(i);
5357
5358    return s;
5359
5360  }
5361
5362}
5363
5364
5365function polynomial( x, coefficients, options={} ) {
5366
5367  var derivative = 'derivative' in options ? options.derivative : false;
5368  var reverse = 'reverse' in options ? options.reverse : false;
5369
5370  if ( reverse ) coefficients = coefficients.slice().reverse();
5371
5372  // Horner's method with highest power coefficient first
5373
5374  var p = coefficients[0];
5375  var q = 0;
5376
5377  for ( var i = 1 ; i < coefficients.length ; i++ ) {
5378    if ( derivative ) q = add( p, mul( q, x ) );
5379    p = add( coefficients[i], mul( p, x ) );
5380  }
5381
5382  if ( derivative ) return { polynomial: p, derivative: q };
5383
5384  return p;
5385
5386}
5387
5388
5389function partialBell( n, k, argumentArray ) {
5390
5391  if ( n === 0 && k === 0 ) return 1;
5392
5393  if ( n === 0 || k === 0 ) return 0;
5394
5395  // evaluate recursively
5396  var s = 0;
5397  var p = 1;
5398
5399  for ( var i = 1 ; i <= n - k + 1 ; i++ ) {
5400
5401    s += p * argumentArray[i-1] * partialBell( n-i, k-1, argumentArray );
5402    p *= ( n - i ) / i;
5403
5404  }
5405
5406  return s;
5407
5408}
5409
5410
5411function findRoot( f, start, options={} ) {
5412
5413  var tolerance = 'tolerance' in options ? options.tolerance : 1e-10;
5414  var maxIter = 100;
5415
5416  if ( Array.isArray(f) ) {
5417
5418    if ( f.length !== start.length )
5419      throw Error( 'Mismatch between equations and starting point for root' );
5420
5421    var root = start.slice();
5422
5423    for ( var i = 0; i < maxIter ; i++ ) {
5424
5425      var J = [], F = [];
5426
5427      for ( var j = 0 ; j < root.length ; j++ ) {
5428        J.push( gradient( f[j], root ) );
5429        F.push( f[j].apply( null, root ) );
5430      }
5431
5432      var delta = luSolve( J, F );
5433
5434      for ( var j = 0 ; j < root.length ; j++ ) root[j] -= delta[j];
5435
5436      if ( delta.every( d => Math.abs(d) < tolerance ) ) return root;
5437
5438    }
5439
5440    throw Error( 'No root found for tolerance ' + tolerance );
5441
5442  }
5443
5444  if ( !Array.isArray(start) && !options.method ) options.method = 'newton';
5445
5446  var method = 'method' in options ? options.method : 'bisect';
5447
5448  switch( method ) {
5449
5450    case 'bisect':
5451
5452      var a = start[0];
5453      var b = start[1];
5454
5455      var fa = f(a);
5456      var fb = f(b);
5457
5458      if ( fa * f(b) >= 0 ) throw Error( 'Change of sign necessary for bisection' );
5459
5460      var root, h;
5461      if ( fa < 0 ) {
5462        root = a;
5463        h = b - a;
5464      } else {
5465        root = b;
5466        h = a - b;
5467      }
5468
5469      for ( var i = 0; i < maxIter ; i++ ) {
5470        h /= 2;
5471        var mid = root + h;
5472        fmid = f(mid);
5473        if ( fmid <= 0 ) root = mid;
5474        if ( fmid === 0 || Math.abs(h) < tolerance ) return root;
5475      }
5476
5477      throw Error( 'No root found for tolerance ' + tolerance );
5478
5479    case 'newton':
5480
5481      var root = start;
5482
5483      if ( isComplex( f(root) ) ) {
5484
5485        for ( var i = 0; i < maxIter ; i++ ) {
5486          var delta = div( f(root), diff( f, root ) );
5487          root = sub( root, delta );
5488          if ( abs(delta) < tolerance ) return root;
5489        }
5490
5491      } else {
5492
5493        for ( var i = 0; i < maxIter ; i++ ) {
5494          var delta = f(root) / diff( f, root );
5495          root -= delta;
5496          if ( Math.abs(delta) < tolerance ) return root;
5497        }
5498
5499      }
5500
5501      throw Error( 'No root found for tolerance ' + tolerance );
5502
5503    default:
5504
5505      throw Error( 'Unsupported root finding method' );
5506
5507  }
5508
5509}
5510
5511
5512function findRoots( f, point, tolerance=1e-10 ) {
5513
5514  console.log( 'Renamed to findRoot' );
5515
5516}
5517
5518
5519function spline( points, value='function', tolerance=1e-10 ) {
5520
5521  // adapted from gsl / cspline.c and reference therein
5522
5523  var a = [], b = [], c = [], d = [];
5524
5525  for ( var i = 0 ; i < points.length ; i++ ) a[i] = points[i][1];
5526
5527  c[0] = 0;
5528  c[ points.length - 1 ] = 0;
5529
5530  var A = matrix( points.length - 2 );
5531  var y = vector( points.length - 2 );
5532
5533  function h( i ) { return points[i+1][0] - points[i][0]; }
5534
5535  for ( var i = 0 ; i < A.length ; i++ ) {
5536    A[i][i] = 2 * ( h(i) + h(i+1) );
5537    y[i] = 3 * ( a[i+2] - a[i+1] ) / h(i+1) - 3 * ( a[i+1] - a[i] ) / h(i);

5538  }
5539  for ( var i = 1 ; i < A.length ; i++ ) {
5540    A[i][i-1] = h(i); 
5541    A[i-1][i] = h(i); 
5542  }
5543
5544  var x = luSolve( A, y );
5545
5546  for ( var i = 0 ; i < x.length ; i++ ) c[i+1] = x[i];
5547
5548  for ( var i = 0 ; i < c.length - 1 ; i++ ) {
5549    b[i] = ( a[i+1] - a[i] ) / h(i) - ( c[i+1] + 2*c[i] ) * h(i) / 3;
5550    d[i] = ( c[i+1] - c[i] ) / 3 / h(i);
5551  }
5552
5553  switch( value ) {
5554
5555    case 'function':
5556
5557      return function( x ) {
5558
5559        if ( x < points[0][0] || x > points[points.length-1][0] )
5560          throw Error( 'Argument outside spline input domain' );
5561
5562        for ( var i = 0 ; i < points.length ; i++ )
5563          if ( x === points[i][0] ) return a[i];
5564
5565        for ( var i = 0 ; i < points.length - 1 ; i++ )
5566          if ( x > points[i][0] && x < points[i+1][0] ) {
5567            var xi = points[i][0];
5568            return a[i] + b[i] * ( x - xi )
5569                   + c[i] * ( x - xi )**2 + d[i] * ( x - xi )**3;
5570          }
5571
5572      }
5573
5574    case 'derivative':
5575
5576      return function( x ) {
5577
5578        if ( x < points[0][0] || x > points[points.length-1][0] )
5579          throw Error( 'Argument outside spline input domain' );
5580
5581        // method does not define b[points.length-1] so fudge endpoint
5582        if ( x === points[points.length-1][0] ) x -= tolerance;
5583
5584        for ( var i = 0 ; i < points.length ; i++ )
5585          if ( x === points[i][0] ) return b[i];
5586
5587        for ( var i = 0 ; i < points.length - 1 ; i++ )
5588          if ( x > points[i][0] && x < points[i+1][0] ) {
5589            var xi = points[i][0];
5590            return b[i] + 2 * c[i] * ( x - xi ) + 3 * d[i] * ( x - xi )**2;
5591          }
5592
5593      }
5594
5595    case 'integral':
5596
5597      return function( x ) {
5598
5599        if ( x < points[0][0] || x > points[points.length-1][0] )
5600          throw Error( 'Argument outside spline input domain' );
5601
5602        var sum = 0;
5603
5604        function F( x, i ) {
5605          var xi = points[i][0];
5606          return a[i] * ( x - xi ) + b[i] * ( x - xi )**2 / 2
5607                   + c[i] * ( x - xi )**3 / 3 + d[i] * ( x - xi )**4 / 4;
5608        }
5609
5610        for ( var i = 0 ; i < points.length - 1 ; i++ )
5611          if ( x < points[i+1][0] ) {
5612            sum += F( x, i ) - F( points[i][0], i );
5613            break;
5614          } else sum += F( points[i+1][0], i ) - F( points[i][0], i );
5615
5616        return sum;
5617
5618      }
5619
5620    default:
5621
5622      throw Error( 'Unsupported spline value' );
5623
5624  }
5625
5626}
5627
5628function padeApproximant( f, n, d, center=0 ) {
5629
5630  if ( !isPositiveInteger(n) || !isPositiveInteger(d) )
5631    throw Error( 'Pade indices must be positive integers' );
5632
5633  var c = [];
5634
5635  if ( typeof f === 'function' ) {
5636    c.push( f(center) );
5637    for ( var i = 1 ; i <= n+d ; i++ ) c.push( diff( f, center, i ) / factorial(i) );
5638  } else {
5639    if ( f.length < n+d+1 ) throw Error( 'Need n+d+1 Pade coefficients' );
5640    c = f;
5641  }
5642
5643  var M = matrix( d ), v = vector( d );
5644
5645  for ( var i = 0 ; i < d ; i++ ) {
5646    v[i] = -c[ n + 1 + i ];
5647    // avoid negative index leading to undefined entry
5648    for ( var j = 0 ; j < d && j <= n+i ; j++ ) M[i][j] = c[ n + i - j ];
5649  }
5650
5651  // need to understand why Mathematica returns results for singular matrices
5652  try { var b = luSolve( M, v ); }
5653  catch { throw Error( 'Singular Pade matrix encountered' ); };
5654  b.unshift( 1 );
5655
5656  var a = [ c[0] ];
5657
5658  for ( var i = 1 ; i <= n ; i++ ) {
5659    var s = c[i];
5660    // cannot exceed i when d > n
5661    for( var j = 1 ; j <= d && j <= i ; j++ ) s += b[j] * c[i-j];
5662    a.push( s );
5663  }
5664
5665  return { N: chop(a), D: chop(b) };
5666
5667}
5668
5669
5670function fourierSinCoefficient( f, n, period ) {
5671
5672  if ( !Number.isInteger(n) ) throw Error( 'Nonintegral Fourier index' );
5673
5674  if ( n === 0 ) return 0;
5675
5676  if ( typeof f === 'function' ) {
5677
5678    var T = period || 2*pi;
5679
5680    return 2/T * integrate( t => f(t) * sin( 2*n*pi/T * t ), [0,T], { method: 'tanh-sinh' } );
5681
5682  }
5683
5684  if ( Array.isArray(f) ) {
5685
5686    var s = 0, N = f.length;
5687
5688    for ( var i = 0 ; i < N ; i++ ) s += f[i][1] * sin( 2*n*pi*i/N );
5689
5690    return 2 * s / N;
5691
5692  }
5693
5694  throw Error( 'Unsupported Fourier input' );
5695
5696}
5697
5698function fourierCosCoefficient( f, n, period ) {
5699
5700  if ( !Number.isInteger(n) ) throw Error( 'Nonintegral Fourier index' );
5701
5702  if ( typeof f === 'function' ) {
5703
5704    var T = period || 2*pi;
5705
5706    if ( n === 0 ) return 1/T * integrate( t => f(t), [0,T], { method: 'tanh-sinh' } );
5707
5708    return 2/T * integrate( t => f(t) * cos( 2*n*pi/T * t ), [0,T], { method: 'tanh-sinh' } );
5709
5710  }
5711
5712  if ( Array.isArray(f) ) {
5713
5714    var s = 0, N = f.length;
5715
5716    if ( n === 0 ) {
5717
5718      for ( var i = 0 ; i < N ; i++ ) s += f[i][1];
5719
5720      return s / N;
5721
5722    }
5723
5724    for ( var i = 0 ; i < N ; i++ ) s += f[i][1] * cos( 2*n*pi*i/N );
5725
5726    return 2 * s / N;
5727
5728  }
5729
5730  throw Error( 'Unsupported Fourier input' );
5731
5732}
5733
5734
5735function eigensystem( A, symmetric=true ) {
5736
5737  if ( symmetric ) return tridiagonalQL( tridiagonalForm(A) );
5738  else throw Error( 'Unsupported eigensystem' );
5739
5740}
5741
5742// sourced from http://math.nist.gov/javanumerics/jama/
5743// no need to reinvent this wheel...
5744
5745function tridiagonalForm( A ) {
5746
5747  var n = A.length;
5748  var V = [];
5749  for ( var i = 0 ; i < n ; i++ ) V[i] = A[i].slice(); // deeper copy
5750
5751  var d = vector( n );
5752  var e = vector( n );
5753
5754  for ( var j = 0 ; j < n ; j++ ) d[j] = V[n-1][j];
5755
5756  // Householder reduction to tridiagonal form
5757   
5758  for ( var i = n - 1 ; i > 0 ; i-- ) {
5759   
5760    // scale to avoid under/overflow
5761   
5762    var scale = 0;
5763    var h = 0;
5764    for ( var k = 0 ; k < i ; k++ ) scale += Math.abs(d[k]);
5765
5766    if ( scale === 0 ) {
5767      e[i] = d[i-1];
5768      for ( var j = 0 ; j < i ; j++ ) {
5769        d[j] = V[i-1][j];
5770        V[i][j] = 0;
5771        V[j][i] = 0;
5772      }
5773    } else {
5774
5775      // generate Householder vector
5776
5777      for ( var k = 0 ; k < i ; k++ ) {
5778        d[k] /= scale;
5779        h += d[k] * d[k];
5780      }
5781      var f = d[i-1];
5782      var g = Math.sqrt(h);
5783      if ( f > 0 ) g = -g;
5784      e[i] = scale * g;
5785      h = h - f * g;
5786      d[i-1] = f - g;
5787      for ( var j = 0; j < i; j++ ) e[j] = 0;
5788
5789      // apply similarity transformation to remaining columns
5790
5791      for ( var j = 0 ; j < i ; j++ ) {
5792        f = d[j];
5793        V[j][i] = f;
5794        g = e[j] + V[j][j] * f;
5795        for ( var k = j + 1 ; k <= i - 1 ; k++ ) {
5796          g += V[k][j] * d[k];
5797          e[k] += V[k][j] * f;
5798        }
5799        e[j] = g;
5800      }
5801      f = 0;
5802      for ( var j = 0 ; j < i ; j++ ) {
5803        e[j] /= h;
5804        f += e[j] * d[j];
5805      }
5806      var hh = f / ( h + h );
5807      for ( var j = 0 ; j < i ; j++ ) {
5808        e[j] -= hh * d[j];
5809      }
5810      for ( var j = 0 ; j < i ; j++ ) {
5811        f = d[j];
5812        g = e[j];
5813        for ( var k = j ; k <= i - 1 ; k++ ) {
5814          V[k][j] -= f * e[k] + g * d[k];
5815        }
5816        d[j] = V[i-1][j];
5817        V[i][j] = 0;
5818      }
5819
5820    }
5821
5822    d[i] = h;
5823
5824  }
5825   
5826  // accumulate transformations
5827   
5828  for ( var i = 0 ; i < n - 1 ; i++ ) {
5829    V[n-1][i] = V[i][i];
5830    V[i][i] = 1;
5831    var h = d[i+1];
5832    if ( h !== 0 ) {
5833      for ( var k = 0 ; k <= i ; k++ ) d[k] = V[k][i+1] / h;
5834      for ( var j = 0 ; j <= i ; j++ ) {
5835        var g = 0;
5836        for ( var k = 0 ; k <= i ; k++ ) g += V[k][i+1] * V[k][j];
5837        for ( var k = 0 ; k <= i ; k++ ) V[k][j] -= g * d[k];
5838      }
5839    }
5840    for ( var k = 0; k <= i; k++) V[k][i+1] = 0;
5841  }
5842  for ( var j = 0 ; j < n ; j++ ) {
5843    d[j] = V[n-1][j];
5844    V[n-1][j] = 0;
5845  }
5846  V[n-1][n-1] = 1;
5847  e[0] = 0;
5848
5849  return { diagonal: d, offDiagonal: e, eigenvectors: V };
5850
5851}
5852
5853
5854function tridiagonalQL( tridiagonalForm ) {
5855
5856  var d = tridiagonalForm.diagonal;
5857  var n = d.length;
5858  var e = tridiagonalForm.offDiagonal;
5859  var V = tridiagonalForm.eigenvectors;
5860
5861  function hypot( a, b) {
5862    var r;
5863    if ( Math.abs(a) > Math.abs(b) ) {
5864      r = b/a;
5865      r = Math.abs(a) * Math.sqrt( 1 + r*r );
5866    } else if (b != 0) {
5867      r = a/b;
5868      r = Math.abs(b) * Math.sqrt( 1 + r*r );
5869    } else r = 0;
5870
5871    return r;
5872
5873  }
5874
5875  for ( var i = 1 ; i < n ; i++ ) e[i-1] = e[i];
5876  e[n-1] = 0;
5877
5878  var f = 0;
5879  var tst1 = 0;
5880  var eps = Math.pow( 2, -52 );
5881
5882  for ( var l = 0 ; l < n ; l++ ) {
5883
5884    // find small subdiagonal element
5885
5886    tst1 = Math.max( tst1, Math.abs(d[l]) + Math.abs(e[l]) );
5887    var m = l;
5888    while ( m < n ) {
5889      if ( Math.abs(e[m]) <= eps*tst1 ) break;
5890      m++;
5891    }
5892
5893    // if m === l, d[l] is an eigenvalue, otherwise iterate
5894
5895    if ( m > l ) {
5896
5897      var iter = 0;
5898      do {
5899
5900        iter = iter + 1;
5901        if ( iter > 1000 ) throw Error( 'Eigenvalues not converging...' );
5902
5903        // compute implicit shift
5904
5905        var g = d[l];
5906        var p = ( d[l+1] - g ) / ( 2 * e[l] );

5907        var r = hypot(p,1);
5908        if ( p < 0 ) r = -r;
5909        d[l] = e[l] / ( p + r );
5910        d[l+1] = e[l] * ( p + r );
5911        var dl1 = d[l+1];
5912        var h = g - d[l];
5913        for ( var i = l + 2 ; i < n ; i++ ) d[i] -= h;
5914        f = f + h;
5915
5916        // implicit QL transformation
5917
5918        p = d[m];
5919        var c = 1;
5920        var c2 = c;
5921        var c3 = c;
5922        var el1 = e[l+1];
5923        var s = 0;
5924        var s2 = 0;
5925        for ( var i = m - 1 ; i >= l ; i-- ) {
5926          c3 = c2;
5927          c2 = c;
5928          s2 = s;
5929          g = c * e[i];
5930          h = c * p;
5931          r = hypot(p,e[i]);
5932          e[i+1] = s * r;
5933          s = e[i] / r;
5934          c = p / r;
5935          p = c * d[i] - s * g;
5936          d[i+1] = h + s * ( c * g + s * d[i] );
5937
5938          // accumulate transformation
5939
5940          for ( var k = 0 ; k < n ; k++ ) {
5941            h = V[k][i+1];
5942            V[k][i+1] = s * V[k][i] + c * h;
5943            V[k][i] = c * V[k][i] - s * h;
5944          }
5945        }
5946        p = -s * s2 * c3 * el1 * e[l] / dl1;
5947        e[l] = s * p;
5948        d[l] = c * p;
5949

5950        // check for convergence
5951   
5952      } while ( Math.abs(e[l]) > eps*tst1 );
5953    }
5954    d[l] = d[l] + f;
5955    e[l] = 0;
5956  }
5957   
5958  // sort eigenvalues and corresponding vectors
5959
5960  for ( var i = 0 ; i < n - 1 ; i++ ) {
5961    var k = i;
5962    var p = d[i];
5963    for ( var j = i + 1 ; j < n ; j++ ) {
5964      if ( d[j] < p ) {
5965        k = j;
5966        p = d[j];
5967      }
5968    }
5969    if ( k != i ) {
5970      d[k] = d[i];
5971      d[i] = p;
5972      for ( var j = 0 ; j < n ; j++ ) {
5973        p = V[j][i];
5974        V[j][i] = V[j][k];
5975        V[j][k] = p;
5976      }
5977    }
5978  }
5979
5980  return { eigenvalues: d, eigenvectors: V };
5981
5982}
5983
5984
5985function hessenbergForm( A ) {
5986
5987}
5988
5989
5990function luDecomposition( A, tolerance=1e-10 ) {
5991
5992  var size = A.length;
5993  var LU = [];
5994  for ( var i = 0 ; i < size ; i++ ) LU[i] = A[i].slice(); // deeper copy
5995
5996  var P = identity( size );
5997  pivots = 0;
5998
5999  for ( var i = 0 ; i < size ; i++ ) {
6000
6001    var maxValue = 0;
6002    var maxIndex = i;
6003
6004    for ( var j = i ; j < size ; j++ ) {
6005      var element = Math.abs( LU[j][i] );
6006      if ( element > maxValue ) {
6007        maxValue = element;
6008        maxIndex = j;
6009      }
6010    }
6011
6012    if ( maxValue < tolerance ) throw Error( 'Matrix is degenerate' );
6013
6014    if ( maxIndex !== i ) {
6015
6016      // pivot matrix rows
6017      var t = LU[i];
6018      LU[i] = LU[maxIndex];
6019      LU[maxIndex] = t;
6020
6021      // pivot permutation rows
6022      var t = P[i];
6023      P[i] = P[maxIndex];
6024      P[maxIndex] = t;
6025
6026      pivots++;
6027
6028    }
6029
6030    for ( var j = i + 1 ; j < size ; j++ ) {
6031      LU[j][i] /= LU[i][i];
6032      for ( var k = i + 1; k < size ; k++ )
6033        LU[j][k] -= LU[j][i] * LU[i][k];
6034    }
6035  }
6036
6037  var L = identity( size );
6038  for ( var i = 1 ; i < size ; i++ )
6039    for ( var j = 0 ; j < i ; j++ ) L[i][j] = LU[i][j];
6040
6041  var U = matrix( size );
6042  for ( var i = 0 ; i < size ; i++ )
6043    for ( var j = i ; j < size ; j++ ) U[i][j] = LU[i][j];
6044
6045  return { L: L, U: U, P: P, pivots: pivots };
6046
6047}
6048
6049function luSolve( A, b ) {
6050
6051  var size = A.length;
6052  var lu = luDecomposition(A);
6053
6054  var x = vector( size );
6055  var y = vector( size );
6056  var pb = vector( size );
6057
6058  for ( var i = 0 ; i < size ; i++ )
6059    for ( var j = 0 ; j < size ; j++ )
6060      pb[i] += lu.P[i][j] * b[j];
6061
6062  // forward solve
6063  for ( var i = 0 ; i < size ; i++ ) {
6064    y[i] = pb[i];
6065    for ( var j = 0 ; j < i ; j++ ) y[i] -= lu.L[i][j] * y[j];
6066    y[i] /= lu.L[i][i];
6067  }
6068
6069  // backward solve
6070  for ( var i = size - 1 ; i >= 0 ; i-- ) {
6071    x[i] = y[i];
6072    for ( var j = i + 1 ; j < size ; j++ ) x[i] -= lu.U[i][j] * x[j];
6073    x[i] /= lu.U[i][i];
6074  }
6075
6076  return x;
6077
6078}
6079
6080function determinant( A ) {
6081
6082  var lu = luDecomposition(A);
6083
6084  var product = 1;
6085  for ( var i = 0 ; i < A.length; i++ ) product *= lu.U[i][i];
6086
6087  return (-1)**lu.pivots * product;
6088
6089}
6090
6091function inverse( A ) {
6092
6093  // calling luSolve for each column is not efficient
6094  //   but avoids code duplication
6095
6096  var I = matrix( A.length );
6097
6098  for ( var i = 0 ; i < A.length ; i++ ) {
6099
6100    var b = vector( A.length );
6101    b[i] = 1;
6102
6103    var x = luSolve( A, b );
6104    for ( var j = 0 ; j < A.length ; j++ ) I[j][i] = x[j];
6105
6106  }
6107
6108  return I;
6109
6110}
6111
6112
6113
6114function vector( size, value=0 ) {
6115
6116  var v = [];
6117  for ( var i = 0 ; i < size ; i++ )
6118    if ( Array.isArray(value) ) v.push( random.apply( null, value ) );
6119    else v.push( value );
6120
6121  return v;
6122
6123}
6124
6125function matrix( rows, columns, value=0 ) {
6126
6127  var columns = columns || rows;
6128
6129  var m = [];
6130  for ( var i = 0 ; i < rows ; i++ ) {
6131    m.push( [] );
6132    for ( var j = 0 ; j < columns ; j++ )
6133      if ( Array.isArray(value) ) m[i].push( random.apply( null, value ) );
6134      else m[i].push( value );
6135  }
6136
6137  return m;
6138
6139}
6140
6141function identity( rows, value=1 ) {
6142
6143  var m = matrix( rows );
6144  for ( var i = 0 ; i < rows ; i++ ) m[i][i] = value;
6145
6146  return m;
6147
6148}
6149
6150function transpose( A ) {
6151
6152  var T = matrix( A[0].length, A.length );
6153
6154  for ( var i = 0 ; i < A.length ; i++ )
6155    for ( var j = 0 ; j < A[0].length ; j++ )
6156      T[j][i] = A[i][j];
6157
6158  return T;
6159
6160}
6161
6162function matrixAdd( A, B ) {
6163
6164  if ( !Array.isArray(A) && !Array.isArray(B) ) throw Error( 'No matrices to add' );
6165  if ( !Array.isArray(A) ) A = matrix( B.length, B[0].length, A );
6166  if ( !Array.isArray(B) ) B = matrix( A.length, A[0].length, B );
6167
6168  var C = matrix( A.length, A[0].length, 0 );
6169
6170  for ( var i = 0 ; i < A.length ; i++ )
6171    for ( var j = 0 ; j < A[0].length ; j++ )
6172      C[i][j] = add( A[i][j], B[i][j] );
6173
6174  return C;

6175
6176}
6177
6178function matrixSub( A, B ) {
6179
6180  if ( !Array.isArray(A) && !Array.isArray(B) ) throw Error( 'No matrices to subtract' );
6181  if ( !Array.isArray(A) ) A = matrix( B.length, B[0].length, A );
6182  if ( !Array.isArray(B) ) B = matrix( A.length, A[0].length, B );
6183
6184  var C = matrix( A.length, A[0].length, 0 );
6185
6186  for ( var i = 0 ; i < A.length ; i++ )
6187    for ( var j = 0 ; j < A[0].length ; j++ )
6188      C[i][j] = sub( A[i][j], B[i][j] );
6189
6190  return C;
6191
6192}
6193
6194function matrixMul( A, B ) {
6195
6196  if ( !Array.isArray(A) && !Array.isArray(B) ) throw Error( 'No matrices to multiply' );
6197  if ( !Array.isArray(A) ) A = identity( B.length, A );
6198  if ( !Array.isArray(B) ) B = identity( A[0].length, B );
6199  if ( A[0].length !== B.length ) throw Error( 'Incompatible matrices to multiply' );
6200
6201  var C = matrix( A.length, B[0].length, 0 );
6202
6203  for ( var i = 0 ; i < A.length ; i++ )
6204    for ( var j = 0 ; j < B[0].length ; j++ )
6205      for ( var k = 0 ; k < A[0].length ; k++ )
6206        C[i][j] = add( C[i][j], mul( A[i][k], B[k][j] ) );
6207
6208  return C;
6209
6210}
6211

Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.