Skip to content

Commit 47d9094

Browse files
unhyperbolicunhyperbolic
committed
Extending the coefficient table for the Lobachevsky function in the kernel so that we can use the kernel's volume function for quad-double precision as well.
1 parent a3c1ba0 commit 47d9094

4 files changed

Lines changed: 142 additions & 66 deletions

File tree

‎src/snappy/extensions/SnapPy/cython_src/core/manifold.pyx‎

Lines changed: 6 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -662,7 +662,7 @@ cdef class Manifold(Triangulation):
662662
manifold_class=self.__class__)
663663
for cover in covers]
664664

665-
cdef _real_volume(self):
665+
cdef _kernel_volume(self):
666666
"""
667667
Return the (real) volume as a Number.
668668
"""
@@ -673,14 +673,8 @@ cdef class Manifold(Triangulation):
673673
solution_type = self.solution_type()
674674
if solution_type in ('not attempted', 'no solution found'):
675675
raise ValueError('Solution type is: %s' % solution_type)
676-
if Number._default_precision > 64:
677-
# must provide a start value to get the correct precision
678-
result = sum(
679-
[z.volume() for z in self._get_tetrahedra_shapes('filled')],
680-
Number(0))
681-
else:
682-
result = Number(<double>volume(self.c_triangulation, &acc))
683-
result.accuracy = acc
676+
result = Real2Number(volume(self.c_triangulation, &acc))
677+
result.accuracy = acc
684678
return result
685679

686680
cdef _cusped_complex_volume(self, Complex *volume, int *accuracy):
@@ -720,7 +714,7 @@ cdef class Manifold(Triangulation):
720714
result = Complex2Number(volume)
721715
result.accuracy = accuracy
722716
else:
723-
result = self._real_volume() + (
717+
result = self._kernel_volume() + (
724718
self._chern_simons() *
725719
Real2Number(PI_SQUARED_BY_2) * Number('I'))
726720
return self._number_(result)
@@ -802,7 +796,7 @@ cdef class Manifold(Triangulation):
802796
number of digits.
803797
804798
>>> vol, accuracy = M.volume(accuracy = True)
805-
>>> accuracy in (10, 63) # Low precision, High precision
799+
>>> accuracy in (10, 59) # Low precision, High precision
806800
True
807801
808802
Inside SageMath, verified computation of the volume of a
@@ -822,7 +816,7 @@ cdef class Manifold(Triangulation):
822816
return verify.compute_volume(
823817
self, verified=verified, bits_prec=bits_prec)
824818

825-
vol = self._real_volume()
819+
vol = self._kernel_volume()
826820
if accuracy:
827821
return (self._number_(vol), vol.accuracy)
828822
else:

‎src/snappy/extensions/SnapPy/kernel/addl_code/dilog.c‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -432,12 +432,4 @@ Complex complex_volume_log(Complex z)
432432
return result;
433433
}
434434

435-
Real Lobachevsky_via_dilog(Real t)
436-
{
437-
Complex z;
438-
z.real = cos(2.0*t);
439-
z.imag = sin(2.0*t);
440-
return 0.5*complex_volume_dilog(z).imag;
441-
}
442-
443435
SNAPPEA_NAMESPACE_END_SCOPE

‎src/snappy/extensions/SnapPy/kernel/addl_code/dilog.h‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#define _dilog_
33

44
#include "kernel.h"
5+
56
SNAPPEA_NAMESPACE_BEGIN_SCOPE
67

