1 // -*- C++ -*-
2 //==============================================================================================
3 //
4 // This file is part of LiDIA --- a library for computational number theory
5 //
6 // Copyright (c) 1994--2002 the LiDIA Group. All rights reserved.
7 //
8 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
9 //
10 //----------------------------------------------------------------------------------------------
11 //
12 //
13 // File : modular_functions.cc
14 // Author : Harald Baier (HB)
15 // Changes : See CVS log
16 //
17 //==============================================================================================
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22
23 # include "LiDIA/modular_functions.h"
24
25 #ifdef VERBOSE
26 # include <cassert>
27 #endif
28
29 #ifdef LIDIA_NAMESPACE
30 namespace LiDIA {
31 #endif
32
33
34 /********************************************************************
35 *
36 * Dedekind's eta-function computed via the Euler sum.
37 * See PhD for details.
38 * The current implementation corresponds to the algorithm
39 * computeEtaViaEulerSum.
40 *
41 *******************************************************************/
42 //
43 // dedekind_eta bases on the algorithm of PhD
44 //
45 bigcomplex
dedekind_eta(const bigcomplex & tau,bool is_fundamental)46 dedekind_eta( const bigcomplex & tau, bool is_fundamental )
47 {
48 // complex variables:
49 // * sum: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/2}+q^{(3n^2+n)/2})
50 // * q: stores e^(2\pi i\tau)
51 // * q_to_n: stores current power of q
52 // * q_to_diff: stores q^diff, where diff stores the difference
53 // between two non-trivial terms
54 bigcomplex sum( 1 ), local_tau( tau ), q, q_to_n, q_to_diff;
55 bigcomplex factor( 1 ), result, tau_inverse, square_root;
56
57 // if is_fundamental = false translate tau to fundamental area
58 // and compute appropriate factor
59 if( is_fundamental == false )
60 {
61 //
62 // ensure that |Re( tau )| <= 1/2 and |tau| >= 1
63 // otherwise shift tau to this stripe
64 //
65 bigint real_part_as_bigint;
66
67 while( true )
68 {
69 // assign to real_part_as_bigint the nearest integer
70 // of Re(local_tau)
71 local_tau.real().bigintify( real_part_as_bigint );
72
73 // if it is not zero, shift it in stripe |Re( tau )| <= 1/2
74 // eta(tau+1) = zeta_24 * eta(tau)
75 if( ! real_part_as_bigint.is_zero() )
76 {
77 subtract( local_tau, local_tau,
78 bigcomplex( bigfloat( real_part_as_bigint ) ) );
79 multiply( factor, factor,
80 exp( 2 * bigcomplex( 0, Pi() ) *
81 bigfloat( real_part_as_bigint, 24 ) ) );
82 }
83
84 // if |local_tau| < 1 invert it to put it out of unit circle
85 // eta(-1/tau) = \sqrt(-i*tau) * eta(tau)
86 if( abs( local_tau ) <= 0.9999 )
87 {
88 invert( tau_inverse, local_tau );
89 multiply( square_root, bigcomplex( 0, 1 ), tau_inverse );
90 sqrt( square_root, square_root );
91
92 if( square_root.real() < 0 )
93 square_root.negate();
94
95 multiply( factor, factor, square_root );
96
97 local_tau.assign( tau_inverse );
98 local_tau.negate();
99 }
100 // otherwise we are done
101 else
102 break;
103 }
104
105 }
106
107 // factor stores the prefactor of Euler sum
108 multiply( factor, factor,
109 exp( 2 * bigcomplex( 0, Pi() ) * local_tau *
110 bigfloat( 1, 24 ) ) );
111
112 exp( q, 2 * bigcomplex( 0, Pi() ) * local_tau );
113
114 // to store the powers of q in use:
115 // * N_minus = ( 3 * n_minus^2 - n_minus ) / 2
116 // * N_plus = ( 3 * n_plus^2 - n_plus ) / 2
117 // * previous_power stores the previous used power of q
118 bigint N_minus( 5 ), N_plus( 7 );
119 int n_minus, n_plus;
120 n_minus = n_plus = 2;
121 bigint previous_power( 2 ), diff;
122
123 // assign the correct values: sum = 1 - q - q^2
124 square( q_to_n, q );
125 subtract( sum, sum, q );
126 subtract( sum, sum, q_to_n );
127
128 while( true )
129 {
130 if( N_minus < N_plus )
131 {
132
133 subtract( diff, N_minus, previous_power );
134 power( q_to_diff, q, diff );
135 multiply( q_to_n, q_to_n, q_to_diff );
136
137 if( n_minus & 1 )
138 subtract( sum, sum, q_to_n );
139 else
140 add( sum, sum, q_to_n );
141
142 if( q_to_n.is_approx_zero() )
143 {
144 multiply( result, factor, sum );
145
146 return( result );
147 }
148
149 previous_power.assign( N_minus );
150 N_minus += 3 * n_minus + 1;
151 n_minus++;
152 }
153
154 else
155 {
156 subtract( diff, N_plus, previous_power );
157 power( q_to_diff, q, diff );
158 multiply( q_to_n, q_to_n, q_to_diff );
159
160 if( n_plus & 1 )
161 subtract( sum, sum, q_to_n );
162 else
163 add( sum, sum, q_to_n );
164
165 if( q_to_n.is_approx_zero() )
166 {
167 multiply( result, factor, sum );
168
169 return( result );
170 }
171
172 previous_power.assign( N_plus );
173 N_plus += 3 * n_plus + 2;
174 n_plus++;
175 }
176 }
177 }
178
179
180 /********************************************************************
181 *
182 * Weber functions computed via the eta-function.
183 * See PhD for details.
184 *
185 *******************************************************************/
186 //
187 // weber_f bases on the algorithm described in PhD
188 //
189 bigcomplex
weber_f(const bigcomplex & tau,bool is_fundamental)190 weber_f( const bigcomplex & tau, bool is_fundamental )
191 {
192 // if tau is not in the fundamental domain, we make use of the
193 // relation weber_f = eta( (tau+1)/2 ) / ( zeta_48 * eta( tau ) )
194 if( is_fundamental == false )
195 {
196 // declare the necessary variables
197 bigcomplex current_tau, numerator, denominator, zeta_48, result;
198
199 // compute zeta_48 = e^{2\pi i / 48}
200 exp( zeta_48, bigcomplex( 0, Pi() ) * bigfloat( 1, 24 ) );
201
202 // compute argument of numerator, numerator, and denominator
203 current_tau = ( tau + 1 ) / 2;
204 numerator = dedekind_eta( current_tau, false );
205 denominator = dedekind_eta( tau, false );
206
207 // compute the result and return it
208 divide( result, numerator, zeta_48 );
209 divide( result, result, denominator );
210
211 return( result );
212 }
213
214 // if is_fundamental is true, use the algorithm of PhD
215
216 // complex variables for computation of eta((tau+1)/2):
217 // * eta_1_2: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/4}+q^{(3n^2+n)/4})
218 // hence the current approximation of \eta(\tau/2)/q^{1/48}
219 // * q_1_2: to store q^{1/2} = e^(\pi i\tau)
220 // * q_1_2_to_n: to store q^{(3n^2+n)/4}
221 // * q_to_diff: stores q^diff, where diff stores the difference
222 // between two non-trivial terms
223 bigcomplex eta_1_2( 1 ), q_1_2, q_1_2_to_n, q_to_diff, result;
224 exp( q_1_2, bigcomplex( 0, Pi() ) * tau );
225
226 // complex variables for computation of eta(tau):
227 // * eta: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/2}+q^{(3n^2+n)/2}),
228 // hence the current approximation of \eta(\tau)/q^{1/24}
229 // * q_to_n: stores current power of q
230 bigcomplex eta( 1 ), q_to_n;
231
232 // to store the powers of q_1_2 in use:
233 // * N_minus = ( 3 * n_minus^2 - n_minus ) / 2
234 // * N_plus = ( 3 * n_plus^2 - n_plus ) / 2
235 // * previous_power stores the previous used power of q
236 bigint N_minus( 5 ), N_plus( 7 );
237 int n_minus, n_plus;
238 n_minus = n_plus = 2;
239 bigint previous_power( 2 ), diff;
240
241 // compute the prefactor q^{-1/48} of the result
242 bigcomplex prefactor;
243 exp( prefactor, bigcomplex( 0, Pi() ) * tau * bigfloat( -1, 24 ) );
244
245 // assign the correct values: eta_1_2 = 1 + q_1_2 - q_1_2^2
246 square( q_1_2_to_n, q_1_2 );
247 add( eta_1_2, eta_1_2, q_1_2 );
248 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
249
250 // assign the correct values: eta = 1 - q - q^2
251 subtract( eta, eta, q_1_2_to_n );
252 subtract( eta, eta, square( q_1_2_to_n ) );
253
254 // eta_false is true iff q_to_n is not zero within precision
255 bool eta_false = true;
256
257 // only evaluate this loop iff q_to_n is not zero within precision
258 while( eta_false )
259 {
260 if( N_minus < N_plus )
261 {
262 // compute the difference of the exponents
263 subtract( diff, N_minus, previous_power );
264 power( q_to_diff, q_1_2, diff );
265 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
266
267 square( q_to_n, q_1_2_to_n );
268
269 // take care of the sign and update eta_1_2 and eta
270 if( n_minus & 1 )
271 {
272 if( N_minus.bit( 0 ) )
273 add( eta_1_2, eta_1_2, q_1_2_to_n );
274 else
275 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
276
277 subtract( eta, eta, q_to_n );
278 }
279 else
280 {
281 if( N_minus.bit( 0 ) )
282 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
283 else
284 add( eta_1_2, eta_1_2, q_1_2_to_n );
285
286 add( eta, eta, q_to_n );
287 }
288
289 // check if q_to_n is zero within precision
290 if( q_to_n.is_approx_zero() )
291 eta_false = false;
292
293 // adapt the exponent N_minus and the counter n_minus
294 previous_power.assign( N_minus );
295 N_minus += 3 * n_minus + 1;
296 n_minus++;
297 }
298 else
299 {
300 // compute the difference of the exponents
301 subtract( diff, N_plus, previous_power );
302 power( q_to_diff, q_1_2, diff );
303 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
304
305 square( q_to_n, q_1_2_to_n );
306
307 // take care of the sign and update eta_1_2 and eta
308 if( n_plus & 1 )
309 {
310 if( N_plus.bit( 0 ) )
311 add( eta_1_2, eta_1_2, q_1_2_to_n );
312 else
313 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
314
315 subtract( eta, eta, q_to_n );
316 }
317 else
318 {
319 if( N_plus.bit( 0 ) )
320 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
321 else
322 add( eta_1_2, eta_1_2, q_1_2_to_n );
323
324 add( eta, eta, q_to_n );
325 }
326
327 // check if q_to_n is zero within precision
328 if( q_to_n.is_approx_zero() )
329 eta_false = false;
330
331 // adapt the exponent N_plus and the counter n_plus
332 previous_power.assign( N_plus );
333 N_plus += 3 * n_plus + 2;
334 n_plus++;
335 }
336 }
337
338 // evaluate this loop until q_to_n is zero within precision
339 while( true )
340 {
341 if( N_minus < N_plus )
342 {
343 // compute the difference of the exponents
344 subtract( diff, N_minus, previous_power );
345 power( q_to_diff, q_1_2, diff );
346 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
347
348 square( q_to_n, q_1_2_to_n );
349
350 // take care of the sign and update eta_1_2
351 if( n_minus & 1 )
352 {
353 if( N_minus.bit( 0 ) )
354 add( eta_1_2, eta_1_2, q_1_2_to_n );
355 else
356 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
357 }
358 else
359 {
360 if( N_minus.bit( 0 ) )
361 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
362 else
363 add( eta_1_2, eta_1_2, q_1_2_to_n );
364 }
365
366 // if q_1_2_to_n is zero within precision return the result:
367 // q^{-1/48} * eta_1_2 / eta = prefactor * eta_1_2 / eta
368 if( q_1_2_to_n.is_approx_zero() )
369 {
370 multiply( result, prefactor, eta_1_2 );
371 divide( result, result, eta );
372
373 return( result );
374 }
375
376 previous_power.assign( N_minus );
377 N_minus += 3 * n_minus + 1;
378 n_minus++;
379 }
380 else
381 {
382 // compute the difference of the exponents
383 subtract( diff, N_plus, previous_power );
384 power( q_to_diff, q_1_2, diff );
385 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
386
387 square( q_to_n, q_1_2_to_n );
388
389 // take care of the sign and update eta_1_2
390 if( n_plus & 1 )
391 {
392 if( N_plus.bit( 0 ) )
393 add( eta_1_2, eta_1_2, q_1_2_to_n );
394 else
395 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
396 }
397 else
398 {
399 if( N_plus.bit( 0 ) )
400 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
401 else
402 add( eta_1_2, eta_1_2, q_1_2_to_n );
403 }
404
405 // if q_1_2_to_n is zero within precision return the result:
406 // q^{-1/48} * eta_1_2 / eta = prefactor * eta_1_2 / eta
407 if( q_1_2_to_n.is_approx_zero() )
408 {
409 multiply( result, prefactor, eta_1_2 );
410 divide( result, result, eta );
411
412 return( result );
413 }
414
415 previous_power.assign( N_plus );
416 N_plus += 3 * n_plus + 2;
417 n_plus++;
418 }
419 }
420 }
421
422 //
423 // weber_f1 bases on the algorithm in PhD
424 //
425 bigcomplex
weber_f1(const bigcomplex & tau,bool is_fundamental)426 weber_f1( const bigcomplex & tau, bool is_fundamental )
427 {
428 // if tau is not in the fundamental domain, we make use of the
429 // relation weber_f1 = eta( tau/2 ) / eta( tau )
430 if( is_fundamental == false )
431 {
432 // declare the necessary variables
433 bigcomplex current_tau, numerator, denominator, result;
434
435 // compute argument of numerator, numerator, and denominator
436 divide( current_tau, tau, 2 );
437 numerator = dedekind_eta( current_tau, false );
438 denominator = dedekind_eta( tau, false );
439
440 // compute the result and return it
441 divide( result, numerator, denominator );
442
443 return( result );
444 }
445
446 // if is_fundamental is true, use the algorithm of PhD
447
448 // complex variables for computation of eta(tau/2):
449 // * eta_1_2: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/4}+q^{(3n^2+n)/4})
450 // hence the current approximation of \eta(\tau/2)/q^{1/48}
451 // * q_1_2: to store q^{1/2} = e^(\pi i\tau)
452 // * q_1_2_to_n: to store q^{(3n^2+n)/4}
453 // * q_to_diff: stores q^diff, where diff stores the difference
454 // between two non-trivial terms
455 bigcomplex eta_1_2( 1 ), q_1_2, q_1_2_to_n, q_to_diff, result;
456 exp( q_1_2, bigcomplex( 0, Pi() ) * tau );
457
458 // complex variables for computation of eta(tau):
459 // * eta: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/2}+q^{(3n^2+n)/2}),
460 // hence the current approximation of \eta(\tau)/q^{1/24}
461 // * q_to_n: stores current power of q
462 bigcomplex eta( 1 ), q_to_n;
463
464 // to store the powers of q_1_2 in use:
465 // * N_minus = ( 3 * n_minus^2 - n_minus ) / 2
466 // * N_plus = ( 3 * n_plus^2 - n_plus ) / 2
467 // * previous_power stores the previous used power of q
468 bigint N_minus( 5 ), N_plus( 7 );
469 int n_minus, n_plus;
470 n_minus = n_plus = 2;
471 bigint previous_power( 2 ), diff;
472
473 // compute the prefactor q^{-1/48} of the result
474 bigcomplex prefactor;
475 exp( prefactor, bigcomplex( 0, Pi() ) * tau * bigfloat( -1, 24 ) );
476
477 // assign the correct values: eta_1_2 = 1 - q_1_2 - q_1_2^2
478 square( q_1_2_to_n, q_1_2 );
479 subtract( eta_1_2, eta_1_2, q_1_2 );
480 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
481
482 // assign the correct values: eta = 1 - q - q^2
483 subtract( eta, eta, q_1_2_to_n );
484 subtract( eta, eta, square( q_1_2_to_n ) );
485
486 // eta_false is true iff q_to_n is not zero within precision
487 bool eta_false = true;
488
489 // only evaluate this loop iff q_to_n is not zero within precision
490 while( eta_false )
491 {
492 if( N_minus < N_plus )
493 {
494 // compute the difference of the exponents
495 subtract( diff, N_minus, previous_power );
496 power( q_to_diff, q_1_2, diff );
497 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
498
499 square( q_to_n, q_1_2_to_n );
500
501 // take care of the sign and update eta_1_2 and eta
502 if( n_minus & 1 )
503 {
504 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
505 subtract( eta, eta, q_to_n );
506 }
507 else
508 {
509 add( eta_1_2, eta_1_2, q_1_2_to_n );
510 add( eta, eta, q_to_n );
511 }
512
513 // check if q_to_n is zero within precision
514 if( q_to_n.is_approx_zero() )
515 eta_false = false;
516
517 // adapt the exponent N_minus and the counter n_minus
518 previous_power.assign( N_minus );
519 N_minus += 3 * n_minus + 1;
520 n_minus++;
521 }
522 else
523 {
524 // compute the difference of the exponents
525 subtract( diff, N_plus, previous_power );
526 power( q_to_diff, q_1_2, diff );
527 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
528
529 square( q_to_n, q_1_2_to_n );
530
531 // take care of the sign and update eta_1_2 and eta
532 if( n_plus & 1 )
533 {
534 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
535 subtract( eta, eta, q_to_n );
536 }
537 else
538 {
539 add( eta_1_2, eta_1_2, q_1_2_to_n );
540 add( eta, eta, q_to_n );
541 }
542
543 // check if q_to_n is zero within precision
544 if( q_to_n.is_approx_zero() )
545 eta_false = false;
546
547 // adapt the exponent N_plus and the counter n_plus
548 previous_power.assign( N_plus );
549 N_plus += 3 * n_plus + 2;
550 n_plus++;
551 }
552 }
553
554 // evaluate this loop until q_to_n is zero within precision
555 while( true )
556 {
557 if( N_minus < N_plus )
558 {
559 // compute the difference of the exponents
560 subtract( diff, N_minus, previous_power );
561 power( q_to_diff, q_1_2, diff );
562 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
563
564 square( q_to_n, q_1_2_to_n );
565
566 // take care of the sign and update eta_1_2
567 if( n_minus & 1 )
568 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
569 else
570 add( eta_1_2, eta_1_2, q_1_2_to_n );
571
572 // if q_1_2_to_n is zero within precision return the result:
573 // q^{-1/48} * eta_1_2 / eta = prefactor * eta_1_2 / eta
574 if( q_1_2_to_n.is_approx_zero() )
575 {
576 multiply( result, prefactor, eta_1_2 );
577 divide( result, result, eta );
578
579 return( result );
580 }
581
582 previous_power.assign( N_minus );
583 N_minus += 3 * n_minus + 1;
584 n_minus++;
585 }
586 else
587 {
588 // compute the difference of the exponents
589 subtract( diff, N_plus, previous_power );
590 power( q_to_diff, q_1_2, diff );
591 multiply( q_1_2_to_n, q_1_2_to_n, q_to_diff );
592
593 square( q_to_n, q_1_2_to_n );
594
595 // take care of the sign and update eta_1_2
596 if( n_plus & 1 )
597 subtract( eta_1_2, eta_1_2, q_1_2_to_n );
598 else
599 add( eta_1_2, eta_1_2, q_1_2_to_n );
600
601 // if q_1_2_to_n is zero within precision return the result:
602 // q^{-1/48} * eta_1_2 / eta = prefactor * eta_1_2 / eta
603 if( q_1_2_to_n.is_approx_zero() )
604 {
605 multiply( result, prefactor, eta_1_2 );
606 divide( result, result, eta );
607
608 return( result );
609 }
610
611 previous_power.assign( N_plus );
612 N_plus += 3 * n_plus + 2;
613 n_plus++;
614 }
615 }
616 }
617
618 //
619 // weber_f2 is described in PhD
620 //
621 bigcomplex
weber_f2(const bigcomplex & tau,bool is_fundamental)622 weber_f2( const bigcomplex & tau, bool is_fundamental )
623 {
624 // if tau is not in the fundamental domain, we make use of the
625 // relation weber_f2 = sqrt( 2 ) * eta( 2*tau ) / eta( tau )
626 if( is_fundamental == false )
627 {
628 // declare the necessary variables
629 bigcomplex current_tau, numerator, denominator, result;
630 bigfloat sqrt_2;
631
632 // compute sqrt_2 = sqrt( 2 )
633 sqrt_2 = sqrt( bigfloat( 2.0 ) );
634
635 // compute argument of numerator, numerator, and denominator
636 multiply( current_tau, tau, 2 );
637 numerator = dedekind_eta( current_tau, false );
638 denominator = dedekind_eta( tau, false );
639
640 // compute the result and return it
641 multiply( result, sqrt_2, numerator );
642 divide( result, result, denominator );
643
644 return( result );
645 }
646
647 // otherwise use algorithm of PhD
648
649 // complex variables for computation of eta(tau):
650 // * eta: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)/2}+q^{(3n^2+n)/2}),
651 // hence the current approximation of \eta(\tau)/q^{1/24}
652 // * q: stores e^(2\pi i\tau)
653 // * q_to_n: stores current power of q
654 // * q_to_diff: stores q^diff, where diff stores the difference
655 // between two non-trivial terms
656 bigcomplex eta( 1 ), q, q_to_n, q_to_diff, result;
657 exp( q, 2 * bigcomplex( 0, Pi() ) * tau );
658
659 // to store the powers of q in use:
660 // * N_minus = ( 3 * n_minus^2 - n_minus ) / 2
661 // * N_plus = ( 3 * n_plus^2 - n_plus ) / 2
662 // * previous_power stores the previous used power of q
663 bigint N_minus( 5 ), N_plus( 7 );
664 int n_minus, n_plus;
665 n_minus = n_plus = 2;
666 bigint previous_power( 2 ), diff;
667
668 // compute the prefactor \sqrt{2} * q^{1/24} of the result
669 bigcomplex prefactor;
670 exp( prefactor, 2 * bigcomplex( 0, Pi() ) * tau * bigfloat( 1, 24 ) );
671 multiply( prefactor, prefactor, sqrt( bigfloat( 2 ) ) );
672
673 // assign the correct values: eta = 1 - q - q^2
674 square( q_to_n, q );
675 subtract( eta, eta, q );
676 subtract( eta, eta, q_to_n );
677
678 // complex variables for computation of eta(2tau):
679 // * eta_2: stores 1 + \sum_{n=1}^oo (-1)^n*(q^{(3n^2-n)}+q^{(3n^2+n)})
680 // hence the current approximation of \eta(2\tau)/q^{1/12}
681 // * q_to_n_square: to store q_to_n^2
682 bigcomplex eta_2( 1 ), q_to_n_square;
683
684 // assign the correct values: eta_2 = 1 - q^2 - q^4
685 subtract( eta_2, eta_2, q_to_n );
686 subtract( eta_2, eta_2, square( q_to_n ) );
687
688 // eta_2_false is true iff q_to_n_square is not zero within precision
689 bool eta_2_false = true;
690
691 // only evaluate this loop iff q_to_n_square is not zero within precision
692 while( eta_2_false )
693 {
694 if( N_minus < N_plus )
695 {
696 // compute the difference of the exponents
697 subtract( diff, N_minus, previous_power );
698 power( q_to_diff, q, diff );
699 multiply( q_to_n, q_to_n, q_to_diff );
700
701 square( q_to_n_square, q_to_n );
702
703 // take care of the sign and update eta and eta_2
704 if( n_minus & 1 )
705 {
706 subtract( eta, eta, q_to_n );
707 subtract( eta_2, eta_2, q_to_n_square );
708 }
709 else
710 {
711 add( eta, eta, q_to_n );
712 add( eta_2, eta_2, q_to_n_square );
713 }
714
715 // check if q_to_n_square is zero within precision
716 if( q_to_n_square.is_approx_zero() )
717 eta_2_false = false;
718
719 // adapt the exponent N_minus and the counter n_minus
720 previous_power.assign( N_minus );
721 N_minus += 3 * n_minus + 1;
722 n_minus++;
723 }
724 else
725 {
726 // compute the difference of the exponents
727 subtract( diff, N_plus, previous_power );
728 power( q_to_diff, q, diff );
729 multiply( q_to_n, q_to_n, q_to_diff );
730
731 square( q_to_n_square, q_to_n );
732
733 // take care of the sign and update eta and eta_2
734 if( n_plus & 1 )
735 {
736 subtract( eta, eta, q_to_n );
737 subtract( eta_2, eta_2, q_to_n_square );
738 }
739 else
740 {
741 add( eta, eta, q_to_n );
742 add( eta_2, eta_2, q_to_n_square );
743 }
744
745 // check if q_to_n_square is zero within precision
746 if( q_to_n_square.is_approx_zero() )
747 eta_2_false = false;
748
749 // adapt the exponent N_plus and the counter n_plus
750 previous_power.assign( N_plus );
751 N_plus += 3 * n_plus + 2;
752 n_plus++;
753 }
754 }
755
756 // evaluate this loop until q_to_n is zero within precision
757 while( true )
758 {
759 if( N_minus < N_plus )
760 {
761 subtract( diff, N_minus, previous_power );
762 power( q_to_diff, q, diff );
763 multiply( q_to_n, q_to_n, q_to_diff );
764
765 if( n_minus & 1 )
766 subtract( eta, eta, q_to_n );
767 else
768 add( eta, eta, q_to_n );
769
770 // if q_to_n is zero within precision return the result:
771 // \sqrt(2) * q^{1/24} * eta_2 / eta = prefactor * eta_2 / eta
772 if( q_to_n.is_approx_zero() )
773 {
774 multiply( result, prefactor, eta_2 );
775 divide( result, result, eta );
776
777 return( result );
778 }
779
780 previous_power.assign( N_minus );
781 N_minus += 3 * n_minus + 1;
782 n_minus++;
783 }
784
785 else
786 {
787 subtract( diff, N_plus, previous_power );
788 power( q_to_diff, q, diff );
789 multiply( q_to_n, q_to_n, q_to_diff );
790
791 if( n_plus & 1 )
792 subtract( eta, eta, q_to_n );
793 else
794 add( eta, eta, q_to_n );
795
796 // if q_to_n is zero within precision return the result:
797 // \sqrt(2) * q^{1/24} * eta_2 / eta = prefactor * eta_2 / eta
798 if( q_to_n.is_approx_zero() )
799 {
800 multiply( result, prefactor, eta_2 );
801 divide( result, result, eta );
802
803 return( result );
804 }
805
806 previous_power.assign( N_plus );
807 N_plus += 3 * n_plus + 2;
808 n_plus++;
809 }
810 }
811 }
812
813
814 /********************************************************************
815 *
816 * Compute \gamma_2 and the modular function j using weber_f2.
817 *
818 *******************************************************************/
819 //
820 // The function gamma2 bases on the direct computation of
821 // weber_f2 via the eta-function.
822 //
823 bigcomplex
gamma_2(const bigcomplex & tau,bool is_fundamental)824 gamma_2( const bigcomplex & tau, bool is_fundamental )
825 {
826 // assign:
827 // f2_to_8 = weber_f2( tau ) ^ 8
828 // assign gamma2 = (f2^24 + 16 ) / f2^8 to result
829 bigcomplex f2_to_8, result, local_tau( tau ), tau_inverse, factor( 1 );
830
831 // if is_fundamental = false translate tau to fundamental area
832 // and compute appropriate factor
833 if( is_fundamental == false )
834 {
835 //
836 // ensure that |Re( tau )| <= 1/2 and |tau| >= 1
837 // otherwise shift tau to this stripe
838 //
839 bigint real_part_as_bigint;
840
841 while( true )
842 {
843 // assign to real_part_as_bigint the nearest integer
844 // of Re(local_tau)
845 local_tau.real().bigintify( real_part_as_bigint );
846
847 // if it is not zero, shift it in stripe |Re( tau )| <= 1/2
848 // gamma_2(tau+1) = 1/zeta_3 * gamma(tau)
849 if( ! real_part_as_bigint.is_zero() )
850 {
851 subtract( local_tau, local_tau,
852 bigcomplex( bigfloat( real_part_as_bigint ) ) );
853 multiply( factor, factor,
854 exp( -2 * bigcomplex( 0, Pi() ) *
855 bigfloat( real_part_as_bigint, 3 ) ) );
856 }
857
858 // if |local_tau| < 1 invert it to put it out of unit circle
859 // gamma_2(-1/tau) = gamma(tau)
860 if( abs( local_tau ) <= 0.9999 )
861 {
862 invert( tau_inverse, local_tau );
863 local_tau.assign( tau_inverse );
864 local_tau.negate();
865 }
866 // otherwise we are done
867 else
868 break;
869 }
870 }
871
872 power( f2_to_8, weber_f2( local_tau, true ), 8 );
873 power( result, f2_to_8, 3 );
874 add( result, result, bigint( 16 ) );
875 divide( result, result, f2_to_8 );
876
877 multiply( result, result, factor );
878
879 return( result );
880 }
881
882
883 bigcomplex
modular_j(const bigcomplex & tau,bool is_fundamental)884 modular_j( const bigcomplex & tau, bool is_fundamental )
885 {
886 // assign:
887 // f2_to_24 = weber_f2( tau ) ^ 24
888 // j = (f2^24 + 16 )^3 / f2^24 to result
889 //
890 bigcomplex f2_to_24, result, local_tau( tau ), tau_inverse;
891
892 if( is_fundamental == false )
893 {
894 //
895 // ensure that |Re( tau )| <= 1/2 and |tau| >= 1
896 // otherwise shift tau to this stripe
897 //
898 bigint real_part_as_bigint;
899
900 while( true )
901 {
902 // assign to real_part_as_bigint the nearest integer
903 // of Re(local_tau)
904 local_tau.real().bigintify( real_part_as_bigint );
905
906 // if it is not zero, shift it in stripe |Re( tau )| <= 1/2
907 // j(tau+1) = j(tau)
908 if( ! real_part_as_bigint.is_zero() )
909 {
910 subtract( local_tau, local_tau,
911 bigcomplex( bigfloat( real_part_as_bigint ) ) );
912 }
913
914 // if |local_tau| < 1 invert it to put it out of unit circle
915 // j(-1/tau) = j(tau)
916 if( abs( local_tau ) <= 0.9999 )
917 {
918 invert( tau_inverse, local_tau );
919 local_tau.assign( tau_inverse );
920 local_tau.negate();
921 }
922 // otherwise we are done
923 else
924 break;
925 }
926
927 }
928
929 power( f2_to_24, weber_f2( local_tau, true ), 24 );
930 add( result, f2_to_24, bigint( 16 ) );
931 power( result, result, 3 );
932 divide( result, result, f2_to_24 );
933
934 return( result );
935 }
936
937
938 #ifdef LIDIA_NAMESPACE
939 } // end of namespace LiDIA
940 # undef IN_NAMESPACE_LIDIA
941 #endif
942