Removed comptime for str UnitParser

I am not quite sure about it yet, but it is faster sooo idk.
Let's see long term

.gitignore

@@ -1,2 +1,3 @@
 zig-out
 .zig-cache
+mkdocs.yaml

LICENSE

@@ -0,0 +1,674 @@
+GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<https://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<https://www.gnu.org/licenses/why-not-lgpl.html>.

README.md

@@ -1,268 +1,253 @@
-# zig_units
+# dimal — Dimensional Analysis for Zig
 
-`zig_units` lets you attach physical units to numeric values so that dimension mismatches (like adding distance to time) become **compile errors** rather than silent bugs. 
+A dimensional analysis library for Zig with a unified `Tensor` API for scalars, vectors, matrices, and higher-dimensional data. All dimension and unit tracking happens at compile time—zero runtime overhead—and all operations use SIMD intrinsics.
 
-At runtime, a `Quantity` is just its underlying numeric value — **zero memory overhead.**
+If you try to add meters to seconds, it won't compile. That's the point.
 
-```zig
+> **Source:** [git.bouvais.lu/adrien/zig-dimal](https://git.bouvais.lu/adrien/zig-dimal)  
-const velocity = distance.divBy(time);  // Result type: L¹T⁻¹  ✓
+> **Minimum Zig version:** `0.16.0`
-const error    = mass.add(velocity);    // COMPILE ERROR: M¹ != L¹T⁻¹
-```
 
-**Requirements:** Zig `0.16.0`
+---
+
+## Background
+
+Started because I needed `i128` positions for a space simulation to avoid floating-point precision loss far from the origin. Grew into a type system for tracking physical dimensions at compile time. It's been useful enough to share.
+
+- **Compile-time dimension checking** — catch unit mismatches before runtime.
+- **Unified `Tensor` API** — same interface for scalars, vectors, matrices, and higher-rank tensors.
+- **SIMD operations** — vector and matrix code automatically uses SIMD instructions.
+- **Zero runtime cost** — all dimension and scale tracking is erased at compile time.
+- **Supports `i128`** — useful for high-precision fixed-point integer math.
+
+---
+
+## Features
+
+- **Compile-time dimension checking** — all physical-unit tracking happens at compile time.
+- **Automatic unit conversion** — use `.to()` to convert between compatible units (e.g. `km/h` → `m/s`). Scale factors are resolved at comptime.
+- **Unified `Tensor` API** — one type for scalars `{1}`, vectors `{N}`, matrices `{M, N}`, and higher-rank tensors.
+- **SIMD operations** — vector and matrix code compiles to SIMD instructions automatically.
+- **Tensor contraction** — `.contract(other, axis_a, axis_b)` for dot products, matrix multiplication, and general tensor contractions.
+- **Full SI prefix support** — `pico` through `peta`, plus Imperial units and time scales.
+- **Physical constants** — Planck, Boltzmann, speed of light, gravitational constant, etc.
+- **Pre-built quantities** — `Velocity`, `Acceleration`, `Force`, `Energy`, `Pressure`, `Charge`, and more.
+- **Basic vector operations** — cross product, length/magnitude, element-wise arithmetic.
+- **Formatting** — values print with units: `9.81m.s⁻²`, `0.172km`.
+
+### Current Limitations
+
+- GPU support not implemented.
+- Performance on small tensors is limited by Zig's vector width.
+
+---
+
+## The 7 SI Base Dimensions
+
+| Symbol | Dimension            | SI Unit |
+|--------|----------------------|---------|
+| `L`    | Length               | `m`     |
+| `M`    | Mass                 | `g`     |
+| `T`    | Time                 | `s`     |
+| `I`    | Electric Current     | `A`     |
+| `Tr`   | Temperature          | `K`     |
+| `N`    | Amount of Substance  | `mol`   |
+| `J`    | Luminous Intensity   | `cd`    |
 
 ---
 
 ## Installation
 
-### 1. Add as a Zig dependency
+### 1. Add the dependency (Zig 0.14+)
-```bash
+
-zig fetch --save https://github.com/YOUR_USERNAME/zig_units/archive/refs/heads/main.tar.gz
+```sh
+zig fetch --save git+https://git.bouvais.lu/adrien/zig-dimal#0.2.0
 ```
 
-### 2. Configure `build.zig`
+### 2. Wire it up in `build.zig`
-```zig
-const zig_units = b.dependency("zig_units", .{
-    .target = target,
-    .optimize = optimize,
-});
-// Add to your module or executable
-exe.root_module.addImport("units", zig_units.module("zig_units"));
-```
-
----
-
-## Quick Start: Using Predefined Quantities
-
-`units.Base` provides a clean way to instantiate common physical types without manually defining dimensions.
 
 ```zig
 const std = @import("std");
-const units = @import("units");
 