78
/* Logarithm of z with the standard branch cut (i.e, negative real axis such
@@ -17,9 +18,6 @@ Complex complex_volume_log(Complex z);
1718

1819
Complex complex_volume_dilog(Complex z);
1920

20-
21-
Real Lobachevsky_via_dilog(Real t);
22-
2321
SNAPPEA_NAMESPACE_END_SCOPE
2422

2523
#endif

‎src/snappy/extensions/SnapPy/kernel/kernel_code/volume.c‎

Lines changed: 135 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -25,10 +25,6 @@
2525
#include "kernel.h"
2626
SNAPPEA_NAMESPACE_BEGIN_SCOPE
2727

28-
#ifdef _QD_REAL_SNAPPY_
29-
#include "dilog.h"
30-
#endif
31-
3228
static Real Lobachevsky(Real theta);
3329

3430

@@ -47,13 +43,13 @@ Real volume(
4743
for (tet = manifold->tet_list_begin.next;
4844
tet != &manifold->tet_list_end;
4945
tet = tet->next)
50-
46+
5147
if (tet->shape[filled] != NULL)
5248

5349
for (i = 0; i < 2; i++) /* i = ultimate, penultimate */
5450

5551
for (j = 0; j < 3; j++)
56-
52+
5753
vol[i] += Lobachevsky(tet->shape[filled]->cwl[i][j].log.imag);
5854

5955

@@ -80,47 +76,143 @@ static Real Lobachevsky(Real theta)
8076
theta_over_pi_squared;
8177
const Real *lobcoefptr;
8278

