1 /++
2 This is a submodule of $(MREF mir,ndslice).
3 
4 Field is a type with `opIndex()(ptrdiff_t index)` primitive.
5 An iterator can be created on top of a field using $(SUBREF iterator, FieldIterator).
6 An ndslice can be created on top of a field using $(SUBREF slice, slicedField).
7 
8 $(BOOKTABLE $(H2 Fields),
9 $(TR $(TH Field Name) $(TH Used By))
10 $(T2 BitField, $(SUBREF topology, bitwise))
11 $(T2 BitpackField, $(SUBREF topology, bitpack))
12 $(T2 CycleField, $(SUBREF topology, cycle) (2 kinds))
13 $(T2 LinspaceField, $(SUBREF topology, linspace))
14 $(T2 MagicField, $(SUBREF topology, magic))
15 $(T2 MapField, $(SUBREF topology, map) and $(SUBREF topology, mapField))
16 $(T2 ndIotaField, $(SUBREF topology, ndiota))
17 $(T2 OrthogonalReduceField, $(SUBREF topology, orthogonalReduceField))
18 $(T2 RepeatField, $(SUBREF topology, repeat))
19 $(T2 SparseField, Used for mutable DOK sparse matrixes )
20 )
21 
22 
23 
24 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
25 Copyright: 2020 Ilya Yaroshenko, Kaleidic Associates Advisory Limited, Symmetry Investments
26 Authors: Ilya Yaroshenko
27 
28 Macros:
29 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, ndslice, $1)$(NBSP)
30 T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
31 +/
32 module mir.ndslice.field;
33 
34 import mir.internal.utility: Iota;
35 import mir.math.common: optmath;
36 import mir.ndslice.internal;
37 import mir.qualifier;
38 
39 @optmath:
40 
41 package template ZeroShiftField(T)
42 {
43     static if (hasZeroShiftFieldMember!T)
44         alias ZeroShiftField = typeof(T.init.assumeFieldsHaveZeroShift());
45     else
46         alias ZeroShiftField = T;
47 }
48 
49 package enum hasZeroShiftFieldMember(T) = __traits(hasMember, T, "assumeFieldsHaveZeroShift");
50 
51 package auto applyAssumeZeroShift(Types...)()
52 {
53     import mir.ndslice.topology;
54     string str;
55     foreach(i, T; Types)
56         static if (hasZeroShiftFieldMember!T)
57             str ~= "_fields[" ~ i.stringof ~ "].assumeFieldsHaveZeroShift, ";
58         else
59             str ~= "_fields[" ~ i.stringof ~ "], ";
60     return str;
61 }
62 
63 auto MapField__map(Field, alias fun, alias fun1)(ref MapField!(Field, fun) f)
64 {
65     import core.lifetime: move;
66     import mir.functional: pipe;
67     return MapField!(Field, pipe!(fun, fun1))(move(f._field));
68 }
69 
70 
71 /++
72 `MapField` is used by $(SUBREF topology, map).
73 +/
74 struct MapField(Field, alias _fun)
75 {
76 @optmath:
77     ///
78     Field _field;
79 
80     ///
81     auto lightConst()() const @property
82     {
83         return MapField!(LightConstOf!Field, _fun)(.lightConst(_field));
84     }
85 
86     ///
87     auto lightImmutable()() immutable @property
88     {
89         return MapField!(LightImmutableOf!Field, _fun)(.lightImmutable(_field));
90     }
91 
92     /++
93     User defined constructor used by $(LREF mapField).
94     +/
95     static alias __map(alias fun1) = MapField__map!(Field, _fun, fun1);
96 
97     auto ref opIndex(T...)(auto ref T index)
98     {
99         import mir.functional: RefTuple, unref;
100         static if (is(typeof(_field[index]) : RefTuple!K, K...))
101         {
102             auto t = _field[index];
103             return mixin("_fun(" ~ _iotaArgs!(K.length, "t.expand[", "].unref, ") ~ ")");
104         }
105         else
106             return _fun(_field[index]);
107     }
108 
109     static if (__traits(hasMember, Field, "length"))
110     auto length() const @property
111     {
112         return _field.length;
113     }
114 
115     static if (__traits(hasMember, Field, "shape"))
116     auto shape() const @property
117     {
118         return _field.shape;
119     }
120 
121     static if (__traits(hasMember, Field, "elementCount"))
122     auto elementCount() const @property
123     {
124         return _field.elementCount;
125     }
126 
127     static if (hasZeroShiftFieldMember!Field)
128     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
129     auto assumeFieldsHaveZeroShift() @property
130     {
131         return _mapField!_fun(_field.assumeFieldsHaveZeroShift);
132     }
133 }
134 
135 /++
136 `VmapField` is used by $(SUBREF topology, map).
137 +/
138 struct VmapField(Field, Fun)
139 {
140 @optmath:
141     ///
142     Field _field;
143     ///
144     Fun _fun;
145 
146     ///
147     auto lightConst()() const @property
148     {
149         return VmapField!(LightConstOf!Field, _fun)(.lightConst(_field));
150     }
151 
152     ///
153     auto lightImmutable()() immutable @property
154     {
155         return VmapField!(LightImmutableOf!Field, _fun)(.lightImmutable(_field));
156     }
157 
158     auto ref opIndex(T...)(auto ref T index)
159     {
160         import mir.functional: RefTuple, unref;
161         static if (is(typeof(_field[index]) : RefTuple!K, K...))
162         {
163             auto t = _field[index];
164             return mixin("_fun(" ~ _iotaArgs!(K.length, "t.expand[", "].unref, ") ~ ")");
165         }
166         else
167             return _fun(_field[index]);
168     }
169 
170     static if (__traits(hasMember, Field, "length"))
171     auto length() const @property
172     {
173         return _field.length;
174     }
175 
176     static if (__traits(hasMember, Field, "shape"))
177     auto shape() const @property
178     {
179         return _field.shape;
180     }
181 
182     static if (__traits(hasMember, Field, "elementCount"))
183     auto elementCount()const @property
184     {
185         return _field.elementCount;
186     }
187 
188     static if (hasZeroShiftFieldMember!Field)
189     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
190     auto assumeFieldsHaveZeroShift() @property
191     {
192         return _vmapField(_field.assumeFieldsHaveZeroShift, _fun);
193     }
194 }
195 
196 /+
197 Creates a mapped field. Uses `__map` if possible.
198 +/
199 auto _mapField(alias fun, Field)(Field field)
200 {
201     import mir.functional: naryFun;
202     static if ((
203         __traits(isSame, fun, naryFun!"a|b") ||
204         __traits(isSame, fun, naryFun!"a^b") ||
205         __traits(isSame, fun, naryFun!"a&b") ||
206         __traits(isSame, fun, naryFun!"a | b") ||
207         __traits(isSame, fun, naryFun!"a ^ b") ||
208         __traits(isSame, fun, naryFun!"a & b")) &&
209         is(Field : ZipField!(BitField!(LeftField, I), BitField!(RightField, I)), LeftField, RightField, I))
210     {
211         import mir.ndslice.topology: bitwiseField;
212         auto f = ZipField!(LeftField, RightField)(field._fields[0]._field, field._fields[1]._field)._mapField!fun;
213         return f.bitwiseField!(typeof(f), I);
214     }
215     else
216     static if (__traits(hasMember, Field, "__map"))
217         return Field.__map!fun(field);
218     else
219         return MapField!(Field, fun)(field);
220 }
221 
222 /+
223 Creates a mapped field. Uses `__vmap` if possible.
224 +/
225 auto _vmapField(Field, Fun)(Field field, Fun fun)
226 {
227     static if (__traits(hasMember, Field, "__vmap"))
228         return Field.__vmap(field, fun);
229     else
230         return VmapField!(Field, Fun)(field, fun);
231 }
232 
233 /++
234 Iterates multiple fields in lockstep.
235 
236 `ZipField` is used by $(SUBREF topology, zipFields).
237 +/
238 struct ZipField(Fields...)
239     if (Fields.length > 1)
240 {
241 @optmath:
242     import mir.functional: RefTuple, Ref, _ref;
243     import std.meta: anySatisfy;
244 
245     ///
246     Fields _fields;
247 
248     ///
249     auto lightConst()() const @property
250     {
251         import std.format;
252         import mir.ndslice.topology: iota;
253         import std.meta: staticMap;
254         return mixin("ZipField!(staticMap!(LightConstOf, Fields))(%(_fields[%s].lightConst,%)].lightConst)".format(_fields.length.iota));
255     }
256 
257     ///
258     auto lightImmutable()() immutable @property
259     {
260         import std.format;
261         import mir.ndslice.topology: iota;
262         import std.meta: staticMap;
263         return mixin("ZipField!(staticMap!(LightImmutableOf, Fields))(%(_fields[%s].lightImmutable,%)].lightImmutable)".format(_fields.length.iota));
264     }
265 
266     auto opIndex()(ptrdiff_t index)
267     {
268         alias Iterators = Fields;
269         alias _iterators = _fields;
270         import mir.ndslice.iterator: _zip_types, _zip_index;
271         return mixin("RefTuple!(_zip_types!Fields)(" ~ _zip_index!Fields ~ ")");
272     }
273 
274     auto opIndexAssign(Types...)(RefTuple!(Types) value, ptrdiff_t index)
275         if (Types.length == Fields.length)
276     {
277         foreach(i, ref val; value.expand)
278         {
279             _fields[i][index] = val;
280         }
281         return opIndex(index);
282     }
283 
284     static if (anySatisfy!(hasZeroShiftFieldMember, Fields))
285     /// Defined if at least one of `Fields` has member `assumeFieldsHaveZeroShift`.
286     auto assumeFieldsHaveZeroShift() @property
287     {
288         import std.meta: staticMap;
289         return mixin("ZipField!(staticMap!(ZeroShiftField, Fields))(" ~ applyAssumeZeroShift!Fields ~ ")");
290     }
291 }
292 
293 /++
294 `RepeatField` is used by $(SUBREF topology, repeat).
295 +/
296 struct RepeatField(T)
297 {
298     import std.traits: Unqual;
299 
300 @optmath:
301     alias UT = Unqual!T;
302 
303     ///
304     UT _value;
305 
306     ///
307     auto lightConst()() const @property @trusted
308     {
309         return RepeatField!(const T)(cast(UT) _value);
310     }
311 
312     ///
313     auto lightImmutable()() immutable @property @trusted
314     {
315         return RepeatField!(immutable T)(cast(UT) _value);
316     }
317 
318     auto ref T opIndex()(ptrdiff_t) @trusted
319     { return cast(T) _value; }
320 }
321 
322 /++
323 `BitField` is used by $(SUBREF topology, bitwise).
324 +/
325 struct BitField(Field, I = typeof(cast()Field.init[size_t.init]))
326     if (__traits(isUnsigned, I))
327 {
328 @optmath:
329     import mir.bitop: ctlz;
330     package(mir) alias E = I;
331     package(mir) enum shift = ctlz(I.sizeof) + 3;
332 
333     ///
334     Field _field;
335 
336     /// optimization for bitwise operations
337     auto __vmap(Fun : LeftOp!(op, bool), string op)(Fun fun)
338         if (op == "|" || op == "&" || op == "^")
339     {
340         import mir.ndslice.topology: bitwiseField;
341         return _vmapField(_field, RightOp!(op, I)(I(0) - fun.value)).bitwiseField;
342     }
343 
344     /// ditto
345     auto __vmap(Fun : RightOp!(op, bool), string op)(Fun fun)
346         if (op == "|" || op == "&" || op == "^")
347     {
348         import mir.ndslice.topology: bitwiseField;
349         return _vmapField(_field, RightOp!(op, I)(I(0) - fun.value)).bitwiseField;
350     }
351 
352     /// ditto
353     auto __vmap(Fun)(Fun fun)
354     {
355         return VmapField!(typeof(this), Fun)(this, fun);
356     }
357 
358     /// ditto
359     alias __map(alias fun) = BitField__map!(Field, I, fun);
360 
361     ///
362     auto lightConst()() const @property
363     {
364         return BitField!(LightConstOf!Field, I)(mir.qualifier.lightConst(_field));
365     }
366 
367     ///
368     auto lightImmutable()() immutable @property
369     {
370         return BitField!(LightImmutableOf!Field, I)(mir.qualifier.lightImmutable(_field));
371     }
372 
373     bool opIndex()(size_t index)
374     {
375         import mir.bitop: bt;
376         return bt!(Field, I)(_field, index) != 0;
377     }
378 
379     bool opIndexAssign()(bool value, size_t index)
380     {
381         import mir.bitop: bta;
382         bta!(Field, I)(_field, index, value);
383         return value;
384     }
385 
386     static if (hasZeroShiftFieldMember!Field)
387     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
388     auto assumeFieldsHaveZeroShift() @property
389     {
390         return BitField!(ZeroShiftField!Field, I)(_field.assumeFieldsHaveZeroShift);
391     }
392 }
393 
394 ///
395 version(mir_test) unittest
396 {
397     import mir.ndslice.iterator: FieldIterator;
398     ushort[10] data;
399     auto f = FieldIterator!(BitField!(ushort*))(0, BitField!(ushort*)(data.ptr));
400     f[123] = true;
401     f++;
402     assert(f[122]);
403 }
404 
405 auto BitField__map(Field, I, alias fun)(BitField!(Field, I) field)
406 {
407     import core.lifetime: move;
408     import mir.functional: naryFun;
409     static if (__traits(isSame, fun, naryFun!"~a") || __traits(isSame, fun, naryFun!"!a"))
410     {
411         import mir.ndslice.topology: bitwiseField;
412         auto f = _mapField!(naryFun!"~a")(move(field._field));
413         return f.bitwiseField!(typeof(f), I);
414     }
415     else
416     {
417         return MapField!(BitField!(Field, I), fun)(move(field));
418     }
419 }
420 
421 /++
422 `BitpackField` is used by $(SUBREF topology, bitpack).
423 +/
424 struct BitpackField(Field, uint pack, I = typeof(cast()Field.init[size_t.init]))
425     if (__traits(isUnsigned, I))
426 {
427     //static assert();
428 @optmath:
429     package(mir) alias E = I;
430     package(mir) enum mask = (I(1) << pack) - 1;
431     package(mir) enum bits = I.sizeof * 8;
432 
433     ///
434     Field _field;
435 
436     ///
437     auto lightConst()() const @property
438     {
439         return BitpackField!(LightConstOf!Field, pack)(.lightConst(_field));
440     }
441 
442     ///
443     auto lightImmutable()() immutable @property
444     {
445         return BitpackField!(LightImmutableOf!Field, pack)(.lightImmutable(_field));
446     }
447 
448     I opIndex()(size_t index)
449     {
450         index *= pack;
451         size_t start = index % bits;
452         index /= bits;
453         auto ret = (_field[index] >>> start) & mask;
454         static if (bits % pack)
455         {
456             sizediff_t end = start - (bits - pack);
457             if (end > 0)
458                 ret ^= cast(I)(_field[index + 1] << (bits - end)) >>> (bits - pack);
459         }
460         return cast(I) ret;
461     }
462 
463     I opIndexAssign()(I value, size_t index)
464     {
465         import std.traits: Unsigned;
466         assert(cast(Unsigned!I)value <= mask);
467         index *= pack;
468         size_t start = index % bits;
469         index /= bits;
470         _field[index] = cast(I)((_field[index] & ~(mask << start)) ^ (value << start));
471         static if (bits % pack)
472         {
473             sizediff_t end = start - (bits - pack);
474             if (end > 0)
475                 _field[index + 1] = cast(I)((_field[index + 1] & ~((I(1) << end) - 1)) ^ (value >>> (pack - end)));
476         }
477         return value;
478     }
479 
480     static if (hasZeroShiftFieldMember!Field)
481     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
482     auto assumeFieldsHaveZeroShift() @property
483     {
484         return BitpackField!(ZeroShiftField!Field, pack, I)(_field.assumeFieldsHaveZeroShift);
485     }
486 }
487 
488 ///
489 unittest
490 {
491     import mir.ndslice.iterator: FieldIterator;
492     ushort[10] data;
493     auto f = FieldIterator!(BitpackField!(ushort*, 6))(0, BitpackField!(ushort*, 6)(data.ptr));
494     f[0] = cast(ushort) 31;
495     f[1] = cast(ushort) 13;
496     f[2] = cast(ushort)  8;
497     f[3] = cast(ushort) 43;
498     f[4] = cast(ushort) 28;
499     f[5] = cast(ushort) 63;
500     f[6] = cast(ushort) 39;
501     f[7] = cast(ushort) 23;
502     f[8] = cast(ushort) 44;
503 
504     assert(f[0] == 31);
505     assert(f[1] == 13);
506     assert(f[2] ==  8);
507     assert(f[3] == 43);
508     assert(f[4] == 28);
509     assert(f[5] == 63);
510     assert(f[6] == 39);
511     assert(f[7] == 23);
512     assert(f[8] == 44);
513     assert(f[9] == 0);
514     assert(f[10] == 0);
515     assert(f[11] == 0);
516 }
517 
518 unittest
519 {
520     import mir.ndslice.slice;
521     import mir.ndslice.topology;
522     import mir.ndslice.sorting;
523     uint[2] data;
524     auto packed = data[].sliced.bitpack!18;
525     assert(packed.length == 3);
526     packed[0] = 5;
527     packed[1] = 3;
528     packed[2] = 2;
529     packed.sort;
530     assert(packed[0] == 2);
531     assert(packed[1] == 3);
532     assert(packed[2] == 5);
533 }
534 
535 ///
536 struct OrthogonalReduceField(FieldsIterator, alias fun, T)
537 {
538     import mir.ndslice.slice: Slice;
539 
540 @optmath:
541     /// non empty slice
542 
543     Slice!FieldsIterator _fields;
544 
545     ///
546     T _initialValue;
547 
548     ///
549     auto lightConst()() const @property
550     {
551         auto fields = _fields.lightConst;
552         return OrthogonalReduceField!(fields.Iterator, fun, T)(fields, _initialValue);
553     }
554 
555     ///
556     auto lightImmutable()() immutable @property
557     {
558         auto fields = _fields.lightImmutable;
559         return OrthogonalReduceField!(fields.Iterator, fun, T)(fields, _initialValue);
560     }
561 
562     /// `r = fun(r, fields[i][index]);` reduction by `i`
563     auto opIndex()(size_t index)
564     {
565         import std.traits: Unqual;
566         auto fields = _fields;
567         T r = _initialValue;
568         if (!fields.empty) do
569         {
570             r = cast(T) fun(r, fields.front[index]);
571             fields.popFront;
572         }
573         while(!fields.empty);
574         return r;
575     }
576 }
577 
578 ///
579 struct CycleField(Field)
580 {
581     import mir.ndslice.slice: Slice;
582 
583 @optmath:
584     /// Cycle length
585     size_t _length;
586     ///
587     Field _field;
588 
589     ///
590     auto lightConst()() const @property
591     {
592         auto field = .lightConst(_field);
593         return CycleField!(typeof(field))(_length, field);
594     }
595 
596     ///
597     auto lightImmutable()() immutable @property
598     {
599         auto field = .lightImmutable(_field);
600         return CycleField!(typeof(field))(_length, field);
601     }
602 
603     ///
604     auto ref opIndex()(size_t index)
605     {
606         return _field[index % _length];
607     }
608 
609     ///
610     static if (!__traits(compiles, &opIndex(size_t.init)))
611     {
612         auto ref opIndexAssign(T)(auto ref T value, size_t index)
613         {
614             return _field[index % _length] = value;
615         }
616     }
617 
618     static if (hasZeroShiftFieldMember!Field)
619     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
620     auto assumeFieldsHaveZeroShift() @property
621     {
622         return CycleField!(ZeroShiftField!Field)(_length, _field.assumeFieldsHaveZeroShift);
623     }
624 }
625 
626 ///
627 struct CycleField(Field, size_t length)
628 {
629     import mir.ndslice.slice: Slice;
630 
631 @optmath:
632     /// Cycle length
633     enum _length = length;
634     ///
635     Field _field;
636 
637     ///
638     auto lightConst()() const @property
639     {
640         auto field = .lightConst(_field);
641         return CycleField!(typeof(field), _length)(field);
642     }
643 
644     ///
645     auto lightImmutable()() immutable @property
646     {
647         auto field = .lightImmutable(_field);
648         return CycleField!(typeof(field), _length)(field);
649     }
650 
651     ///
652     auto ref opIndex()(size_t index)
653     {
654         return _field[index % _length];
655     }
656 
657     ///
658     static if (!__traits(compiles, &opIndex(size_t.init)))
659     {
660         auto ref opIndexAssign(T)(auto ref T value, size_t index)
661         {
662             return _field[index % _length] = value;
663         }
664     }
665 
666     static if (hasZeroShiftFieldMember!Field)
667     /// Defined if `Field` has member `assumeFieldsHaveZeroShift`.
668     auto assumeFieldsHaveZeroShift() @property
669     {
670         return CycleField!(ZeroShiftField!Field, _length)(_field.assumeFieldsHaveZeroShift);
671     }
672 }
673 
674 /++
675 `ndIotaField` is used by $(SUBREF topology, ndiota).
676 +/
677 struct ndIotaField(size_t N)
678     if (N)
679 {
680 @optmath:
681     ///
682     size_t[N - 1] _lengths;
683 
684     ///
685     auto lightConst()() const @property
686     {
687         return ndIotaField!N(_lengths);
688     }
689 
690     ///
691     auto lightImmutable()() const @property
692     {
693         return ndIotaField!N(_lengths);
694     }
695 
696     ///
697     size_t[N] opIndex()(size_t index) const
698     {
699         size_t[N] indices;
700         foreach_reverse (i; Iota!(N - 1))
701         {
702             indices[i + 1] = index % _lengths[i];
703             index /= _lengths[i];
704         }
705         indices[0] = index;
706         return indices;
707     }
708 }
709 
710 /++
711 `LinspaceField` is used by $(SUBREF topology, linspace).
712 +/
713 struct LinspaceField(T)
714 {
715     ///
716     size_t _length;
717 
718     ///
719     T _start = cast(T) 0, _stop = cast(T) 0;
720 
721     ///
722     auto lightConst()() scope const @property
723     {
724         return LinspaceField!T(_length, _start, _stop);
725     }
726 
727     ///
728     auto lightImmutable()() scope const @property
729     {
730         return LinspaceField!T(_length, _start, _stop);
731     }
732 
733     // no fastmath
734     ///
735     T opIndex()(sizediff_t index) scope const
736     {
737         sizediff_t d = _length - 1;
738         auto v = typeof(T.init.re)(d - index);
739         auto w = typeof(T.init.re)(index);
740         v /= d;
741         w /= d;
742         auto a = v * _start;
743         auto b = w * _stop;
744         return a + b;
745     }
746 
747 @optmath:
748 
749     ///
750     size_t length(size_t dimension = 0)() scope const @property
751         if (dimension == 0)
752     {
753         return _length;
754     }
755 
756     ///
757     size_t[1] shape()() scope const @property @nogc
758     {
759         return [_length];
760     }
761 }
762 
763 /++
764 Magic square field.
765 +/
766 struct MagicField
767 {
768 @optmath:
769 @safe pure nothrow @nogc:
770 
771     /++
772     Magic Square size.
773     +/
774     size_t _n;
775 
776 scope const:
777 
778     ///
779     MagicField lightConst()() @property
780     {
781         return this;
782     }
783 
784     ///
785     MagicField lightImmutable()() @property
786     {
787         return this;
788     }
789 
790     ///
791     size_t length(size_t dimension = 0)() @property
792         if(dimension <= 2)
793     {
794         return _n * _n;
795     }
796 
797     ///
798     size_t[1] shape() @property
799     {
800         return [_n * _n];
801     }
802 
803     ///
804     size_t opIndex(size_t index)
805     {
806         pragma(inline, false);
807         auto d = index / _n;
808         auto m = index % _n;
809         if (_n & 1)
810         {
811             //d = _n - 1 - d;     // MATLAB synchronization
812             //index = d * _n + m; // ditto
813             auto r = (index + 1 - d + (_n - 3) / 2) % _n;
814             auto c = (_n * _n - index + 2 * d) % _n;
815             return r * _n + c + 1;
816         }
817         else
818         if ((_n & 2) == 0)
819         {
820             auto a = (d + 1) & 2;
821             auto b = (m + 1) & 2;
822             return a != b ? index + 1: _n * _n - index;
823         }
824         else
825         {
826             auto n = _n / 2 ;
827             size_t shift;
828             ptrdiff_t q;
829             ptrdiff_t p = m - n;
830             if (p >= 0)
831             {
832                 m = p;
833                 shift = n * n;
834                 auto mul = m <= n / 2 + 1;
835                 q = d - n;
836                 if (q >= 0)
837                 {
838                     d = q;
839                     mul = !mul;
840                 }
841                 if (mul)
842                 {
843                     shift *= 2;
844                 }
845             }
846             else
847             {
848                 auto mul = m < n / 2;
849                 q = d - n;
850                 if (q >= 0)
851                 {
852                     d = q;
853                     mul = !mul;
854                 }
855                 if (d == n / 2 && (m == 0 || m == n / 2))
856                 {
857                     mul = !mul;
858                 }
859                 if (mul)
860                 {
861                     shift = n * n * 3; 
862                 }
863             }
864             index = d * n + m;
865             auto r = (index + 1 - d + (n - 3) / 2) % n;
866             auto c = (n * n - index + 2 * d) % n;
867             return r * n + c + 1 + shift;
868         }
869     }
870 }
871 
872 /++
873 `SparseField` is used to represent Sparse ndarrays in mutable DOK format.
874 +/
875 struct SparseField(T)
876 {
877     ///
878     T[size_t] _table;
879 
880     ///
881     auto lightConst()() const @trusted
882     {
883         return SparseField!(const T)(cast(const(T)[size_t])_table);
884     }
885 
886     ///
887     auto lightImmutable()() immutable @trusted
888     {
889         return SparseField!(immutable T)(cast(immutable(T)[size_t])_table);
890     }
891 
892     ///
893     T opIndex()(size_t index)
894     {
895         import std.traits: isScalarType;
896         static if (isScalarType!T)
897             return _table.get(index, cast(T)0);
898         else
899             return _table.get(index, null);
900     }
901 
902     ///
903     T opIndexAssign()(T value, size_t index)
904     {
905         import std.traits: isScalarType;
906         static if (isScalarType!T)
907         {
908             if (value != 0)
909                 _table[index] = value;
910             else
911                 _table.remove(index);
912         }
913         else
914         {
915             if (value !is null)
916                 _table[index] = value;
917             else
918                 _table.remove(index);
919         }
920         return value;
921     }
922 
923     ///
924     T opIndexUnary(string op)(size_t index)
925         if (op == `++` || op == `--`)
926     {
927         import std.traits: isScalarType;
928         mixin (`auto value = ` ~ op ~ `_table[index];`);
929         static if (isScalarType!T)
930         {
931             if (value == 0)
932                 _table.remove(index);
933         }
934         else
935         {
936             if (value is null)
937                 _table.remove(index);
938         }
939         return value;
940     }
941 
942     ///
943     T opIndexOpAssign(string op)(T value, size_t index)
944         if (op == `+` || op == `-`)
945     {
946         import std.traits: isScalarType;
947         mixin (`value = _table[index] ` ~ op ~ `= value;`); // this works
948         static if (isScalarType!T)
949         {
950             if (value == 0)
951                 _table.remove(index);
952         }
953         else
954         {
955             if (value is null)
956                 _table.remove(index);
957         }
958         return value;
959     }
960 }