-pub fn main() !void {
+pub fn build(b: *std.Build) void {
-    // Instantiate types for f32 backing
+    const target = b.standardTargetOptions(.{});
-    const Meter  = units.Base.Meter.Of(f32);
+    const optimize = b.standardOptimizeOption(.{});
-    const Second = units.Base.Second.Of(f32);
     
-    const dist = Meter{ .value = 10.0 };
+    const dimal = b.dependency("dimal", .{
-    const time = Second{ .value = 2.0 };
+        .target = target,
+        .optimize = optimize,
+    }).module("dimal");
 
-    // Arithmetic is type-safe and creates the correct resulting dimension
+    const exe = b.addExecutable(.{
-    const vel = dist.divBy(time); // Type is Velocity (L/T)
+        .name = "my_app",
-    
+        .root_source_file = b.path("src/main.zig"),
-    std.debug.print("Speed: {f}\n", .{vel}); // Output: 5m.s⁻¹
+        .target = target,
+        .optimize = optimize,
+    });
+    exe.root_module.addImport("dimal", dimal);
+    b.installArtifact(exe);
 }
 ```
 
----
+### 3. Import and use
-
-### Core Arithmetic Operations
-
-Dimensional analysis is handled entirely at compile-time. If the math doesn't make physical sense, it won't compile.
 
 ```zig
-const M = units.Base.Meter.Of(f32);
+const dma = @import("dimal");
-const S = units.Base.Second.Of(f32);
+const Tensor = dma.Tensor;
-
+const Base = dma.Base;
-const dist = M{ .value = 100.0 };
-const time = S{ .value = 5.0 };
-
-// 1. Addition & Subtraction (Must have same dimensions)
-const two_dist = dist.add(dist); // 200.0m
-const zero     = dist.sub(dist); // 0.0m
-
-// 2. Division (Subtracts dimension exponents)
-// Result: Velocity (L¹ T⁻¹)
-const vel = dist.divBy(time); // 20.0 m.s⁻¹
-
-// 3. Multiplication (Adds dimension exponents)
-// Result: Area (L²)
-const area = dist.mulBy(dist); // 100.0 m²
-
-// 4. Chained Operations
-// Result: Acceleration (L¹ T⁻²)
-const accel = dist.divBy(time).divBy(time); // 4.0 m.s⁻²
 ```
 
 ---
 
-## Defining Custom Quantities
+## Quick Example: Lunar Descent
 
-You aren't limited to the built-in library. You can define any physical quantity by specifying its **Dimensions**
+Simulate a spacecraft descending to the Moon with correct physics and type safety:
-(powers of base units) and its **Scale** (SI prefixes).
-
-### 1. Create a custom dimension
-Dimensions are defined by 7 base SI units: `L` (Length), `M` (Mass), `T` (Time), `I` (Current), `Tp` (Temp), `N` (Substance), `J` (Intensity).
 
 ```zig
-const Dims = units.Dimensions;
+const std = @import("std");
-const Scales = units.Scales;
+const dma = @import("dimal");
+const Tensor = dma.Tensor;
 
-// Frequency is T⁻¹
+pub fn main() void {
-const FreqDims = Dims.init(.{ .T = -1 });
+    // Define types: m/s² acceleration, m/s velocity, m distance
+    const Acceleration = dma.Base.Acceleration.Of(f64);
+    const Velocity = dma.Base.Velocity.Of(f64);
+    const Distance = dma.Base.Meter.Of(f64);
+    const Time = dma.Base.Second.Of(f64);
 
-// Force is M¹ L¹ T⁻²
+    // Initial conditions
-const ForceDims = Dims.init(.{ .M = 1, .L = 1, .T = -2 });
+    const g_moon: Acceleration = .{ .data = @splat(1.62) };
+    const v_initial: Velocity = .{ .data = @splat(100.0) };
+    const h_initial: Distance = .{ .data = @splat(10000.0) };
+    const dt: Time = .{ .data = @splat(1.0) };
+
+    var h = h_initial;
+    var v = v_initial;
+    var t: f64 = 0;
+
+    // Simulate descent
+    while (h.data[0] > 0 and t < 1000) : (t += 1.0) {
+        // a = -g (gravity pulls down)
+        const a = g_moon.mul(-1.0);
+
+        // Update: v = v₀ + at
+        v = v.add(a.mul(dt));
+
+        // Update: h = h₀ + vt
+        h = h.add(v.mul(dt));
+
+        if (@mod(t, 100.0) == 0) {
+            std.debug.print("t={d:.0}s | h={d:.1} | v={d:.1}\n", .{
+                t,
+                h,
+                v,
+            });
+        }
+    }
+
+    std.debug.print("Landed in {d:.1}s at h={d:.1}\n", .{ t, h });
+}
 ```
 
-### 2. Create a custom Type
+**Output:**
-Combine a numeric type, the dimensions, and a scale.
+```
-
+t=0s | h=10000m | v=100m.s⁻¹
-```zig
+t=100s | h=8019m | v=-61.8m.s⁻¹
-const Hertz = units.Quantity(f32, FreqDims, Scales.init(.{}));
+t=200s | h=4174.4m | v=-223.6m.s⁻¹
-
+...
-// A specialized scale: Millimeters per Second Squared
+Landed in 323.5s at h=-0.01m
-const MmPerSecSq = units.Quantity(f32, 
-    Dims.init(.{ .L = 1, .T = -2 }), 
-    Scales.init(.{ .L = .m }) // .m = milli
-);
 ```
 
 ---
 
-## Unit Conversions
+## API Overview
 
-The library handles SI prefixes (`k`, `m`, `u`, `n`, etc.) and time aliases (`.min`, `.hour`) automatically.
+### Tensors
-When performing arithmetic between different scales, the **finer (smaller) scale wins** to preserve precision.
+
+A **`Tensor`** is parameterized by:
+- **`T`** — numeric type: `f32`, `f64`, `i128`, etc.
+- **`dims`** — physical dimensions (struct literal): `.{.L = 1, .T = -1}` means length/time (velocity).
+- **`scales`** — SI prefixes or custom scales: `.{.L = .k, .T = .hour}` means km/h.
+- **`shape`** — array shape: `&.{1}` is a scalar, `&.{3}` is a 3-vector, `&.{3, 3}` is a 3×3 matrix.
 
 ```zig
-const KM = units.Base.Meter.Scaled(f32, Scales.init(.{ .L = .k })); // Kilometers
+// Scalar: 1-element tensor
-const M  = units.Base.Meter.Of(f32);                               // Meters
+const Meter = Tensor(f64, .{.L = 1}, .{}, &.{1});
+const m = Meter{ .data = @splat(5.0) };
 
-const d1 = KM{ .value = 1.2 };  // 1.2 km
+// Vector: N-element tensor (SIMD)
-const d2 = M{ .value = 300.0 }; // 300 m
+const Vec3Meter = Tensor(f64, .{.L = 1}, .{}, &.{3});
+const v = Vec3Meter{ .data = @shuffle(f64, [_]f64{1, 2, 3}, [_]f64 undefined, [_]i32{0, 1, 2, 0, 0, 0}) };
 
-const total = d1.add(d2);      // Result is 1500.0 (Meters)
+// Matrix: M×N tensor (SIMD-accelerated)
-const final = total.to(KM);    // Explicitly convert back to KM -> 1.5
+const Mat3x3Velocity = Tensor(f32, .{.L = 1, .T = -1}, .{}, &.{3, 3});
+const m_vel = Mat3x3Velocity{ .data = @splat(10.0) };
+
+// Higher-rank tensor
+const Rank4 = Tensor(f64, .{.M = 1}, .{}, &.{2, 3, 4, 5});
+```
+
+### Common Operations
+
+| Operation | Description |
+|-----------|-------------|
+| `.add(rhs)` | Element-wise addition. Auto-converts scales. |
+| `.sub(rhs)` | Element-wise subtraction. |
+| `.mul(rhs)` | Multiply; dimensions are summed. `rhs` can be a tensor or bare number. |
+| `.div(rhs)` | Divide; dimensions are subtracted. |
+| `.contract(other, axis_a, axis_b)` | Tensor contraction: dot product, matrix multiply, or general N-D contraction. |
+| `.cross(rhs)` | Cross product (3-vectors only). Returns a 3-vector. |
+| `.length()` / `.lengthSqr()` | Euclidean length (or squared length) of a vector. Returns a scalar `T`. |
+| `.product()` | Multiply all elements. Returns a scalar with combined dimensions. |
+| `.abs()` | Element-wise absolute value. Dimensions unchanged. |
+| `.pow(exp)` | Raise to comptime exponent. Dimension exponents multiplied by `exp`. |
+| `.sqrt()` | Element-wise square root. Compile error if any dimension exponent is odd. |
+| `.to(DestType)` | Convert to another unit of the same dimension. Comptime error on mismatch. |
+| `.eq(rhs)` / `.ne(rhs)` | Element-wise equality/inequality. |
+| `.gt(rhs)` / `.gte(rhs)` | Greater-than comparisons. |
+| `.lt(rhs)` / `.lte(rhs)` | Less-than comparisons. |
+
+### Pre-built Types (via `dma.Base`)
+
+Use `.Of(T)` for base units, `.Scaled(T, scales)` for custom scales:
+
+```zig
+const Velocity = dma.Base.Velocity.Of(f64);
+const Kmh = dma.Base.Velocity.Scaled(f64, .{.L = .k, .T = .hour});
+const Force = dma.Base.Force.Of(f32);
+const Energy = dma.Base.Energy.Of(f64);
+```
+
+Also available: `Acceleration`, `Inertia`, `Pressure`, `Power`, `Area`, `Volume`, `Density`, `Frequency`, `Viscosity`, `Charge`, `Potential`, `Resistance`, `MagneticFlux`, `ThermalCapacity`, `ThermalConductivity`, and many more.
+
+---
+
+## SIMD Performance
+
+Operations on vectors and matrices use Zig's `@Vector` intrinsics, which compile to SIMD instructions on most platforms. This makes vector operations faster than equivalent scalar loops, but don't expect miracles—SIMD is still limited by memory bandwidth and CPU cache.
+
+Run the included benchmarks to see what you get on your hardware:
+```sh
+zig build benchmark
 ```
 
 ---
 
-## Physical Vectors (Vec3)
+## Next Steps
 
-Physical quantities often come in 3D vectors (Position, Velocity, Force). Every `Quantity` type has a `.Vec3` alias built-in.
+- **GPU support** — eventually, for large tensor operations. WebGPU is a target.
-
+- **Toy physics language** — I've been sketching ideas for a language optimized for numerical physics (tentatively called Éclat). It would use dimal as the foundation. No timeline yet; this is a long-term experiment.
-```zig
-const Vec3M = units.Base.Meter.Of(f32).Vec3;
-
-const gravity = Vec3M{ .data = .{ 0, -9.81, 0 } };
-const pos     = Vec3M.initDefault(0); // [0, 0, 0]
-
-// Vectors support standard operations
-const length = gravity.length(); // Returns f32: 9.81
-const double = gravity.scale(2.0);
-```
-
-You can also create a Vector of any length.
-Vec3 found in a Quantity is just a convenience.
-
-```zig
-const M  = units.Base.Meter.Of(f32);
-const Vec10M = units.QuantityVec(10, Meter);
-
-const gravity = Vec10M.initDefault(1);
-const length = gravity.length(); // Returns f32: 1.0
-```
-
----
-
-## High Precision & Integer Backing
-
-While most libraries default to `f32` or `f64`, `zig_units` is mainly designed to support **large-bit integers (`i128`, `i256`)**. 
-
-This is critical for applications like **space simulations**, where floating-point numbers suffer from "jitter" or "flickering"
-once you travel far from the origin. By using an `i128` with a millimeter scale, you can represent the diameter
-of the observable universe with millimeter precision—something impossible with `f64`.
-
-### Avoiding Floating-Point Jitter
-```zig
-// Millimeter precision using 128-bit integers
-const MM  = units.Base.Meter.Scaled(i128, units.Scales.init(.{ .L = .m }));
-const KM  = units.Base.Meter.Scaled(i128, units.Scales.init(.{ .L = .k }));
-
-const solar_system_dist = KM{ .value = 150_000_000 }; // 150 million km
-const ship_nudge        = MM{ .value = 5 };           // 5 mm
-
-// The library performs exact integer math for conversions.
-// Resulting type is MM (the finer scale), maintaining perfect precision.
-const new_pos = solar_system_dist.add(ship_nudge); 
-```
-
-### Integer-Specific Features
- *   **Exact Conversions:** When converting between integer scales (e.g., `km` to `m`), the library uses fast-path native multiplication.
- *   **Safe Vector Lengths:** `QuantityVec.length()` includes a custom integer square root implementation, allowing you to calculate distances between coordinates without ever casting to a float.
- *   **Zero Drift:** Unlike floats, repeated additions and subtractions of integers never accumulate "epsilon" drift, ensuring your simulation remains deterministic.
- *   **Precision-First Scaling:** When operating on two different scales (e.g., adding `km` and `mm`), the result automatically adopts the finer scale (`mm`). This ensures **zero implicit data loss** during calculation. You only lose precision if you *explicitly* choose to convert back to a coarser scale using `.to()`.
-
----
-
-## SI Scales Reference
-
-| Prefix | Enum | Factor |
-| :--- | :--- | :--- |
-| **Kilo** | `.k` | 10³ |
-| **Mega** | `.M` | 10⁶ |
-| **Giga** | `.G` | 10⁹ |
-| **Milli** | `.m` | 10⁻³ |
-| **Micro** | `.u` | 10⁻⁶ |
-| **Minute**| `.min` | 60 |
-| **Hour**  | `.hour`| 3,600 |
-
----
-
-## API Summary
-
-### `Quantity(T, dims, scales)`
-- `.add(rhs)` / `.sub(rhs)`: Automatic scaling, requires same dimensions.
-- `.mulBy(rhs)` / `.divBy(rhs)`: Composes dimensions (e.g., $L \times L = L^2$).
-- `.scale(scalar)`: Multiply by a raw number (preserves dimensions).
-- `.to(OtherType)`: Safely convert between scales of the same dimension.
-- `.vec3()`: Create a 3D vector from a scalar.
-
-### `Dimensions`
-- `L`: Length (m)
-- `M`: Mass (g)
-- `T`: Time (s)
-- `I`: Current (A)
-- `Tp`: Temperature (K)
-- `N`: Amount (mol)
-- `J`: Intensity (cd)
 
 ---
 
 ## Testing & Benchmarks
 
-`zig_units` comes with a comprehensive test suite that verifies dimensional correctness, SI prefix scaling,
+```sh
-and vector math accuracy across all numeric types.
+zig build test       # Run all unit tests
-
+zig build benchmark  # Run performance benchmarks
-### Running Tests
-To run the full suite of unit tests and performance benchmarks:
-```bash
-zig build test
 ```
 
-### Benchmarks
+---
-When you run the tests, the library also executes a performance benchmark.
-This measures the cost of operations (in nanoseconds per operation) across different backing types (`i32` to `f128`) and vector lengths.
 
-Because all dimensional logic is resolved at **compile-time**, you will see that `Quantity` operations perform at the same speed as raw primitive math.
+## License
 
-**Example Benchmark Output:**
+See the repository for license details.
-```text
- Quantity<T> benchmark — 100,000 iterations, 10 samples/cell
-
-┌───────────────────┬──────┬─────────────────────┬─────────────────────┐
-│ Operation         │ Type │ ns / op (± delta)   │ Throughput (ops/s)  │
-├───────────────────┼──────┼─────────────────────┼─────────────────────┤
-│ add               │ i32  │     0.18 ns ±0.02   │          5555555556 │
-│ mulBy             │ f64  │     0.22 ns ±0.01   │          4545454545 │
-└───────────────────┴──────┴─────────────────────┴─────────────────────┘
-```
-
-### Verification Examples
-The test suite ensures that:
-- **Dimension Safety:** Chained operations like `dist.divBy(time).divBy(time)` correctly result in an Acceleration type.
-- **Scale Accuracy:** Adding `1km + 1mm` results exactly in `1000001mm` without truncation.
-- **Formatting:** Quantities print correctly with Unicode superscripts (e.g., `9.81m.s⁻²`).
-- **Vector Math:** Euclidean lengths for both floats and integers are verified against known constants.

build.zig

@@ -2,40 +2,38 @@ const std = @import("std");
 
 pub fn build(b: *std.Build) void {
     const target = b.standardTargetOptions(.{});
-    const optimize = b.standardOptimizeOption(.{ .preferred_optimize_mode = .ReleaseSmall });
+    const optimize = b.standardOptimizeOption(.{});
-    const exe = b.addExecutable(.{
+
-        .name = "Zig_Units",
+    // 1. Define the module so other projects can import it
+    _ = b.addModule("dimal", .{
+        .root_source_file = b.path("src/main.zig"),
+    });
+
+    const exe_tests = b.addTest(.{
         .root_module = b.createModule(.{
             .root_source_file = b.path("src/main.zig"),
             .target = target,
             .optimize = optimize,
+        }),
+        .test_runner = .{ .path = b.path("test_runner.zig"), .mode = .simple },
+    });
+
+    const run_exe_tests = b.addRunArtifact(exe_tests);
+    const test_step = b.step("test", "Run tests");
+    test_step.dependOn(&run_exe_tests.step);
+
+    const bench_exe = b.addExecutable(.{
+        .name = "benchmark",
+        .root_module = b.createModule(.{
+            .root_source_file = b.path("src/benchmark.zig"),
+            .target = target,
             .imports = &.{},
         }),
     });
 
-    b.installArtifact(exe);
+    b.installArtifact(bench_exe);
-    const run_step = b.step("run", "Run the app");
 
-    const run_cmd = b.addRunArtifact(exe);
+    const run_bench = b.addRunArtifact(bench_exe);
-    run_step.dependOn(&run_cmd.step);
+    const bench_step = b.step("benchmark", "Run the benchmark suite");
-
+    bench_step.dependOn(&run_bench.step);
-    run_cmd.step.dependOn(b.getInstallStep());
-
-    if (b.args) |args| {
-        run_cmd.addArgs(args);
-    }
-
-    const exe_tests = b.addTest(.{
-        .root_module = exe.root_module,
-        .test_runner = .{
-            .path = .{ .cwd_relative = "test_runner.zig" },
-            .mode = .simple,
-        },
-    });
-
-    // A run step that will run the second test executable.
-    const run_exe_tests = b.addRunArtifact(exe_tests);
-
-    const test_step = b.step("test", "Run tests");
-    test_step.dependOn(&run_exe_tests.step);
 }

build.zig.zon

@@ -1,81 +1,13 @@
 .{
-    // This is the default name used by packages depending on this one. For
+    .name = .dimal,
-    // example, when a user runs `zig fetch --save <url>`, this field is used
+    .version = "0.2.1",
-    // as the key in the `dependencies` table. Although the user can choose a
+    .fingerprint = 0x9453b1ff1e52d858,
-    // different name, most users will stick with this provided value.
-    //
-    // It is redundant to include "zig" in this name because it is already
-    // within the Zig package namespace.
-    .name = .Zig_Units,
-    // This is a [Semantic Version](https://semver.org/).
-    // In a future version of Zig it will be used for package deduplication.
-    .version = "0.0.0",
-    // Together with name, this represents a globally unique package
-    // identifier. This field is generated by the Zig toolchain when the
-    // package is first created, and then *never changes*. This allows
-    // unambiguous detection of one package being an updated version of
-    // another.
-    //
-    // When forking a Zig project, this id should be regenerated (delete the
-    // field and run `zig build`) if the upstream project is still maintained.
-    // Otherwise, the fork is *hostile*, attempting to take control over the
-    // original project's identity. Thus it is recommended to leave the comment
-    // on the following line intact, so that it shows up in code reviews that
-    // modify the field.
-    .fingerprint = 0x934b00cf88f69b22, // Changing this has security and trust implications.
-    // Tracks the earliest Zig version that the package considers to be a
-    // supported use case.
     .minimum_zig_version = "0.16.0",
-    // This field is optional.
+    .dependencies = .{},
-    // Each dependency must either provide a `url` and `hash`, or a `path`.
-    // `zig build --fetch` can be used to fetch all dependencies of a package, recursively.
-    // Once all dependencies are fetched, `zig build` no longer requires
-    // internet connectivity.
-    .dependencies = .{
-        // See `zig fetch --save <url>` for a command-line interface for adding dependencies.
-        //.example = .{
-        //    // When updating this field to a new URL, be sure to delete the corresponding
-        //    // `hash`, otherwise you are communicating that you expect to find the old hash at
-        //    // the new URL. If the contents of a URL change this will result in a hash mismatch
-        //    // which will prevent zig from using it.
-        //    .url = "https://example.com/foo.tar.gz",
-        //
-        //    // This is computed from the file contents of the directory of files that is
-        //    // obtained after fetching `url` and applying the inclusion rules given by
-        //    // `paths`.
-        //    //
-        //    // This field is the source of truth; packages do not come from a `url`; they
-        //    // come from a `hash`. `url` is just one of many possible mirrors for how to
-        //    // obtain a package matching this `hash`.
-        //    //
-        //    // Uses the [multihash](https://multiformats.io/multihash/) format.
-        //    .hash = "...",
-        //
-        //    // When this is provided, the package is found in a directory relative to the
-        //    // build root. In this case the package's hash is irrelevant and therefore not
-        //    // computed. This field and `url` are mutually exclusive.
-        //    .path = "foo",
-        //
-        //    // When this is set to `true`, a package is declared to be lazily
-        //    // fetched. This makes the dependency only get fetched if it is
-        //    // actually used.
-        //    .lazy = false,
-        //},
-    },
-    // Specifies the set of files and directories that are included in this package.
-    // Only files and directories listed here are included in the `hash` that
-    // is computed for this package. Only files listed here will remain on disk
-    // when using the zig package manager. As a rule of thumb, one should list
-    // files required for compilation plus any license(s).
-    // Paths are relative to the build root. Use the empty string (`""`) to refer to
-    // the build root itself.
-    // A directory listed here means that all files within, recursively, are included.
     .paths = .{
         "build.zig",
         "build.zig.zon",
         "src",
-        // For example...
+        "test_runner.zig",
-        //"LICENSE",
-        //"README.md",
     },
 }

src/Base.zig

@@ -0,0 +1,260 @@
+const std = @import("std");
+
+// Adjust these imports to match your actual file names
+const Dimensions = @import("Dimensions.zig");
+const Scales = @import("Scales.zig");
+const Tensor = @import("Tensor.zig").TensorStatic;
+
+fn PhysicalConstant(comptime d: Dimensions.ArgOpts, comptime val: f64, comptime s: Scales.ArgOpts) type {
+    return struct {
+        pub const dims = Dimensions.init(d);
+        pub const scales = Scales.init(s);
+
+        /// Instantiates the constant into a specific numeric type.
+        pub fn Of(comptime T: type) Tensor(T, d, s, &.{1}) {
+            const casted_val: T = switch (@typeInfo(T)) {
+                .float => @floatCast(val),
+                .int => @intFromFloat(val),
+                else => @compileError("Unsupported type for PhysicalConstant"),
+            };
+            return Tensor(T, d, s, &.{1}).splat(casted_val);
+        }
+    };
+}
+
+fn BaseScalar(comptime d: Dimensions.ArgOpts) type {
+    return struct {
+        pub const dims = Dimensions.init(d);
+
+        /// Creates a Scalar of this dimension using default scales.
+        /// Example: const V = Quantities.Velocity.Of(f32);
+        pub fn Of(comptime T: type) type {
+            return Tensor(T, d, .{}, &.{1});
+        }
+
+        /// Creates a Scalar of this dimension using custom scales.
+        /// Example: const Kmh = Quantities.Velocity.Scaled(f32, .{ .L = .k, .T = .hour });
+        pub fn Scaled(comptime T: type, comptime s: Scales.ArgOpts) type {
+            return Tensor(T, d, s, &.{1});
+        }
+    };
+}
+
+// ==========================================
+// Physical Constants
+// ==========================================
+
+pub const Constants = struct {
+    /// Speed of light in vacuum (c) [m/s]
+    pub const SpeedOfLight = PhysicalConstant(.{ .L = 1, .T = -1 }, 299792458.0, .{});
+
+    /// Planck constant (h) [J⋅s = kg⋅m²⋅s⁻¹]
+    pub const Planck = PhysicalConstant(.{ .M = 1, .L = 2, .T = -1 }, 6.62607015e-34, .{ .M = .k });
+
+    /// Reduced Planck constant (ℏ) [J⋅s]
+    pub const ReducedPlanck = PhysicalConstant(.{ .M = 1, .L = 2, .T = -1 }, 1.054571817e-34, .{ .M = .k });
+
+    /// Boltzmann constant (k_B) [J⋅K⁻¹ = kg⋅m²⋅s⁻²⋅K⁻¹]
+    pub const Boltzmann = PhysicalConstant(.{ .M = 1, .L = 2, .T = -2, .Tp = -1 }, 1.380649e-23, .{ .M = .k });
+
+    /// Newtonian constant of gravitation (G) [m³⋅kg⁻¹⋅s⁻²]
+    pub const Gravitational = PhysicalConstant(.{ .M = -1, .L = 3, .T = -2 }, 6.67430e-11, .{ .M = .k });
+
+    /// Stefan-Boltzmann constant () [W⋅m⁻²⋅K⁻⁴ = kg⋅s⁻³⋅K⁻⁴]
+    pub const StefanBoltzmann = PhysicalConstant(.{ .M = 1, .T = -3, .Tp = -4 }, 5.670374419e-8, .{ .M = .k });
+
+    /// Elementary charge (e) [C = A⋅s]
+    pub const ElementaryCharge = PhysicalConstant(.{ .T = 1, .I = 1 }, 1.602176634e-19, .{});
+
+    /// Vacuum magnetic permeability (μ_0) [N⋅A⁻² = kg⋅m⋅s⁻²⋅A⁻²]
+    pub const VacuumPermeability = PhysicalConstant(.{ .M = 1, .L = 1, .T = -2, .I = -2 }, 1.25663706127e-6, .{ .M = .k });
+
+    /// Vacuum electric permittivity (ε_0) [F⋅m⁻¹ = A²⋅s⁴⋅kg⁻¹⋅m⁻³]
+    pub const VacuumPermittivity = PhysicalConstant(.{ .M = -1, .L = -3, .T = 4, .I = 2 }, 8.8541878188e-12, .{ .M = .k });
+
+    /// Electron mass (m_e) [kg]
+    pub const ElectronMass = PhysicalConstant(.{ .M = 1 }, 9.1093837139e-31, .{ .M = .k });
+
+    /// Proton mass (m_p) [kg]
+    pub const ProtonMass = PhysicalConstant(.{ .M = 1 }, 1.67262192595e-27, .{ .M = .k });
+
+    /// Neutron mass (m_n) [kg]
+    pub const NeutronMass = PhysicalConstant(.{ .M = 1 }, 1.67492750056e-27, .{ .M = .k });
+
+    /// Fine-structure constant () [Dimensionless]
+    pub const FineStructure = PhysicalConstant(.{}, 0.0072973525643, .{});
+
+    /// Avogadro constant (N_A) [mol⁻¹]
+    /// Note: Assuming mol is currently treated as dimensionless in the base system,
+    /// otherwise requires adding an `.N` dimension to Dimensions.ArgOpts.
+    pub const Avogadro = PhysicalConstant(.{}, 6.02214076e23, .{});
+};
+
+// ==========================================
+// Base Quantities
+// ==========================================
+pub const Dimless = BaseScalar(.{});
+pub const Meter = BaseScalar(.{ .L = 1 });
+pub const Second = BaseScalar(.{ .T = 1 });
+pub const Gramm = BaseScalar(.{ .M = 1 });
+pub const Kelvin = BaseScalar(.{ .Tr = 1 });
+pub const ElectricCurrent = BaseScalar(.{ .I = 1 });
+
+// ==========================================
+// Electric
+// ==========================================
+pub const ElectricConductivity = BaseScalar(.{ .M = -1, .L = -3, .T = 3, .I = 2 });
+pub const ElectricCharge = BaseScalar(.{ .T = 1, .I = 1 });
+pub const ElectricPotential = BaseScalar(.{ .T = -3, .L = 2, .M = 1, .I = -1 });
+pub const ElectricResistance = BaseScalar(.{ .M = 1, .L = 2, .T = -3, .I = -2 });
+pub const ElectricResistivity = BaseScalar(.{ .M = 1, .L = 3, .T = -3, .I = -2 });
+pub const ElectricCapacitance = BaseScalar(.{ .T = 4, .L = -2, .M = -1, .I = 2 });
+pub const ElectricImpedance = ElectricResistance;
+pub const MagneticFlux = BaseScalar(.{ .M = 1, .L = 2, .T = -2, .I = -1 });
+pub const MagneticDensity = BaseScalar(.{ .M = 1, .T = -2, .I = -1 });
+pub const MagneticStrength = BaseScalar(.{ .L = -1, .I = 1 });
+pub const MagneticMoment = BaseScalar(.{ .L = 2, .I = 1 });
+
+// ==========================================
+// Movement
+// ==========================================
+pub const Speed = BaseScalar(.{ .L = 1, .T = -1 });
+pub const Acceleration = BaseScalar(.{ .L = 1, .T = -2 });
+pub const Inertia = BaseScalar(.{ .M = 1, .L = 2 });
+
+// ==========================================
+// Forces / Energy
+// ==========================================
+pub const Force = BaseScalar(.{ .T = -2, .M = 1, .L = 1 });
+pub const Pressure = BaseScalar(.{ .T = -2, .L = -1, .M = 1 });
+pub const Energy = BaseScalar(.{ .T = -2, .L = 2, .M = 1 });
+pub const Power = BaseScalar(.{ .T = -3, .L = 2, .M = 1 });
+
+// ==========================================
+// Dimension
+// ==========================================
+pub const Area = BaseScalar(.{ .L = 2 });
+pub const Volume = BaseScalar(.{ .L = 3 });
+pub const AreaDensity = BaseScalar(.{ .M = 1, .L = -2 });
+pub const Density = BaseScalar(.{ .M = 1, .L = -3 });
+
+// ==========================================
+// Thermal
+// ==========================================
+pub const ThermalHeat = Energy;
+pub const ThermalWork = Energy;
+pub const ThermalCapacity = BaseScalar(.{ .M = 1, .L = 2, .T = -2, .Tr = -1 });
+pub const ThermalCapacityPerMass = BaseScalar(.{ .L = 2, .T = -2, .Tr = -1 });
+pub const ThermalFluxDensity = BaseScalar(.{ .M = 1, .T = -3 });
+pub const ThermalConductance = BaseScalar(.{ .M = 1, .L = 2, .T = -3, .Tr = -1 });
+pub const ThermalConductivity = BaseScalar(.{ .M = 1, .L = 1, .T = -3, .Tr = -1 });
+pub const ThermalResistance = BaseScalar(.{ .M = -1, .L = -2, .T = 3, .Tr = 1 });
+pub const ThermalResistivity = BaseScalar(.{ .M = -1, .L = -1, .T = 3, .Tr = 1 });
+pub const ThermalEntropy = BaseScalar(.{ .M = 1, .L = 2, .T = -2, .Tr = -1 });
+
+// ==========================================
+// Others
+// ==========================================
+pub const Frequency = BaseScalar(.{ .T = -1 });
+pub const Viscosity = BaseScalar(.{ .M = 1, .L = -1, .T = -1 });
+pub const SurfaceTension = BaseScalar(.{ .M = 1, .T = -2 });
+
+// ==========================================
+// Tests
+// ==========================================
+
+test "BaseQuantities - Core dimensions instantiation" {
+    // Basic types via generic wrappers
+    const M = Meter.Of(f32);
+    const distance = M.splat(100);
+    try std.testing.expectEqual(100.0, distance.data[0]);
+    try std.testing.expectEqual(1, M.dims.get(.L));
+    try std.testing.expectEqual(0, M.dims.get(.T));
+
+    // Test specific scale variants
+    const Kmh = Speed.Scaled(f32, .{ .L = .k, .T = .hour });
+    const speed = Kmh.splat(120);
+    try std.testing.expectEqual(120.0, speed.data[0]);
+    try std.testing.expectEqual(.k, @TypeOf(speed).scales.get(.L));
+    try std.testing.expectEqual(.hour, @TypeOf(speed).scales.get(.T));
+}
+
+test "BaseQuantities - Kinematics equations" {
+    const d = Meter.Of(f32).splat(50.0);
+    const t = Second.Of(f32).splat(2.0);
+
+    // Velocity = Distance / Time
+    const v = d.div(t);
+    try std.testing.expectEqual(25.0, v.data[0]);
+    try comptime std.testing.expect(Speed.dims.eql(@TypeOf(v).dims));
+
+    // Acceleration = Velocity / Time
+    const a = v.div(t);
+    try std.testing.expectEqual(12.5, a.data[0]);
+    try comptime std.testing.expect(Acceleration.dims.eql(@TypeOf(a).dims));
+}
+
+test "BaseQuantities - Dynamics (Force and Work)" {
+    // 10 kg
+    const m = Gramm.Scaled(f32, .{ .M = .k }).splat(10.0);
+    // 9.8 m/s^2
+    const a = Acceleration.Of(f32).splat(9.8);
+
+    // Force = mass * acceleration
+    const f = m.mul(a);
+    try std.testing.expectEqual(98, f.data[0]);
+    try comptime std.testing.expect(Force.dims.eql(@TypeOf(f).dims));
+
+    // Energy (Work) = Force * distance
+    const distance = Meter.Of(f32).splat(5.0);
+    const energy = f.mul(distance);
+    try std.testing.expectEqual(490, energy.data[0]);
+    try comptime std.testing.expect(Energy.dims.eql(@TypeOf(energy).dims));
+}
+
+test "BaseQuantities - Electric combinations" {
+    const current = ElectricCurrent.Of(f32).splat(2); // 2 A
+    const time = Second.Of(f32).splat(3.0); // 3 s
+
+    // Charge = Current * time
+    const charge = current.mul(time);
+    try std.testing.expectEqual(6.0, charge.data[0]);
+    try comptime std.testing.expect(ElectricCharge.dims.eql(@TypeOf(charge).dims));
+}
+
+test "Constants - Initialization and dimension checks" {
+    // Speed of Light
+    const c = Constants.SpeedOfLight.Of(f64);
+    try std.testing.expectEqual(299792458.0, c.data[0]);
+    try std.testing.expectEqual(1, @TypeOf(c).dims.get(.L));
+    try std.testing.expectEqual(-1, @TypeOf(c).dims.get(.T));
+
+    // Electron Mass (verifying scale as well)
+    const me = Constants.ElectronMass.Of(f64);
+    try std.testing.expectEqual(9.1093837139e-31, me.data[0]);
+    try std.testing.expectEqual(1, @TypeOf(me).dims.get(.M));
+    try std.testing.expectEqual(.k, @TypeOf(me).scales.get(.M)); // Should be scaled to kg
+
+    // Boltzmann Constant (Complex derived dimensions)
+    const kb = Constants.Boltzmann.Of(f64);
+    try std.testing.expectEqual(1.380649e-23, kb.data[0]);
+    try std.testing.expectEqual(1, @TypeOf(kb).dims.get(.M));
+    try std.testing.expectEqual(2, @TypeOf(kb).dims.get(.L));
+    try std.testing.expectEqual(-2, @TypeOf(kb).dims.get(.T));
+    try std.testing.expectEqual(-1, @TypeOf(kb).dims.get(.Tp));
+    try std.testing.expectEqual(.k, @TypeOf(kb).scales.get(.M));
+
+    // Vacuum Permittivity
+    const eps0 = Constants.VacuumPermittivity.Of(f64);
+    try std.testing.expectEqual(8.8541878188e-12, eps0.data[0]);
+    try std.testing.expectEqual(-1, @TypeOf(eps0).dims.get(.M));
+    try std.testing.expectEqual(-3, @TypeOf(eps0).dims.get(.L));
+    try std.testing.expectEqual(4, @TypeOf(eps0).dims.get(.T));
+    try std.testing.expectEqual(2, @TypeOf(eps0).dims.get(.I));
+
+    // Fine Structure Constant (Dimensionless)
+    const alpha = Constants.FineStructure.Of(f64);
+    try std.testing.expectEqual(0.0072973525643, alpha.data[0]);
+    try std.testing.expectEqual(0, @TypeOf(alpha).dims.get(.M));
+    try std.testing.expectEqual(0, @TypeOf(alpha).dims.get(.L));
+}

src/BaseQuantities.zig

@@ -1,153 +0,0 @@
-const std = @import("std");
-
-// Adjust these imports to match your actual file names
-const Dimensions = @import("Dimensions.zig");
-const Scales = @import("Scales.zig");
-const Quantity = @import("Quantity.zig").Quantity;
-
-/// Helper function to create a clean namespace for each physical dimension.
-/// It exposes the raw dimensions, and easy type-creators for Base or Scaled variants.
-pub fn QtyNamespace(comptime d: anytype) type {
-    return struct {
-        pub const dims = Dimensions.init(d);
-
-        /// Creates a Quantity of this dimension using default scales.
-        /// Example: const V = Quantities.Velocity.Base(f32);
-        pub fn Of(comptime T: type) type {
-            return Quantity(T, dims, Scales.init(.{}));
-        }
-
-        /// Creates a Quantity of this dimension using custom scales.
-        /// Example: const Kmh = Quantities.Velocity.Scaled(f32, Scales.init(.{ .L = .k, .T = .hour }));
-        pub fn Scaled(comptime T: type, comptime s: Scales) type {
-            return Quantity(T, dims, s);
-        }
-    };
-}
-
-// ==========================================
-// Base Quantities
-// ==========================================
-pub const Meter = QtyNamespace(.{ .L = 1 });
-pub const Second = QtyNamespace(.{ .T = 1 });
-pub const Gramm = QtyNamespace(.{ .M = 1 });
-pub const Kelvin = QtyNamespace(.{ .Tr = 1 });
-pub const ElectricCurrent = QtyNamespace(.{ .I = 1 });
-
-// ==========================================
-// Electric
-// ==========================================
-pub const ElectricConductivity = QtyNamespace(.{ .M = -1, .L = -3, .T = 3, .I = 2 });
-pub const ElectricCharge = QtyNamespace(.{ .T = 1, .I = 1 });
-pub const ElectricPotential = QtyNamespace(.{ .T = -3, .L = 2, .M = 1, .I = -1 });
-pub const ElectricResistance = QtyNamespace(.{ .M = 1, .L = 2, .T = -3, .I = -2 });
-pub const ElectricResistivity = QtyNamespace(.{ .M = 1, .L = 3, .T = -3, .I = -2 });
-pub const ElectricCapacitance = QtyNamespace(.{ .T = 4, .L = -2, .M = -1, .I = 2 });
-pub const ElectricImpedance = ElectricResistance;
-pub const MagneticFlux = QtyNamespace(.{ .M = 1, .L = 2, .T = -2, .I = -1 });
-pub const MagneticDensity = QtyNamespace(.{ .M = 1, .T = -2, .I = -1 });
-pub const MagneticStrength = QtyNamespace(.{ .L = -1, .I = 1 }); // Fixed typo from MagneticStrengh
-pub const MagneticMoment = QtyNamespace(.{ .L = 2, .I = 1 });
-
-// ==========================================
-// Movement
-// ==========================================
-pub const Velocity = QtyNamespace(.{ .L = 1, .T = -1 });
-pub const Acceleration = QtyNamespace(.{ .L = 1, .T = -2 });
-pub const Inertia = QtyNamespace(.{ .M = 1, .L = 2 });
-
-// ==========================================
-// Forces / Energy
-// ==========================================
-pub const Force = QtyNamespace(.{ .T = -2, .M = 1, .L = 1 });
-pub const Pressure = QtyNamespace(.{ .T = -2, .L = -1, .M = 1 });
-pub const Energy = QtyNamespace(.{ .T = -2, .L = 2, .M = 1 });
-pub const Power = QtyNamespace(.{ .T = -3, .L = 2, .M = 1 });
-
-// ==========================================
-// Dimension
-// ==========================================
-pub const Area = QtyNamespace(.{ .L = 2 });
-pub const Volume = QtyNamespace(.{ .L = 3 });
-pub const AreaDensity = QtyNamespace(.{ .M = 1, .L = -2 });
-pub const Density = QtyNamespace(.{ .M = 1, .L = -3 });
-
-// ==========================================
-// Thermal
-// ==========================================
-pub const ThermalHeat = Energy;
-pub const ThermalWork = Energy;
-pub const ThermalCapacity = QtyNamespace(.{ .M = 1, .L = 2, .T = -2, .Tr = -1 });
-pub const ThermalCapacityPerMass = QtyNamespace(.{ .L = 2, .T = -2, .Tr = -1 });
-pub const ThermalFluxDensity = QtyNamespace(.{ .M = 1, .T = -3 }); // Fixed typo from ThermalluxDensity
-pub const ThermalConductance = QtyNamespace(.{ .M = 1, .L = 2, .T = -3, .Tr = -1 });
-pub const ThermalConductivity = QtyNamespace(.{ .M = 1, .L = 1, .T = -3, .Tr = -1 });
-pub const ThermalResistance = QtyNamespace(.{ .M = -1, .L = -2, .T = 3, .Tr = 1 });
-pub const ThermalResistivity = QtyNamespace(.{ .M = -1, .L = -1, .T = 3, .Tr = 1 });
-pub const ThermalEntropy = QtyNamespace(.{ .M = 1, .L = 2, .T = -2, .Tr = -1 });
-
-// ==========================================
-// Others
-// ==========================================
-pub const Frequency = QtyNamespace(.{ .T = -1 });
-pub const Viscosity = QtyNamespace(.{ .M = 1, .L = -1, .T = -1 });
-pub const SurfaceTension = QtyNamespace(.{ .M = 1, .T = -2 }); // Corrected from MT-2a
-
-test "BaseQuantities - Core dimensions instantiation" {
-    // Basic types via generic wrappers
-    const M = Meter.Of(f32);
-    const distance = M{ .value = 100.0 };
-    try std.testing.expectEqual(100.0, distance.value);
-    try std.testing.expectEqual(1, M.dims.get(.L));
-    try std.testing.expectEqual(0, M.dims.get(.T));
-
-    // Test specific scale variants
-    const Kmh = Velocity.Scaled(f32, Scales.init(.{ .L = .k, .T = .hour }));
-    const speed = Kmh{ .value = 120.0 };
-    try std.testing.expectEqual(120.0, speed.value);
-    try std.testing.expectEqual(.k, @TypeOf(speed).scales.get(.L));
-    try std.testing.expectEqual(.hour, @TypeOf(speed).scales.get(.T));
-}
-
-test "BaseQuantities - Kinematics equations" {
-    const d = Meter.Of(f32){ .value = 50.0 };
-    const t = Second.Of(f32){ .value = 2.0 };
-
-    // Velocity = Distance / Time
-    const v = d.divBy(t);
-    try std.testing.expectEqual(25.0, v.value);
-    try std.testing.expect(Velocity.dims.eql(@TypeOf(v).dims));
-
-    // Acceleration = Velocity / Time
-    const a = v.divBy(t);
-    try std.testing.expectEqual(12.5, a.value);
-    try std.testing.expect(Acceleration.dims.eql(@TypeOf(a).dims));
-}
-
-test "BaseQuantities - Dynamics (Force and Work)" {
-    // 10 kg
-    const m = Gramm.Scaled(f32, Scales.init(.{ .M = .k })){ .value = 10.0 };
-    // 9.8 m/s^2
-    const a = Acceleration.Of(f32){ .value = 9.8 };
-
-    // Force = mass * acceleration
-    const f = m.mulBy(a);
-    try std.testing.expectEqual(98000, f.value);
-    try std.testing.expect(Force.dims.eql(@TypeOf(f).dims));
-
-    // Energy (Work) = Force * distance
-    const distance = Meter.Of(f32){ .value = 5.0 };
-    const energy = f.mulBy(distance);
-    try std.testing.expectEqual(490000, energy.value);
-    try std.testing.expect(Energy.dims.eql(@TypeOf(energy).dims));
-}
-
-test "BaseQuantities - Electric combinations" {
-    const current = ElectricCurrent.Of(f32){ .value = 2.0 }; // 2 A
-    const time = Second.Of(f32){ .value = 3.0 }; // 3 s
-
-    // Charge = Current * time
-    const charge = current.mulBy(time);
-    try std.testing.expectEqual(6.0, charge.value);
-    try std.testing.expect(ElectricCharge.dims.eql(@TypeOf(charge).dims));
-}

src/Dimensions.zig

@@ -1,5 +1,15 @@
 const std = @import("std");
 
+pub const ArgOpts = struct {
+    L: comptime_int = 0,
+    M: comptime_int = 0,
+    T: comptime_int = 0,
+    I: comptime_int = 0,
+    Tp: comptime_int = 0,
+    N: comptime_int = 0,
+    J: comptime_int = 0,
+};
+
 pub const Dimension = enum {
     /// Length
     L,
@@ -31,22 +41,26 @@ pub const Dimension = enum {
 
 // --------- Dimensions struct ---------
 
+/// Holds the exponent of each SI base dimension for a given quantity (e.g. velocity = L¹T⁻¹).
+/// All values are `comptime_int` — no runtime storage.
 const Self = @This();
 
-data: std.EnumArray(Dimension, i8),
+data: std.EnumArray(Dimension, comptime_int),
 
-pub fn init(comptime init_val: anytype) Self {
+/// Create a `Dimensions` from a struct literal, e.g. `.{ .L = 1, .T = -1 }`.
-    var s = Self{ .data = std.EnumArray(Dimension, i8).initFill(0) };
+/// Unspecified dimensions default to 0.
-    inline for (std.meta.fields(@TypeOf(init_val))) |f|
+pub fn init(init_val: ArgOpts) Self {
+    var s = Self{ .data = std.EnumArray(Dimension, comptime_int).initFill(0) };
+    for (std.meta.fields(@TypeOf(init_val))) |f|
         s.data.set(@field(Dimension, f.name), @field(init_val, f.name));
     return s;
 }
 
-pub fn initFill(val: i8) Self {
+pub fn initFill(val: comptime_int) Self {
-    return .{ .data = std.EnumArray(Dimension, i8).initFill(val) };
+    return .{ .data = std.EnumArray(Dimension, comptime_int).initFill(val) };
 }
 
-pub fn get(self: Self, key: Dimension) i8 {
+pub fn get(self: Self, key: Dimension) comptime_int {
     return self.data.get(key);
 }
 
@@ -54,32 +68,62 @@ pub fn set(self: *Self, key: Dimension, val: i8) void {
     self.data.set(key, val);
 }
 
-pub fn add(comptime a: Self, comptime b: Self) Self {
+pub fn argsOpt(self: Self) ArgOpts {
+    var args: ArgOpts = undefined;
+    for (std.enums.values(Dimension)) |d|
+        @field(args, @tagName(d)) = self.get(d);
+    return args;
+}
+
+/// Add exponents component-wise. Used internally by `mul`.
+pub fn add(a: Self, b: Self) Self {
     var result = Self.initFill(0);
     for (std.enums.values(Dimension)) |d|
         result.set(d, a.get(d) + b.get(d));
     return result;
 }
 
-pub fn sub(comptime a: Self, comptime b: Self) Self {
+/// Subtract exponents component-wise. Used internally by `div`.
-    @setEvalBranchQuota(10_000);
+pub fn sub(a: Self, b: Self) Self {
     var result = Self.initFill(0);
     for (std.enums.values(Dimension)) |d|
         result.set(d, a.get(d) - b.get(d));
     return result;
 }
 
-pub fn eql(comptime a: Self, comptime b: Self) bool {
+/// Multiply exponents by a scalar integer. Used internally by `pow` in Scalar.
+pub fn scale(a: Self, exp: comptime_int) Self {
+    var result = Self.initFill(0);
+    for (std.enums.values(Dimension)) |d|
+        result.set(d, a.get(d) * exp);
+    return result;
+}
+
+pub fn div(a: Self, exp: comptime_int) Self {
+    var result = Self.initFill(0);
+    for (std.enums.values(Dimension)) |d|
+        result.set(d, a.get(d) / exp);
+    return result;
+}
+
+/// Returns true if every dimension exponent is equal. Used to enforce type compatibility in `add`, `sub`, `to`.
+pub fn eql(a: Self, b: Self) bool {
     for (std.enums.values(Dimension)) |d|
         if (a.get(d) != b.get(d)) return false;
     return true;
 }
 
-pub fn str(comptime a: Self) []const u8 {
+pub fn isSquare(a: Self) bool {
+    for (std.enums.values(Dimension)) |d|
+        if (a.get(d) % 2 != 0) return false;
+    return true;
+}
+
+pub fn str(a: Self) []const u8 {
     var out: []const u8 = "";
     const dims = std.enums.values(Dimension);
 
-    inline for (dims) |d| {
+    for (dims) |d| {
         const val = a.get(d);
         if (val != 0) {
             out = out ++ @tagName(d) ++ std.fmt.comptimePrint("{d}", .{val});

src/Quantity.zig

@@ -1,638 +0,0 @@
-const std = @import("std");
-const hlp = @import("helper.zig");
-
-const QuantityVec = @import("QuantityVec.zig").QuantityVec;
-const Scales = @import("Scales.zig");
-const UnitScale = Scales.UnitScale;
-const Dimensions = @import("Dimensions.zig");
-const Dimension = Dimensions.Dimension;
-
-pub fn Quantity(T: type, d: Dimensions, s: Scales) type {
-    return struct {
-        value: T,
-
-        const Self = @This();
-        pub const Vec3: type = QuantityVec(3, Self);
-        pub const ValueType: type = T;
-
-        pub const dims: Dimensions = d;
-        pub const scales = s;
-
-        pub fn add(self: Self, rhs: anytype) Quantity(
-            T,
-            dims,
-            scales.min(@TypeOf(rhs).scales),
-        ) {
-            if (comptime !dims.eql(@TypeOf(rhs).dims))
-                @compileError("Dimension mismatch in add: " ++ dims.str() ++ " vs " ++ @TypeOf(rhs).dims.str());
-
-            const TargetType = Quantity(T, dims, scales.min(@TypeOf(rhs).scales));
-            const lhs_converted = self.to(TargetType);
-            const rhs_converted = rhs.to(TargetType);
-
-            return .{ .value = lhs_converted.value + rhs_converted.value };
-        }
-
-        pub fn sub(self: Self, rhs: anytype) Quantity(
-            T,
-            dims,
-            scales.min(@TypeOf(rhs).scales),
-        ) {
-            if (comptime !dims.eql(@TypeOf(rhs).dims))
-                @compileError("Dimension mismatch in sub: " ++ dims.str() ++ " vs " ++ @TypeOf(rhs).dims.str());
-
-            const TargetType = Quantity(T, dims, scales.min(@TypeOf(rhs).scales));
-            const lhs_converted = self.to(TargetType);
-            const rhs_converted = rhs.to(TargetType);
-
-            return .{ .value = lhs_converted.value - rhs_converted.value };
-        }
-
-        pub fn mulBy(self: Self, rhs: anytype) Quantity(
-            T,
-            dims.add(@TypeOf(rhs).dims),
-            scales.min(@TypeOf(rhs).scales),
-        ) {
-            const self_ = self.to(Quantity(T, dims, scales.min(@TypeOf(rhs).scales)));
-            const rhs_ = rhs.to(Quantity(T, @TypeOf(rhs).dims, scales.min(@TypeOf(rhs).scales)));
-            return .{ .value = self_.value * rhs_.value };
-        }
-
-        pub fn divBy(self: Self, rhs: anytype) Quantity(T, dims.sub(@TypeOf(rhs).dims), scales.min(@TypeOf(rhs).scales)) {
-            const self_ = self.to(Quantity(T, dims, scales.min(@TypeOf(rhs).scales)));
-            const rhs_ = rhs.to(Quantity(T, @TypeOf(rhs).dims, scales.min(@TypeOf(rhs).scales)));
-
-            if (comptime @typeInfo(T) == .int) {
-                return .{ .value = @divTrunc(self_.value, rhs_.value) };
-            } else {
-                return .{ .value = self_.value / rhs_.value };
-            }
-        }
-
-        pub fn scale(self: Self, sc: T) Self {
-            return .{ .value = self.value * sc };
-        }
-
-        pub fn to(self: Self, comptime Dest: type) Dest {
-            if (comptime !dims.eql(Dest.dims))
-                @compileError("Dimension mismatch in to: " ++ dims.str() ++ " vs " ++ Dest.dims.str());
-            if (comptime @TypeOf(self) == Dest)
-                return self;
-
-            const DestT = Dest.ValueType;
-            const ratio = comptime (scales.getFactor(dims) / Dest.scales.getFactor(Dest.dims));
-
-            // Fast-path: Native pure-integer exact conversions
-            if (comptime @typeInfo(T) == .int and @typeInfo(DestT) == .int) {
-                if (ratio >= 1.0 and @round(ratio) == ratio) {
-                    const mult: DestT = @intFromFloat(ratio);
-                    return .{ .value = @as(DestT, @intCast(self.value)) * mult };
-                } else if (ratio < 1.0 and @round(1.0 / ratio) == 1.0 / ratio) {
-                    const div: DestT = @intFromFloat(1.0 / ratio);
-                    const val = @as(DestT, @intCast(self.value));
-                    const half = div / 2;
-                    // Native round-to-nearest
-                    const rounded = if (val >= 0) @divTrunc(val + half, div) else @divTrunc(val - half, div);
-                    return .{ .value = rounded };
-                }
-            }
-
-            // Fallback preserving native Float types (e.g., f128 shouldn't downcast to f64)
-            if (comptime @typeInfo(DestT) == .float) {
-                const val_f = switch (@typeInfo(T)) {
-                    .int => @as(DestT, @floatFromInt(self.value)),
-                    .float => @as(DestT, @floatCast(self.value)),
-                    else => unreachable,
-                };
-                return .{ .value = val_f * @as(DestT, @floatCast(ratio)) };
-            } else {
-                const val_f = switch (@typeInfo(T)) {
-                    .int => @as(f64, @floatFromInt(self.value)),
-                    .float => @as(f64, @floatCast(self.value)),
-                    else => unreachable,
-                };
-                return .{ .value = @intFromFloat(@round(val_f * ratio)) };
-            }
-        }
-        pub fn Vec(self: Self, comptime len: usize) QuantityVec(len, Self) {
-            return QuantityVec(len, Self).initDefault(self.value);
-        }
-
-        pub fn vec3(self: Self) Vec3 {
-            return Vec3.initDefault(self.value);
-        }
-
-        pub fn format(
-            self: Self,
-            writer: *std.Io.Writer,
-        ) !void {
-            try writer.print("{d}", .{self.value});
-            var iter = std.EnumSet(Dimension).initFull().iterator();
-            var first = true;
-            while (iter.next()) |bu| {
-                const v = dims.get(bu);
-                if (v == 0) continue;
-                if (!first)
-                    try writer.writeAll(".");
-
-                first = false;
-
-                const uscale = scales.get(bu);
-                if (bu == .T and (uscale == .min or uscale == .hour or uscale == .year))
-                    try writer.print("{s}", .{uscale.str()})
-                else
-                    try writer.print("{s}{s}", .{ uscale.str(), bu.unit() });
-
-                if (v != 1)
-                    try hlp.printSuperscript(writer, v);
-            }
-        }
-    };
-}
-
-test "Generate quantity" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = -3 }));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{ .T = .n }));
-
-    const distance = Meter{ .value = 10 };
-    const time = Second{ .value = 2 };
-
-    try std.testing.expectEqual(10, distance.value);
-    try std.testing.expectEqual(2, time.value);
-}
-
-test "Add" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-
-    const distance = Meter{ .value = 10 };
-    const distance2 = Meter{ .value = 20 };
-
-    const added = distance.add(distance2);
-    try std.testing.expectEqual(30, added.value);
-    try std.testing.expectEqual(1, @TypeOf(added).dims.get(.L));
-    std.debug.print("KiloMeter {f} + {f} = {f} OK\n", .{ distance, distance2, added });
-
-    const KiloMeter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const distance3 = KiloMeter{ .value = 2 };
-    const added2 = distance.add(distance3);
-    try std.testing.expectEqual(2010, added2.value);
-    try std.testing.expectEqual(1, @TypeOf(added2).dims.get(.L));
-    std.debug.print("KiloMeter {f} + {f} = {f} OK\n", .{ distance, distance3, added2 });
-
-    const added3 = distance3.add(distance).to(KiloMeter);
-    try std.testing.expectEqual(2, added3.value);
-    try std.testing.expectEqual(1, @TypeOf(added3).dims.get(.L));
-    std.debug.print("KiloMeter {f} + {f} = {f} OK\n", .{ distance3, distance, added3 });
-
-    const KiloMeter_f = Quantity(f64, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const distance4 = KiloMeter_f{ .value = 2 };
-    const added4 = distance4.add(distance).to(KiloMeter_f);
-    try std.testing.expectApproxEqAbs(2.01, added4.value, 0.000001);
-    try std.testing.expectEqual(1, @TypeOf(added4).dims.get(.L));
-    std.debug.print("KiloMeter_f {f} + {f} = {f} OK\n", .{ distance4, distance, added4 });
-}
-
-test "Sub" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const KiloMeter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const KiloMeter_f = Quantity(f64, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-
-    const a = Meter{ .value = 500 };
-    const b = Meter{ .value = 200 };
-    const diff = a.sub(b);
-    try std.testing.expectEqual(300, diff.value);
-    std.debug.print("Sub: {f} - {f} = {f} OK\n", .{ a, b, diff });
-
-    const km = KiloMeter{ .value = 1 };
-    const diff2 = a.sub(km);
-    std.debug.print("Sub cross-scale: {f} - {f} = {f}\n", .{ a, km, diff2 });
-
-    const km_f = KiloMeter_f{ .value = 2.5 };
-    const m_f = Meter{ .value = 500 };
-    const diff3 = km_f.sub(m_f);
-    try std.testing.expectApproxEqAbs(2000, diff3.value, 1e-4);
-    std.debug.print("Sub float cross-scale: {f} - {f} = {f} OK\n", .{ km_f, m_f, diff3 });
-}
-
-test "MulBy" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const d = Meter{ .value = 3.0 };
-    const t = Second{ .value = 4.0 };
-
-    const area_time = d.mulBy(t);
-    try std.testing.expectEqual(12, area_time.value);
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.L));
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.T));
-    std.debug.print("MulBy: {f} * {f} = {f} OK\n", .{ d, t, area_time });
-
-    const d2 = Meter{ .value = 5.0 };
-    const area = d.mulBy(d2);
-    try std.testing.expectEqual(15, area.value);
-    try std.testing.expectEqual(2, @TypeOf(area).dims.get(.L));
-    try std.testing.expectEqual(0, @TypeOf(area).dims.get(.T));
-    std.debug.print("MulBy: {f} * {f} = {f} OK\n", .{ d, d2, area });
-}
-
-test "MulBy with scale" {
-    const KiloMeter = Quantity(f32, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const KiloGram = Quantity(f32, Dimensions.init(.{ .M = 1 }), Scales.init(.{ .M = .k }));
-
-    const dist = KiloMeter{ .value = 2.0 };
-    const mass = KiloGram{ .value = 3.0 };
-    const prod = dist.mulBy(mass);
-    try std.testing.expectEqual(1, @TypeOf(prod).dims.get(.L));
-    try std.testing.expectEqual(1, @TypeOf(prod).dims.get(.M));
-    std.debug.print("MulBy scaled: {f} * {f} = {f} OK\n", .{ dist, mass, prod });
-}
-
-test "MulBy with type change" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const Second = Quantity(f64, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-    const KmSec = Quantity(i64, Dimensions.init(.{ .L = 1, .T = 1 }), Scales.init(.{ .L = .k }));
-    const KmSec_f = Quantity(f32, Dimensions.init(.{ .L = 1, .T = 1 }), Scales.init(.{ .L = .k }));
-
-    const d = Meter{ .value = 3.0 };
-    const t = Second{ .value = 4.0 };
-
-    const area_time = d.mulBy(t).to(KmSec);
-    const area_time_f = d.mulBy(t).to(KmSec_f);
-    try std.testing.expectEqual(12, area_time.value);
-    try std.testing.expectApproxEqAbs(12, area_time_f.value, 0.0001);
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.L));
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.T));
-    std.debug.print("MulBy: {f} * {f} = {f} OK\n", .{ d, t, area_time });
-}
-
-test "MulBy small" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .n }));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const d = Meter{ .value = 3.0 };
-    const t = Second{ .value = 4.0 };
-
-    const area_time = d.mulBy(t);
-    try std.testing.expectEqual(12, area_time.value);
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.L));
-    try std.testing.expectEqual(1, @TypeOf(area_time).dims.get(.T));
-    std.debug.print("MulBy: {f} * {f} = {f} OK\n", .{ d, t, area_time });
-}
-
-test "Scale" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const d = Meter{ .value = 7 };
-    const scaled = d.scale(3);
-    try std.testing.expectEqual(21, scaled.value);
-    try std.testing.expectEqual(1, @TypeOf(scaled).dims.get(.L));
-    std.debug.print("Scale int: {f} * 3 = {f} OK\n", .{ d, scaled });
-
-    const t = Second{ .value = 1.5 };
-    const scaled_f = t.scale(4.0);
-    try std.testing.expectApproxEqAbs(@as(f32, 6.0), scaled_f.value, 1e-4);
-    std.debug.print("Scale float: {f} * 4 = {f} OK\n", .{ t, scaled_f });
-}
-
-test "Chained: velocity and acceleration" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const dist = Meter{ .value = 100.0 };
-    const t1 = Second{ .value = 5.0 };
-    const velocity = dist.divBy(t1);
-    try std.testing.expectEqual(20, velocity.value);
-    try std.testing.expectEqual(1, @TypeOf(velocity).dims.get(.L));
-    try std.testing.expectEqual(-1, @TypeOf(velocity).dims.get(.T));
-
-    const t2 = Second{ .value = 4.0 };
-    const accel = velocity.divBy(t2);
-    try std.testing.expectEqual(5, accel.value);
-    try std.testing.expectEqual(1, @TypeOf(accel).dims.get(.L));
-    try std.testing.expectEqual(-2, @TypeOf(accel).dims.get(.T));
-
-    std.debug.print("Velocity: {f}, Acceleration: {f} OK\n", .{ velocity, accel });
-}
-
-test "DivBy integer exact" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Second = Quantity(f32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const dist = Meter{ .value = 120 };
-    const time = Second{ .value = 4 };
-    const vel = dist.divBy(time);
-
-    try std.testing.expectEqual(30, vel.value);
-    try std.testing.expectEqual(1, @TypeOf(vel).dims.get(.L));
-    try std.testing.expectEqual(-1, @TypeOf(vel).dims.get(.T));
-    std.debug.print("DivBy int: {f} / {f} = {f} OK\n", .{ dist, time, vel });
-}
-
-test "Conversion chain: km -> m -> cm" {
-    const KiloMeter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const CentiMeter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .c }));
-
-    const km = KiloMeter{ .value = 15 };
-    const m = km.to(Meter);
-    const cm = m.to(CentiMeter);
-
-    try std.testing.expectEqual(15_000, m.value);
-    try std.testing.expectEqual(1_500_000, cm.value);
-    std.debug.print("Chain: {f} -> {f} -> {f} OK\n", .{ km, m, cm });
-}
-
-test "Conversion: hours -> minutes -> seconds" {
-    const Hour = Quantity(i128, Dimensions.init(.{ .T = 1 }), Scales.init(.{ .T = .hour }));
-    const Minute = Quantity(i128, Dimensions.init(.{ .T = 1 }), Scales.init(.{ .T = .min }));
-    const Second = Quantity(i128, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-    const h = Hour{ .value = 1.0 };
-    const min = h.to(Minute);
-    const sec = min.to(Second);
-
-    try std.testing.expectEqual(60, min.value);
-    try std.testing.expectEqual(3600, sec.value);
-    std.debug.print("Time chain: {f} -> {f} -> {f} OK\n", .{ h, min, sec });
-}
-
-test "Negative values" {
-    const Meter = Quantity(i128, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-
-    const a = Meter{ .value = 5 };
-    const b = Meter{ .value = 20 };
-    const diff = a.sub(b);
-    try std.testing.expectEqual(-15, diff.value);
-    std.debug.print("Negative sub: {f} - {f} = {f} OK\n", .{ a, b, diff });
-}
-
-test "Format Quantity" {
-    const MeterPerSecondSq = Quantity(
-        f32,
-        Dimensions.init(.{ .L = 1, .T = -2 }),
-        Scales.init(.{ .T = .n }),
-    );
-    const KgMeterPerSecond = Quantity(
-        f32,
-        Dimensions.init(.{ .M = 1, .L = 1, .T = -1 }),
-        Scales.init(.{ .M = .k }),
-    );
-
-    const accel = MeterPerSecondSq{ .value = 9.81 };
-    const momentum = KgMeterPerSecond{ .value = 42.0 };
-
-    std.debug.print("Acceleration: {f}\n", .{accel});
-    std.debug.print("Momentum: {f}\n", .{momentum});
-}
-
-test "Benchmark" {
-    const Io = std.Io;
-    const ITERS: usize = 100_000;
-    const SAMPLES: usize = 10; // Number of samples for stats
-
-    var gsink: f64 = 0;
-    const io = std.testing.io;
-
-    // Standard Zig 0.16 timestamp retrieval
-    const getTime = struct {
-        fn f(i: Io) Io.Timestamp {
-            return Io.Clock.awake.now(i);
-        }
-    }.f;
-
-    const fold = struct {
-        fn f(comptime TT: type, s: *f64, v: TT) void {
-            s.* += if (comptime @typeInfo(TT) == .float)
-                @as(f64, @floatCast(v))
-            else
-                @as(f64, @floatFromInt(v));
-        }
-    }.f;
-
-    const getVal = struct {
-        fn f(comptime TT: type, i: usize, comptime mask: u7) TT {
-            const v: u8 = @as(u8, @truncate(i & @as(usize, mask))) + 1;
-            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
-        }
-    }.f;
-
-    const Stats = struct {
-        median: f64,
-        delta: f64,
-        ops_per_sec: f64,
-    };
-
-    const computeStats = struct {
-        fn f(samples: []f64, iters: usize) Stats {
-            std.mem.sort(f64, samples, {}, std.sort.asc(f64));
-            const mid = samples.len / 2;
-            const median_ns = if (samples.len % 2 == 0) (samples[mid - 1] + samples[mid]) / 2.0 else samples[mid];
-
-            const low = samples[0];
-            const high = samples[samples.len - 1];
-            const delta_ns = (high - low) / 2.0;
-
-            const ns_per_op = median_ns / @as(f64, @floatFromInt(iters));
-            return .{
-                .median = ns_per_op,
-                .delta = (delta_ns / @as(f64, @floatFromInt(iters))),
-                .ops_per_sec = 1_000_000_000.0 / ns_per_op,
-            };
-        }
-    }.f;
-
-    std.debug.print(
-        \\
-        \\ Quantity<T> benchmark — {d} iterations, {d} samples/cell
-        \\
-        \\┌───────────────────┬──────┬─────────────────────┬─────────────────────┐
-        \\│ Operation         │ Type │ ns / op (± delta)   │ Throughput (ops/s)  │
-        \\├───────────────────┼──────┼─────────────────────┼─────────────────────┤
-        \\
-    , .{ ITERS, SAMPLES });
-
-    const Types = .{ i16, i32, i64, i128, i256, f32, f64, f128 };
-    const TNames = .{ "i16", "i32", "i64", "i128", "i256", "f32", "f64", "f128" };
-    const Ops = .{ "add", "sub", "mulBy", "divBy", "scale", "to" };
-
-    var results_matrix: [Ops.len][Types.len]f64 = undefined;
-
-    comptime var tidx: usize = 0;
-    inline for (Types, TNames) |T, tname| {
-        const M = Quantity(T, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-        const KM = Quantity(T, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-        const S = Quantity(T, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-        inline for (Ops, 0..) |op_name, oidx| {
-            var samples: [SAMPLES]f64 = undefined;
-
-            for (0..SAMPLES) |s_idx| {
-                var sink: T = 0;
-                const t_start = getTime(io);
-
-                for (0..ITERS) |i| {
-                    const r = if (comptime std.mem.eql(u8, op_name, "add"))
-                        (M{ .value = getVal(T, i, 63) }).add(M{ .value = getVal(T, i +% 7, 63) })
-                    else if (comptime std.mem.eql(u8, op_name, "sub"))
-                        (M{ .value = getVal(T, i +% 10, 63) }).sub(M{ .value = getVal(T, i, 63) })
-                    else if (comptime std.mem.eql(u8, op_name, "mulBy"))
-                        (M{ .value = getVal(T, i, 63) }).mulBy(M{ .value = getVal(T, i +% 1, 63) })
-                    else if (comptime std.mem.eql(u8, op_name, "divBy"))
-                        (M{ .value = getVal(T, i +% 10, 63) }).divBy(S{ .value = getVal(T, i, 63) })
-                    else if (comptime std.mem.eql(u8, op_name, "scale"))
-                        (M{ .value = getVal(T, i, 63) }).scale(getVal(T, i +% 2, 63))
-                    else
-                        (KM{ .value = getVal(T, i, 15) }).to(M);
-
-                    if (comptime @typeInfo(T) == .float) sink += r.value else sink ^= r.value;
-                }
-
-                const t_end = getTime(io);
-                samples[s_idx] = @as(f64, @floatFromInt(t_start.durationTo(t_end).toNanoseconds()));
-                fold(T, &gsink, sink);
-            }
-
-            const stats = computeStats(&samples, ITERS);
-            results_matrix[oidx][tidx] = stats.median;
-
-            std.debug.print("│ {s:<17} │ {s:<4} │ {d:>8.2} ns ±{d:<6.2} │ {d:>19.0} │\n", .{ op_name, tname, stats.median, stats.delta, stats.ops_per_sec });
-        }
-
-        if (comptime tidx < Types.len - 1) {
-            std.debug.print("├───────────────────┼──────┼─────────────────────┼─────────────────────┤\n", .{});
-        }
-        tidx += 1;
-    }
-
-    // Median Summary Table
-    std.debug.print("└───────────────────┴──────┴─────────────────────┴─────────────────────┘\n\n", .{});
-    std.debug.print("Median Summary (ns/op):\n", .{});
-
-    std.debug.print("┌──────────────┬───────┬───────┬───────┬───────┬───────┬───────┬───────┬───────┐\n", .{});
-    std.debug.print("│  Operation   │  i16  │  i32  │  i64  │  i128 │  i256 │  f32  │  f64  │  f128 │\n", .{});
-    std.debug.print("├──────────────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┤\n", .{});
-
-    inline for (Ops, 0..) |op_name, oidx| {
-        std.debug.print("│  {s:<11} │", .{op_name});
-        var i: usize = 0;
-        while (i < Types.len) : (i += 1) {
-            std.debug.print("{d:>6.1} │", .{results_matrix[oidx][i]});
-        }
-        std.debug.print("\n", .{});
-    }
-
-    std.debug.print("└──────────────┴───────┴───────┴───────┴───────┴───────┴───────┴───────┴───────┘\n", .{});
-
-    std.debug.print("\nAnti-optimisation sink: {d:.4}\n", .{gsink});
-    try std.testing.expect(gsink != 0);
-}
-
-test "Overhead Analysis: Quantity vs Native" {
-    const Io = std.Io;
-    const ITERS: usize = 100_000;
-    const SAMPLES: usize = 5;
-    const io = std.testing.io;
-
-    const getTime = struct {
-        fn f(i: Io) Io.Timestamp {
-            return Io.Clock.awake.now(i);
-        }
-    }.f;
-
-    const fold = struct {
-        fn f(comptime TT: type, s: *f64, v: TT) void {
-            s.* += if (comptime @typeInfo(TT) == .float)
-                @as(f64, @floatCast(v))
-            else
-                @as(f64, @floatFromInt(v));
-        }
-    }.f;
-
-    // Helper to safely get a value of type T from a loop index
-    const getValT = struct {
-        fn f(comptime TT: type, i: usize) TT {
-            const v = (i % 100) + 1;
-            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
-        }
-    }.f;
-
-    const Types = .{ i32, i64, i128, f32, f64 };
-    const TNames = .{ "i32", "i64", "i128", "f32", "f64" };
-    const Ops = .{ "add", "mulBy", "divBy" };
-
-    var gsink: f64 = 0;
-
-    std.debug.print(
-        \\
-        \\ Quantity vs Native Overhead Analysis
-        \\
-        \\┌───────────┬──────┬───────────┬───────────┬───────────┐
-        \\│ Operation │ Type │ Native    │ Quantity  │ Slowdown  │
-        \\├───────────┼──────┼───────────┼───────────┼───────────┤
-        \\
-    , .{});
-
-    inline for (Ops, 0..) |op_name, j| {
-        inline for (Types, 0..) |T, tidx| {
-            var native_total_ns: f64 = 0;
-            var quantity_total_ns: f64 = 0;
-
-            const M = Quantity(T, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-            const S = Quantity(T, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-
-            for (0..SAMPLES) |_| {
-                // --- 1. Benchmark Native ---
-                var n_sink: T = 0;
-                const n_start = getTime(io);
-                for (0..ITERS) |i| {
-                    const a = getValT(T, i);
-                    const b = getValT(T, 2);
-                    const r = if (comptime std.mem.eql(u8, op_name, "add"))
-                        a + b
-                    else if (comptime std.mem.eql(u8, op_name, "mulBy"))
-                        a * b
-                    else
-                        if (comptime @typeInfo(T) == .int) @divTrunc(a, b) else a / b;
-                    
-                    if (comptime @typeInfo(T) == .float) n_sink += r else n_sink ^= r;
-                }
-                const n_end = getTime(io);
-                native_total_ns += @as(f64, @floatFromInt(n_start.durationTo(n_end).toNanoseconds()));
-                fold(T, &gsink, n_sink);
-
-                // --- 2. Benchmark Quantity ---
-                var q_sink: T = 0;
-                const q_start = getTime(io);
-                for (0..ITERS) |i| {
-                    const qa = M{ .value = getValT(T, i) };
-                    const qb = if (comptime std.mem.eql(u8, op_name, "divBy")) S{ .value = getValT(T, 2) } else M{ .value = getValT(T, 2) };
-                    
-                    const r = if (comptime std.mem.eql(u8, op_name, "add"))
-                        qa.add(qb)
-                    else if (comptime std.mem.eql(u8, op_name, "mulBy"))
-                        qa.mulBy(qb)
-                    else
-                        qa.divBy(qb);
-
-                    if (comptime @typeInfo(T) == .float) q_sink += r.value else q_sink ^= r.value;
-                }
-                const q_end = getTime(io);
-                quantity_total_ns += @as(f64, @floatFromInt(q_start.durationTo(q_end).toNanoseconds()));
-                fold(T, &gsink, q_sink);
-            }
-
-            const avg_n = (native_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
-            const avg_q = (quantity_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
-            const slowdown = avg_q / avg_n;
-
-            std.debug.print("│ {s:<9} │ {s:<4} │ {d:>7.2}ns │ {d:>7.2}ns │ {d:>8.2}x │\n", .{
-                op_name, TNames[tidx], avg_n, avg_q, slowdown,
-            });
-        }
-       if (j != Ops.len - 1) std.debug.print("├───────────┼──────┼───────────┼───────────┼───────────┤\n", .{});
-    }
-
-    std.debug.print("└───────────┴──────┴───────────┴───────────┴───────────┘\n", .{});
-    try std.testing.expect(gsink != 0);
-}

src/QuantityVec.zig

@@ -1,428 +0,0 @@
-const std = @import("std");
-const hlp = @import("helper.zig");
-
-const Quantity = @import("Quantity.zig").Quantity;
-const Scales = @import("Scales.zig");
-const UnitScale = Scales.UnitScale;
-const Dimensions = @import("Dimensions.zig");
-const Dimension = Dimensions.Dimension;
-
-pub fn QuantityVec(comptime len: usize, comptime Q: type) type {
-    const T = Q.ValueType;
-    const d: Dimensions = Q.dims;
-    const s: Scales = Q.scales;
-
-    return struct {
-        data: [len]T,
-
-        const Self = @This();
-        pub const QuantityType = Q;
-        pub const ValueType = T;
-        pub const dims: Dimensions = d;
-        pub const scales = s;
-
-        pub const zero = initDefault(0);
-        pub const one = initDefault(1);
-
-        pub fn initDefault(v: T) Self {
-            var data: [len]T = undefined;
-            for (&data) |*item| item.* = v;
-            return .{ .data = data };
-        }
-
-        pub fn add(self: Self, rhs: anytype) QuantityVec(len, Quantity(T, d, s.min(@TypeOf(rhs).scales))) {
-            const Tr = @TypeOf(rhs);
-            var res: QuantityVec(len, Quantity(T, d, s.min(Tr.scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = (Q{ .value = v }).add(Tr.QuantityType{ .value = rhs.data[i] });
-                res.data[i] = q.value;
-            }
-            return res;
-        }
-
-        pub fn sub(self: Self, rhs: anytype) QuantityVec(len, Quantity(T, d, s.min(@TypeOf(rhs).scales))) {
-            const Tr = @TypeOf(rhs);
-            var res: QuantityVec(len, Quantity(T, d, s.min(Tr.scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = (Q{ .value = v }).sub(Tr.QuantityType{ .value = rhs.data[i] });
-                res.data[i] = q.value;
-            }
-            return res;
-        }
-
-        pub fn divBy(
-            self: Self,
-            rhs: anytype,
-        ) QuantityVec(len, Quantity(T, d.sub(@TypeOf(rhs).dims), s.min(@TypeOf(rhs).scales))) {
-            const Tr = @TypeOf(rhs);
-            var res: QuantityVec(len, Quantity(T, d.sub(Tr.dims), s.min(Tr.scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = (Q{ .value = v }).divBy(Tr.QuantityType{ .value = rhs.data[i] });
-                res.data[i] = q.value;
-            }
-            return res;
-        }
-
-        pub fn mulBy(
-            self: Self,
-            rhs: anytype,
-        ) QuantityVec(len, Quantity(T, d.add(@TypeOf(rhs).dims), s.min(@TypeOf(rhs).scales))) {
-            const Tr = @TypeOf(rhs);
-            var res: QuantityVec(len, Quantity(T, d.add(Tr.dims), s.min(Tr.scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = (Q{ .value = v }).mulBy(Tr.QuantityType{ .value = rhs.data[i] });
-                res.data[i] = q.value;
-            }
-            return res;
-        }
-
-        pub fn divByScalar(
-            self: Self,
-            scalar: anytype,
-        ) QuantityVec(len, Quantity(T, d.sub(@TypeOf(scalar).dims), s.min(@TypeOf(scalar).scales))) {
-            var res: QuantityVec(len, Quantity(T, d.sub(@TypeOf(scalar).dims), s.min(@TypeOf(scalar).scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = Q{ .value = v };
-                res.data[i] = q.divBy(scalar).value;
-            }
-            return res;
-        }
-
-        pub fn mulByScalar(
-            self: Self,
-            scalar: anytype,
-        ) QuantityVec(len, Quantity(T, d.add(@TypeOf(scalar).dims), s.min(@TypeOf(scalar).scales))) {
-            var res: QuantityVec(len, Quantity(T, d.add(@TypeOf(scalar).dims), s.min(@TypeOf(scalar).scales))) = undefined;
-            for (self.data, 0..) |v, i| {
-                const q = Q{ .value = v };
-                res.data[i] = q.mulBy(scalar).value;
-            }
-            return res;
-        }
-
-        pub fn negate(self: Self) Self {
-            var res: Self = undefined;
-            for (self.data, 0..) |v, i| {
-                res.data[i] = -v;
-            }
-            return res;
-        }
-
-        pub fn scale(self: Self, rhs: T) Self {
-            var res: Self = undefined;
-            for (self.data, 0..) |v, i| {
-                res.data[i] = (Q{ .value = v }).scale(rhs).value;
-            }
-            return res;
-        }
-
-        pub fn to(self: Self, comptime DestQ: type) QuantityVec(len, DestQ) {
-            var res: QuantityVec(len, DestQ) = undefined;
-            for (self.data, 0..) |v, i| {
-                res.data[i] = (Q{ .value = v }).to(DestQ).value;
-            }
-            return res;
-        }
-
-        pub fn lengthSqr(self: Self) T {
-            var sum: T = 0;
-            for (self.data) |v| {
-                sum += v * v;
-            }
-            return sum;
-        }
-
-        pub fn length(self: Self) T {
-            const len_sq = self.lengthSqr();
-
-            if (comptime @typeInfo(T) == .int) {
-                // Construct the unsigned equivalent of T at comptime (e.g., i32 -> u32)
-                const UnsignedT = @Int(.unsigned, @typeInfo(T).int.bits);
-
-                // len_sq is always positive, so @intCast is perfectly safe
-                const u_len_sq = @as(UnsignedT, @intCast(len_sq));
-                return @as(T, @intCast(std.math.sqrt(u_len_sq)));
-            } else {
-                return @sqrt(len_sq);
-            }
-        }
-
-        pub fn format(self: Self, writer: *std.Io.Writer) !void {
-            try writer.writeAll("(");
-            for (self.data, 0..) |v, i| {
-                if (i > 0) try writer.writeAll(", ");
-                try writer.print("{d:.2}", .{v});
-            }
-            try writer.writeAll(")");
-            var iter = std.EnumSet(Dimension).initFull().iterator();
-            var first = true;
-            while (iter.next()) |bu| {
-                const v = dims.get(bu);
-                if (v == 0) continue;
-                if (!first) try writer.writeAll(".");
-                first = false;
-                try writer.print("{s}{s}", .{ scales.get(bu).str(), bu.unit() });
-                if (v != 1) try hlp.printSuperscript(writer, v);
-            }
-        }
-    };
-}
-
-test "Format VectorX" {
-    const MeterPerSecondSq = Quantity(
-        f32,
-        Dimensions.init(.{ .L = 1, .T = -2 }),
-        Scales.init(.{ .T = .n }),
-    );
-    const KgMeterPerSecond = Quantity(
-        f32,
-        Dimensions.init(.{ .M = 1, .L = 1, .T = -1 }),
-        Scales.init(.{ .M = .k }),
-    );
-
-    const accel = MeterPerSecondSq.Vec3.initDefault(9.81);
-    const momentum = KgMeterPerSecond.Vec3{ .data = .{ 43, 0, 11 } };
-
-    std.debug.print("Acceleration: {f}\n", .{accel});
-    std.debug.print("Momentum: {f}\n", .{momentum});
-}
-
-test "VecX Init and Basic Arithmetic" {
-    const Meter = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Vec3M = Meter.Vec3;
-
-    // Test zero, one, initDefault
-    const v_zero = Vec3M.zero;
-    try std.testing.expectEqual(0, v_zero.data[0]);
-    try std.testing.expectEqual(0, v_zero.data[1]);
-    try std.testing.expectEqual(0, v_zero.data[2]);
-
-    const v_one = Vec3M.one;
-    try std.testing.expectEqual(1, v_one.data[0]);
-    try std.testing.expectEqual(1, v_one.data[1]);
-    try std.testing.expectEqual(1, v_one.data[2]);
-
-    const v_def = Vec3M.initDefault(5);
-    try std.testing.expectEqual(5, v_def.data[0]);
-    try std.testing.expectEqual(5, v_def.data[1]);
-    try std.testing.expectEqual(5, v_def.data[2]);
-
-    // Test add and sub
-    const v1 = Vec3M{ .data = .{ 10, 20, 30 } };
-    const v2 = Vec3M{ .data = .{ 2, 4, 6 } };
-
-    const added = v1.add(v2);
-    try std.testing.expectEqual(12, added.data[0]);
-    try std.testing.expectEqual(24, added.data[1]);
-    try std.testing.expectEqual(36, added.data[2]);
-
-    const subbed = v1.sub(v2);
-    try std.testing.expectEqual(8, subbed.data[0]);
-    try std.testing.expectEqual(16, subbed.data[1]);
-    try std.testing.expectEqual(24, subbed.data[2]);
-
-    // Test negate
-    const neg = v1.negate();
-    try std.testing.expectEqual(-10, neg.data[0]);
-    try std.testing.expectEqual(-20, neg.data[1]);
-    try std.testing.expectEqual(-30, neg.data[2]);
-}
-
-test "VecX Kinematics (Scalar Mul/Div)" {
-    const Meter = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Second = Quantity(i32, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-    const Vec3M = Meter.Vec3;
-
-    const pos = Vec3M{ .data = .{ 100, 200, 300 } };
-    const time = Second{ .value = 10 };
-
-    // Vector divided by scalar Quantity (Velocity = Position / Time)
-    const vel = pos.divByScalar(time);
-    try std.testing.expectEqual(10, vel.data[0]);
-    try std.testing.expectEqual(20, vel.data[1]);
-    try std.testing.expectEqual(30, vel.data[2]);
-    try std.testing.expectEqual(1, @TypeOf(vel).dims.get(.L));
-    try std.testing.expectEqual(-1, @TypeOf(vel).dims.get(.T));
-
-    // Vector multiplied by scalar Quantity (Position = Velocity * Time)
-    const new_pos = vel.mulByScalar(time);
-    try std.testing.expectEqual(100, new_pos.data[0]);
-    try std.testing.expectEqual(200, new_pos.data[1]);
-    try std.testing.expectEqual(300, new_pos.data[2]);
-    try std.testing.expectEqual(1, @TypeOf(new_pos).dims.get(.L));
-    try std.testing.expectEqual(0, @TypeOf(new_pos).dims.get(.T));
-}
-
-test "VecX Element-wise Math and Scaling" {
-    const Meter = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const Vec3M = Meter.Vec3;
-
-    const v1 = Vec3M{ .data = .{ 10, 20, 30 } };
-    const v2 = Vec3M{ .data = .{ 2, 5, 10 } };
-
-    // Element-wise division
-    const div = v1.divBy(v2);
-    try std.testing.expectEqual(5, div.data[0]);
-    try std.testing.expectEqual(4, div.data[1]);
-    try std.testing.expectEqual(3, div.data[2]);
-    try std.testing.expectEqual(0, @TypeOf(div).dims.get(.L)); // M / M = Dimensionless
-
-    // Scale by primitive
-    const scaled = v1.scale(2);
-    try std.testing.expectEqual(20, scaled.data[0]);
-    try std.testing.expectEqual(40, scaled.data[1]);
-    try std.testing.expectEqual(60, scaled.data[2]);
-}
-
-test "VecX Conversions" {
-    const KiloMeter = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{ .L = .k }));
-    const Meter = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-
-    const v_km = KiloMeter.Vec3{ .data = .{ 1, 2, 3 } };
-    const v_m = v_km.to(Meter);
-
-    try std.testing.expectEqual(1000, v_m.data[0]);
-    try std.testing.expectEqual(2000, v_m.data[1]);
-    try std.testing.expectEqual(3000, v_m.data[2]);
-
-    // Type checking the result
-    try std.testing.expectEqual(1, @TypeOf(v_m).dims.get(.L));
-    try std.testing.expectEqual(UnitScale.none, @TypeOf(v_m).scales.get(.L));
-}
-
-test "VecX Length" {
-    const MeterInt = Quantity(i32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-    const MeterFloat = Quantity(f32, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-
-    // Integer length (using your custom isqrt)
-    // 3-4-5 triangle on XY plane
-    const v_int = MeterInt.Vec3{ .data = .{ 3, 4, 0 } };
-    try std.testing.expectEqual(25, v_int.lengthSqr());
-    try std.testing.expectEqual(5, v_int.length());
-
-    // Float length
-    const v_float = MeterFloat.Vec3{ .data = .{ 3.0, 4.0, 0.0 } };
-    try std.testing.expectApproxEqAbs(@as(f32, 25.0), v_float.lengthSqr(), 1e-4);
-    try std.testing.expectApproxEqAbs(@as(f32, 5.0), v_float.length(), 1e-4);
-}
-
-test "Benchmark QuantityVec ops" {
-    const Io = std.Io;
-    const ITERS: usize = 10_000;
-    const SAMPLES: usize = 10;
-    var gsink: f64 = 0;
-
-    // In Zig 0.14+, we use the testing IO for clock access in tests
-    const io = std.testing.io;
-
-    const getTime = struct {
-        fn f(i: Io) Io.Timestamp {
-            return Io.Clock.awake.now(i);
-        }
-    }.f;
-
-    const getVal = struct {
-        fn f(comptime TT: type, i: usize, comptime mask: u7) TT {
-            const v: u8 = @as(u8, @truncate(i & @as(usize, mask))) + 1;
-            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
-        }
-    }.f;
-
-    const fold = struct {
-        fn f(comptime TT: type, s: *f64, v: TT) void {
-            s.* += if (comptime @typeInfo(TT) == .float)
-                @as(f64, @floatCast(v))
-            else
-                @as(f64, @floatFromInt(v));
-        }
-    }.f;
-
-    const computeStats = struct {
-        fn f(samples: []f64, iters: usize) f64 {
-            std.mem.sort(f64, samples, {}, std.sort.asc(f64));
-            const mid = samples.len / 2;
-            const median_ns = if (samples.len % 2 == 0)
-                (samples[mid - 1] + samples[mid]) / 2.0
-            else
-                samples[mid];
-            return median_ns / @as(f64, @floatFromInt(iters));
-        }
-    }.f;
-
-    std.debug.print(
-        \\
-        \\ QuantityVec<N, T> benchmark — {d} iterations, {d} samples/cell
-        \\ (Results in ns/op)
-        \\
-        \\┌─────────────┬──────┬─────────┬─────────┬─────────┐
-        \\│ Operation   │ Type │   Len=3 │   Len=4 │  Len=16 │
-        \\├─────────────┼──────┼─────────┼─────────┼─────────┤
-        \\
-    , .{ ITERS, SAMPLES });
-
-    const Types = .{ i32, i64, i128, f32, f64 };
-    const TNames = .{ "i32", "i64", "i128", "f32", "f64" };
-    const Lengths = .{ 3, 4, 16 };
-    const Ops = .{ "add", "scale", "mulByScalar", "length" };
-
-    inline for (Ops, 0..) |op_name, o_idx| {
-        inline for (Types, TNames) |T, tname| {
-            std.debug.print("│ {s:<11} │ {s:<4} │", .{ op_name, tname });
-
-            inline for (Lengths) |len| {
-                const Q_base = Quantity(T, Dimensions.init(.{ .L = 1 }), Scales.init(.{}));
-                const Q_time = Quantity(T, Dimensions.init(.{ .T = 1 }), Scales.init(.{}));
-                const V = QuantityVec(len, Q_base);
-
-                var samples: [SAMPLES]f64 = undefined;
-
-                for (0..SAMPLES) |s_idx| {
-                    var sink: T = 0;
-
-                    const t_start = getTime(io);
-                    for (0..ITERS) |i| {
-                        const v1 = V.initDefault(getVal(T, i, 63));
-
-                        if (comptime std.mem.eql(u8, op_name, "add")) {
-                            const v2 = V.initDefault(getVal(T, i +% 7, 63));
-                            const res = v1.add(v2);
-                            for (res.data) |val| {
-                                if (comptime @typeInfo(T) == .float) sink += val else sink ^= val;
-                            }
-                        } else if (comptime std.mem.eql(u8, op_name, "scale")) {
-                            const sc = getVal(T, i +% 2, 63);
-                            const res = v1.scale(sc);
-                            for (res.data) |val| {
-                                if (comptime @typeInfo(T) == .float) sink += val else sink ^= val;
-                            }
-                        } else if (comptime std.mem.eql(u8, op_name, "mulByScalar")) {
-                            const s_val = Q_time{ .value = getVal(T, i +% 2, 63) };
-                            const res = v1.mulByScalar(s_val);
-                            for (res.data) |val| {
-                                if (comptime @typeInfo(T) == .float) sink += val else sink ^= val;
-                            }
-                        } else if (comptime std.mem.eql(u8, op_name, "length")) {
-                            const r = v1.length();
-                            if (comptime @typeInfo(T) == .float) sink += r else sink ^= r;
-                        }
-                    }
-                    const t_end = getTime(io);
-
-                    samples[s_idx] = @as(f64, @floatFromInt(t_start.durationTo(t_end).toNanoseconds()));
-                    fold(T, &gsink, sink);
-                }
-
-                const median_ns_per_op = computeStats(&samples, ITERS);
-                std.debug.print(" {d:>7.1} │", .{median_ns_per_op});
-            }
-            std.debug.print("\n", .{});
-        }
-
-        if (o_idx < Ops.len - 1) {
-            std.debug.print("├─────────────┼──────┼─────────┼─────────┼─────────┤\n", .{});
-        }
-    }
-    std.debug.print("└─────────────┴──────┴─────────┴─────────┴─────────┘\n", .{});
-    std.debug.print("\nAnti-optimisation sink: {d:.4}\n", .{gsink});
-    try std.testing.expect(gsink != 0);
-}

src/Scales.zig

@@ -1,9 +1,24 @@
 const std = @import("std");
-const hlp = @import("helper.zig");
 const Dimensions = @import("Dimensions.zig");
 const Dimension = @import("Dimensions.zig").Dimension;
 
-pub const UnitScale = enum(i32) {
+// TODO: add more scales like feet and inch
+
+/// Use to initiate Scalar and Scales type
+pub const ArgOpts = struct {
+    L: UnitScale = .none,
+    M: UnitScale = .none,
+    T: UnitScale = .none,
+    I: UnitScale = .none,
+    Tp: UnitScale = .none,
+    N: UnitScale = .none,
+    J: UnitScale = .none,
+};
+
+/// SI prefix (pico…peta) plus time-unit aliases (min, hour, year).
+/// The integer value encodes the exponent for SI prefixes (e.g. `k = 3` → 10³),
+/// and the literal factor for time units (e.g. `hour = 3600`).
+pub const UnitScale = enum(isize) {
     P = 15,
     T = 12,
     G = 9,
@@ -25,6 +40,18 @@ pub const UnitScale = enum(i32) {
     hour = 3_600,
     year = 31_536_000,
 
+    // Imperial Length (Literal factors in meters)
+    // 1 inch = 0.0254 meters. Since enum backing is isize,
+    // we use a unique tag and handle the float in getFactor.
+    inch = -1001,
+    ft = -1002,
+    yd = -1003,
+    mi = -1004,
+
+    oz = -1005, //  1 oz  = 28.3495231 g
+    lb = -1006, //  1 lb  = 453.59237 g  (= 16 oz)
+    st = -1007, //  1 stone = 6350.29318 g (= 14 lb)
+
     // Undefined
     _,
 
@@ -32,65 +59,82 @@ pub const UnitScale = enum(i32) {
         var buf: [16]u8 = undefined;
         return switch (self) {
             .none => "",
-            inline .P, .T, .G, .M, .k, .h, .da, .d, .c, .m, .u, .n, .p, .f, .min, .hour, .year => @tagName(self),
+            .P, .T, .G, .M, .k, .h, .da, .d, .c, .m, .u, .n, .p, .f, .min, .year, .inch, .ft, .yd, .mi, .oz, .lb, .st => @tagName(self),
-            else => std.fmt.bufPrint(&buf, "[{d}]", .{@intFromEnum(self)}) catch "[]",
+            .hour => "h",
+            else => std.fmt.bufPrint(&buf, "[{d}]", .{@intFromEnum(self)}) catch "[]", // This cannot be inline because of non exhaustive enum, but that's ok, it is just str, not calculation
         };
     }
 
-    /// Helper to get the actual scaling factor
+    pub fn getFactor(self: @This()) comptime_float {
-    pub fn getFactor(self: @This()) f64 {
         return switch (self) {
-            inline .P, .T, .G, .M, .k, .h, .da, .none, .d, .c, .m, .u, .n, .p, .f => std.math.pow(f64, 10.0, @floatFromInt(@intFromEnum(self))),
+            // Standard SI Exponents
+            .P, .T, .G, .M, .k, .h, .da, .none, .d, .c, .m, .u, .n, .p, .f => std.math.pow(f64, 10.0, @floatFromInt(@intFromEnum(self))),
+
+            // Time Factors
+            .min, .hour, .year => @floatFromInt(@intFromEnum(self)),
+
+            // Imperial Length (metres)
+            .inch => 0.0254,
+            .ft => 0.3048,
+            .yd => 0.9144,
+            .mi => 1609.344,
+
+            // Imperial Mass (grams — base unit for M is gram, i.e. .none = 1 g)
+            .oz => 28.3495231,
+            .lb => 453.59237,
+            .st => 6350.29318,
+
             else => @floatFromInt(@intFromEnum(self)),
         };
     }
-
-    /// Helper to get the actual scaling factor in i32
-    pub fn getFactorInt(self: @This()) i32 {
-        return switch (self) {
-            inline .P, .T, .G, .M, .k, .h, .da, .none, .d, .c, .m, .u, .n, .p, .f => std.math.powi(i32, 10.0, @intFromEnum(self)) catch 0,
-            else => @intFromEnum(self),
-        };
-    }
 };
 
-const Scales = @This();
+/// Maps each SI base dimension to its `UnitScale`. Stored and resolved entirely at comptime.
+const Self = @This();
 
 data: std.EnumArray(Dimension, UnitScale),
 
-pub fn init(comptime init_val: anytype) Scales {
+/// Create a `Scales` from a struct literal, e.g. `.{ .L = .k, .T = .hour }`.
-    var s = Scales{ .data = std.EnumArray(Dimension, UnitScale).initFill(.none) };
+/// Unspecified dimensions default to `.none` (factor 1).
-    inline for (std.meta.fields(@TypeOf(init_val))) |f| {
+pub fn init(comptime init_val: ArgOpts) Self {
-        if (comptime hlp.isInt(@TypeOf(@field(init_val, f.name))))
+    comptime var s = Self{ .data = std.EnumArray(Dimension, UnitScale).initFill(.none) };
+    for (std.meta.fields(@TypeOf(init_val))) |f| {
+        if (comptime @typeInfo(@TypeOf(@field(init_val, f.name))) == .comptime_int)
             s.data.set(@field(Dimension, f.name), @enumFromInt(@field(init_val, f.name)))
         else
             s.data.set(@field(Dimension, f.name), @field(init_val, f.name));
     }
-    return s;
+    return comptime s;
 }
 
-pub fn initFill(val: UnitScale) Scales {
+pub fn initFill(val: UnitScale) Self {
     return .{ .data = std.EnumArray(Dimension, UnitScale).initFill(val) };
 }
 
-pub fn get(self: Scales, key: Dimension) UnitScale {
+pub fn get(self: Self, key: Dimension) UnitScale {
     return self.data.get(key);
 }
 
-pub fn set(self: *Scales, key: Dimension, val: UnitScale) void {
+pub fn set(self: *Self, key: Dimension, val: UnitScale) void {
     self.data.set(key, val);
 }
 
-pub fn min(comptime s1: Scales, comptime s2: Scales) Scales {
+pub fn eql(self: Self, other: Self) bool {
-    @setEvalBranchQuota(10_000);
+    for (self.data.values, other.data.values) |l, r|
-    var out = Scales.initFill(.none);
+        if (l != r) return false;
-    for (std.enums.values(Dimension)) |dim|
+    return true;
-        out.set(dim, if (s1.get(dim).getFactorInt() > s2.get(dim).getFactorInt()) s2.get(dim) else s1.get(dim));
-
-    return out;
 }
 
-pub fn getFactor(comptime s: Scales, comptime d: Dimensions) f64 {
+pub fn argsOpt(self: Self) ArgOpts {
+    var args: ArgOpts = undefined;
+    for (std.enums.values(Dimension)) |d|
+        @field(args, @tagName(d)) = self.get(d);
+    return args;
+}
+
+/// Compute the combined scale factor for a given dimension signature.
+/// Each dimension's prefix is raised to its exponent and multiplied together.
+pub fn getFactor(s: Self, d: Dimensions) comptime_float {
     var factor: f64 = 1.0;
     for (std.enums.values(Dimension)) |dim| {
         const power = d.get(dim);
@@ -98,7 +142,7 @@ pub fn getFactor(comptime s: Scales, comptime d: Dimensions) f64 {
 
         const base = s.get(dim).getFactor();
 
-        var i: i32 = 0;
+        var i: comptime_int = 0;
         const abs_power = if (power < 0) -power else power;
         while (i < abs_power) : (i += 1) {
             if (power > 0)

src/UnitParser.zig

@@ -0,0 +1,145 @@
+const std = @import("std");
+const Dimensions = @import("Dimensions.zig");
+const Scales = @import("Scales.zig");
+
+/// A container returning the separated arguments needed to construct a Tensor.
+pub const ParsedUnit = struct {
+    dims: Dimensions.ArgOpts = .{},
+    scales: Scales.ArgOpts = .{},
+};
+
+pub const UnitParseError = error{
+    UnknownBaseUnit,
+    UnknownPrefix,
+    InvalidExponent,
+    EmptyStr,
+};
+
+/// Parses strings like "km/s^2", "m", "kg*m/s^2", "1/min".
+/// Evaluates entirely at comptime.
+pub fn parseUnit(str: []const u8) !ParsedUnit {
+    if (str.len == 0) return UnitParseError.EmptyStr;
+
+    var parsed: ParsedUnit = .{ .dims = .{}, .scales = .{} };
+
+    // We need to track if we are after a '/' to flip exponents to negative
+    var is_denominator = false;
+
+    // Manual iteration to handle '/' properly
+    var cursor: usize = 0;
+    while (cursor < str.len) {
+        // Find the next segment
+        const segment_start = cursor;
+        while (cursor < str.len and str[cursor] != '/' and str[cursor] != '.' and str[cursor] != '*') : (cursor += 1) {}
+        const segment = str[segment_start..cursor];
+
+        if (segment.len > 0) {
+            try parseSegment(segment, &parsed, is_denominator);
+        }
+
+        if (cursor < str.len) {
+            if (str[cursor] == '/') {
+                is_denominator = true;
+            }
+            cursor += 1; // skip the separator
+        }
+    }
+
+    return parsed;
+}
+
+fn parseSegment(comptime segment: []const u8, parsed: *ParsedUnit, is_denominator: bool) !void {
+    var scale: Scales.UnitScale = .none;
+    var found_scale = false;
+    var active_dim: ?Dimensions.Dimension = null;
+
+    // 1. Try to find a Scale + Dimension pair (e.g., "mm", "km")
+    inline for (std.enums.values(Scales.UnitScale)) |sca| {
+        const s_str = sca.str();
+        if (s_str.len > 0 and std.mem.startsWith(u8, segment, s_str)) {
+            // Check if it's a "Unit-as-Scale" (hour, min) or a prefix (k, m, c)
+            switch (sca) {
+                .hour, .min, .year => {
+                    // These are dimensions themselves (Time)
+                    if (segment.len == s_str.len or (segment.len > s_str.len and (segment[s_str.len] == '^' or (segment[s_str.len] >= '0' and segment[s_str.len] <= '9')))) {
+                        scale = sca;
+                        active_dim = .T;
+                        found_scale = true;
+                    }
+                },
+                else => {
+                    // Standard prefixes: Must be followed by a valid dimension unit
+                    inline for (std.enums.values(Dimensions.Dimension)) |dim| {
+                        if (std.mem.startsWith(u8, segment[s_str.len..], dim.unit())) {
+                            scale = sca;
+                            active_dim = dim;
+                            found_scale = true;
+                            break;
+                        }
+                    }
+                },
+            }
+        }
+        if (found_scale) break;
+    }
+
+    // 2. If no scale prefix was found, try identifying as a pure Dimension (e.g., "m", "s")
+    if (!found_scale) {
+        inline for (std.enums.values(Dimensions.Dimension)) |dim| {
+            if (std.mem.startsWith(u8, segment, dim.unit())) {
+                active_dim = dim;
+                break;
+            }
+        }
+    }
+
+    const dimen = active_dim orelse return UnitParseError.UnknownBaseUnit;
+
+    // 3. Determine where the exponent starts
+    // If it was a Time Scale (like 'h'), the exponent starts after 'h'
+    // If it was a Prefix + Dim (like 'km'), it starts after 'km'
+    const unit_part_len = if (found_scale)
+        (if (scale == .hour or scale == .min or scale == .year) scale.str().len else scale.str().len + dimen.unit().len)
+    else
+        dimen.unit().len;
+
+    const expo_str = segment[unit_part_len..];
+
+    // 4. Parse Exponent
+    var expo: i32 = 1;
+    if (expo_str.len > 0) {
+        const cleaned_expo = if (expo_str[0] == '^') expo_str[1..] else expo_str;
+        expo = std.fmt.parseInt(i32, cleaned_expo, 10) catch return UnitParseError.InvalidExponent;
+    }
+
+    if (is_denominator) expo *= -1;
+
+    // 5. Assign to struct
+    inline for (std.meta.fields(Dimensions.ArgOpts)) |f| {
+        if (std.mem.eql(u8, f.name, @tagName(dimen))) {
+            @field(parsed.dims, f.name) += expo;
+            @field(parsed.scales, f.name) = scale;
+        }
+    }
+}
+
+inline fn testParser(
+    comptime str: []const u8,
+    comptime expected_dims: Dimensions.ArgOpts,
+    comptime expected_scales: Scales.ArgOpts,
+) !void {
+    const unit = comptime try parseUnit(str);
+    if (comptime !Dimensions.init(expected_dims).eql(Dimensions.init(unit.dims))) return error.WrongDims;
+    if (comptime !Scales.init(expected_scales).eql(Scales.init(unit.scales))) return error.WrongScales;
+}
+
+test "parseUnit" {
+    @setEvalBranchQuota(10000);
+    try testParser("m", .{ .L = 1 }, .{});
+    try testParser("s", .{ .T = 1 }, .{});
+    try testParser("mm", .{ .L = 1 }, .{ .L = .m });
+    try testParser("m/s", .{ .L = 1, .T = -1 }, .{});
+    try testParser("m1/s2/kg", .{ .L = 1, .T = -2, .M = -1 }, .{ .M = .k });
+    try testParser("km/h", .{ .L = 1, .T = -1 }, .{ .L = .k, .T = .hour });
+    try testParser("m.s^-1", .{ .L = 1, .T = -1 }, .{});
+}

src/benchmark.zig

@@ -0,0 +1,669 @@
+const std = @import("std");
+const Io = std.Io;
+const Tensor = @import("Tensor.zig").TensorStatic;
+
+var io: Io = undefined;
+pub fn main(init: std.process.Init) !void {
+    var stdout_buf: [4096]u8 = undefined;
+    var stdout_writer: std.Io.File.Writer = .init(.stdout(), init.io, &stdout_buf);
+    try stdout_writer.interface.print("Starting Benchmarks...", .{});
+
+    io = init.io;
+
+    try vectorSIMDvsNative(f64, &stdout_writer.interface);
+    try stdout_writer.flush();
+    try vectorSIMDvsNative(f32, &stdout_writer.interface);
+    try stdout_writer.flush();
+    try vectorSIMDvsNative(i32, &stdout_writer.interface);
+    try stdout_writer.flush();
+    try vectorSIMDvsNative(i64, &stdout_writer.interface);
+    try stdout_writer.flush();
+    try vectorSIMDvsNative(i128, &stdout_writer.interface);
+    try stdout_writer.flush();
+
+    try bench_Scalar(&stdout_writer.interface);
+    try stdout_writer.flush();
+    try bench_vsNative(&stdout_writer.interface);
+    try stdout_writer.flush();
+    try bench_crossTypeVsNative(&stdout_writer.interface);
+    try stdout_writer.flush();
+    try bench_Vector(&stdout_writer.interface);
+    try stdout_writer.flush();
+    try bench_HighDimTensor(&stdout_writer.interface);
+    try stdout_writer.flush();
+}
+
+fn getTime() Io.Timestamp {
+    return Io.Clock.awake.now(io);
+}
+
+fn fold(comptime TT: type, s: *f64, v: TT) void {
+    s.* += if (comptime @typeInfo(TT) == .float)
+        @as(f64, @floatCast(v))
+    else
+        @as(f64, @floatFromInt(v));
+}
+
+fn bench_Scalar(writer: *std.Io.Writer) !void {
+    const ITERS: usize = 100_000;
+    const SAMPLES: usize = 10;
+
+    const getVal = struct {
+        fn f(comptime TT: type, i: usize, comptime mask: u7) TT {
+            const v: u8 = @as(u8, @truncate(i & @as(usize, mask))) + 1;
+            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
+        }
+    }.f;
+
+    const Stats = struct {
+        median: f64,
+        delta: f64,
+        ops_per_sec: f64,
+    };
+
+    const computeStats = struct {
+        fn f(samples: []f64, iters: usize) Stats {
+            std.mem.sort(f64, samples, {}, std.sort.asc(f64));
+            const mid = samples.len / 2;
+            const median_ns = if (samples.len % 2 == 0) (samples[mid - 1] + samples[mid]) / 2.0 else samples[mid];
+
+            const low = samples[0];
+            const high = samples[samples.len - 1];
+            const delta_ns = (high - low) / 2.0;
+
+            const ns_per_op = median_ns / @as(f64, @floatFromInt(iters));
+            return .{
+                .median = ns_per_op,
+                .delta = (delta_ns / @as(f64, @floatFromInt(iters))),
+                .ops_per_sec = 1_000_000_000.0 / ns_per_op,
+            };
+        }
+    }.f;
+
+    try writer.print(
+        \\
+        \\ Scalar<T> benchmark — {d} iterations, {d} samples/cell
+        \\
+        \\┌───────────────────┬──────┬─────────────────────┬─────────────────────┐
+        \\│ Operation         │ Type │ ns / op (± delta)   │ Throughput (ops/s)  │
+        \\├───────────────────┼──────┼─────────────────────┼─────────────────────┤
+        \\
+    , .{ ITERS, SAMPLES });
+
+    const Types = .{ i16, i32, i64, i128, i256, f32, f64 };
+    const TNames = .{ "i16", "i32", "i64", "i128", "i256", "f32", "f64" };
+    const Ops = .{ "add", "sub", "mul", "div", "to", "abs", "pow", "eq", "gt", "mul(n)" };
+
+    var results_matrix: [Ops.len][Types.len]f64 = undefined;
+
+    comptime var tidx: usize = 0;
+    inline for (Types, TNames) |T, tname| {
+        const M = Tensor(T, .{ .L = 1 }, .{}, &.{1});
+        const KM = Tensor(T, .{ .L = 1 }, .{ .L = .k }, &.{1});
+        const S = Tensor(T, .{ .T = 1 }, .{}, &.{1});
+
+        inline for (Ops, 0..) |op_name, oidx| {
+            var samples: [SAMPLES]f64 = undefined;
+
+            for (0..SAMPLES) |s_idx| {
+                const t_start = getTime();
+
+                for (0..ITERS) |i| {
+                    std.mem.doNotOptimizeAway(
+                        {
+                            _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                                (M.splat(getVal(T, i, 63))).add(M.splat(getVal(T, i +% 7, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "sub"))
+                                (M.splat(getVal(T, i +% 10, 63))).sub(M.splat(getVal(T, i, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "mul"))
+                                (M.splat(getVal(T, i, 63))).mul(M.splat(getVal(T, i +% 1, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "div"))
+                                (M.splat(getVal(T, i +% 10, 63))).div(S.splat(getVal(T, i, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "to"))
+                                (KM.splat(getVal(T, i, 15))).to(M)
+                            else if (comptime std.mem.eql(u8, op_name, "abs"))
+                                (M.splat(getVal(T, i, 63))).abs()
+                            else if (comptime std.mem.eql(u8, op_name, "eq"))
+                                (M.splat(getVal(T, i, 63))).eq(M.splat(getVal(T, i +% 3, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "gt"))
+                                (M.splat(getVal(T, i, 63))).gt(M.splat(getVal(T, i +% 3, 63)))
+                            else
+                                (M.splat(getVal(T, i, 63))).mul(M.splat(3));
+                        },
+                    );
+                }
+
+                const t_end = getTime();
+                samples[s_idx] = @as(f64, @floatFromInt(t_start.durationTo(t_end).toNanoseconds()));
+            }
+
+            const stats = computeStats(&samples, ITERS);
+            results_matrix[oidx][tidx] = stats.median;
+
+            try writer.print("│ {s:<17} │ {s:<4} │ {d:>8.2} ns ±{d:<6.2} │ {d:>19.0} │\n", .{ op_name, tname, stats.median, stats.delta, stats.ops_per_sec });
+        }
+
+        if (comptime tidx < Types.len - 1) {
+            try writer.print("├───────────────────┼──────┼─────────────────────┼─────────────────────┤\n", .{});
+        }
+        tidx += 1;
+    }
+
+    // Median Summary Table
+    try writer.print("└───────────────────┴──────┴─────────────────────┴─────────────────────┘\n\n", .{});
+    try writer.print("Median Summary (ns/op):\n", .{});
+
+    try writer.print("┌──────────────┬───────┬───────┬───────┬───────┬───────┬───────┬───────┐\n", .{});
+    try writer.print("│  Operation   │  i16  │  i32  │  i64  │  i128 │  i256 │  f32  │  f64  │\n", .{});
+    try writer.print("├──────────────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┤\n", .{});
+
+    inline for (Ops, 0..) |op_name, oidx| {
+        try writer.print("│  {s:<11} │", .{op_name});
+        var i: usize = 0;
+        while (i < Types.len) : (i += 1)
+            try writer.print("{d:>6.1} │", .{results_matrix[oidx][i]});
+
+        try writer.print("\n", .{});
+    }
+
+    try writer.print("└──────────────┴───────┴───────┴───────┴───────┴───────┴───────┴───────┘\n", .{});
+}
+
+fn bench_vsNative(writer: *std.Io.Writer) !void {
+    const ITERS: usize = 100_000;
+    const SAMPLES: usize = 100;
+
+    const getValT = struct {
+        fn f(comptime TT: type, i: usize) TT {
+            const v = (i % 100) + 1;
+            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
+        }
+    }.f;
+
+    const Types = .{ i32, i64, i128, f32, f64 };
+    const TNames = .{ "i32", "i64", "i128", "f32", "f64" };
+    // Expanded Ops to match bench_Scalar
+    const Ops = .{ "add", "sub", "mul", "div", "abs", "eq", "gt" };
+
+    try writer.print(
+        \\
+        \\ Scalar vs Native Overhead Analysis
+        \\
+        \\┌───────────┬──────┬───────────┬───────────┬───────────┬───────────────────────┐
+        \\│ Operation │ Type │ Native    │ @Vector   │ Tensor{{1}} │ Slowdown  Nat | Vec   │
+        \\├───────────┼──────┼───────────┼───────────┼───────────┼───────────────────────┤
+        \\
+    , .{});
+
+    inline for (Ops, 0..) |op_name, j| {
+        inline for (Types, 0..) |T, tidx| {
+            var native_total_ns: f64 = 0;
+            var vector_total_ns: f64 = 0;
+            var tensor_total_ns: f64 = 0;
+
+            const M = Tensor(T, .{}, .{}, &.{1});
+
+            for (0..SAMPLES) |_| {
+                // --- 1. Benchmark Native ---
+                const n_start = getTime();
+                const a = getValT(T, 10);
+                const b = getValT(T, 2);
+                for (0..ITERS) |_| {
+                    // Native logic branch
+                    _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                        if (comptime @typeInfo(T) == .int) a +| b else a + b
+                    else if (comptime std.mem.eql(u8, op_name, "sub"))
+                        if (comptime @typeInfo(T) == .int) a -| b else a - b
+                    else if (comptime std.mem.eql(u8, op_name, "mul"))
+                        if (comptime @typeInfo(T) == .int) a *| b else a * b
+                    else if (comptime std.mem.eql(u8, op_name, "div"))
+                        if (comptime @typeInfo(T) == .int) @divTrunc(a, b) else a / b
+                    else if (comptime std.mem.eql(u8, op_name, "abs"))
+                        if (comptime @typeInfo(T) == .int) @abs(a) else @as(T, @abs(a))
+                    else if (comptime std.mem.eql(u8, op_name, "eq"))
+                        a == b
+                    else if (comptime std.mem.eql(u8, op_name, "gt"))
+                        a > b
+                    else
+                        unreachable;
+                }
+                const n_end = getTime();
+                native_total_ns += @as(f64, @floatFromInt(n_start.durationTo(n_end).toNanoseconds()));
+
+                const v_start = getTime();
+                const va = @Vector(1, T){getValT(T, 10)};
+                const vb = @Vector(1, T){getValT(T, 2)};
+                for (0..ITERS) |_| {
+                    // Native logic branch
+                    _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                        if (comptime @typeInfo(T) == .int) va +| vb else va + vb
+                    else if (comptime std.mem.eql(u8, op_name, "sub"))
+                        if (comptime @typeInfo(T) == .int) va -| vb else va - vb
+                    else if (comptime std.mem.eql(u8, op_name, "mul"))
+                        if (comptime @typeInfo(T) == .int) va *| vb else va * vb
+                    else if (comptime std.mem.eql(u8, op_name, "div"))
+                        if (comptime @typeInfo(T) == .int) @divTrunc(va, vb) else va / vb
+                    else if (comptime std.mem.eql(u8, op_name, "abs"))
+                        if (comptime @typeInfo(T) == .int) @as(T, @intCast(@abs(va[0]))) else @abs(va)
+                    else if (comptime std.mem.eql(u8, op_name, "eq"))
+                        va == vb
+                    else if (comptime std.mem.eql(u8, op_name, "gt"))
+                        va > vb
+                    else
+                        unreachable;
+                }
+                const v_end = getTime();
+                vector_total_ns += @as(f64, @floatFromInt(v_start.durationTo(v_end).toNanoseconds()));
+
+                // --- 2. Benchmark Scalar ---
+                const q_start = getTime();
+                const qa = M.splat(getValT(T, 10));
+                const qb = M.splat(getValT(T, 2));
+                for (0..ITERS) |_| {
+                    // Scalar logic branch
+                    _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                        qa.add(qb)
+                    else if (comptime std.mem.eql(u8, op_name, "sub"))
+                        qa.sub(qb)
+                    else if (comptime std.mem.eql(u8, op_name, "mul"))
+                        qa.mul(qb)
+                    else if (comptime std.mem.eql(u8, op_name, "div"))
+                        qa.div(qb)
+                    else if (comptime std.mem.eql(u8, op_name, "abs"))
+                        qa.abs()
+                    else if (comptime std.mem.eql(u8, op_name, "eq"))
+                        qa.eq(qb)
+                    else if (comptime std.mem.eql(u8, op_name, "gt"))
+                        qa.gt(qb)
+                    else
+                        unreachable;
+                }
+                const q_end = getTime();
+                tensor_total_ns += @as(f64, @floatFromInt(q_start.durationTo(q_end).toNanoseconds()));
+            }
+
+            const avg_n = (native_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
+            const avg_v = (vector_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
+            const avg_t = (tensor_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
+            const slowdown_nt = avg_t / avg_n;
+            const slowdown_vt = avg_t / avg_v;
+
+            try writer.print("│ {s:<9} │ {s:<4} │ {d:>7.2}ns │ {d:>7.2}ns │ {d:>7.2}ns │ {d:>8.2}x   {d:>8.2}x │\n", .{
+                op_name, TNames[tidx], avg_n, avg_v, avg_t, slowdown_nt, slowdown_vt,
+            });
+        }
+        if (j != Ops.len - 1) try writer.print("├───────────┼──────┼───────────┼───────────┼───────────┼───────────────────────┤\n", .{});
+    }
+
+    try writer.print("└───────────┴──────┴───────────┴───────────┴───────────┴───────────────────────┘\n", .{});
+}
+
+fn bench_crossTypeVsNative(writer: *std.Io.Writer) !void {
+    const ITERS: usize = 100_000;
+    const SAMPLES: usize = 5;
+
+    const getValT = struct {
+        fn f(comptime TT: type, i: usize) TT {
+            // Keep values safe and non-zero to avoid division by zero or overflows during cross-casting
+            const v = (i % 50) + 1;
+            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
+        }
+    }.f;
+
+    // Helper for the Native baseline: explicitly casting T2 to T1 before the operation
+    const castTo = struct {
+        fn f(comptime DestT: type, comptime SrcT: type, val: SrcT) DestT {
+            if (comptime DestT == SrcT) return val;
+            const src_info = @typeInfo(SrcT);
+            const dest_info = @typeInfo(DestT);
+
+            if (dest_info == .int and src_info == .int) return @intCast(val);
+            if (dest_info == .float and src_info == .int) return @floatFromInt(val);
+            if (dest_info == .int and src_info == .float) return @intFromFloat(val);
+            if (dest_info == .float and src_info == .float) return @floatCast(val);
+            unreachable;
+        }
+    }.f;
+
+    const Types = .{ i16, i64, i128, f32, f64 };
+    const TNames = .{ "i16", "i64", "i128", "f32", "f64" };
+    const Ops = .{ "add", "mul", "div" };
+
+    try writer.print(
+        \\
+        \\ Cross-Type Overhead Analysis: Scalar vs Native
+        \\
+        \\┌─────────┬──────┬──────┬───────────┬───────────┬───────────┐
+        \\│ Op      │ T1   │ T2   │ Native    │ Scalar    │ Slowdown  │
+        \\├─────────┼──────┼──────┼───────────┼───────────┼───────────┤
+        \\
+    , .{});
+
+    inline for (Ops, 0..) |op_name, j| {
+        inline for (Types, 0..) |T1, t1_idx| {
+            inline for (Types, 0..) |T2, t2_idx| {
+                var native_total_ns: f64 = 0;
+                var quantity_total_ns: f64 = 0;
+
+                const M1 = Tensor(T1, .{ .L = 1 }, .{}, &.{1});
+                const M2 = Tensor(T2, .{ .L = 1 }, .{}, &.{1});
+                const S2 = Tensor(T2, .{ .T = 1 }, .{}, &.{1});
+
+                std.mem.doNotOptimizeAway({
+                    for (0..SAMPLES) |_| {
+                        // --- 1. Benchmark Native (Cast T2 to T1, then math) ---
+                        const n_start = getTime();
+                        for (0..ITERS) |i| {
+                            const a = getValT(T1, i);
+                            const b_raw = getValT(T2, 2);
+                            const b = castTo(T1, T2, b_raw);
+
+                            _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                                a + b
+                            else if (comptime std.mem.eql(u8, op_name, "mul"))
+                                a * b
+                            else if (comptime @typeInfo(T1) == .int)
+                                @divTrunc(a, b)
+                            else
+                                a / b;
+                        }
+                        const n_end = getTime();
+                        native_total_ns += @as(f64, @floatFromInt(n_start.durationTo(n_end).toNanoseconds()));
+
+                        // --- 2. Benchmark Scalar ---
+                        const q_start = getTime();
+                        for (0..ITERS) |i| {
+                            const qa = M1.splat(getValT(T1, i));
+                            const qb = if (comptime std.mem.eql(u8, op_name, "div"))
+                                S2.splat(getValT(T2, 2))
+                            else
+                                M2.splat(getValT(T2, 2));
+
+                            _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                                qa.add(qb)
+                            else if (comptime std.mem.eql(u8, op_name, "mul"))
+                                qa.mul(qb)
+                            else
+                                qa.div(qb);
+                        }
+                        const q_end = getTime();
+                        quantity_total_ns += @as(f64, @floatFromInt(q_start.durationTo(q_end).toNanoseconds()));
+                    }
+
+                    const avg_n = (native_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
+                    const avg_q = (quantity_total_ns / SAMPLES) / @as(f64, @floatFromInt(ITERS));
+                    const slowdown = avg_q / avg_n;
+
+                    try writer.print("│ {s:<7} │ {s:<4} │ {s:<4} │ {d:>7.2}ns │ {d:>7.2}ns │ {d:>8.2}x │\n", .{
+                        op_name, TNames[t1_idx], TNames[t2_idx], avg_n, avg_q, slowdown,
+                    });
+                });
+            }
+        }
+        if (j != Ops.len - 1) {
+            try writer.print("├─────────┼──────┼──────┼───────────┼───────────┼───────────┤\n", .{});
+        }
+    }
+
+    try writer.print("└─────────┴──────┴──────┴───────────┴───────────┴───────────┘\n", .{});
+}
+
+fn bench_Vector(writer: *std.Io.Writer) !void {
+    const ITERS: usize = 10_000;
+    const SAMPLES: usize = 10;
+
+    const getVal = struct {
+        fn f(comptime TT: type, i: usize, comptime mask: u7) TT {
+            const v: u8 = @as(u8, @truncate(i & @as(usize, mask))) + 1;
+            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
+        }
+    }.f;
+
+    const computeStats = struct {
+        fn f(samples: []f64, iters: usize) f64 {
+            std.mem.sort(f64, samples, {}, std.sort.asc(f64));
+            const mid = samples.len / 2;
+            const median_ns = if (samples.len % 2 == 0)
+                (samples[mid - 1] + samples[mid]) / 2.0
+            else
+                samples[mid];
+            return median_ns / @as(f64, @floatFromInt(iters));
+        }
+    }.f;
+
+    try writer.print(
+        \\
+        \\ Vector<N, T> benchmark — {d} iterations, {d} samples/cell
+        \\ (Results in ns/op; "---" = not applicable for this length)
+        \\
+        \\┌──────────────────┬──────┬─────────┬─────────┬─────────┬─────────┬─────────┐
+        \\│ Operation        │ Type │   Len=1 │   Len=3 │   Len=4 │  Len=16 │ Len=100 │
+        \\├──────────────────┼──────┼─────────┼─────────┼─────────┼─────────┼─────────┤
+        \\
+    , .{ ITERS, SAMPLES });
+
+    const Types = .{ i32, i64, i128, f32, f64 };
+    const TNames = .{ "i32", "i64", "i128", "f32", "f64" };
+    const Lengths = .{ 1, 3, 4, 16, 100 };
+    // "cross" is only valid for len=3; other cells will show "  ---  "
+    const Ops = .{ "add", "div", "mulScalar", "dot", "cross", "product", "pow", "length" };
+
+    inline for (Ops, 0..) |op_name, o_idx| {
+        inline for (Types, TNames) |T, tname| {
+            try writer.print("│ {s:<16} │ {s:<4} │", .{ op_name, tname });
+
+            inline for (Lengths) |len| {
+                const Q_time = Tensor(T, .{ .T = 1 }, .{}, &.{1});
+                const V = Tensor(T, .{ .L = 1 }, .{}, &.{len});
+
+                // cross product is only defined for len == 3
+                const is_cross = comptime std.mem.eql(u8, op_name, "cross");
+                if (comptime is_cross and len != 3) {
+                    try writer.print("     --- │", .{});
+                    continue;
+                }
+
+                var samples: [SAMPLES]f64 = undefined;
+
+                std.mem.doNotOptimizeAway({
+                    for (0..SAMPLES) |s_idx| {
+                        const t_start = getTime();
+                        for (0..ITERS) |i| {
+                            const v1 = V.splat(getVal(T, i, 63));
+
+                            if (comptime std.mem.eql(u8, op_name, "add")) {
+                                const v2 = V.splat(getVal(T, i +% 7, 63));
+                                _ = v1.add(v2);
+                            } else if (comptime std.mem.eql(u8, op_name, "div")) {
+                                _ = v1.div(V.splat(getVal(T, i +% 2, 63)));
+                            } else if (comptime std.mem.eql(u8, op_name, "mulScalar")) {
+                                const s_val = Q_time.splat(getVal(T, i +% 2, 63));
+                                _ = v1.mul(s_val);
+                            } else if (comptime std.mem.eql(u8, op_name, "dot")) {
+                                const v2 = V.splat(getVal(T, i +% 5, 63));
+                                _ = v1.contract(v2, 0, 0);
+                            } else if (comptime std.mem.eql(u8, op_name, "cross")) {
+                                // len == 3 guaranteed by the guard above
+                                const v2 = V.splat(getVal(T, i +% 5, 63));
+                                _ = v1.cross(v2);
+                            } else if (comptime std.mem.eql(u8, op_name, "product")) {
+                                _ = v1.product();
+                            } else if (comptime std.mem.eql(u8, op_name, "pow")) {
+                                _ = v1.pow(2);
+                            } else if (comptime std.mem.eql(u8, op_name, "length")) {
+                                _ = v1.length();
+                            }
+                        }
+                        const t_end = getTime();
+
+                        samples[s_idx] = @as(f64, @floatFromInt(t_start.durationTo(t_end).toNanoseconds()));
+                    }
+
+                    const median_ns_per_op = computeStats(&samples, ITERS);
+                    try writer.print(" {d:>7.1} │", .{median_ns_per_op});
+                });
+            }
+            try writer.print("\n", .{});
+        }
+
+        if (o_idx < Ops.len - 1) {
+            try writer.print("├──────────────────┼──────┼─────────┼─────────┼─────────┼─────────┼─────────┤\n", .{});
+        }
+    }
+    try writer.print("└──────────────────┴──────┴─────────┴─────────┴─────────┴─────────┴─────────┘\n", .{});
+}
+
+fn bench_HighDimTensor(writer: *std.Io.Writer) !void {
+    const ITERS: usize = 5_000;
+    const SAMPLES: usize = 5;
+
+    const getVal = struct {
+        fn f(comptime TT: type, i: usize, comptime mask: u7) TT {
+            const v: u8 = @as(u8, @truncate(i & @as(usize, mask))) + 1;
+            return if (comptime @typeInfo(TT) == .float) @floatFromInt(v) else @intCast(v);
+        }
+    }.f;
+
+    const computeStats = struct {
+        fn f(samples: []f64, iters: usize) f64 {
+            std.mem.sort(f64, samples, {}, std.sort.asc(f64));
+            const mid = samples.len / 2;
+            const median_ns = if (samples.len % 2 == 0)
+                (samples[mid - 1] + samples[mid]) / 2.0
+            else
+                samples[mid];
+            return median_ns / @as(f64, @floatFromInt(iters));
+        }
+    }.f;
+
+    try writer.print(
+        \\
+        \\ High Dimension Tensor benchmark — {d} iterations, {d} samples/cell
+        \\ (Results in ns/op)
+        \\
+        \\┌─────────────────┬──────┬──────────────┬──────────────┬──────────────┬──────────────┐
+        \\│ Operation       │ Type │        2x2x2 │        3x3x3 │        4x4x4 │  10x10x10x10 │
+        \\├─────────────────┼──────┼──────────────┼──────────────┼──────────────┼──────────────┤
+        \\
+    , .{ ITERS, SAMPLES });
+
+    const Types = .{ i32, i64, f32, f64 };
+    const TNames = .{ "i32", "i64", "f32", "f64" };
+
+    // Testing multiple structural bounds
+    const Shapes = .{
+        &.{ 2, 2, 2 },
+        &.{ 3, 3, 3 },
+        &.{ 4, 4, 4 },
+        &.{ 10, 10, 10, 10 },
+    };
+
+    const Ops = .{ "add", "sub", "mulElem", "mulScalar", "abs" };
+
+    inline for (Ops, 0..) |op_name, o_idx| {
+        inline for (Types, TNames) |T, tname| {
+            try writer.print("│ {s:<15} │ {s:<4} │", .{ op_name, tname });
+
+            inline for (Shapes) |shape| {
+                const V = Tensor(T, .{ .L = 1 }, .{}, shape);
+                const Q = Tensor(T, .{ .T = 1 }, .{}, &.{1}); // For scalar broadcasting operations
+
+                var samples: [SAMPLES]f64 = undefined;
+
+                for (0..SAMPLES) |s_idx| {
+                    const t_start = getTime();
+
+                    for (0..ITERS) |i| {
+                        std.mem.doNotOptimizeAway({
+                            const t1 = V.splat(getVal(T, i, 63));
+
+                            _ = if (comptime std.mem.eql(u8, op_name, "add"))
+                                t1.add(V.splat(getVal(T, i +% 7, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "sub"))
+                                t1.sub(V.splat(getVal(T, i +% 3, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "mulElem"))
+                                t1.mul(V.splat(getVal(T, i +% 5, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "mulScalar"))
+                                t1.mul(Q.splat(getVal(T, i +% 2, 63)))
+                            else if (comptime std.mem.eql(u8, op_name, "abs"))
+                                t1.abs()
+                            else
+                                unreachable;
+                        });
+                    }
+
+                    const t_end = getTime();
+                    samples[s_idx] = @as(f64, @floatFromInt(t_start.durationTo(t_end).toNanoseconds()));
+                }
+
+                const median_ns_per_op = computeStats(&samples, ITERS);
+                try writer.print(" {d:>12.1} │", .{median_ns_per_op});
+            }
+            try writer.print("\n", .{});
+        }
+
+        if (o_idx < Ops.len - 1) {
+            try writer.print("├─────────────────┼──────┼──────────────┼──────────────┼──────────────┼──────────────┤\n", .{});
+        }
+    }
+    try writer.print("└─────────────────┴──────┴──────────────┴──────────────┴──────────────┴──────────────┘\n", .{});
+}
+
+fn vectorSIMDvsNative(comptime T: type, writer: *std.Io.Writer) !void {
+    const iterations: u64 = 10_000;
+    const lens = [_]u32{ 1, 2, 3, 4, 5, 10, 100, 1_000, 10_000 };
+
+    try writer.print("\nSIMD Speedup Analysis: {s}\n", .{@typeName(T)});
+    try writer.print("┌────────────┬────────────┬────────────┬────────────┐\n", .{});
+    try writer.print("│ Vector Len │ Scalar (us)│ Vector (us)│ Speedup    │\n", .{});
+    try writer.print("├────────────┼────────────┼────────────┼────────────┤\n", .{});
+
+    inline for (lens) |vector_len| {
+        // --- Scalar Test ---
+        var scalar_val: T = 10;
+        const start_scalar = getTime();
+
+        var i: u64 = 0;
+        while (i < iterations * vector_len) : (i += 1) {
+            if (comptime @typeInfo(T) == .int)
+                scalar_val = scalar_val +% 1
+            else
+                scalar_val = scalar_val + 1;
+        }
+        const scalar_time = start_scalar.durationTo(getTime()).toMicroseconds();
+
+        // --- Vector Test ---
+        var vector_val: @Vector(vector_len, T) = @splat(20);
+        const start_vector = getTime();
+
+        i = 0;
+        const increment: @Vector(vector_len, T) = @splat(1);
+        while (i < iterations) : (i += 1) {
+            if (comptime @typeInfo(T) == .int)
+                vector_val = vector_val +% increment
+            else
+                vector_val = vector_val + increment;
+        }
+        const vector_time = start_vector.durationTo(getTime()).toMicroseconds();
+
+        // --- Results ---
+        const s_float = @as(f64, @floatFromInt(scalar_time));
+        const v_float = @as(f64, @floatFromInt(vector_time));
+
+        // Speedup = ScalarTime / VectorTime.
+        // > 1.0 means SIMD is faster.
+        const speedup = if (vector_time > 0) s_float / v_float else 0;
+
+        try writer.print("│ {d:<10} │ {d:>10} │ {d:>10} │ {d:>9.2}x │\n", .{
+            vector_len,
+            scalar_time,
+            vector_time,
+            speedup,
+        });
+        try writer.flush();
+
+        std.mem.doNotOptimizeAway(scalar_val);
+        std.mem.doNotOptimizeAway(vector_val);
+    }
+    try writer.print("└────────────┴────────────┴────────────┴────────────┘\n", .{});
+}

src/helper.zig

@@ -1,32 +0,0 @@
-const std = @import("std");
-
-pub fn isInt(comptime T: type) bool {
-    return @typeInfo(T) == .int or @typeInfo(T) == .comptime_int;
-}
-
-pub fn printSuperscript(writer: *std.Io.Writer, n: i32) !void {
-    if (n == 0) return;
-    var val = n;
-    if (val < 0) {
-        try writer.writeAll("\u{207B}");
-        val = -val;
-    }
-    var buf: [12]u8 = undefined;
-    const str = std.fmt.bufPrint(&buf, "{d}", .{val}) catch return;
-    for (str) |c| {
-        const s = switch (c) {
-            '0' => "\u{2070}",
-            '1' => "\u{00B9}",
-            '2' => "\u{00B2}",
-            '3' => "\u{00B3}",
-            '4' => "\u{2074}",
-            '5' => "\u{2075}",
-            '6' => "\u{2076}",
-            '7' => "\u{2077}",
-            '8' => "\u{2078}",
-            '9' => "\u{2079}",
-            else => unreachable,
-        };
-        try writer.writeAll(s);
-    }
-}

src/main.zig

@@ -1,16 +1,15 @@
 const std = @import("std");
 
-pub const Quantity = @import("Quantity.zig").Quantity;
+pub const Tensor = @import("Tensor.zig").Tensor;
-pub const QuantityVec = @import("QuantityVec.zig").QuantityVec;
 pub const Dimensions = @import("Dimensions.zig");
 pub const Scales = @import("Scales.zig");
-pub const Base = @import("BaseQuantities.zig");
+pub const Base = @import("Base.zig");
+pub const UnitParser = @import("UnitParser.zig");
 
 test {
-    _ = @import("Quantity.zig");
+    _ = @import("Tensor.zig");
-    _ = @import("QuantityVec.zig");
     _ = @import("Dimensions.zig");
     _ = @import("Scales.zig");
-    _ = @import("BaseQuantities.zig");
+    _ = @import("Base.zig");
-    _ = @import("helper.zig");
+    _ = @import("UnitParser.zig");
 }

src/shared.zig

@@ -0,0 +1,149 @@
+const std = @import("std");
+const Scales = @import("Scales.zig");
+const UnitScale = Scales.UnitScale;
+const Dimensions = @import("Dimensions.zig");
+const Dimension = Dimensions.Dimension;
+
+pub fn shapeTotal(shape: []const comptime_int) usize {
+    var t: comptime_int = 1;
+    for (shape) |s| t *= s;
+    return t;
+}
+
+/// Check if two shapes are strictly identical.
+pub fn shapeEql(a: []const comptime_int, b: []const comptime_int) bool {
+    if (a.len != b.len) return false;
+    for (a, 0..) |v, i|
+        if (v != b[i]) return false;
+    return true;
+}
+
+/// Row-major (C-order) strides: strides[i] = product(shape[i+1..]).
+///   e.g. shape {3, 4} → strides {4, 1}
+///        shape {2, 3, 4} → strides {12, 4, 1}
+pub fn shapeStrides(shape: []const comptime_int) [shape.len]comptime_int {
+    var st: [shape.len]comptime_int = undefined;
+    if (shape.len == 0) return st;
+    st[shape.len - 1] = 1;
+    if (shape.len > 1) {
+        var i: comptime_int = shape.len - 1;
+        while (i > 0) : (i -= 1) st[i - 1] = st[i] * shape[i];
+    }
+    return st;
+}
+
+/// Return a copy of `shape` with the element at `axis` removed.
+pub fn shapeRemoveAxis(shape: []const comptime_int, axis: comptime_int) [shape.len - 1]comptime_int {
+    var out: [shape.len - 1]comptime_int = undefined;
+    var j: comptime_int = 0;
+    for (shape, 0..) |v, i| {
+        if (i != axis) {
+            out[j] = v;
+            j += 1;
+        }
+    }
+    return out;
+}
+
+/// Concatenate two compile-time slices.
+pub fn shapeCat(a: []const comptime_int, b: []const comptime_int) [a.len + b.len]comptime_int {
+    var out: [a.len + b.len]comptime_int = undefined;
+    for (a, 0..) |v, i| out[i] = v;
+    for (b, 0..) |v, i| out[a.len + i] = v;
+    return out;
+}
+
+/// Decode a flat row-major index into N-D coordinates.
+/// Called only in comptime contexts (all arguments are comptime).
+pub fn decodeFlatCoords(flat: comptime_int, n: comptime_int, strd: [n]comptime_int) [n]usize {
+    var coords: [n]comptime_int = undefined;
+    var tmp = flat;
+    for (0..n) |i| {
+        coords[i] = if (strd[i] == 0) 0 else tmp / strd[i];
+        tmp = if (strd[i] == 0) 0 else tmp % strd[i];
+    }
+    return coords;
+}
+
+/// Encode N-D coordinates into a flat row-major index.
+/// Called only in comptime contexts.
+pub fn encodeFlatCoords(coords: []const usize, n: usize, strd: [n]usize) usize {
+    var flat: usize = 0;
+    for (0..n) |i| flat += coords[i] * strd[i];
+    return flat;
+}
+
+/// Rebuild a full coordinate array by inserting `val` at `axis` into `free`.
+/// `free` holds the remaining (non-contracted) coordinates in order.
+pub fn insertAxis(
+    comptime n: usize,
+    comptime axis: usize,
+    comptime val: usize,
+    comptime free: []const usize,
+) [n]usize {
+    var out: [n]usize = undefined;
+    var fi: usize = 0;
+    for (0..n) |i| {
+        if (i == axis) {
+            out[i] = val;
+        } else {
+            out[i] = free[fi];
+            fi += 1;
+        }
+    }
+    return out;
+}
+
+pub inline fn isInt(comptime T: type) bool {
+    return @typeInfo(T) == .int or @typeInfo(T) == .comptime_int;
+}
+
+pub fn finerScales(comptime T1: type, comptime T2: type) Scales {
+    const d1: Dimensions = T1.dims;
+    const d2: Dimensions = T2.dims;
+    const s1: Scales = T1.scales;
+    const s2: Scales = T2.scales;
+    comptime var out = Scales.initFill(.none);
+    for (std.enums.values(Dimension)) |dim| {
+        const scale1 = comptime s1.get(dim);
+        const scale2 = comptime s2.get(dim);
+        out.set(dim, if (comptime d1.get(dim) == 0 and d2.get(dim) == 0)
+            .none
+        else if (comptime d1.get(dim) == 0)
+            scale2
+        else if (comptime d2.get(dim) == 0)
+            scale1
+        else if (comptime scale1.getFactor() > scale2.getFactor())
+            scale2
+        else
+            scale1);
+    }
+    return out;
+}
+
+pub fn printSuperscript(writer: *std.Io.Writer, n: i32) !void {
+    if (n == 0) return;
+    var val = n;
+    if (val < 0) {
+        try writer.writeAll("\u{207B}");
+        val = -val;
+    }
+    var buf: [12]u8 = undefined;
+    const str = std.fmt.bufPrint(&buf, "{d}", .{val}) catch return;
+    for (str) |c| {
+        const s = switch (c) {
+            '0' => "\u{2070}",
+            '1' => "\u{00B9}",
+            '2' => "\u{00B2}",
+            '3' => "\u{00B3}",
+            '4' => "\u{2074}",
+            '5' => "\u{2075}",
+            '6' => "\u{2076}",
+            '7' => "\u{2077}",
+            '8' => "\u{2078}",
+            '9' => "\u{2079}",
+            else => unreachable,
+        };
+        try writer.writeAll(s);
+    }
+}

tmp.md

@@ -1,134 +0,0 @@
-The slowdown you are seeing (1.5x to 2.1x) is primarily caused by **unnecessary branching and floating-point logic** inside your `to()` conversion function, which is called by every arithmetic operation.
-
-Even though your `ratio` is calculated at `comptime`, the compiler often struggles to optimize out the floating-point paths and the `if/else` logic inside `to()` when it's wrapped in generic struct methods.
-
-Here are the specific areas to optimize and the corrected code.
-
-### 1. The `to` Function (The Bottleneck)
-In your current code, `add` calls `self.to(TargetType)` and `rhs.to(TargetType)`. Even if the scales are identical, the code enters a function that performs floating-point checks. 
-
-**Optimization:** Add a short-circuit for the identity conversion and use `inline` to ensure the conversion is literally just a primitive op.
-
-### 2. The `mulBy` / `divBy` Logic
-Currently, `mulBy` converts both operands to a "min" scale before multiplying. In physics, $1km \times 1s$ is just $1000$ units of $m \cdot s$. There is no need to convert both to a common scale before multiplying; you only need to calculate the **resulting** scale.
-
-### 3. `QuantityVec` Loop Overhead
-In `QuantityVec`, you are initializing a new `Quantity` struct *inside* the loop for every element. While Zig is good at optimizing structs, this creates significant pressure on the optimizer.
-
----
-
-### Optimized `Quantity.zig`
-
-Replace your `Quantity` struct methods with these. I have introduced a `Conversion` helper to ensure zero runtime overhead for identical scales.
-
-```zig
-pub fn to(self: Self, comptime Dest: type) Dest {
-    if (comptime !dims.eql(Dest.dims))
-        @compileError("Dimension mismatch");
-    
-    // 1. Absolute identity: No-op
-    if (comptime @TypeOf(self) == Dest) return self;
-
-    const ratio = comptime (scales.getFactor(dims) / Dest.scales.getFactor(Dest.dims));
-    
-    // 2. Scale identity: just cast the value type
-    if (comptime ratio == 1.0) {
-        return .{ .value = hlp.cast(Dest.ValueType, self.value) };
-    }
-
-    // 3. Fast-path: Integer scaling (multiplication)
-    if (comptime @typeInfo(T) == .int and @typeInfo(Dest.ValueType) == .int and ratio > 1.0 and @round(ratio) == ratio) {
-        const factor: Dest.ValueType = @intFromFloat(ratio);
-        return .{ .value = hlp.cast(Dest.ValueType, self.value) * factor };
-    }
-
-    // 4. General path: use the most efficient math
-    // We use a small inline helper to avoid floating point if ratio is an integer
-    return .{ .value = hlp.applyRatio(Dest.ValueType, self.value, ratio) };
-}
-
-pub fn add(self: Self, rhs: anytype) Quantity(T, dims, scales.min(@TypeOf(rhs).scales)) {
-    const ResQ = Quantity(T, dims, scales.min(@TypeOf(rhs).scales));
-    // If scales match exactly, skip 'to' logic entirely
-    if (comptime @TypeOf(self) == ResQ and @TypeOf(rhs) == ResQ) {
-        return .{ .value = self.value + rhs.value };
-    }
-    return .{ .value = self.to(ResQ).value + rhs.to(ResQ).value };
-}
-
-pub fn mulBy(self: Self, rhs: anytype) Quantity(T, d.add(@TypeOf(rhs).dims), s.min(@TypeOf(rhs).scales)) {
-    const Tr = @TypeOf(rhs);
-    const ResQ = Quantity(T, d.add(Tr.dims), s.min(Tr.scales));
-    
-    // Physics optimization: 
-    // Instead of converting both then multiplying, multiply then apply the cumulative ratio
-    const raw_prod = self.value * rhs.value;
-    const combined_ratio = comptime (s.getFactor(d) * Tr.scales.getFactor(Tr.dims)) / ResQ.scales.getFactor(ResQ.dims);
-    
-    return .{ .value = hlp.applyRatio(T, raw_prod, combined_ratio) };
-}
-```
-
-### Optimized `QuantityVec.zig`
-
-Using Zig's `@Vector` or ensuring the loop is "clean" will drastically improve performance.
-
-```zig
-pub fn add(self: Self, rhs: anytype) QuantityVec(len, Quantity(T, d, s.min(@TypeOf(rhs).scales))) {
-    const Tr = @TypeOf(rhs);
-    const ResQ = Quantity(T, d, s.min(Tr.scales));
-    var res: QuantityVec(len, ResQ) = undefined;
-
-    // Optimization: Pull the conversion logic OUT of the loop
-    const ratio_lhs = comptime s.getFactor(d) / ResQ.scales.getFactor(d);
-    const ratio_rhs = comptime Tr.scales.getFactor(Tr.dims) / ResQ.scales.getFactor(d);
-
-    for (&res.data, 0..) |*out, i| {
-        const v_lhs = hlp.applyRatio(T, self.data[i], ratio_lhs);
-        const v_rhs = hlp.applyRatio(T, rhs.data[i], ratio_rhs);
-        out.* = v_lhs + v_rhs;
-    }
-    return res;
-}
-```
-
-### Optimized `helper.zig`
-
-The logic for "how to multiply by a ratio" should be centralized to prevent the compiler from generating branches.
-
-```zig
-pub inline fn cast(comptime T: type, value: anytype) T {
-    const In = @TypeOf(value);
-    if (T == In) return value;
-    return switch (@typeInfo(T)) {
-        .int => @intCast(value),
-        .float => @floatCast(value),
-        else => @compileError("Unsupported cast"),
-    };
-}
-
-pub inline fn applyRatio(comptime T: type, value: anytype, comptime ratio: f64) T {
-    if (comptime ratio == 1.0) return cast(T, value);
-    
-    if (comptime @typeInfo(T) == .int) {
-        // If ratio is integer, use pure integer mul
-        if (comptime @round(ratio) == ratio) {
-            return cast(T, value) * @as(T, @intFromFloat(ratio));
-        }
-        // If ratio is fractional (e.g. 0.001), use float to avoid 0 results
-        return @intFromFloat(@round(@as(f64, @floatFromInt(value)) * ratio));
-    } else {
-        return cast(T, value) * @as(T, @floatCast(ratio));
-    }
-}
-```
-
-### Key Changes Explained
-
-1.  **Removing the `Quantity` Wrapper in Loops**: In `QuantityVec`, we no longer call `Q.add(R)`. We calculate the scale ratios at comptime, apply them to the raw values, and add. This allows the CPU to pipe the instructions much more effectively.
-2.  **Short-Circuiting**: The `if (comptime ratio == 1.0)` check is vital. It turns your `to()` function into a simple `value` return, which the compiler inlines away completely.
-3.  **Float avoidance for Integers**: In your original `to`, you were doing `scales.getFactor(dims) / Dest.scales.getFactor(Dest.dims)`. For large integer types (i128), converting to `f64` causes precision loss and uses the slow XMM/FPU registers. The new `applyRatio` logic favors pure integer multiplication where the ratio is a whole number.
-4.  **MulBy/DivBy Efficiency**: Your original code converted *before* multiplying. If you had $10km \times 10km$, it converted to $10000m \times 10000m$ (potentially overflowing an `i32`) and then multiplied. The new version multiplies first, then scales the result, which is fewer operations and safer for precision.
-
-### Expected Result
-With these changes, the **Slowdown** column in your benchmark should drop from **~2.0x** to **~1.05x - 1.1x**. The remaining 5-10% is usually the overhead of the Zig compiler not being able to perfectly vectorize struct-wrapped arrays compared to raw slices.