83-
const static Real lobcoef[30] = {
84-
5.4831135561607547882413838888201e-1,
85-
1.0823232337111381915160036965412e-1,
86-
4.8444907713545197129262758561472e-2,
87-
2.7891037672165120538296812180796e-2,
88-
1.8199901365960328824311744707278e-2,
89-
1.2823667776324462157674846128817e-2,
90-
9.5243928393815114745643670962394e-3,
91-
7.3530535460250636167039159668208e-3,
92-
5.8479755397266959054962366178048e-3,
93-
4.7619093045811136799814912001842e-3,
94-
3.9525701124525799515137944808229e-3,
95-
3.3333335320272968375315987081340e-3,
96-
2.8490028914574211634331659107080e-3,
97-
2.4630541963678177950454607261557e-3,
98-
2.1505376364114568439132649094016e-3,
99-
1.8939393943803620897114593734851e-3,
100-
1.6806722690053911275277764855335e-3,
101-
1.5015015015233512340706336099638e-3,
102-
1.3495276653220485553945730785293e-3,
103-
1.2195121951230603594927151084491e-3,
104-
1.1074197120711266596728487987281e-3,
105-
1.0101010101010675186059356321779e-3,
106-
9.2506938020352840967144128733280e-4,
107-
8.5034013605442478972252664720549e-4,
108-
7.8431372549019677504189889653457e-4,
109-
7.2568940493468811469129541350086e-4,
110-
6.7340067340067343805464730535677e-4,
111-
6.2656641604010025932192218654462e-4,
112-
5.8445353594389246257711686275038e-4,
113-
5.4644808743169398954500641421420e-4};
79+
/*
80+
* Coefficients when using (theta/pi)^2 (rather than 2theta as in Milnor's
81+
* article are: lobcoef[n-1] = B_n / (2n * (2n+1)!).
82+
*/
83+
84+
const static Real lobcoef[] = {
85+
Real_literal(0.54831135561607547882413838888200839640631663373559947924518607645666916),
86+
Real_literal(0.10823232337111381915160036965411679027747509519187269076829762154441206),
87+
Real_literal(0.048444907713545197129262758561472406090562737620612074373448031619253913),
88+
Real_literal(0.027891037672165120538296812180795901812748910851384722786919727851460638),
89+
Real_literal(0.018199901365960328824311744707278527581927627846626863950141618084296208),
90+
Real_literal(0.012823667776324462157674846128817175268723283185359018006833858360682718),
91+
Real_literal(0.0095243928393815114745643670962393841309283971063729004712643352595510756),
92+
Real_literal(0.0073530535460250636167039159668208582501708631543417024349493029675646109),
93+
Real_literal(0.0058479755397266959054962366178048066650903783463096686031764586548984213),
94+
Real_literal(0.0047619093045811136799814912001842069016473514008924314083336023326244644),
95+
Real_literal(0.0039525701124525799515137944808228725414819566870426073449401170083655117),
96+
Real_literal(0.0033333335320272968375315987081340264526707583463961243426528808232410724),
97+
Real_literal(0.0028490028914574211634331659107079960988564924922739825809417305868056694),
98+
Real_literal(0.0024630541963678177950454607261556749221734562386528554084027580525485374),
99+
Real_literal(0.0021505376364114568439132649094016445914434834700386131974541613088459364),
100+
Real_literal(0.0018939393943803620897114593734851242378675179971268603873680134045497318),
101+
Real_literal(0.0016806722690053911275277764855335073275480051330640716043986822774314086),
102+
Real_literal(0.0015015015015233512340706336099638583066746522837594728140000571608701835),
103+
Real_literal(0.0013495276653220485553945730785293350018297884932694278385761099459486526),
104+
Real_literal(0.0012195121951230603594927151084490894046813352545778780499888766201286338),
105+
Real_literal(0.0011074197120711266596728487987281535793552578838553521397551670792696940),
106+
Real_literal(0.0010101010101010675186059356321778714866966416572581116884737311954715217),
107+
Real_literal(0.00092506938020352840967144128733281205763842899124669629843259385009940028),
108+
Real_literal(0.00085034013605442478972252664720550450550130447223043549800638015528487459),
109+
Real_literal(0.00078431372549019677504189889653458063498755971697147651834285997213070194),
110+
Real_literal(0.00072568940493468811469129541350086966613353368683991032196330402314473094),
111+
Real_literal(0.00067340067340067343805464730535677605201169421218471570149734048250963608),
112+
Real_literal(0.00062656641604010025932192218654463205691672613189236961636505306150304972),
113+
Real_literal(0.00058445353594389246257711686275039311791015283232573619050290510052940166),
114+
Real_literal(0.00054644808743169398954500641421420402887115055012803013595120991570526388),
115+
Real_literal(0.00051203277009728622642951379134650382923778341875683484772951362214662668),
116+
Real_literal(0.00048076923076923076925683178299258002601403221348263766209165106917992326),
117+
Real_literal(0.00045228403437358661239871213349439220469538002442190272651695245467437249),
118+
Real_literal(0.00042625745950554134697501625395951440614581658273676506258910109058253172),
119+
Real_literal(0.00040241448692152917505064266919406627759952704195610391016186252974479011),
120+
Real_literal(0.00038051750380517503805183095823318630877975054069405700162050767736825230),
121+
Real_literal(0.00036036036036036036036037943767899221574891464472901952614625815337966503),
122+
Real_literal(0.00034176349965823650034176802286049242959357590379707824504090373336509638),
123+
Real_literal(0.00032456994482310938007140646177294716396351475071612603518494296542986160),
124+
Real_literal(0.00030864197530864197530864223061130017068964208861437088716272851881731828),
125+
Real_literal(0.00029385836027034969144872177690130859190645619735433204075507035958070047),
126+
Real_literal(0.00028011204481792717086834735341702753069028381067264206340490557243926157),
127+
Real_literal(0.00026730820636193531141406041510951538923845674963408877795559083381442895),
128+
Real_literal(0.00025536261491317671092951991910908433804526695253698011185563263977549274),
129+
Real_literal(0.00024420024420024420024420024439750758655587697272670911321955799120557960),
130+
Real_literal(0.00023375409069658719027582982706917914720350141068627308664252699279781143),
131+
Real_literal(0.00022396416573348264277715565510649207120585548868104092087916833244365033),
132+
Real_literal(0.00021477663230240549828178694158346687596295880348877877999587138073681975),
133+
Real_literal(0.00020614306328592042877757163471514233016303526020183670196216961395205878),
134+
Real_literal(0.00019801980198019801980198019801995819027826158649650535718470614680570092),
135+
Real_literal(0.00019036740909956215495907100704363168003546644830602865483193454009135441),
136+
Real_literal(0.00018315018315018315018315018315019218018435463612415058382191296281604803),
137+
Real_literal(0.00017633574325515782049021336624934091446863764385636787513395261454155277),
138+
Real_literal(0.00016989466530750934420659191301393188608357305803217213477759251626113235),
139+
Real_literal(0.00016380016380016380016380016380016392635085630710800020561027927419353004),
140+
Real_literal(0.00015802781289506953223767383059418460692051113127743130677131138961607846),
141+
Real_literal(0.00015255530129672006102212051868802441619347804876901545794128487662694200),
142+
Real_literal(0.00014736221632773356911287945770704391571427382145891896281239861597518719),
143+
Real_literal(0.00014242985329725110383136305369605469349227539354864777387907333055519734),
144+
Real_literal(0.00013774104683195592286501377410468319569591134773089035882927656524241241),
145+
Real_literal(0.00013328002132480341196854591496734639480049035001087111838931889604806538),
146+
Real_literal(0.00012903225806451612903225806451612903226413158374617954842822662343359260),
147+
Real_literal(0.00012498437695288088988876390451193600800046930927486850082235822131221729),
148+
Real_literal(0.00012112403100775193798449612403100775193834044765952709771922503950573763),
149+
Real_literal(0.00011743981209630064591896652965355255431599937568634925774427251553256592),
150+
Real_literal(0.00011392116655274550011392116655274550011394209056645203718650263341661608),
151+
Real_literal(0.00011055831951354339414040906578220011055832459013245760039864898052927931),
152+
Real_literal(0.00010734220695577501073422069557750107342207078997821165187864373198858670),
153+
Real_literal(0.00010426441455531227192159316025440517151496224271288238319522064975493931),
154+
Real_literal(0.00010131712259371833839918946301925025329280655698731648777166568520260448),
155+
Real_literal(0.000098493056239535112774549394267704126859056454187550815874774458607052536),
156+
Real_literal(0.000095785440613026819923371647509578544061302686287504871493692171413730712),
157+
Real_literal(0.000093187960115553070543285807473674401267356258616194922616581545745022541),
158+
Real_literal(0.000090694721567204788681298748412842372573916198331452650884240353876968909),
159+
Real_literal(0.000088300220750551876379690949227373068432671081739571676129130996508773933),
160+
Real_literal(0.000085999312005503955968352253181974544203646370844097205713846072743848142),
161+
Real_literal(0.000083787180561374109761206535400083787180561374113430295518236324546825852),
162+
Real_literal(0.000081659317328107137024334476563775926833251674016899203352322197727142874),
163+
Real_literal(0.000079611495900007961149590000796114959000079611496117897545031001380594531),
164+
Real_literal(0.000077639751552795031055900621118012422360248447205022067373119844882407762),
165+
Real_literal(0.000075740362038930546088010300689237294554267969400906692186733765092884420),
166+
Real_literal(0.000073909830007390983000739098300073909830007390983003899795666036594999131),
167+
Real_literal(0.000072144866892720582930524493182310078637904913065436165576714549367475980),
168+
Real_literal(0.000070442378134685826993519301211608903916596224288532169115575520723840951),
169+
Real_literal(0.000068799449604403164774681802545579635362917096663226740157631752531744860),
170+
Real_literal(0.000067213335125688936685038311601021642693910471837612593564174914690296269),
171+
Real_literal(0.000065681444991789819376026272577996715927750410509031201429356229569166303),
172+
Real_literal(0.000064201335387776065742167437082691319979455572675911659632799637185261893),
173+
Real_literal(0.000062770698637875839558094281589354089511016257610947210006284714531379688),
174+
Real_literal(0.000061387354205033763044812768569674647022713321055862492366637842038051102),
175+
Real_literal(0.000060049240377109229568245961688584639404311535459076442692796049462557084),
176+
Real_literal(0.000058754406580493537015276145710928319623971797884841363104628857398648107),
177+
Real_literal(0.000057501006267609683169455465470645736300385256741992984877821619701641521),
178+
Real_literal(0.000056287290329843521332883035010694585162670269053247776652175444503985811),
179+
Real_literal(0.000055111600992008817856158721410856985395425737117663268117973945261198842),
180+
Real_literal(0.000053972366148531951640759930915371329879101899827288428324706353043561812),
181+
Real_literal(0.000052868094105207507269362939466032249537404176579434311393076385265544716),
182+
Real_literal(0.000051797368693670361545633481819123588521703097482647881487620944619519558),
183+
Real_literal(0.000050758844728693974925130704025176386985432211562864829196487614294000870),
184+
Real_literal(0.000049751243781094527363184079601990049751243781094527363184079632950324765),
185+
Real_literal(0.000048773350241428083695069014290591620738428522655221187143344883947489731),
186+
Real_literal(0.000047824007651841224294595887135341941654710664753706360593017696742890745),
187+
Real_literal(0.000046902115285399371511655175648421743820646311148632803339430608776486158),
188+
Real_literal(0.000046006624953993375046006624953993375046006624953993375046006625065829285),
189+
Real_literal(0.000045136538027533288196795305800045136538027533288196795305800045163968239),
190+
Real_literal(0.000044290902648595978386039507485162547612720347240676764992470546556467760),
191+
Real_literal(0.000043468811128015648772006085633557922190828080851988698109106715932970323),
192+
Real_literal(0.000042669397508107185526540365250042669397508107185526540365250042669802678),
193+
Real_literal(0.000041891835281303673913954170332202253780738134137656570734363872481352850),
194+
Real_literal(0.000041135335252982311805841217605923488276429452900041135335252982311830254),
195+
Real_literal(0.000040399143538156991071789278067304973134569547125600937260130085242198859),
196+
Real_literal(0.000039682539682539682539682539682539682539682539682539682539682539682541154),
197+
Real_literal(0.000038984834899224201785505438384468441776149078008654633347627772796382569),
198+
Real_literal(0.000038305370412931893051405807094154600474986593120355473837432007967517135),
199+
Real_literal(0.000037643515904385469602860907208733295689817428947863730472426124600037665),
200+
Real_literal(0.000036998668047950273790143554832026047062305756992748261062601746337131869),
201+
Real_literal(0.000036370249136206583015093653391525731951263866157483178759774504455355521),
202+
Real_literal(0.000035757705785596796109561610527068583279696774654938139168990917542730459),
203+
Real_literal(0.000035160507717731444042051967230406807074294152807566541260855806757849583),
204+
Real_literal(0.000034578146611341632088520055325034578146611341632088520055325034578146611)
205+
};
114206

115207
/*
116-
* As long as DBL_EPSILON > 5e-22 there will be enough lobcoefs
117-
* for the series. If DBL_EPSILON is smaller than this you'll
208+
* As long as REAL_EPSILON > 2e-77 there will be enough lobcoefs
209+
* for the series. If REAL_EPSILON is smaller than this you'll
118210
* need to add more coefficients to the list.
119211
*/
120212

121-
#ifdef _QD_REAL_SNAPPY_
122-
return Lobachevsky_via_dilog(theta);
123-
#endif
213+
#if REAL_DIG > 64
214+
#error "No support for precision beyond quad-double."
215+
#endif
124216

125217
/*
126218
* Milnor (Lemma 1, p. 17) shows that the Lobachevsky function is
@@ -159,7 +251,7 @@ static Real Lobachevsky(Real theta)
159251
term = *lobcoefptr++ * product;
160252
sum += term;
161253
}
162-
while (term > DBL_EPSILON);
254+
while (term > REAL_EPSILON);
163255

164256
return theta*(1.0 - log(2.0*theta) + sum);
165257
}

0 commit comments

Comments
 (0)